peclet.coupling — CFD-DEM coupling, unresolved and resolved

pip install peclet[cfd-dem]        # the `peclet-coupling` distribution; an sdist, built against
                                   # your MPI and Kokkos prefix (see docs/DEPLOYMENT.md)

Two-way coupling of peclet.flow (Eulerian fluid) and peclet.dem (Lagrangian particles). Two drivers, one vocabulary:

• CfdDem — the unresolved point-particle driver for dilute-to-dense suspensions and packed beds: a grain is a point with a drag closure (Multiphysics Phase 6, ../docs/MULTIPHYSICS_PLAN.md).
• ResolvedCfdDem — the resolved driver: each grain is an analytic SDF instance in the flow solver's scene, the fluid resolves its surface, and there is no drag correlation anywhere (see Resolved coupling below).

Both live in the caller's own PHYSICAL units: the drivers take the cell size and the lower corner from the flow solver itself (flow.spacing, flow.origin), so a solver built with Solver(cells, extent=...) couples to a DEM stated in metres with no conversion anywhere, and a cell-unit solver (no extent) keeps spacing 1. Nothing in the coupling API is stated in cells.

Design (CfdDem)

Physics-free glue. The compute kernels (particle↔grid deposition, drag laws, momentum feedback) live in the _coupling nanobind extension and run in place on the arrays the two solvers already expose — the fluid grid fields zero-copy through flow.diagnostics.field_view(...), the particle drag round-tripped through the dem host API. There is no C++ link between flow and dem: the Python CfdDem driver (python/peclet_coupling/driver.py) composes them. This mirrors the suite's architecture (Python is the composition layer).

Per fluid step (CfdDem.step()):

1. Void fraction — scatter each particle's volume onto the grid (trilinear, wall-aware: near an immersed solid the weights re-normalise over the fluid corners so no hold-up leaks into walls), fold the ghost deposits (periodic wrap on periodic axes; same-side fold onto the boundary cell at a non-periodic domain face — a grain resting on the distributor scatters part of its volume below z=0, and that hold-up belongs to the bottom cell, not to a ghost the fluid never owns), and ε = clamp(1 − Vsolid/Vcell, eps_min, 1). The floor eps_min defaults to 0.25, a physical regularisation rather than a guard: real voidage bottoms out near random close packing (~0.36 monodisperse, ~0.25 for wide bidisperse mixes), and anything lower can only come from interpenetrated particles or deposit artefacts and must not reach the volume-averaged fluid (whose projection amplifies the interstitial velocity by 1/ε). The older 0.4 floor under-predicted dense-bed drag ~3x; the interim 0.05 guard let interpenetration artefacts detonate a bed. The same 0.25 is the kernel default (_coupling.compute_void_fraction); the fixed-bed tests pass eps_min=0.05 explicitly because their uniform lattice never clamps. A particle whose trilinear stencil falls outside the domain by more than one ghost layer (e.g. pushed through a DEM wall by a violent contact solve) is dropped from the exchange entirely — no deposit, zero drag — so a runaway escapee can never feed a diverging β·u_p source into the boundary row. The optional volume filter is a PHYSICAL width: smooth_length (MFIX's DES_DIFFUSE_WIDTH, Capecelatro & Desjardins' δ_f) is a length in the caller's units, set from the particle diameter — the point-particle approximation is what needs δ_f ≫ d_p, and it must not move when the mesh does. It is realised as n explicit diffusion sweeps with a per-axis coefficient α_a = C/h_a², C = 1/(2 Σ_a 1/h_a²), n = round(σ² Σ_a 1/h_a²), so the Gaussian is a ball in space and not in index — on a cubic mesh that reduces to α = 1/6, n = round(3 (σ/h)²), term for term and bit for bit the cubic formula it generalises. Gate: tests/test_smoothing_isotropy.py (three physical σ equal to 8.9e-10 on a (1, 2, 0.5) cell, the single-α ablation off by exactly the spacing ratios, the cubic mesh bitwise).
2. Drag + feedback — gather the fluid velocity and ε at each particle, evaluate the drag law ("stokes", "schiller_naumann", "ergun", "di_felice", "wen_yu", "gidaspow", "beetstra" — Beetstra–van der Hoef–Kuipers 2007, "tang" — Tang et al. 2015; one literature name each), write the drag force to the particles and deposit the reaction onto the fluid momentum source.
3. Advance — apply the drag to the particles and sub-step dem dem_substeps times (drag held constant), then advance the fluid one step (its RHS/operator now carry the feedback).

Implicit drag (the key stability piece)

An explicit reaction force −β(u−u_p) in the fluid RHS diverges for the stiff drag coefficient β of a dense bed (β·dt/ρ ≫ 1; local β reaches ~10³). So the default feedback is semi-implicit: the coupling deposits the linear-drag coefficient density onto flow's drag_beta field (added to the momentum diagonal by flow.enable_drag()) and the target β·u_p onto force_* (the RHS), so the fluid solve becomes (ρ/dt + β)u = … + β u_p — unconditionally stable for any β. The particle side stays explicit (fine for moving particles at moderate β). implicit_drag=False selects the explicit −F/Vcell feedback (dilute only).

Two fluid modes

• porous=True — volume-averaged (use this for beds). The fluid solves the full volume-averaged continuity ∂ε/∂t + ∇·(εu) = 0 (u = the interstitial gas velocity) with a SIMPLE-like eps- and drag-weighted pressure projection — scheme, defaults and validation in flow/doc/porous_drag_scheme.md. The pressure-force split is Model B: the gas carries the full −∇p, the particles get drag + gravity, and the literature (Model-A) drag closures are converted once inside the kernel, β_B = β_A/ε (model_b flag). Gas convection (implicit FOU + explicit deferred-correction TVD) is enabled by the driver by default (advection=True).
• porous=False — dilute simplification. The fluid stays incompressible (div u = 0); ε enters the drag correlation only. Cheap and validated dilute→moderate. Note this is not "Model B": Models A and B both use the full continuity and differ only in the −ε∇p vs −∇p split.

Other scope notes: deposition uses atomic_add ⇒ results are tolerance-, not bit-exact; the "ergun" drag kind is the superficial-velocity form built for the incompressible mode — for porous beds use "gidaspow" (its dense branch is the classic interstitial Ergun form).

Resolved coupling (ResolvedCfdDem)

ResolvedCfdDem(flow, dem, *, radius, mu, rho, fluid_dt, dem_substeps=20, periodic=(bx, by, bz), gravity=(0, 0, 0), rho_p=None, move_particles=True, buoyancy=True, apply_torque=False, force_method="reaction") is Layer 4 of ../docs/ANALYTIC_SDF_GEOMETRY.md: one kSphere scene node, one instance per grain, installed in the flow solver's scene (set_scene). Per coupling step it pushes dem's positions / quaternions / velocities / angular velocities into the scene as instance transforms + rigid-body motion, flow.rebuild_geometry() re-derives the SDF, cut-cell overlay, apertures and pressure operator, the fluid steps, and the hydrodynamic load comes back — by default the discrete reaction (the momentum the fluid actually lost to each grain, exactly conservative; flow's hydro_force_torque_reaction()), or the reconstructed traction integral with force_method="traction" (flow's diagnostics.hydro_force_torque_traction(), a diagnostic; it under-reads the drag by a resolution-independent ~30 %). Gravity/buoyancy is added (rho_p, buoyancy), the force (and, with apply_torque=True, the torque) is handed to dem, and dem sub-steps at the DEM timestep with the load held constant (weak, explicit coupling: the fluid load is lagged by one fluid step). Pure Python — no compiled kernels of its own, so it imports even without the _coupling extension.

The bridge is a pure identity in the caller's units (dem state in, force and torque out; no scale factor anywhere). The hydrodynamic torque is validated (a spinning sphere reproduces the Stokes torque 8πμa³Ω to ~2–3 %) but apply_torque is off by default: dem's default inverse inertia is not the grain's, so set the physical principal inertia ((2/5) m R² for a sphere) before enabling it. The per-grain results of the last step are the last_force / last_torque (N,3) arrays (traction mode also fills last_force_pressure / last_force_viscous). periodic is per-axis like everywhere else in peclet, but the flow scene's periodic images are all-or-nothing, so a mixed triple is refused. Gallery: peclet-examples/examples/rotating-sphere-torque.

Backends

CfdDem runs on whatever Kokkos backend peclet.flow was built for. On a CUDA/HIP build the coupling kernels run on-device, so the driver array-programs through CuPy and takes the grid fields (flow.diagnostics.field_view) and particle state (dem.get_*_view) zero-copy via DLPack; on a host build it uses NumPy over the same buffers. Detected automatically from peclet.flow.execution_space.

Validation (tests/)

Both cases pass identically on host-openmp and CUDA (RTX 5080):

• test_terminal_velocity.py — single settling sphere: the slip velocity matches Stokes to 0.1–0.2 % and Schiller–Naumann to 1.4–1.6 % (the lab-frame speed is ~2× the slip because the particle drags its own Stokeslet flow, so the physical comparison is the slip).
• test_fixed_bed_ergun.py — uniform fixed bed (one particle per cell, ε = 0.6): the measured (f_drive, U) pair lands on the Ergun curve to 0.0 % across the viscous, transition and inertial (Re_p ≈ 6) regimes — validating ε deposition, both Ergun terms, and the stable two-way feedback.
• test_fixed_bed_ergun_porous.py — the same bed on the volume-averaged (porous, Model B) path with the Gidaspow closure: (f_drive, U = ε·u_interstitial) lands on the Ergun curve to ~3 % across all three regimes with no fitted factors — validating the eps-weighted projection, the interstitial kinematics and the β_B = β_A/ε conversion together.
• test_mpi_fixed_bed_ergun.py — the fixed-bed Ergun benchmark run distributed (flow init_mpi, each rank couples its ORB block; particle deposits fold across ranks + periodically via the reverse/add-reduce halo diagnostics.exchange_field_add, deposit origin shifted by the block origin). The superficial velocity U (reduced over ranks) lands on the Ergun curve to 0.0 % and is bit-identical at np=1/2/4 — the distributed deposition + fold + solve reproduce the coupled physics exactly.
• P2G/G2P conservation + the gather/scatter adjoint identity: core test_particle_grid (host + CUDA).

Multi-rank coupling

CfdDem runs distributed when the flow solver is decomposed (flow.init_mpi(...), world size > 1): each rank couples its local block, the deposit grid map is shifted by the block origin (so particles in global coordinates land locally), and cross-rank + periodic ghost deposits (void fraction + drag reaction) fold onto their owner with the reverse halo (diagnostics.exchange_field_add) instead of the single-rank NumPy fold.

One partition from construction on. A coupled run owns ONE decomposition, chosen by the combined cost of the two codes: the weight field 1 + gamma * particle_count over the global grid (fluid work 1 per cell, gamma per particle; gamma = 1 until its calibration). A moving multi-rank CfdDem builds it at construction and moves BOTH codes onto its weighted ORB before the first step; CfdDem.rebalance(gamma) does the same again later, through the same code path. flow's diagnostics.rebalance_by_weights(w) builds it ALIGNED for its pressure multigrid — every split plane on a multiple of 2^a cells, the largest a within a 1.05 weight imbalance — and returns 2^a; CfdDem passes that to dem.migrate_to_weights(w, align=2^a) and keeps it for the moving step's per-step migration, so the two codes own identical blocks. Give the flow + dem the same global grid / domain (flow.init_mpi before its geometry, dem.init_mpi + enable_mpi_step) before constructing CfdDem; where each rank's particles start does not matter, the construction migrates them. The construction co-rebalance is what closes the start-up trap: init_mpi alone leaves flow on its default partition (splits snapped to powers of two) and dem on the equal-cell ORB, which differ on some grids — 48³ at np = 4 gives flow 32|16 and dem 24|24, where the deposit used to couple ~7.5k particles into cells their rank does not own, silently. A fixed bed (move_particles=False) is not co-rebalanced at construction: dem never migrates there, so its particles stay in flow's init_mpi block where the caller put them; rebalance() is still available when its dem is distributed.

Co-location is asserted, never assumed. After every co-rebalance (construction included), and once more at the first moving step as a safety net, every rank checks that each particle dem owns lies in flow's block (and that flow's block is aligned as dem was told); a mismatch raises RuntimeError on every rank. test_mpi_construction_partition.py (np = 1 first, then 2 and 4): the 48³ heap, no rebalance() call, co-location holds and the velocities reproduce np = 1. test_mpi_rebalance.py (np = 1 first, then 2 and 4): a 32³ heap, rebalance(4) → alignment 2, the mean particle and fluid velocities reproduce np = 1 (0 and 1.6e-16 relative), and moving dem back onto the equal-cell ORB makes the assertion fire. It needs flow 413e75a: a size-changing redistribute used to zero the porous eps^n, 5e-2 off after three steps.

Moving particles (move_particles=True): each fluid step CfdDem first migrates dem onto flow's grid partition (dem.migrate_to_weights) so every owned particle sits in its rank's block, then runs the DISTRIBUTED DEM substeps (dem.step_mpi, requires dem.init_mpi + dem.enable_mpi_step). A rank that momentarily owns no particles still runs the halo collectives (the per-particle kernels are skipped). Validated test_mpi_fixed_bed_ergun.py (static, bit-identical np 1/2/4) and test_mpi_moving_suspension.py (drifting cloud crossing rank boundaries: the distributed migrate + step + deposit-fold + gather reproduce single-rank to ~4e-8, np 1/2/4) and test_mpi_rebalance.py (above). Call flow's init_mpi BEFORE its geometry (set_solid / set_pressure_geometry); flow raises otherwise.

Polydisperse, multi-rank: a per-particle radius must follow its particle through every migration (construction, rebalance(), the per-step migrate_to_weights). dem carries each particle's size as its scale (world radius = scale × global_scale × base radius), so the driver keeps only the constant radius/scale and re-derives the owned radii from dem.get_scales() after each migration. The radius array must therefore be dem's own radii (set_scales after set_positions); one that is not proportional to the scales raises ValueError at the first co-rebalance. Until 2026-09-25 the driver kept the construction-time array, so every migration handed particles another particle's radius (or an array of the wrong length): test_mpi_polydisperse_moving.py (48³, radii 0.5/0.8) missed np = 1 by 5.5 (np 2) and 7.7 (np 4) relative; it now reproduces it to 5e-8.

np=4 of the moving test used to fail twice over (fixed 2026-09-24). Up to dem 8abcfd2 it deadlocked: the cloud's y = 16 row sits on the 2x2 block face, so the upper ranks owe ghosts to the lower ones and receive none, and dem skipped every exchange on a rank without ghosts (fixed in dem 6d30722). It then missed np=1 by 1.1e-2 because the test called flow's set_pressure_geometry BEFORE init_mpi: the pressure multigrid is built with the geometry, so every rank solved its own block-periodic pressure problem (np=2 hid it: the cloud is symmetric within each x-half). flow now raises on that order, the test calls init_mpi first, and np=4 reproduces np=1 to 2e-8.

Known limitation (not the coupling — every distributed coupling op is bit-identical to single-rank in isolation): a sustained dilute settling suspension in a triply-periodic box with no buoyancy is an ill-posed, numerically unstable configuration — at np>1 the flow solve's reduction-floor non-determinism seeds that instability. Well-posed cases (bounded / driven flow, denser beds) are unaffected.

Note: dem.get_velocities() (host copy getter) has a pre-existing failure after a periodic DEM step on CUDA (a Kokkos strided-subview-after-resize limitation, unrelated to the coupling); the driver uses the zero-copy device views throughout and exposes last_slip for inspection.

Build

cmake -S . -B build -DCMAKE_PREFIX_PATH="$PWD/../extern/install/host-openmp"
cmake --build build -j        # -> build/peclet/coupling/_coupling.*.so (+ the staged .py files)
# run the tests (all three build trees on PYTHONPATH); each test is a pytest function AND a script:
export PYTHONPATH="$PWD/build:$PWD/../flow/build:$PWD/../dem/build"
OMP_NUM_THREADS=4 OMP_PROC_BIND=false pytest tests -q -k "terminal or ergun or isotropy"
python tests/test_fixed_bed_ergun.py
# the test_mpi_*.py tests need an MPI build of flow + dem and mpi4py (they SKIP under pytest else):
mpirun -np 2 python tests/test_mpi_fixed_bed_ergun.py

CI (.github/workflows/ci.yml) builds Kokkos (OpenMP), flow, dem and coupling from source and runs the single-rank pytest battery on every push.

Follow-ups

Kernel-width (vs trilinear) deposition smoothing; the ρε volume-averaged inertia and ∇·[εμ(∇u+∇uᵀ)] viscous forms in the gas momentum (accuracy — see flow/doc/porous_drag_scheme.md §6); a PEA-style implicit particle-drag substep for very stiff moving beds (m_p/β < Δt — the fluid side is already implicit).