# mod-PATH3DU

> 3D particle tracking engine for groundwater flow models.
> Rust core with Python bindings via PyO3.

## Import

    import mp3du  # canonical import

## Key Types

- SimulationConfig: simulation parameters (from JSON or Python)
- WaterlooConfig: velocity interpolation configuration
- ParticleStart: initial particle position and cell
- TrajectoryResult: particle path and final status
- GridHandle: loaded unstructured grid
- WaterlooFieldHandle: fitted velocity field
- CellProperties: per-cell physical properties
- CellFlows: per-cell flow budget arrays
- SspaConfig: SSP&A (kriging-based) velocity fitting configuration
- SspaInputs: hydrated SSP&A inputs (heads, porosity, conductivity, well mask)
- VelocityFieldHandle: alias for WaterlooFieldHandle (returned by both fit methods)

## Key Functions

- hydrate_sspa_inputs(): build SspaInputs from NumPy arrays
- fit_sspa(): fit SSP&A velocity field (consumes GridHandle)
- hydrate_cell_properties(): build CellProperties from NumPy arrays
- hydrate_cell_flows(): build CellFlows from NumPy arrays
- hydrate_waterloo_inputs(): build WaterlooInputs from NumPy arrays
- build_grid(): build GridHandle from vertices and centers
- fit_waterloo(): fit Waterloo velocity field (consumes GridHandle)
- run_simulation(): run particle tracking on a fitted field

## Links

- Full reference: https://sspa-inc.github.io/mp3du-rs-docs/reference/
- Schema: https://sspa-inc.github.io/mp3du-rs-docs/reference/schema-reference/
- Python API: https://sspa-inc.github.io/mp3du-rs-docs/reference/python-api/
- Quickstart: https://sspa-inc.github.io/mp3du-rs-docs/getting-started/quickstart/
- Examples: https://sspa-inc.github.io/mp3du-rs-docs/examples/
- Full machine-readable reference: https://sspa-inc.github.io/mp3du-rs-docs/llms-full.txt

## ⚠️ Critical Pitfall: Water Table vs Layer Type

The `water_table` array in `hydrate_cell_flows()` controls saturated thickness
which directly divides velocity. It MUST match the MODFLOW layer type (LAYTYP):

    LAYTYP==0 (confined)    → water_table = top   (NOT head!)
    LAYTYP==1 (unconfined)  → water_table = head
    LAYTYP>0  (convertible) → water_table = min(head, top)

Setting water_table=head for confined layers silently produces 2×–17× velocity
errors. See llms-full.txt and docs/reference/units-and-conventions.md.

## ⚠️ Critical Pitfall: Face Flow Signs Depend on MODFLOW Version

Both hydrate_cell_flows() and hydrate_waterloo_inputs() expect face_flow
with positive = INTO cell (flow entering the cell through that face).
Convert once, then pass the SAME array to both.

    MODFLOW-USG / MF6 (FLOW-JA-FACE): raw positive = INTO cell
        face_flow = flowja                   # pass directly — already positive = INTO

    MODFLOW-2005 / NWT: raw stores DIRECTIONAL flows (positive = in direction
    of increasing row/column), NOT per-cell flows.  You must assemble into
    per-cell per-face flows first (flip sign for reverse-direction neighbors).
    After assembly the result is positive = OUT — negate once:
        face_flow = -assembled               # negate once → positive = INTO

    # Pass face_flow to BOTH functions — no per-function negation.

Never negate q_well. See llms-full.txt for the full conversion reference.

## ⚠️ Critical Pitfall: Vertex Winding and Face Index Alignment

All cell polygons passed to `build_grid()` MUST be wound Clockwise (CW)
when viewed from above. There is NO runtime validation — wrong
winding silently inverts face normals and produces a 180° reversed
velocity field. If your source uses CCW winding, reverse the vertex list
and negate all face_flow values.

CRITICAL: face indexing is defined by vertex order. If you reorder
vertices, you MUST reorder every face-level array the same way
(`face_flow`, `face_neighbor`, `face_vx1/vy1/vx2/vy2`, `face_length`,
`noflow_mask`, and any source-specific face/neighbour IDs).
Reversing only the vertex list while still iterating neighbour/face IDs
in the original order misaligns flows with geometric edges and can flip
velocity directions even when the scalar flow signs are otherwise correct.

## ⚠️ Critical Pitfall: Z Coordinate is Local [0, 1]

ParticleStart.z and trajectory output z are LOCAL normalized coordinates
(0 = cell bottom, 1 = cell top, 0.5 = midpoint), NOT physical elevations.
Passing a physical elevation (e.g. z=5.0 when top=10, bot=0) causes
immediate ExitedDomain.

    z_local = (z_physical - bot) / (top - bot)
    z_physical = bot + z_local * (top - bot)

## ⚠️ Critical Pitfall: SSP&A Conductivity Argument Exclusivity

hydrate_sspa_inputs() accepts exactly one of `hydraulic_conductivity=` or `hhk=`.
Passing both raises ValueError. Omitting both raises ValueError.
`hhk` is the legacy/domain shorthand for horizontal hydraulic conductivity.

## ⚠️ Critical Pitfall: SSP&A Array Shapes

All arrays passed to hydrate_sspa_inputs() must have shape (n_cells,):
  heads, porosity, well_mask, and conductivity (hhk or hydraulic_conductivity).
Mismatched lengths raise ValueError at hydrate time or fit time.

## ⚠️ Critical Pitfall: well_mask Is Cell-Based

well_mask is a boolean array of length n_cells (one entry per grid cell),
NOT a list of well indices or a list matching len(drifts).
Set True for cells that contain a well. The mask is used to exclude
well cells from kriging conditioning data.

## ⚠️ Critical Pitfall: fit_sspa() Consumes the GridHandle

After calling fit_sspa(), the GridHandle is invalidated and cannot be
reused. If you need the grid for multiple fits, rebuild it.

## ⚠️ Critical Pitfall: SSP&A Supported Drift Types

Only three drift types are supported: well, linesink, noflow.
Do not invent other type values. Unsupported types raise ValueError.
  - well: requires x, y
  - linesink: requires x1, y1, x2, y2
  - noflow: requires x1, y1, x2, y2
All drifts require: type, event, term, name, value.

## ⚠️ Critical Pitfall: SSP&A Performance

fit_sspa() neighborhood precomputation is currently O(n²) in cell count.
For a 201×201 grid (40,401 cells), fitting takes ~350 seconds.
Plan accordingly for large grids.

## ⚠️ Critical Pitfall: SSP&A Line/NoFlow Capture Not Yet Implemented

Runtime capture (particle interception) is implemented for well drifts only.
Line-sink and no-flow drifts affect the velocity field but do NOT capture
particles at runtime. This is a known limitation.