# mod-PATH3DU — Full Machine-Readable Reference

## Canonical Import

    import mp3du
    mp3du.version()  # returns version string

## Python Type Stubs

```python
"""Type stubs for the mp3du Rust particle tracker."""

from typing import Any, Dict, List, Optional, Tuple

import numpy as np
import numpy.typing as npt


def version() -> str:
    """Return the version string of the installed mp3du library."""
    ...


class WaterlooConfig:
    """Configuration for the Waterloo polynomial velocity-field fitting method.

    The Waterloo method fits a polynomial velocity field to the discrete
    face-by-face flows from an unstructured MODFLOW grid by minimising a
    least-squares system at a set of interior control points.

    Args:
        order_of_approx: Order of the polynomial approximation (total number
            of basis functions). Higher values capture more spatial detail but
            increase memory and solve time. Default is ``35``.
        n_control_points: Number of control points distributed inside the
            fitting domain. Must be ≥ ``order_of_approx``. Default is ``122``.
    """

    def __init__(self, order_of_approx: int = 35, n_control_points: int = 122) -> None: ...

    @property
    def order_of_approx(self) -> int:
        """Order of the polynomial approximation (number of basis functions)."""
        ...

    @property
    def n_control_points(self) -> int:
        """Number of control points used in the polynomial least-squares fit."""
        ...


class SimulationConfig:
    """Top-level simulation configuration object.

    Typically constructed from a JSON string (or file) via
    `SimulationConfig.from_json()`.  The JSON must conform to the
    [Schema Reference](../../reference/schema-reference.md).
    """

    @staticmethod
    def from_json(json_str: str) -> "SimulationConfig":
        """Deserialize a ``SimulationConfig`` from a JSON string.

        Args:
            json_str: A JSON string that conforms to the mp3du configuration schema.

        Returns:
            A validated ``SimulationConfig`` instance.

        Raises:
            SchemaError: If the JSON is malformed or fails schema validation.
        """
        ...

    def to_json(self) -> str:
        """Serialize this configuration to a JSON string.

        Returns:
            A compact JSON representation of the configuration.
        """
        ...

    def validate(self) -> None:
        """Validate the configuration against the schema.

        Raises:
            SchemaError: If any field value is out of range or inconsistent.
        """
        ...


class ParticleStart:
    """Defines a single particle's starting position and initial tracking step.

    Attributes:
        id: User-supplied integer identifier for this particle.
        x: Starting X coordinate (model units, global).
        y: Starting Y coordinate (model units, global).
        z: Starting Z coordinate — **local normalized value in [0, 1]**.
            ``0.0`` = cell bottom, ``1.0`` = cell top, ``0.5`` = layer
            midpoint.  This is **NOT** a physical elevation.  To convert:
            ``z_local = (z_physical - bot) / (top - bot)``.
            Passing a physical elevation (e.g. ``z=5.0`` when ``top=10,
            bot=0``) causes immediate ``ExitedDomain`` because ``5.0 > 1.0``.
        cell_id: Zero-based index of the MODFLOW cell that contains the
            particle's starting position.
        initial_dt: Initial time-step size for the adaptive tracker (model
            time units).  Set to a small positive value; the solver adjusts
            automatically afterwards.
    """

    id: int
    x: float
    y: float
    z: float
    cell_id: int
    initial_dt: float

    def __init__(
        self,
        id: int,
        x: float,
        y: float,
        z: float,
        cell_id: int,
        initial_dt: float,
    ) -> None: ...


class TrajectoryResult:
    """Contains the complete tracking trajectory for a single particle.

    Attributes:
        particle_id: The ``id`` supplied in the corresponding
            `ParticleStart`.
        final_status: Short string tag for the terminal state, e.g.
            ``"captured_well"``, ``"exited_domain"``, ``"max_steps"``.
        termination_reason: Human-readable description of why tracking
            stopped.
    """

    particle_id: int
    final_status: str
    termination_reason: str

    def to_records(self) -> List[Dict[str, Any]]:
        """Convert the trajectory to a list of per-step record dicts.

        Each record contains the keys ``t``, ``x``, ``y``, ``z``,
        ``cell_id``, and ``dt``.  The list can be passed directly to
        ``pandas.DataFrame(result.to_records())``.

        Returns:
            A list of dicts, one per recorded trajectory step.
        """
        ...

    def __len__(self) -> int:
        """Return the number of recorded trajectory steps."""
        ...


class GridHandle:
    """Opaque handle to a compiled unstructured grid.

    Produced by `build_grid()`.  The handle is *consumed* (invalidated) by
    `fit_waterloo()` or `fit_sspa()` — it cannot be reused after fitting.
    """

    def is_loaded(self) -> bool:
        """Return ``True`` if the grid data is still resident in memory.

        Returns ``False`` after the handle has been consumed by a fitting call.
        """
        ...

    def n_cells(self) -> int:
        """Return the number of cells in the grid."""
        ...


class WaterlooFieldHandle:
    """Opaque handle to a fitted velocity field.

    Returned by both `fit_waterloo()` and `fit_sspa()`.  Pass this handle
    to `run_simulation()` to perform particle tracking.
    """

    def is_loaded(self) -> bool:
        """Return ``True`` if the fitted field data is resident in memory."""
        ...

    def method_name(self) -> str:
        """Return the name of the fitting method used (e.g. ``"waterloo"`` or ``"sspa"``)."""
        ...

    def n_cells(self) -> int:
        """Return the number of cells in the underlying grid."""
        ...


# Alias — fit_sspa() and fit_waterloo() both return the same handle type.
VelocityFieldHandle = WaterlooFieldHandle


class SspaConfig:
    """Configuration for the SSP&A (kriging-based) velocity-field fitting method.

    Args:
        search_radius: Neighbourhood search radius used when selecting
            conditioning data for each kriging estimate (model units).
        krig_offset: Small offset added to the kriging variogram nugget to
            improve numerical stability. Default is ``0.1``.
    """

    def __init__(self, search_radius: float, krig_offset: float = 0.1) -> None: ...

    @property
    def search_radius(self) -> float:
        """Neighbourhood search radius for kriging (model units)."""
        ...

    @property
    def krig_offset(self) -> float:
        """Variogram nugget offset for numerical stability."""
        ...


class SspaInputs:
    """Hydrated inputs for the SSP&A fitting method.

    Produced by `hydrate_sspa_inputs()`.  Pass to `fit_sspa()`.
    """

    def n_cells(self) -> int:
        """Return the number of cells represented in these inputs."""
        ...


def hydrate_sspa_inputs(
    heads: npt.NDArray[np.float64],
    porosity: npt.NDArray[np.float64],
    well_mask: npt.NDArray[np.bool_],
    hydraulic_conductivity: Optional[npt.NDArray[np.float64]] = None,
    hhk: Optional[npt.NDArray[np.float64]] = None,
) -> SspaInputs:
    """Hydrate SSP&A fitting inputs from NumPy arrays.

    Args:
        heads: Simulated hydraulic heads, shape ``(n_cells,)``.
        porosity: Effective porosity per cell, shape ``(n_cells,)``.
        well_mask: Boolean array marking cells that contain a well,
            shape ``(n_cells,)``.
        hydraulic_conductivity: Horizontal hydraulic conductivity per cell,
            shape ``(n_cells,)``.  Provide exactly one of
            ``hydraulic_conductivity`` or ``hhk``.
        hhk: Alias for ``hydraulic_conductivity`` (legacy shorthand).
            Provide exactly one of ``hydraulic_conductivity`` or ``hhk``.

    Returns:
        An `SspaInputs` object ready to pass to `fit_sspa()`.
    """
    ...


def fit_sspa(
    config: SspaConfig,
    grid: GridHandle,
    inputs: SspaInputs,
    drifts: List[Dict[str, Any]],
) -> VelocityFieldHandle:
    """Fit the SSP&A (kriging-based) velocity field.

    Consumes the GridHandle (it will no longer be loaded after the call).
    """
    ...


class CellProperties:
    """Per-cell physical properties used during velocity-field fitting and tracking.

    All arrays are indexed by zero-based cell ID.  Produced by
    `hydrate_cell_properties()`.
    """

    @property
    def top(self) -> List[float]:
        """Top elevation of each cell (model units)."""
        ...

    @property
    def bot(self) -> List[float]:
        """Bottom elevation of each cell (model units)."""
        ...

    @property
    def porosity(self) -> List[float]:
        """Effective porosity of each cell (dimensionless, 0–1)."""
        ...

    @property
    def retardation(self) -> List[float]:
        """Retardation factor of each cell (dimensionless, ≥ 1)."""
        ...

    @property
    def hhk(self) -> List[float]:
        """Horizontal hydraulic conductivity of each cell (L/T)."""
        ...

    @property
    def vhk(self) -> List[float]:
        """Vertical hydraulic conductivity of each cell (L/T)."""
        ...

    @property
    def disp_long(self) -> List[float]:
        """Longitudinal dispersivity of each cell (L)."""
        ...

    @property
    def disp_trans_h(self) -> List[float]:
        """Horizontal transverse dispersivity of each cell (L)."""
        ...

    @property
    def disp_trans_v(self) -> List[float]:
        """Vertical transverse dispersivity of each cell (L)."""
        ...

    def n_cells(self) -> int:
        """Return the number of cells."""
        ...


class CellFlows:
    """Per-cell flow data (heads, budgets, face flows) used during particle tracking.

    All arrays are indexed by zero-based cell ID.  Produced by
    `hydrate_cell_flows()`.

    !!! note "Sign conventions"
        Pass raw MODFLOW values — **do not negate**.
        ``face_flow`` positive = into cell (MODFLOW-USG / MF6 convention).
        ``q_well`` negative = extraction, positive = injection.
    """

    @property
    def head(self) -> List[float]:
        """Simulated hydraulic head in each cell (model units)."""
        ...

    @property
    def water_table(self) -> List[float]:
        """Water-table elevation in each cell (model units)."""
        ...

    @property
    def q_top(self) -> List[float]:
        """Volumetric flux entering through the cell top (L³/T, raw MODFLOW sign)."""
        ...

    @property
    def q_bot(self) -> List[float]:
        """Volumetric flux entering through the cell bottom (L³/T, raw MODFLOW sign)."""
        ...

    @property
    def q_vert(self) -> List[float]:
        """Net vertical flow term (L³/T, raw MODFLOW sign)."""
        ...

    @property
    def q_well(self) -> List[float]:
        """Well volumetric flux (L³/T). Negative = extraction, positive = injection."""
        ...

    @property
    def q_other(self) -> List[float]:
        """Sum of all other (non-well, non-face) boundary fluxes (L³/T)."""
        ...

    @property
    def q_storage(self) -> List[float]:
        """Storage term (L³/T, positive = release from storage)."""
        ...

    @property
    def has_well(self) -> List[bool]:
        """``True`` for each cell that contains an active well."""
        ...

    @property
    def face_offset(self) -> List[int]:
        """Start index into ``face_flow`` / ``face_neighbor`` for each cell."""
        ...

    @property
    def face_flow(self) -> List[float]:
        """Face-by-face volumetric fluxes (L³/T).  Positive = into cell."""
        ...

    @property
    def face_neighbor(self) -> List[Optional[int]]:
        """Zero-based cell index of the neighbouring cell across each face,
        or ``None`` for domain-boundary faces."""
        ...

    def n_cells(self) -> int:
        """Return the number of cells."""
        ...


class WaterlooInputs:
    """Hydrated inputs for the Waterloo polynomial velocity-field fitting method.

    Produced by `hydrate_waterloo_inputs()`.  Pass to `fit_waterloo()`.

    !!! note "Sign convention"
        ``face_flow`` here uses the same convention as ``CellFlows``
        (positive = INTO cell).  Pass the **same** ``face_flow`` array
        to both ``hydrate_cell_flows()`` and ``hydrate_waterloo_inputs()``.
    """

    def n_cells(self) -> int:
        """Return the number of cells represented in these inputs."""
        ...

    def __len__(self) -> int:
        """Return the number of cells (same as ``n_cells()``)."""
        ...


def hydrate_cell_properties(
    top: npt.NDArray[np.float64],
    bot: npt.NDArray[np.float64],
    porosity: npt.NDArray[np.float64],
    retardation: npt.NDArray[np.float64],
    hhk: npt.NDArray[np.float64],
    vhk: npt.NDArray[np.float64],
    disp_long: npt.NDArray[np.float64],
    disp_trans_h: npt.NDArray[np.float64],
    disp_trans_v: npt.NDArray[np.float64],
) -> CellProperties:
    """Hydrate per-cell physical properties from NumPy arrays.

    All arrays must have shape ``(n_cells,)`` and be in cell-ID order.

    Args:
        top: Top elevation of each cell (model units).
        bot: Bottom elevation of each cell (model units).
        porosity: Effective porosity (dimensionless, 0–1).
        retardation: Retardation factor (dimensionless, ≥ 1).
        hhk: Horizontal hydraulic conductivity (L/T).
        vhk: Vertical hydraulic conductivity (L/T).
        disp_long: Longitudinal dispersivity (L).
        disp_trans_h: Horizontal transverse dispersivity (L).
        disp_trans_v: Vertical transverse dispersivity (L).

    Returns:
        A `CellProperties` object ready for use in `fit_waterloo()`.
    """
    ...


def hydrate_cell_flows(
    head: npt.NDArray[np.float64],
    water_table: npt.NDArray[np.float64],
    q_top: npt.NDArray[np.float64],
    q_bot: npt.NDArray[np.float64],
    q_vert: npt.NDArray[np.float64],
    q_well: npt.NDArray[np.float64],
    q_other: npt.NDArray[np.float64],
    q_storage: npt.NDArray[np.float64],
    has_well: npt.NDArray[np.bool_],
    face_offset: npt.NDArray[np.uint64],
    face_flow: npt.NDArray[np.float64],
    face_neighbor: npt.NDArray[np.int64],
    # Optional boundary condition metadata (IFACE-based capture)
    bc_cell_ids: Optional[npt.NDArray[np.int64]] = None,
    bc_iface: Optional[npt.NDArray[np.int32]] = None,
    bc_flow: Optional[npt.NDArray[np.float64]] = None,
    bc_type_id: Optional[npt.NDArray[np.int32]] = None,
    bc_type_names: Optional[List[str]] = None,
    is_domain_boundary: Optional[npt.NDArray[np.bool_]] = None,
    has_water_table: Optional[npt.NDArray[np.bool_]] = None,
) -> CellFlows:
    """Hydrate cell flow data from NumPy arrays.

    mp3du API convention (target):

    - ``face_flow``: positive = **into** cell.
    - ``q_well``: negative = extraction, positive = injection (raw MODFLOW sign).

    The transformation from raw MODFLOW output depends on version:

    - **MODFLOW-USG / MF6** (``FLOW-JA-FACE``): raw is positive = IN.
      **Pass directly**: ``face_flow = flowja``.
    - **MODFLOW-2005 / NWT** (after directional→per-face assembly):
      result is positive = OUT. **Negate**: ``face_flow = -assembled``.

    Pass the **same** ``face_flow`` array to both ``hydrate_cell_flows()``
    and ``hydrate_waterloo_inputs()``.  No per-function negation needed.

    Never negate ``q_well``.
    See ``docs/reference/units-and-conventions.md`` for the full reference.
    """
    ...


def hydrate_waterloo_inputs(
    centers_xy: npt.NDArray[np.float64],
    radii: npt.NDArray[np.float64],
    perimeters: npt.NDArray[np.float64],
    areas: npt.NDArray[np.float64],
    q_vert: npt.NDArray[np.float64],
    q_well: npt.NDArray[np.float64],
    q_other: npt.NDArray[np.float64],
    face_offset: npt.NDArray[np.uint64],
    face_vx1: npt.NDArray[np.float64],
    face_vy1: npt.NDArray[np.float64],
    face_vx2: npt.NDArray[np.float64],
    face_vy2: npt.NDArray[np.float64],
    face_length: npt.NDArray[np.float64],
    face_flow: npt.NDArray[np.float64],
    noflow_mask: npt.NDArray[np.bool_],
) -> WaterlooInputs:
    """Hydrate Waterloo velocity-fitting inputs from NumPy arrays.

    mp3du API convention (target):

    - ``face_flow``: positive = **INTO** cell — same as ``hydrate_cell_flows()``.
    - ``q_well``: raw MODFLOW sign (negative = extraction). **Do NOT negate.**

    Pass the **same** ``face_flow`` array used for ``hydrate_cell_flows()``.
    No per-function negation or sign flip is needed.

    The Waterloo method subtracts the analytic well singularity during fitting
    and adds it back during evaluation; both must use the same ``q_well`` sign.
    See ``docs/reference/units-and-conventions.md`` for the full reference.
    """
    ...


def build_grid(
    vertices: List[List[Tuple[float, float]]],
    centers: List[Tuple[float, float, float]],
) -> GridHandle:
    """Compile an unstructured polygon grid into a `GridHandle`.

    Args:
        vertices: For each cell, an ordered list of ``(x, y)`` vertex
            coordinates defining the cell polygon.
        centers: For each cell, a ``(x, y, z)`` tuple giving the cell
            centroid (x, y) and representative elevation (z).

    Returns:
        A `GridHandle` that can be passed to `fit_waterloo()` or
        `fit_sspa()`.  The handle is consumed (invalidated) by the
        fitting call.
    """
    ...


def fit_waterloo(
    config: WaterlooConfig,
    grid: GridHandle,
    fit_inputs: WaterlooInputs,
    cell_properties: CellProperties,
    cell_flows: CellFlows,
) -> WaterlooFieldHandle:
    """Fit the Waterloo velocity field.

    Consumes the GridHandle (it will no longer be loaded after the call).
    Well locations are derived automatically from cell_flows.has_well
    and the grid cell centres.

    Sign conventions: fit_inputs must use face_flow positive = INTO cell
    (same array as hydrate_cell_flows; q_well = raw MODFLOW sign).
    See hydrate_waterloo_inputs() and docs/reference/units-and-conventions.md.
    """
    ...


def run_simulation(
    config: SimulationConfig,
    velocity_field: VelocityFieldHandle,
    particles: List[ParticleStart],
    parallel: bool = True,
) -> List[TrajectoryResult]:
    """Run particle tracking on a fitted velocity field.

    Well capture is controlled by capture.capture_radius in the
    SimulationConfig:
      - Omitted / null: capture immediately on cell entry (strong sink).
      - A positive number (e.g. 0.5): capture only when within that
        distance of the cell centre (weak sink).
    """
    ...

```

## JSON Schema

```json
{
  "$schema": "http://json-schema.org/draft-07/schema#",
  "title": "MP3DU SimulationConfig",
  "type": "object",
  "required": [
    "velocity_method",
    "solver",
    "adaptive",
    "dispersion",
    "retardation_enabled",
    "capture",
    "initial_dt",
    "max_dt",
    "direction"
  ],
  "properties": {
    "velocity_method": {
      "type": "string",
      "enum": ["Waterloo"]
    },
    "solver": {
      "type": "string",
      "enum": [
        "Euler",
        "Rk4StepDoubling",
        "DormandPrince",
        "CashKarp",
        "VernerRobust",
        "VernerEfficient"
      ]
    },
    "adaptive": {
      "type": "object",
      "required": [
        "tolerance",
        "safety",
        "alpha",
        "min_scale",
        "max_scale",
        "max_rejects",
        "min_dt",
        "euler_dt"
      ],
      "properties": {
        "tolerance": { "type": "number", "exclusiveMinimum": 0 },
        "safety": { "type": "number", "exclusiveMinimum": 0 },
        "alpha": { "type": "number", "exclusiveMinimum": 0 },
        "min_scale": { "type": "number", "exclusiveMinimum": 0 },
        "max_scale": { "type": "number", "exclusiveMinimum": 0 },
        "max_rejects": { "type": "integer", "minimum": 0 },
        "min_dt": { "type": "number", "exclusiveMinimum": 0 },
        "euler_dt": { "type": "number", "exclusiveMinimum": 0 }
      },
      "additionalProperties": false
    },
    "dispersion": {
      "oneOf": [
        {
          "type": "object",
          "required": ["method"],
          "properties": {
            "method": { "const": "None" }
          },
          "additionalProperties": false
        },
        {
          "type": "object",
          "required": ["method", "alpha_l", "alpha_th", "alpha_tv"],
          "properties": {
            "method": { "const": "Gsde" },
            "alpha_l": { "type": "number", "minimum": 0 },
            "alpha_th": { "type": "number", "minimum": 0 },
            "alpha_tv": { "type": "number", "minimum": 0 }
          },
          "additionalProperties": false
        },
        {
          "type": "object",
          "required": ["method", "alpha_l", "alpha_th", "alpha_tv"],
          "properties": {
            "method": { "const": "Ito" },
            "alpha_l": { "type": "number", "minimum": 0 },
            "alpha_th": { "type": "number", "minimum": 0 },
            "alpha_tv": { "type": "number", "minimum": 0 }
          },
          "additionalProperties": false
        }
      ]
    },
    "retardation_enabled": {
      "type": "boolean"
    },
    "capture": {
      "type": "object",
      "required": [
        "max_time",
        "max_steps",
        "stagnation_velocity",
        "stagnation_limit"
      ],
      "properties": {
        "max_time": { "type": "number", "exclusiveMinimum": 0 },
        "max_steps": { "type": "integer", "minimum": 1 },
        "stagnation_velocity": { "type": "number", "minimum": 0 },
        "stagnation_limit": { "type": "integer", "minimum": 1 },
        "capture_radius": {
          "type": "number",
          "exclusiveMinimum": 0,
          "description": "Distance from well centre (cell centre) within which a particle is captured. Omit or set to null for immediate capture on cell entry."
        },
        "face_epsilon": {
          "type": "number",
          "exclusiveMinimum": 0,
          "default": 1e-6,
          "description": "Tolerance for top/bottom face proximity checks during IFACE-based capture. Particle is considered at the face when z > 1 - face_epsilon (top) or z < face_epsilon (bottom). Default: 1e-6."
        }
      },
      "additionalProperties": false
    },
    "initial_dt": { "type": "number", "exclusiveMinimum": 0 },
    "max_dt": { "type": "number", "exclusiveMinimum": 0 },
    "direction": { "type": "number", "enum": [1.0, -1.0] }
  },
  "additionalProperties": false
}

```

## Minimal Valid Configuration

```json
{
  "velocity_method": "Waterloo",
  "solver": "DormandPrince",
  "adaptive": {
    "tolerance": 1e-6,
    "safety": 0.9,
    "alpha": 0.2,
    "min_scale": 0.2,
    "max_scale": 5.0,
    "max_rejects": 10,
    "min_dt": 1e-10,
    "euler_dt": 1.0
  },
  "dispersion": {
    "method": "None"
  },
  "retardation_enabled": false,
  "capture": {
    "max_time": 365250.0,
    "max_steps": 1000000,
    "stagnation_velocity": 1e-12,
    "stagnation_limit": 100
  },
  "initial_dt": 1.0,
  "max_dt": 100.0,
  "direction": 1.0
}

```

## Complete Configuration

```json
{
  "velocity_method": "Waterloo",
  "solver": "DormandPrince",
  "adaptive": {
    "tolerance": 1e-6,
    "safety": 0.9,
    "alpha": 0.2,
    "min_scale": 0.2,
    "max_scale": 5.0,
    "max_rejects": 10,
    "min_dt": 1e-10,
    "euler_dt": 1.0
  },
  "dispersion": {
    "method": "Gsde",
    "alpha_l": 10.0,
    "alpha_th": 1.0,
    "alpha_tv": 0.1
  },
  "retardation_enabled": true,
  "capture": {
    "max_time": 365250.0,
    "max_steps": 1000000,
    "stagnation_velocity": 1e-12,
    "stagnation_limit": 100,
    "face_epsilon": 1e-6
  },
  "initial_dt": 1.0,
  "max_dt": 100.0,
  "direction": 1.0
}

```

## SSP&A Minimal Workflow

```python
import json
import numpy as np
import mp3du

# 1. Build grid (example: single square cell)
vertices = [[(0,0), (100,0), (100,100), (0,100)]]
centers = [(50.0, 50.0, 0.0)]
grid = mp3du.build_grid(vertices, centers)

# 2. Hydrate SSP&A inputs (all arrays shape (n_cells,))
n = grid.n_cells()
inputs = mp3du.hydrate_sspa_inputs(
    heads=np.array([10.0]),
    porosity=np.array([0.3]),
    well_mask=np.array([True]),
    hhk=np.array([100.0]),          # OR hydraulic_conductivity=
)

# 3. Define drifts (list of dicts)
drifts = [
    {"type": "well", "event": "SS", "term": 1, "name": "W1",
      "value": -500.0, "x": 50.0, "y": 50.0},
]

# 4. Fit velocity field (consumes grid!)
cfg = mp3du.SspaConfig(search_radius=300.0, krig_offset=0.1)
field = mp3du.fit_sspa(cfg, grid, inputs, drifts)
# grid is now invalid — do not reuse

# 5. Run simulation
sim_cfg = mp3du.SimulationConfig.from_json(json.dumps({
    "velocity_method": "Waterloo",
    "solver": "DormandPrince",
    "adaptive": {"tolerance": 1e-6},
    "capture": {"max_time": 1000.0, "max_steps": 100000},
    "initial_dt": 0.1,
    "direction": 1.0,
}))
particles = [mp3du.ParticleStart(id=0, x=25.0, y=25.0, z=0.5, cell_id=0, initial_dt=0.1)]
results = mp3du.run_simulation(sim_cfg, field, particles)
for r in results:
    print(r.final_status, len(r.to_records()), "steps")
```

## SSP&A Drift Type Reference

Supported drift types for fit_sspa() drifts parameter:

| Type | Aliases | Required Keys | Capture |
|------|---------|--------------|---------|
| well | — | type, event, term, name, value, x, y | ✅ Implemented |
| linesink | line_sink, line-sink | type, event, term, name, value, x1, y1, x2, y2 | ❌ Deferred |
| noflow | no_flow, no-flow | type, event, term, name, value, x1, y1, x2, y2 | ❌ Deferred |

Line segments with the same (type, event, term, name) are grouped into one drift.
Mixed values within a group are rejected. Zero-length segments are rejected.

## Error Patterns

Placeholder — populated after error catalog is built.