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 — the same source runs on CUDA, HIP, and OpenMP backends, selected at build time by the install prefix. Part of the peclet suite; split out of peclet-flow (its former peclet.flow.pnm module — the repo's original "pnm_from_sdf" feature).

Install / build

# 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

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 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:

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.

# 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.