Metadata-Version: 2.4
Name: peclet-pnm
Version: 1.0.2
Summary: peclet.pnm — Kokkos pore-network extraction from SDF geometry (pores, watershed segmentation, throat topology)
Author-Email: Frank Peters <e.a.j.f.peters@gmail.com>
License-Expression: MIT
License-File: LICENSE
Classifier: Development Status :: 5 - Production/Stable
Classifier: Intended Audience :: Science/Research
Classifier: Operating System :: POSIX :: Linux
Classifier: Programming Language :: C++
Classifier: Programming Language :: Python :: 3
Classifier: Topic :: Scientific/Engineering :: Physics
Project-URL: Homepage, https://github.com/computational-chemical-engineering/peclet
Project-URL: Documentation, https://github.com/computational-chemical-engineering/peclet
Project-URL: Source, https://github.com/computational-chemical-engineering/peclet-pnm
Requires-Python: >=3.10
Requires-Dist: numpy>=1.20
Description-Content-Type: text/markdown

# peclet-pnm

**`peclet.pnm` — GPU pore-network extraction from SDF geometry.**

Given a signed-distance-field (SDF) description of a porous solid (negative inside the solid,
positive in the pore space), `peclet.pnm` extracts the pore network:

- **`SDFReader`** — pure-C++ VTI (VTK ImageData) reader for SDF volumes.
- **`extract_pores`** — pore detection: local maxima of the SDF + weighted centroids and radii.
- **`segment_volume`** — marker-controlled watershed segmentation of the pore space
  (marker init → union-find connected-component labelling → flood fill).
- **`extract_topology`** — label adjacency (pore-to-pore throats and pore-solid contacts) from boundary pairs between basins.
- **`extract_pore_network`** — the fused pipeline (SDF uploaded once, segmentation device-resident
  across all three stages): returns `(pores, segmentation, connections)` in one call.

The compute is [Kokkos](https://github.com/kokkos/kokkos) — the same source runs on **CUDA, HIP, and
OpenMP** backends, selected at build time by the install prefix. Part of the
[peclet](https://github.com/computational-chemical-engineering/peclet) suite; split out of
[peclet-flow](https://github.com/computational-chemical-engineering/peclet-flow) (its former
`peclet.flow.pnm` module — the repo's original "pnm_from_sdf" feature).

## Install / build

```bash
# From the peclet suite checkout (Kokkos prefix bootstrapped once by ../tools/bootstrap_deps.sh):
CMAKE_PREFIX_PATH="$PWD/../extern/install/nvidia-cuda" pip install .

# Or a dev cmake build (nanobind found via the active interpreter):
cmake -S . -B build -DCMAKE_PREFIX_PATH="$PWD/../extern/install/nvidia-cuda"
cmake --build build -j            # -> build/peclet/pnm/_pnm.*.so ; PYTHONPATH=$PWD/build to import

# Tests (one tree per backend): the single-rank C++ contract on synthetic SDFs, the Python binding
# smoke, the 7199-pore packing_ring gate (SKIPPED by ctest when ../flow/data/packing_ring.vti is
# absent) and, with -DPECLET_PNM_MPI=ON, the distributed np=1,2,4 suite.
cmake -S . -B build_dev -DCMAKE_PREFIX_PATH="$PWD/../extern/install/nvidia-cuda" \
  -DPECLET_PNM_MPI=ON -DPECLET_PNM_BUILD_TESTS=ON -DMPIEXEC_EXECUTABLE=/usr/bin/mpirun
cmake --build build_dev -j && OMP_PROC_BIND=false ctest --test-dir build_dev --output-on-failure
```

Without a Kokkos prefix on `CMAKE_PREFIX_PATH`, the build vendors Kokkos (OpenMP+Serial) via
FetchContent, so `pip install .` works standalone on any Linux with a C++20 toolchain.

## Usage

```python
import peclet.pnm as pnm

sdf_3d, origin_zyx, spacing_zyx = pnm.SDFReader.read_vti("packing.vti")  # (Nz,Ny,Nx) C-order
pores = pnm.extract_pores(sdf_3d, origin_zyx, spacing_zyx)               # Pore(x,y,z,radius) list
seg = pnm.segment_volume(sdf_3d, spacing_zyx)                            # int32 (Nz,Ny,Nx) labels
conns = pnm.extract_topology(seg)                                        # (M,2) int32 label pairs
throats = conns[(conns[:, 0] > 0) & (conns[:, 1] > 0)]                   # pore-pore pairs only

# or fused (SDF uploaded once, segmentation stays device-resident across stages):
pores, seg, conns = pnm.extract_pore_network(sdf_3d, origin_zyx, spacing_zyx)
```

Conventions: the SDF array is `(Nz, Ny, Nx)` C-order (x fastest) and every triple that describes
it (`origin_zyx`, `spacing_zyx`, `shape_zyx`, `grad_p_zyx`) is z-y-x, marked by the suffix; SDF
sign is negative inside the solid — see the suite's `docs/CONVENTIONS.md` and `docs/NAMING.md` §1.7.
(`Pore.x/y/z` are three self-named scalars, not a triple, so they carry no suffix.) Arrays in,
arrays out: `segment_volume` returns the labels as an int32 array of the SDF's shape (the kernels'
flat x-fastest vector re-shaped in place, no copy), `extract_topology` reads that array back
without a copy and returns the pairs as an `(M, 2)` int32 array, and the network-flow dict holds
NumPy arrays; only the pores are a Python `list[Pore]`.
Labels from `segment_volume`: pores `1, 2, …`, solid grains `-1, -2, …`, `0` = unreached solid
debris. Precision: the SDF is float32 and the geometry kernels compute in float32 (`origin_zyx` /
`spacing_zyx` are narrowed to float32, so pore centres and radii are float32 in the input unit
system); the network-flow MAC fields are float64.

Smoke tests: `python scripts/test_extraction.py <sdf.vti>` and
`python scripts/verify_segmentation.py <sdf.vti>` (writes a labelled `.vti` + a pore-pair edge list).

## Network flow: throat flow rates + pore pressures from a DNS

`extract_network_flow` turns a converged [peclet-flow](https://github.com/computational-chemical-engineering/peclet-flow)
velocity/pressure field on the same grid into pore-network flow data — the method carried over
from the Voronoi-tessellation PNM of sphere packings (`pnm_voronoi`), where the throat flow was
∫u·n over the Voronoi facet and the pore pressure a trilinear sample at the pore center:

```python
s = peclet.flow.Solver(nx, ny, nz)
...; s.set_body_force(fx, 0, 0); s.set_solid(sdf_xyz, cutcell_pressure=True); ...steps...
net = pnm.extract_network_flow(
    sdf_zyx, origin_zyx, spacing_zyx,
    s.get_uf().T, s.get_vf().T, s.get_wf().T, s.get_p().T,   # zero-copy transposes to zyx
    s.get_ox().T, s.get_oy().T, s.get_oz().T,                # cut-cell face openness
    grad_p_zyx=[0, 0, -fx])                                  # body force f == -grad p_macro
net["pores"]           # list[Pore] in label order (pores[k] is label k+1)
net["throats"]         # (M,2) int32 label pairs a < b, one row per interface PATCH
net["throat_flow"]     # Q through each pore-pore interface (o·u·A summed over MAC faces)
net["pore_pressure"]   # periodic p interpolated at each pore center (basin SDF peak)
net["throat_dp"]       # total-pressure drop P_i - P_j (periodic parts + macro gradient
                       # along the throat-anchored min-image path)
net["pore_residual"]   # signed flux over each pore's whole boundary — ~ solver tolerance
```

The openness arrays must be the ones the velocity field was projected with: with peclet.flow's
cut-cell IBM that requires `set_solid(..., cutcell_pressure=True)` — without it every openness
flow reports is 0 and every throat flux comes back 0 (the binding cannot check this; the
precondition lives in the docstring).

On the voxel network the throat integral is exact: a throat is a set of grid-aligned MAC faces
and the openness-weighted face velocity is the discrete flux carrier, so per-pore mass balance
holds to the pressure-solve tolerance (`pore_residual` is the built-in check). Fluxes are
accumulated on **flow basins** (gradient-ascent assignment of *every* cell, including cut cells
whose center is inside the solid) — keyed on the segmentation labels alone, the near-wall
staircase flux would bypass the interface (measured 6% on a tube).

**Both IBM variants are supported.** With the cut-cell IBM the bookkeeping is machine-exact
(pass `get_ox()...`). With the **ghost-cell IBM** (`set_ghost_projection(True)`) pass flow's
`get_ox_proj()/get_oy_proj()/get_oz_proj()` — the binary (COUPLED) openness the ghost projection
conserves. Ghost-cell IBM is pointwise 2nd-order but not locally mass-conserving at the wall, so
there the network data is truncation-accurate: `pore_residual` becomes the per-pore wall leak
(measured 3.2e-2·F at a 4-cell tube radius, converging at order ~2.7 under refinement).

**MPI:** `extract_network_flow_mpi(sdf_local, global_shape_zyx, ..., u_local, ...)` runs the whole
pipeline distributed on the core ORB blocks (fields from a distributed peclet.flow run on the
same decomposition); every rank returns the identical global network. Matches the single-rank
result to accumulation-order tolerance (`tests/kokkos_mpi/test_pnm_flow_mpi`, np = 1, 2, 4).

Validated in `scripts/verify_network_flow.py` (chamber-tube chain + asymmetric tube lattice +
the ghost-IBM chain, DNS by peclet.flow): every throat carries the DNS flux to ~1e-11 relative
(cut-cell), residuals ~1e-12·F, g = Q/dp > 0 on all throats, and the dp sum around each loop
equals the macroscopic drop.
`scripts/demo_network_flow_packing.py` runs the pipeline on a real sphere packing.

**Throats are per-patch:** a throat is a *connected patch* of interface faces (CCL over the
interface, core faces = both cells fluid-centered, wall-film faces attached by propagation), so
two disjoint interfaces between the same two pores — e.g. two parallel tubes, or a direct contact
plus one through the periodic wrap — are separate parallel throats and the throat list can repeat
a pore pair (validated: two capsules of different radii report two (1,2) throats whose fluxes sum
to the DNS flux exactly). Remaining caveat: on loose packings (porosity ≳ 0.6) intra-pore
pressure variation is comparable to throat drops, so per-throat g = Q/dp scatters — a property of
the point-pressure PNM abstraction, not of the extraction.

## Distributed (MPI) extraction

Built with `-DPECLET_PNM_MPI=ON`, the module also runs the whole pipeline **multi-rank**: the SDF
is decomposed over ranks by the shared peclet-core ORB (the same deterministic partition flow/dem
use), every stage runs per-rank on a 1-cell ghost layer (core `GridHalo` exchange), and the result
is **bit-exact to the single-rank pipeline** — labels are global voxel ids, so the CCL fixpoint,
the watershed flood (Jacobi), the gradient-path pore basins, and the renumbering are all
decomposition-independent.

```python
# mpirun -np 4 python extract.py
import peclet.pnm as pnm
(oz, oy, ox), (sz, sy, sx) = pnm.mpi_block(global_shape_zyx)  # this rank's ORB block, in VOXELS
local = sdf[oz:oz + sz, oy:oy + sy, ox:ox + sx]
pores, seg, conns = pnm.extract_pore_network_mpi(local, global_shape_zyx, origin_zyx, spacing_zyx)
# pores: the pores whose peak this rank owns; seg: this rank's block (int32, local.shape, global
# label ids); conns: global (M,2), identical on every rank. origin_zyx stays the GLOBAL grid's
# physical origin — mpi_block's offset_zyx is an integer voxel offset, a different thing.
# Rank and size come from mpi4py (MPI.COMM_WORLD.rank / .size); the module has no mpi_rank().
```

Validated by `tests/kokkos_mpi` (ctest, np = 1, 2, 4, OpenMP + CUDA): per-voxel segmentation ids,
the pore set, and the connection list all match the single-rank oracle exactly (pore centroid
positions to 1e-5·spacing on GPU — FMA contraction noise; radii and everything integer bitwise).

## License

MIT.
