peclet.flow — Eulerian Navier–Stokes solver¶

The incompressible cut-cell IBM Navier–Stokes solver on a staggered MAC grid (staggered Solver, collocated SolverColocated), with geometric VoF two-phase flow, analytic moving geometry (scenes) and the distributed (MPI) solve. execution_space reports the compiled-in Kokkos backend; has_mpi whether this build carries the multi-rank API. Regenerate from an MPI-enabled build.

Note

Auto-generated from the installed module docstrings. Drive simulations from Python; the full C++ API is on each repo's Doxygen site.

`peclet.flow`¶

peclet.flow — the Eulerian incompressible Navier–Stokes solver.

A Kokkos cut-cell Immersed-Boundary-Method solver on a staggered MAC grid (grid-agnostic by design: Cartesian cut-cell today, able to consume an unstructured Voronoi grid from :mod:peclet.voro). The compiled backend (Serial / OpenMP / CUDA / HIP) is chosen at build time — peclet.flow.execution_space reports which one this build has.

:class:peclet.flow.Solver — the staggered MAC solver.
:class:peclet.flow.SolverColocated — the collocated/cell-centered variant.
:func:peclet.flow.march_to_steady — march a solver to its certified steady state, optionally Anderson-accelerated (returns a :class:peclet.flow.MarchResult).

Pore-network extraction lives in the companion :mod:peclet.pnm package (peclet-pnm; it was peclet.flow.pnm before 2026-07).

peclet is an implicit (PEP 420) namespace shared with the other peclet-* packages, so it has no top-level __init__.py.

`Solver`¶

One incompressible Navier-Stokes solver on a block of the grid (staggered MAC, the reference) -- or, for SolverColocated, the same API with cell-centered velocities and the ABC approximate projection. Construct on a (physical) domain; set_rho / set_mu / set_dt and the domain BCs; give it a geometry (set_solid, set_pressure_geometry or a scene); step(); read out with the get_* arrays. Developer instruments and tuning live on .diagnostics.

Method / property	Description
`add_field`	add_field(self, name: str) -> None Register a new zero-initialised cell-centred field on the grid (for transported scalars or material properties). Idempotent.
`add_scalar`	add_scalar(self, name: str, diffusivity: float = 0.0, scheme: str = 'koren', iters: int = 50) -> None Register a transported scalar (temperature/concentration/...): constant diffusivity (grid units), advection scheme 'fou' (first-order upwind), 'koren' (Koren TVD, the default) or 'sou' (second-order upwind), and RB-GS diffusion sweeps. The scalar is a registered field (get_field/set_field/diagnostics.field_view). Requires geometry (set_solid/set_pressure_geometry) for the openness-weighted operators.
`advance_scalars`	advance_scalars(self) -> None Advance all registered scalars one dt with the current velocity (also done by step()).
`advect_vof`	advect_vof(self, dt: float) -> None KINEMATIC colour advection: advance C ONCE with the solver's CURRENT face velocity over `dt` (seconds), with NO Navier-Stokes step. This is the advection-benchmark entry point (Zalesak, LeVeque, and the cut-cell conservation gates): a frozen projected velocity advecting a colour field is a pure statement about the transport scheme, with the momentum and pressure solves out of the picture.\n\nIt THROWS if the current velocity is not discretely divergence-free to 1e-10, measured with max_open_divergence_projected() (the NON-mutating diagnostic). That is not a nicety: Weymouth-Yue's exact conservation is CONDITIONAL on sum_f o_f u_f = 0 per cell, because the dilation term adds H(C-1/2) times that residual to every full cell's volume budget. Run step() to a steady state (or project()) and advect with the solver's own output; never with an analytic sample, which is solenoidal only to O(h^2). It ALSO throws when no cut-cell pressure operator exists (WO-R2): without one the divergence guard measures nothing at all — max_open_divergence() returns 0.0 — and a cell-centre-sampled LeVeque field (true max
`advect_vof_blocks`	advect_vof_blocks(self, dt: float) -> None Advance every block by dt with the CURRENT (projected) face velocity and union the result into 'C' — the block twin of advect_vof(dt), and it carries the same precondition: it RAISES unless the face field is discretely divergence-free to 1e-10, because Weymouth-Yue's exact conservation is conditional on it. Per step: gather the face velocity from the ranks that own each block's cells to its master (plain Isend/Irecv — the block table and the decomposition are both replicated, so every rank knows every message size and an NBX handshake would only rediscover it), run the three sweeps on the dense block, re-centre the box, replicate the table, and scatter the inner colour back with UNPACK_MAX. NOTE on the union: max() from an empty union CLIPS the negative Weymouth-Yue round-off residue to an exact 0 — measured up to 6.2e-17 on the LeVeque gate, and every measured union/global-field difference was exactly that and nothing else.
`apply_phase_change`	apply_phase_change(self, dt: float) -> None KINEMATIC phase-change step: build mdot / the PLIC areas / the normals, deposit the divergence source (census only — nothing is projected) and apply the interface regression over `dt`. No Navier-Stokes step, no advection. This is the P0a / P1 driver: an interface regressing under a prescribed or thermally-computed flux, with the momentum and pressure solves out of the picture.
`balanced_force_projection`	Whether the balanced-force projection runs (see set_balanced_force_projection).
`block_origin`	block_origin(self) -> list[int] This rank's inner-block origin in GLOBAL cells ([0,0,0] single-rank). Shift the coupling deposit origin by this so particles in global coordinates land in the local block.
`capillary_dt`	capillary_dt(self) -> float The Brackbill (1992) capillary time-step limit sqrt((rho_1+rho_2) h^3 / (4 pi sigma)), +inf when surface tension is off. Denner & van Wachem (2015) verified this IS the stability boundary of an explicit CSF — the prefactor 1/(4 pi), the h^{3/2} scaling and the SUM of the phase densities (both phases oscillate). The density sum is the declared phase pair under enable_vof_momentum, otherwise min(rho)+max(rho) over the field.
`cell_centers`	cell_centers(self) -> tuple The physical cell-center coordinates of THIS rank's inner block as three 1-D float64 arrays (x, y, z), i.e. origin + (i + 1/2) * spacing per axis. np.meshgrid(s.cell_centers(), indexing='ij') is the grid an SDF for set_solid is sampled on. Under MPI these are the LOCAL block's centers in GLOBAL coordinates. CELL-UNIT TRAP: without an extent this is i + 1/2, while set_scene puts cell i's centre at i -- the same body given as a scene and as an SDF sampled here lands half a cell apart. With an extent both use origin + (i + 1/2) spacing. Unifying the cell-unit convention changes results; deferred to 2.0.0.
`cells`	This rank's inner block resolution (nx, ny, nz) -- the shape of every field array (get_u, get_p, the SDF for set_solid). Under MPI it is the LOCAL block from mpi_block().
`clear_divergence_sink`	clear_divergence_sink(self) -> None Drop the auto-balanced sink region.
`clear_divergence_source`	clear_divergence_source(self) -> None Drop the prescribed divergence source (the field stays registered).
`compute_vof_curvature`	compute_vof_curvature(self) -> dict VoF rung V3: compute the interface curvature from the CURRENT colour field and store it in the registered cell fields 'kappa' and 'kappa_branch'. Returns THIS RANK's branch census as a dict. 'kappa' is the RAW registered field, kappa = 2H in 1/h (cell units); vof_curvature() returns the same thing in 1/length. POSITIVE for a convex blob of liquid: a sphere of liquid of radius R cells reads +2/R, an infinite cylinder +1/R, a plane 0. The cascade is Popinet's (2009): a standard height function on 7-cell column sums of C over the 3x3 transverse patch in the direction of the largest
`contact_angle`	contact_angle(self) -> float The prescribed static contact angle in degrees (90 if none was set).
`decomposition_levels`	This solver's decomposition mode (0 = aligned ORB, >= 2 = coarse-first with that depth).
`decomposition_max_imbalance`	The load-imbalance budget the coarse-first search stays within.
`diagnostics`	The developer tier of this solver (instruments, tuning, ablations, internals).
`enable_cell_force`	enable_cell_force(self) -> None Allocate + register the per-cell body-force fields force_x/force_y/force_z and route them into the momentum RHS, for an external writer (e.g. CFD-DEM drag feedback) to fill directly via diagnostics.field_view('force_z'). They persist across steps until overwritten.
`enable_drag`	enable_drag(self) -> None Enable implicit (semi-implicit) linear drag for CFD-DEM: allocate the per-cell 'drag_beta' field (added to the momentum diagonal so a -beta(u-u_p) source is treated implicitly -> unconditionally stable for the stiff beta of a dense bed) plus force_x/y/z (which carry betau_p, the RHS target). Fill 'drag_beta' and 'force_*' via diagnostics.field_view each step.
`enable_phase_change`	enable_phase_change(self, rho_gas: float, rho_liquid: float, h_lv: float = 1.0) -> None Enable VoF phase change (VOF_PLAN section 9, rungs P0/P1; the Boyd & Ling 2023 / Malan 2021 pattern). Registers 'mdot' (the interfacial mass flux, kg m^-2 s^-1 in solver units, POSITIVE = evaporation) and 'pc_source' (the deposited divergence source, 1/s). What it does per step, at the HEAD of step(): (1) reconstruct the PLIC plane of every interfacial cell and take its POLYGON AREA A analytically (A =
`enable_vof`	enable_vof(self) -> None Enable the geometric (PLIC + Weymouth-Yue) VoF colour field 'C'. Registers 'C' as an ordinary cell field (get_field/set_field/diagnostics.field_view/closure input) and allocates the colour field's OWN g=3 working block + halo; the solver's G=2 is untouched. C is advected at the end of every step() with the just-projected face velocities. RUNG V2a SCOPE — read before using: there is NO surface tension (rung V4) and NO momentum-consistent transport (rung V2b). Mass and momentum are advected by different fluxes, so a mixed cell carries a spurious interfacial momentum source of order delta-rho; the literature (Rudman 1998, Arrufat 2021) is unambiguous that this breaks down around density ratio 1000 unless the resolution is absurd. V2a is valid only at MODEST density ratios for cases with motion. A high-ratio case AT REST (the hydrostatic acid test) is exact, because there is no momentum to mis-advect. COLLOCATED (SolverColocated) is supported since rung V8 (2026-09-02), ALL-FLUID only: the colour is advected by the PROJECTED face field uf_/vf_/wf_ (the field the ABC approximate projection makes exactly divergence-free, which is what Weymouth-Yue's conservation proof needs), and the CSF enters the implicit momentum predictor as the mass-adjoint face integral (doc/collocated_varrho_forces.md). An immersed solid on that grid THROWS (the cut-cell colour transport of rung V5a is staggered-only), and so does enable_vof_momentum. On the collocated grid the SPURIOUS-CURRENT number to read is the FACE field (get_uf/get_vf/get_wf): the CELL field additionally carries the approximate projection's invisible odd-even mode, which centerToFace annihilates and the projection therefore cannot remove. An IMMERSED SOLID is supported since rung V5a: the geometric fluxes are openness-weighted and the update is done in fluid-volume units, so sum eps_eff*C is conserved exactly against the projection's own openness-weighted divergence; it needs set_solid(..., cutcell_pressure=True) (the staircase operator has no face openness to weight with) and it approximates the solid-clipped flux polygon by (whole-cell PLIC slab) x (open area) — see diagnostics.vof_diagnostics()['clipped_volume'].
`enable_vof_block_csf`	enable_vof_block_csf(self) -> None Form the surface-tension force PER BLOCK (rung W2): each marker runs its own curvature cascade on its own dense box and forms the V4 balanced-force face force there (sigma kappa_f (C(i) - C(i-s))/h with the same selective face-curvature rule addCsfRhs uses), and the three face fields are scattered into the global RHS with UNPACK_SUM — TBFsolver's VOF.f90::computeSurfaceTension structure on the suite's kernels. WHY THE FORCE AND NOT THE CURVATURE IS SCATTERED: kappa is not additive and the union colour is a max, so a face between two OVERLAPPING markers has no single (kappa, dC) pair to build a force from. The force is the additive quantity, and forming it where each marker's own colour still exists is the only place the balanced-force pairing (the same face difference the projection's gradient uses) is available per marker. Requires set_surface_tension(sigma) and the block container; STAGGERED only. Once on, step() drives the whole two-phase stage through the blocks: the union C feeds the closures exactly as before, and the colour is advected by the blocks in the advect_vof slot.
`enable_vof_blocks`	enable_vof_blocks(self, seeds: collections.abc.Sequence[collections.abc.Sequence[float]]
`enable_vof_blocks_from_colors`	enable_vof_blocks_from_colors(self, boxes: collections.abc.Sequence[collections.abc.Sequence[int]], colors: list) -> None Restart the block container from a vof_block_color() checkpoint: one marker per (lo_x, lo_y, lo_z, hi_x, hi_y, hi_z) INNER box (exactly the 'lo'/'hi' of diagnostics.vof_block_stats(), NOT grown by the margin again) with its colour written directly. Exact whatever the markers are doing -- nothing is gathered out of the union field.
`enable_vof_blocks_from_field`	enable_vof_blocks_from_field(self, boxes: collections.abc.Sequence[collections.abc.Sequence[int]]) -> None Seed the block container from the CURRENT colour field: one marker per given global index box (lo_x, lo_y, lo_z, hi_x, hi_y, hi_z; half-open), its colour copied out of the field set_vof installed. The general seeding path — enable_vof_blocks(spheres) is a convenience over it — and the way a marker of any shape enters (a quasi-2-D Hysing cylinder, a scanned bubble, an arbitrary blob found with a host connected-component labelling). The boxes are the BUBBLE extents; the container grows each by its own 3-cell margin. They must not share a cell: the seeding gather is a copy, not a union, so a cell inside two boxes would be handed to both markers.
`enable_vof_momentum`	enable_vof_momentum(self, rho_gas: float, rho_liquid: float) -> None Advect rhou on the HALF-SHIFTED momentum control volumes with the SAME geometric fluxes, the same sweep order and one frozen dilation flag as the colour field, then recover u = (rho^c u)/rho^c. This is what makes the scheme usable above density ratio ~100 (Rudman 1998; Arrufat 2021: a raindrop accurate within 15% at 15 cells/diameter instead of ~200). The two phase densities are REQUIRED, and must be the ones the rho closure produces at C=0 and C=1 — the momentum flux is rho_gas(a-F) + rho_liquidF with F the geometric liquid flux, so the scheme has to know which density each phase carries. Turning this on MOVES the VoF advection from the end of step() to its head (the momentum advection must precede the predictor that consumes it), so the advecting field is u^n — the previous step's projected output. u at step 0 must therefore be discretely divergence-free (rest or a uniform field both are). Requires the variable-density path, staggered layout, explicit advection, no immersed solid and no porous continuity, and is refused together with the per-bubble block container (enable_vof_blocks, in either order): the momentum rides the structured colour's fluxes, not the markers'.
`extent`	The GLOBAL domain side lengths. Equal to the cell counts without a physical domain.
`field_names`	field_names(self) -> list[str] Names of all registered fields (velocity u/v/w, p, sdf, plus any added), sorted.
`get_field`	get_field(self, name: str) -> numpy.ndarray[dtype=float64] Return a registered field's inner region as a Fortran-order (nx,ny,nz) float64 array.
`get_ox`	get_ox(self) -> numpy.ndarray[dtype=float64] The -x face openness field of the cut-cell operator: the fluid area fraction (0..1, dimensionless) of the face between cells (i-1,j,k) and (i,j,k), as a Fortran-order (nx,ny,nz) float64 array over THIS rank's inner cells. Computed by set_solid from the SDF (1 everywhere without a solid); it is what the momentum operator and the geometric cut-cell projection see. get_ox_proj is the openness the PROJECTION actually conserves fluxes through, which differs under diagnostics.set_ghost_projection (binary/coupled faces).
`get_ox_proj`	get_ox_proj(self) -> numpy.ndarray[dtype=float64] -x face openness whose fluxes the projection CONSERVES (binary/COUPLED under diagnostics.set_ghost_projection, geometric cut-cell otherwise). Use for flux bookkeeping (peclet.pnm extract_network_flow).
`get_oy`	get_oy(self) -> numpy.ndarray[dtype=float64] The geometric -y face openness per inner cell, (nx,ny,nz); see get_ox.
`get_oy_proj`	get_oy_proj(self) -> numpy.ndarray[dtype=float64] -y face openness the projection conserves (see get_ox_proj).
`get_oz`	get_oz(self) -> numpy.ndarray[dtype=float64] The geometric -z face openness per inner cell, (nx,ny,nz); see get_ox.
`get_oz_proj`	get_oz_proj(self) -> numpy.ndarray[dtype=float64] -z face openness the projection conserves (see get_ox_proj).
`get_p`	get_p(self) -> numpy.ndarray[dtype=float64] Return the physical pressure as a Fortran-order (nx,ny,nz) float64 array (index [x,y,z]).
`get_u`	get_u(self) -> numpy.ndarray[dtype=float64] Return the x-velocity component as a Fortran-order (nx,ny,nz) float64 array (index [x,y,z]).
`get_uf`	get_uf(self) -> numpy.ndarray[dtype=float64] Return the divergence-free FACE x-velocity (collocated: projected MAC field; staggered: == get_u).
`get_v`	get_v(self) -> numpy.ndarray[dtype=float64] Return the y-velocity component as a Fortran-order (nx,ny,nz) float64 array (index [x,y,z]).
`get_vf`	get_vf(self) -> numpy.ndarray[dtype=float64] Return the divergence-free FACE y-velocity (collocated: projected MAC field; staggered: == get_v).
`get_vof`	get_vof(self) -> numpy.ndarray[dtype=float64] The colour field's inner region as a Fortran-order (nx,ny,nz) float64 array (== get_field('C')).
`get_w`	get_w(self) -> numpy.ndarray[dtype=float64] Return the z-velocity component as a Fortran-order (nx,ny,nz) float64 array (index [x,y,z]).
`get_wf`	get_wf(self) -> numpy.ndarray[dtype=float64] Return the divergence-free FACE z-velocity (collocated: projected MAC field; staggered: == get_w).
`ghost_width`	ghost_width(self) -> int Ghost-layer width g of the velocity block (diagnostics.field_view returns an (n+2g) buffer).
`global_cells`	The GLOBAL grid resolution (gnx, gny, gnz); == `cells` single-rank, the numbers init_mpi gets under MPI (the CFD-DEM co-decomposition weight field is this shape).
`has_cutcell_pressure`	has_cutcell_pressure(self) -> bool True once the cut-cell pressure operator exists (set_solid or set_pressure_geometry was called). The porous continuity requires it; the coupling driver auto-sets an all-fluid geometry when absent.
`has_field`	has_field(self, name: str) -> bool Whether a field of this name is registered.
`has_moving_instance`	has_moving_instance(self) -> bool True when at least one scene instance carries a nonzero velocity.
`has_scalar`	has_scalar(self, name: str) -> bool Whether a transported scalar of this name is registered.
`has_scene`	has_scene(self) -> bool True once set_scene has installed an analytic CSG scene on this solver. False for geometry set any other way (set_solid with a sampled SDF, set_pressure_geometry), and false again if no scene has been given yet. set_solid_from_scene and the scene-derived per-body hydrodynamic force/motion bookkeeping all require this to be True first.
`hydro_force_torque`	hydro_force_torque(self) -> numpy.ndarray[dtype=float64] DEPRECATED (warns) -- the canonical names are hydro_force_torque_reaction() for the force and torque on the bodies, and diagnostics.hydro_force_torque_traction() for what this returns. It is NOT the force on the bodies: it is the reconstructed surface-traction integral, which under-reads the drag by ~30 % at any resolution (the viscous part; traction / reaction 0.685-0.730 over phi 0.008-0.45, N 24-128). Returns exactly what diagnostics.hydro_force_torque_traction() returns (shape (4, n_instances, 3): force, torque, pressure part, viscous part); see that method for the details. Where the reaction refuses (SolverColocated, porous, variable properties, domain BCs, implicit advection) the traction is the only estimate there is -- call it by its diagnostic name and allow for the under-read. Kept working under the suite's rule for a shipped name (docs/NAMING.md section 0); removal needs a major release.
`hydro_force_torque_reaction`	hydro_force_torque_reaction(self) -> numpy.ndarray[dtype=float64] THE hydrodynamic force and torque on each scene instance: the DISCRETE REACTION -- the momentum the fluid actually lost to each body, assembled from the composed step's budget (time term, body force, and the viscous fluxes at u; the pressure telescopes inside each owner region and needs no field). Shape (2, n_instances, 3): [0] force, [1] torque about the instance centre. The force is exactly conservative (sum = f V_fluid at steady state, to the solver residual) and as accurate as the flow solution it sustains: Zick & Homsy's drag within 0.7 % at N = 24 (tests/python/test_hydro_force_units.py). The torque carries the closed-form transposed-stress wall term; measured 2026-09-30 on a steadily spinning sphere it reads 1.019-1.027 of 8 pi mu a^3 Omega / (1 - phi) at R/h 6-12 (over-reads ~2 %, falling with resolution), and on a translating sphere whose true torque is zero it is zero to 1e-15 on a symmetric placement and
`init_mpi`	init_mpi(self, gnx: int, gny: int, gnz: int) -> None Wire the multi-rank step: pass the GLOBAL grid dims (gnx,gny,gnz). The Solver must have been constructed with this rank's LOCAL block dims (from mpi_block). MPI_Init is called if needed. Call it BEFORE the geometry (set_solid / set_pressure_geometry / set_solid_from_scene): the geometry builds the pressure and velocity multigrids for the partition in force when it runs, so init_mpi raises if a geometry already exists (built first, every rank would solve its own block-periodic pressure problem).
`instance_center`	instance_center(self, i: int) -> list[float] The resolved centre of rotation of instance i (world coordinates).
`max_open_divergence`	max_open_divergence(self) -> float Return the max cut-cell velocity-flux divergence max
`max_open_divergence_projected`	max_open_divergence_projected(self) -> float max
`max_porous_residual`	max_porous_residual(self) -> float Residual of the volume-averaged continuity max
`num_scene_instances`	num_scene_instances(self) -> int Number of instances in the scene set by set_scene (0 before it). Instance indices run 0..num_scene_instances()-1 in set_instance_transform / set_instance_motion, scene_owner() and the hydrodynamic force/torque accessors.
`origin`	The physical lower corner of the GLOBAL domain (0,0,0 by default).
`physical_units`	True when the solver was given an extent, i.e. lengths are physical rather than cells.
`pressure_bottom_extent`	The agglomerated-bottom extent threshold (see set_pressure_bottom_extent).
`pressure_telescope`	pressure_telescope(self) -> bool Whether the pressure multigrid's coarse levels telescope (merge ORB siblings onto fewer ranks once a level's block turns odd) rather than stop coarsening in place. Reads back set_pressure_telescope's current setting; ON by default since 2026-09-02 (the FoxBerry ladder), so this is normally True even though it was never called explicitly -- check it rather than assuming the legacy in-place-only hierarchy.
`rank`	rank(self) -> int This rank's index in MPI_COMM_WORLD.
`rebuild_geometry`	rebuild_geometry(self) -> None Re-derive all geometry from the current instance transforms. Velocity and pressure are PRESERVED across the rebuild (set_solid zeroes them by design). Full rebuild: the measured cost is ~65% momentum/IBM stencils, ~35% pressure/MG, scene sampling in the noise. Cells uncovered by the motion inherit zero, not an extrapolated fluid value.
`refresh_wall_velocity`	refresh_wall_velocity(self) -> None Re-derive ONLY the wall-velocity fields and the momentum operator that folds them in, for a driver that changes an instance's VELOCITY every step while its transform stays put -- the linearised oscillating body, a shear cell driven by counter-moving plates. set_instance_motion alone does not reach those fields (they are built in set_solid_from_scene), so the alternative is a full rebuild_geometry() per step to update a boundary condition. Refuses if a transform has changed since the last geometry build, because it does NOT re-sample the SDF, the apertures, the ownership field or the pressure operator. Velocity and pressure are untouched.
`reset_vof_bc_volumes`	reset_vof_bc_volumes(self) -> None Zero the per-face boundary liquid ledger (both the per-step and the running totals).
`set_advection`	set_advection(self, on: bool) -> None Enable/disable explicit high-order momentum advection (default scheme SOU). Off ⇒ Stokes.
`set_advection_scheme`	set_advection_scheme(self, scheme: str) -> None High-order momentum advection scheme: 'sou' (second-order upwind, the default) or 'koren' (Koren TVD).
`set_backflow_stabilization`	set_backflow_stabilization(self, beta: float) -> None Outflow backflow-stabilization coefficient (Bazilevs 2009 / Esmaily-Moghadam 2011): beta in [0,1] scales the dissipative outflow term that prevents backflow divergence when flow reverses at the outlet (e.g. a separated wake / BFS recirculation). Default 0.2; 0 = off. Inert where the outlet is purely outgoing.
`set_balanced_force_projection`	set_balanced_force_projection(self, enabled: bool) -> None Balanced-force projection (flow doc/collocated_varrho_forces.md §4.6), on both grids. Once per step, before the momentum predictor, the gradient part of the forces (set_body_force, the per-cell force fields, the CSF) is moved into the pressure with one extra solve of the projection's own operator, so static balances -- a hydrostatic column, a constant-curvature drop -- are exact from the first step at every viscosity, dt, density ratio and immersed-solid aperture. The solve does not depend on the flow state, so stability and the converged steady state are identical on and off; off leaves a decaying transient balance residue. get_p() stays the total physical pressure; its balanced part is the registered field 'p_balanced', which is RESTART STATE: registered as soon as the option is active, restoring it with 'p' (set_field) keeps the split, and omitting it (or switching the option on after steps with it off) re-splits P -- never a second copy of the balanced pressure. Cost: one extra pressure solve per step while the forces change, solved for the increment of the previous split and skipped when that is already converged -- a static interface costs no iterations. Default: ON on SolverColocated with variable density or surface tension, OFF otherwise; with an inflow/outflow face that default resolves OFF (a one-time stderr notice says so); an explicit call always wins. Refused (named error) with porous continuity, the ghost projection, set_fluid_only_constraint(2), the block CSF, inflow/outflow faces, the non-incremental pressure, and on the constant-density SolverColocated path (nothing to balance there).
`set_body_force`	set_body_force(self, force: collections.abc.Sequence[float]) -> None Set the body force per unit volume (fx, fy, fz) as one 3-sequence — e.g. a mean pressure gradient.
`set_capillary_cfl`	set_capillary_cfl(self, factor: float) -> None Safety factor on the capillary limit: step() raises when dt > factor * capillary_dt(). Default 1.0 — the Brackbill formula was measured to be the boundary itself, so it carries no built-in margin. Set it huge to disable the check, the same escape hatch set_vof_cfl_limit is for the Weymouth-Yue cap.
`set_collocated_scheme`	set_collocated_scheme(self, name: str) -> None Collocated cut-cell projection scheme (no effect on the staggered solver). DEFAULT since 2026-08-25: AUTO = 'ghost' where supported, falling back to 'gauge-exact' with a stderr notice on porous/variable-rho/domain-BC/Chebyshev configurations (ghost v1 limits); any explicit selection here disables AUTO. 'ghost' -- the fluid-only constraint scheme (binary-openness divergence + directional closures + the same gauge-exact gradient). Family-free, unconditionally stable (dt 60..1e20, no stabilizer), protocol-independent (C2), both-bed clean ladder -1.4% (R=8) -> +0.22% asymptote, Z&H anchor -0.018% at N=128. Costs: nonsymmetric BiCGStab pressure solve (~2.3-2.7x), ~1.6 KB/cell overlay (single-GPU size cap), fragmentation guard. MPI-capable (np=1,2,4 ctests); at-scale np>=16 hardening in progress. Equivalent to diagnostics.set_ghost_projection(True, 2, 2); the (1, 2) mixed mode stays quarantined there. 'gauge-exact' (the AUTO fallback) -- the aperture constraint with the directional gauge-exact pressure gradient. Cheapest scheme measured (symmetric MG-PCG, no fragmentation guard). CAVEATS from the 2026-08 campaign (doc/collocated_invisible_subspace.md): possesses an attractor FAMILY of steady states (support-inconsistent gradient/constraint pair) selected by march protocol, and a rotational-update instability whose wall-blend stopgap (diagnostics.set_rotational_wall_weight) has a resolution-dependent margin -- UNSTABLE at (R>=16 cells/radius, dt>=600); use fixed dt<=60 protocols at high resolution. Clean-protocol bias vs the staggered reference: -2.5% (R=8) -> ~+0.2% asymptote. 'embed' -- the complete Basilisk embed.h port (doc/history/collocated_embed_port_plan.md): the finite-volume momentum operator with the true-normal dirichlet_gradient wall drag applied as a defect correction on the IBM matrix, the openness-weighted -grad(P) predictor and cell correction, the wall-aware face map and the solid-cut-cell sliver mask. The live candidate for removing the collocated accuracy ceiling; its two intermediate rungs are diagnostics.set_face_interp(5
`set_contact_angle`	set_contact_angle(self, theta: float) -> None Rung V5b (WO-S). Prescribe a STATIC contact angle, in DEGREES, measured THROUGH THE LIQUID, on every immersed SDF wall. Needs set_solid(..., cutcell_pressure=True) + enable_vof; with no call the wall stays at the neutral 90-degree fill of rung V5a and every V5a number is byte-identical. What it changes is PASS 1 of the solid-band fill and nothing else: the band cells receive the volume fractions of the plane that continues the fluid-side interface into the solid at angle theta (m . n_w = cos theta, with n_w = grad(sdf)/
`set_contact_angle_dynamic`	set_contact_angle_dynamic(self, enabled: bool, theta_e: float = 90.0, slip_length_cells: float = 0.0, mu_liquid: float = 0.0, sigma: float = 0.0) -> None Rung V6 (WO-V6). The DYNAMIC contact angle -- `enabled=False` turns it (and the hysteresis) off so the static angle stands again; `enabled=True` imposes at the grid scale the Cox-Voinov apparent angle of a contact line moving at speed U_cl, with the outer cut-off at the CELL SIZE and an EXPLICIT slip length lambda (Afkhami, Zaleski & Bussmann, JCP 228:5370, 2009): theta_Delta^3 = theta_e^3 + 9 Ca_cl ln(Delta/lambda), Ca_cl = mu_l U_cl / sigma angles in radians internally, arguments in DEGREES; Ca_cl > 0 advancing, < 0 receding; theta_Delta is clamped into [1, 179] degrees (diagnostics.set_contact_angle_clamp). `slip_length_cells` is lambda/Delta and must lie in (0, 1). It is not a tuning knob you may omit from a report: a VoF contact line's NUMERICAL slip is proportional to the cell size, so without an explicit lambda the imposed angle is silently grid-dependent (VOF_PLAN section 6). ALWAYS state lambda alongside a dynamic-wetting result. `mu_liquid` is the liquid dynamic viscosity entering Ca_cl (solver units); `sigma` defaults to whatever set_surface_tension holds. theta_e also becomes the static base angle, so this call replaces set_contact_angle. U_cl is measured as -u.t_hat at the anchor fluid cell of each band cell, with t_hat the IN-WALL direction of the fluid-side PLIC normal (m points into the gas, so the liquid advances along -t_hat), then smoothed with a 3-point mean along t_hat. Needs set_solid(..., cutcell_pressure=True) + enable_vof; with no call every V5b number is byte-identical.
`set_contact_angle_field`	set_contact_angle_field(self, theta: ndarray[dtype=float64, order='F']) -> None Per-cell static contact angle in DEGREES, as a Fortran-order (nx,ny,nz) array. Only the value AT the solid band cell being filled is read, so cells away from a wall are irrelevant. theta is a field so the dynamic-angle rung (V6) changes only what fills it.
`set_contact_angle_hysteresis`	set_contact_angle_hysteresis(self, theta_a: float, theta_r: float) -> None Rung V6 (WO-V6). Advancing / receding contact-angle hysteresis, in DEGREES. While the measured apparent angle lies in [theta_r, theta_a] the contact line is PINNED and the fill imposes the APPARENT angle itself — which reproduces the current interface exactly, because the V5b fill is idempotent (WO-S finding 1), so no Young force appears and the line does not move. Above theta_a the advancing angle is imposed and below theta_r the receding one, each with the Cox-Voinov correction of set_contact_angle_dynamic when that is also configured. Sets the static base to (theta_a+theta_r)/2 if no angle was set yet.
`set_decomposition`	set_decomposition(self, levels: int, max_imbalance: float = 1.05) -> None Choose how this solver's shared MPI decomposition is built. levels=0 (DEFAULT) = the aligned ORB: split positions on the FINE grid are snapped to a power of two (capped at 16, i.e. 5 nested multigrid levels). levels >= 2 = COARSE-FIRST: build the ORB on the grid coarsened levels-1 times and refine the partition upward, so every block is a multiple of the coarsening factor BY CONSTRUCTION and the hierarchy nests for the full requested depth. Each candidate depth is built and MEASURED, and the deepest one whose max-block/min-block ratio stays within max_imbalance (default 1.05) is taken, else the aligned ORB. Call BEFORE init_mpi() (a later call raises), and pass the SAME levels / max_imbalance to flow.mpi_block() -- both derive the same partition and must agree.
`set_divergence_sink`	set_divergence_sink(self, weights: ndarray[dtype=float64, order='F']) -> None WO-P23: an AUTO-BALANCED sink region for the phase-change divergence source. Give a non-negative weight per cell (Fortran-order (nx,ny,nz)); after each deposit the solver subtracts (global sum of the phase-change source) * w / (global sum of w), so a CLOSED domain's Poisson RHS is exactly compatible every step with no user bookkeeping. Put the weights in the LIQUID far from the interface, where the exact solution simply has the liquid leaving. This is what makes a closed-box phase-change run possible without the variable-density outflow operator (whose density-ratio inconsistency is WO-R2's subject). Registered as the field 'div_sink'.
`set_divergence_source`	set_divergence_source(self, array: ndarray[dtype=float64, order='F']) -> None A PRESCRIBED extra divergence source (1/s) on the inner cells, added to the Poisson RHS beside the phase-change deposit, so the projection solves div(open u) = S. Registered as the field 'div_source'. Its use is to make a CLOSED (periodic) domain compatible with a net vapour production by carrying a balancing sink; with an outflow face the outflow does that job and this is unnecessary.
`set_domain_bc`	set_domain_bc(self, face: str, type: str, velocity: collections.abc.Sequence[float] = [0.0, 0.0, 0.0]) -> None Set a per-face domain BC. `face` is one of '-x', '+x', '-y', '+y', '-z', '+z'; `type` is 'periodic' (the default on every face), 'wall' (no-slip), 'inflow' (Dirichlet `velocity` -- an inlet, or a lid with a tangential velocity), 'outflow' (zero-gradient velocity, p = 0) or 'slip' (free-slip / symmetry plane: zero normal velocity, zero normal derivative of the tangential components, pressure Neumann like a wall; `velocity` is ignored). The TYPE must be set BEFORE the geometry (set_solid / set_pressure_geometry / set_solid_from_scene): it is folded into the operators when the geometry is built, and a call that would CHANGE it raises afterwards. A VALUE update on a face whose type is unchanged (ramping an inflow jet, a lid's tangential speed) is allowed at any time and takes effect the next step. A half-domain closed by a 'slip' face reproduces the full symmetric domain pointwise (tests/kokkos/test_freeslip). Both grids (Solver and SolverColocated).
`set_domain_bc_profile`	set_domain_bc_profile(self, face: str, profile: ndarray[dtype=float64, order='C']) -> None Prescribe a per-position inlet velocity profile (Nb,Nc,3) over a face (one of '-x', '+x', '-y', '+y', '-z', '+z'; sets it to 'inflow'). Like set_domain_bc, the TYPE (making a face inflow for the first time) must be set BEFORE the geometry; updating the profile on a face that is ALREADY inflow is allowed after the geometry.
`set_dt`	set_dt(self, dt: float) -> None Set the time step dt (the caller's time unit); the momentum solve is scaled by 1/dt (well-conditioned at large dt). Under a physical domain the FIRST call pins the reference time the internal scales are built on; any later call (step_adaptive makes them) rebuilds the momentum operator at the next step.
`set_field`	set_field(self, name: str, array: ndarray[dtype=float64, order='F']) -> None Write a Fortran-order (nx,ny,nz) float64 array into a registered field's inner region (ghosts refilled on the next diagnostics.exchange_field/step).
`set_implicit_advection`	set_implicit_advection(self, on: bool) -> None Use implicit-FOU advection with deferred-correction TVD.
`set_instance_motion`	set_instance_motion(self, instance: int, lin_vel: collections.abc.Sequence[float] = [0.0, 0.0, 0.0], ang_vel: collections.abc.Sequence[float] = [0.0, 0.0, 0.0], center: collections.abc.Sequence[float]
`set_instance_transform`	set_instance_transform(self, instance: int, translation: collections.abc.Sequence[float], quat: collections.abc.Sequence[float] = [0.0, 0.0, 0.0, 1.0]) -> None Move one scene instance (quat as (x,y,z,w)). Takes effect at the next rebuild_geometry(): the SDF, the cut-cell overlay, the apertures and the pressure operator are all derived from the transforms, so moving one without rebuilding would run the solver on stale geometry.
`set_mass_flux`	set_mass_flux(self, array: ndarray[dtype=float64, order='F']) -> None P0: prescribe a per-cell interfacial mass flux from a Fortran-order (nx,ny,nz) float64 array (== set_field('mdot') plus the ghost fill). Turns the thermal mass flux OFF.
`set_mass_flux_uniform`	set_mass_flux_uniform(self, mdot: float) -> None P0: prescribe a UNIFORM interfacial mass flux (kg m^-2 s^-1, solver units; positive = evaporation, liquid consumed). Turns the thermal mass flux OFF.
`set_mu`	set_mu(self, mu: float) -> None Set the dynamic viscosity mu (the caller's units); a change after the geometry rebuilds the momentum operator at the next step.
`set_phase_change_energy`	set_phase_change_energy(self, enabled: bool, rho_cp_gas: float = 0.0, rho_cp_liquid: float = 0.0) -> None CONSISTENT rhoc_pT transport (VOF_PLAN section 9 item 6) for the scalar named by set_phase_change_thermal; `enabled=False` returns to the constant-diffusivity scalar operator and the Koren TVD advective term (the rung P0/P1 energy path). Two things change together: (1) TRANSPORT: H = (rho c_p) T is advected with the colour advection's OWN geometric fluxes, sweep order and frozen dilation flag (vof/energy_advect.hpp), and T is recovered as H / (rho c_p)(C^{n+1}) -- so a uniform temperature is preserved EXACTLY at any heat capacity ratio. With two different fluxes the heat carried into a mixed cell is divided by a capacity built from another flux: an error of order d(rho c_p), i.e. ~2000x at water/steam, which is the artificial interfacial heating this removes. (2) DIFFUSION: the implicit operator becomes A_C = rho c_p(C)/dt + sum_f k_f open_f with k_f the arithmetic mean of the cells' k(C) (the k_gas/k_liquid of set_phase_change_thermal) EXCEPT at a face touching a per-cell Dirichlet (interfacial) cell, where the pure neighbour's own k is used -- a Dirichlet row is an identity row, so that coefficient's only job is the conductance with which the pure cell reaches a boundary condition that already sits at the interface. Units: rho*c_p in J/(cell^3 K), k in W/(cell K). Requires set_phase_change_thermal.
`set_phase_change_thermal`	set_phase_change_thermal(self, enabled: bool, scalar: str = '', t_sat: float = 0.0, k_gas: float = 0.0, k_liquid: float = 0.0, r_int: float = 0.0) -> None P1: compute mdot each step from a registered scalar (the temperature) instead of prescribing it; `enabled=False` stops that (back to the prescribed mdot field): mdot = ( k_g dT_g/dn - k_l dT_l/dn ) / h_lv , n = m/
`set_porous_continuity`	set_porous_continuity(self, on: bool = True) -> None Enable the volume-averaged (porous) continuity for unresolved CFD-DEM (staggered only): the projection enforces d(eps)/dt + div(eps u) = 0 instead of div(u)=0, so the fluid velocity is NOT solenoidal where the void fraction changes (bubbling/expansion). Binds the 'eps' field (void fraction from the particle deposition, created seeded to 1 if absent; the coupling writes it each step BEFORE step()). eps=1 everywhere reduces exactly to div(u)=0. Pair with max_porous_residual() for the meaningful convergence check.
`set_pressure_bottom`	set_pressure_bottom(self, mode: str) -> None Coarse-level (bottom) solve of the pressure multigrid. A V-cycle converges at a domain-independent rate only if its COARSEST level is effectively solved, and a geometric hierarchy cannot always get small enough: an axis stops coarsening once it turns odd, and under MPI once ANY rank's block turns odd -- so at fixed cells/rank the coarsest GLOBAL grid grows with the rank count and the bottom is progressively under-solved. 'auto' (DEFAULT) agglomerates the coarsest level onto a global operator and solves it exactly whenever it exceeds set_pressure_bottom_extent (4) cells on any axis, and uses the cheap smoothed bottom otherwise (byte-identical to 'smoother' then). 'smoother' = never agglomerate (legacy). 'agglomerated' = always. Measured on one GPU (2048x64x64 channel): a smoothed bottom needs 13.5 pressure iterations/step at 4 levels and 6.0 at 6, against 4.4 at full geometric depth; agglomerated it is 4.0 at BOTH depths, and faster in wall-clock than the deep hierarchy. Works on the cut-cell IBM and ghost-projection paths (per-fluid-component null-space projection, 2026-08-13).
`set_pressure_bottom_extent`	set_pressure_bottom_extent(self, cells: int) -> None The threshold set_pressure_bottom('auto') applies: agglomerate the coarsest level onto a global operator once the coarsest GLOBAL grid exceeds this many cells on any AXIS (not cell count -- what a few smoother sweeps cannot fix is a mode spanning many cells along one axis). DEFAULT 4.
`set_pressure_chebyshev`	set_pressure_chebyshev(self, on: bool, max_iter: int = 120, rtol: float = 1e-09) -> None Use the communication-light Chebyshev pressure accelerator. Mutually exclusive with the two CG drivers (on=True clears an FCG selection); on=False deselects Chebyshev and falls back to FCG if that is selected, otherwise to MG-PCG.
`set_pressure_fcg`	set_pressure_fcg(self, on: bool, max_iter: int = 200, rtol: float = 1e-08) -> None Use the FLEXIBLE MG-CG pressure accelerator: the same V-cycle-preconditioned CG as set_pressure_pcg, the same stopping estimate/mean removal, and the same cap+tolerance (they are shared -- it is the same solve), differing only in the beta recurrence (Polak-Ribiere r^T(z_{k+1}-z_k)/r^T z_k instead of Fletcher-Reeves). Costs one extra level-0 vector and one extra global dot per iteration; buys tolerance of a preconditioner that is not symmetric with respect to the fine operator. On a symmetric preconditioner the two betas are algebraically identical, so FCG reproduces PCG's iteration count. Unlike set_pressure_pcg this flag GENUINELY selects: on=True clears the Chebyshev selection (so it works after diagnostics.set_density_mode/set_porous), on=False returns to MG-PCG, and a later set_pressure_chebyshev(True, ...) wins over it.
`set_pressure_geometry`	set_pressure_geometry(self, sdf: ndarray[dtype=float64, order='F']) -> None Set an all-fluid SDF for the cut-cell pressure operator without an immersed solid (the channel/BFS domain-BC path). For a no-slip immersed BODY in an inflow/outflow domain, call set_solid(sdf, cutcell_pressure=True) instead -- do NOT also call this (a second geometry setter overwrites the SDF and wipes the solid).
`set_pressure_multigrid`	set_pressure_multigrid(self, on: bool, levels: int = 4) -> None Set the pressure multigrid depth (levels=1 => pure RB-GS, no coarse grid).
`set_pressure_pcg`	set_pressure_pcg(self, on: bool, max_iter: int = 200, rtol: float = 1e-08) -> None Use the MG-PCG pressure accelerator (single-GPU default) and set its iteration cap and relative tolerance. on=True GENUINELY selects: it clears both competing selections (Chebyshev and FCG), so it works after diagnostics.set_density_mode / set_porous -- last set wins. on=False RAISES: MG-PCG is the terminal fallback of the driver dispatch, so it cannot be deselected on its own -- select the driver you want instead (set_pressure_chebyshev(True, ...) or set_pressure_fcg(True, ...)). Under diagnostics.set_ghost_projection the nonsymmetric operator is solved by BiCGStab; this call still sets that solve's cap/tolerance (they are shared) and does not change its method.
`set_pressure_solver_params`	set_pressure_solver_params(self, iters: int) -> None Set the pressure smoother iteration count.
`set_pressure_telescope`	set_pressure_telescope(self, on: bool) -> None Coarse-level TELESCOPING of the pressure multigrid (multi-rank). The geometric hierarchy coarsens a level in place, which needs every rank's block even on that axis; once a block turns odd the hierarchy used to STOP, leaving a coarsest grid far too large to solve (384^3 on 1536 ranks: blocks die at 3x6x4 with the coarsest extent still 48, and the pressure iteration count grows 16.6 -> 38.7 across the ladder). With this on, a blocked level instead merges ORB siblings onto fewer ranks (the merged block's origin is its parent's split value, so parity is restored), moves the residual down / correction up within rank groups, and the roots keep coarsening on a sub-communicator to a tiny bottom. Idle ranks below the telescope point cost nothing that matters; the iteration count becomes a property of the problem rather than the rank count. ON by default since 2026-09-02 (the FoxBerry ladder); False restores the in-place-only hierarchy, byte-identical to before telescoping existed. See suite/docs/archive/MG_TELESCOPING_PLAN.md.
`set_property_model`	set_property_model(self, target: str, kind: str, field: str, params: collections.abc.Sequence[float] = [], field2: str = '') -> None Register a device closure writing a property/body-force field from input field(s). target: a registered field (a property 'mu'/'rho'/… or a body-force component 'force_x'/'force_y'/'force_z'). kind: 'linear' (params [p0,p1,p2]: p0+p1field+p2field2), 'boussinesq' (params [rho0,g,beta,T0]: rho0gbeta(field-T0) buoyancy), 'arrhenius' (params [mu_ref,B,Tref]: mu_refexp(B*(1/field-1/Tref))). Applied at the top of step().
`set_property_table`	set_property_table(self, target: str, field: str, x: collections.abc.Sequence[float], y: collections.abc.Sequence[float]) -> None Register a tabulated property: target = piecewise-linear interpolation of (x, y) at the input field value (x ascending, clamped at the ends).
`set_rho`	set_rho(self, rho: float) -> None Set the fluid density rho (the caller's units). Under a physical domain the FIRST call pins the reference density the internal scales are built on; a later change rebuilds the momentum operator at the next step.
`set_scalar_bc`	set_scalar_bc(self, name: str, face: str, type: str, value: float = 0.0) -> None Scalar boundary condition on a domain face ('-x', '+x', '-y', '+y', '-z', '+z'): type 'periodic', 'neumann' (zero flux, adiabatic) or 'dirichlet' (the given value). Single-rank.
`set_scene`	set_scene(self, node_ints: ndarray[dtype=int32, order='C'], node_reals: ndarray[dtype=float64, order='C'], instance_ints: ndarray[dtype=int32, order='C'], instance_reals: ndarray[dtype=float64, order='C'], periodic: bool = False) -> None Install an analytic geometry scene from core's flat encoding (3 ints + 16 reals per node, 2 ints + 18 reals per instance -- reals[17] is the centre-pinned flag; legacy 17-real records are accepted, reading an all-zero centre as 'follows the body'). Coordinates are the caller's own PHYSICAL ones when the solver was given an extent -- so ONE scene serves flow, dem and voro unchanged -- and CELL UNITS on the global inner grid otherwise (cell (i,j,k)'s centre at (i,j,k) -- half a cell from cell_centers(), which returns i + 1/2 in cell units). The scene is replicated on every rank, so scene-derived geometry needs no communication and -- unlike diagnostics.set_exact_crossings -- is NOT single-rank only. periodic=True treats the scene as min-image periodic over the global grid (one instance per body, no images); periodic=False leaves images to the caller.
`set_solid`	set_solid(self, sdf: ndarray[dtype=float64, order='F'], cutcell_pressure: bool = False) -> None Set the solid SDF as a Fortran-order (nx,ny,nz) float64 array sampled at cell_centers() (negative inside the solid, positive in fluid, in the solver's length unit). cutcell_pressure=True enables the open-face-weighted cut-cell pressure operator (proper no-slip); it composes with domain BCs, so this is the single call for a no-slip immersed body in an inflow/outflow domain. The multigrid coarse operators are always the rediscretized per-level cut-cell operators (the former `pressure_coarse=` mode selector is gone).
`set_solid_from_scene`	set_solid_from_scene(self, cutcell_pressure: bool = True) -> None Sample the installed scene onto this rank's inner grid and install it as the solid, entirely on device (no nxnynz float64 host round trip).
`set_state`	set_state(self, u: ndarray[dtype=float64, order='F'], v: ndarray[dtype=float64, order='F'], w: ndarray[dtype=float64, order='F']) -> None Upload an initial velocity field (u,v,w each a Fortran-order (nx,ny,nz) float64 array).
`set_superficial_velocity`	set_superficial_velocity(self, enabled: bool, axis: str = 'x', velocity: float = 0.0) -> None Hold the SUPERFICIAL velocity along a periodic axis at `velocity` (physical units): the volume mean of that velocity component, i.e. the volume flux of the MIXTURE through any cross-section divided by the cross-section's full area. In multiphase flow 'velocity' alone is ambiguous between superficial and interstitial (per-phase); this is the superficial one. At the end of every step() one uniform shift is added to that component on every face; on an axis with no boundary a uniform shift of the normal velocity leaves every cell divergence unchanged, so the field stays discretely divergence-free. `velocity = 0` is a closed (batch) column: the flow can circulate but not pass through (TBFsolver's constant flow rate at 0). The mean is a device reduction, global under MPI; single-rank the step makes no host round trip for it. Staggered Solver, all-fluid domain (set_pressure_geometry), periodic `axis` only -- step() raises otherwise. set_superficial_velocity(False) turns it off. diagnostics.last_superficial_velocity_shift() is the shift the last step applied.
`set_surface_tension`	set_surface_tension(self, sigma: float) -> None VoF rung V4: turn on the BALANCED-FORCE continuum surface force with coefficient sigma (0 = off, the default). Enables VoF and the curvature cascade, which then runs once per step at the head, from the same colour field the density closure sees. The force at the staggered velocity unknown u_c(i) is F = sigma * kappa_f * (C(i) - C(i - s_c)) / h i.e. the colour difference is the projection's OWN face difference — the identical discrete gradient operator the pressure uses (Francois et al. 2006; Popinet 2009). That is the whole content of 'balanced force': with a constant kappa the force is exactly a discrete gradient, so the projection annihilates it and a static drop stays at machine zero. Face-interpolating a cell-centred sigmakappagrad(C) instead — the obvious way to reuse the per-cell body force — is not a discrete gradient of anything and leaves spurious currents of order sigmakappa/mu that no curvature accuracy removes. UNITS: sigma is a force per unit LENGTH in the caller's own units -- the system rho, mu and set_body_force are stated in (so per CELL without a physical domain, where the cell size is 1). SIGN: kappa is positive for a convex blob of liquid and C is the liquid fraction, so the equilibrium pressure is P = sigmakappa*C + const — the Young-Laplace overpressure INSIDE the drop. Surface tension is EXPLICIT, so step() enforces the Brackbill capillary time step (see capillary_dt()).
`set_velocity`	set_velocity(self, c: int, array: ndarray[dtype=float64, order='F']) -> None Write component c's inner velocity from a Fortran-order (nx,ny,nz) float64 array (solid rows re-masked). Initial-condition hook: a uniform stream around a fixed body is the Galilean twin of the same body towed through fluid at rest.
`set_velocity_multigrid`	set_velocity_multigrid(self, on: bool, levels: int = 4, vcycles: int = 8) -> None Enable velocity (momentum) multigrid for the implicit diffusion solve.
`set_velocity_residual_tolerance`	set_velocity_residual_tolerance(self, rtol: float) -> None Residual-based momentum stop: a component's implicit solve ends once max
`set_vof`	set_vof(self, array: ndarray[dtype=float64, order='F']) -> None Set the colour field from a Fortran-order (nx,ny,nz) float64 array: C = the LIQUID volume fraction of the cell, in [0,1]. Enables VoF if it is not on yet. Initialise it SHARP (0/1 with the exact fraction in the interface cells) — geometric VoF keeps a sharp interface and a diffuse initial profile is simply a different problem.
`set_vof_backflow`	set_vof_backflow(self, face: str, value: float = 0.0) -> None inletOutlet colour on OUTFLOW face 'face' ('-x' ... '+z'; default 0 = gas): where the boundary face velocity points back INTO the domain, the colour ghost carries this value instead of the zero-gradient copy. This is OpenFOAM's inletOutletFvPatchField (Rusche 2002 thesis section 4), the standard VoF outlet — an outlet is a place where you know what leaves (whatever is inside) but must STATE what comes back. Zero-gradient on a reversing face re-injects whatever happens to be sitting at the outlet, which for a draining film feeds the film back into the domain. Where the fluid leaves, zero-gradient is kept.
`set_vof_block_assign`	set_vof_block_assign(self, mode: str, every: int = 0) -> None Master assignment of the block container (rung W1). 'round-robin' by block id (rung W0's, blind to block SIZE), 'lpt' = LONGEST-PROCESSING-TIME greedy on the block cell counts, 'orb' = weighted ORB over a 1-D block space (core's BlockDecomposer<1> with the cell counts as weights). 'every' > 0 re-runs the assignment every 'every' block steps and MIGRATES the colour of any block that changed master (nothing else in a block is state, so the migration is exact and the bitwise gates hold across it). All three are pure functions of the REPLICATED block table, so every rank computes the same assignment without an exchange -- which is what lets a re-assignment happen mid-run without breaking a bitwise gate. LPT is the measured winner (see the WO-W12 findings): the ORB's blocks must be CONTIGUOUS in block id, which LPT is free of.
`set_vof_cfl_limit`	set_vof_cfl_limit(self, value: float) -> None Weymouth-Yue boundedness cap on the interface-local Courant number. Default 0.25 — the PROVEN 3D bound 1/(2(N-1)) from Weymouth's thesis eq. A.33; the widely-quoted 0.5 is the 2D value. step() throws when it is exceeded. Conservation is INDEPENDENT of boundedness (the telescoping is algebraic), so an over-CFL run loses 0<=C<=1, not volume.
`set_vof_inflow`	set_vof_inflow(self, face: str, value: float) -> None Colour of the fluid ENTERING through inflow face 'face' ('-x', '+x', '-y', '+y', '-z', '+z'); 1 = liquid, 0 = gas. The face must already be an inflow (set_domain_bc(face, 'inflow', ...)) or this raises. The value may be FRACTIONAL, and it then means 'this fraction of the incoming flux is liquid' — a flux statement, not a sub-cell interface position. That is exactly how it is implemented: on a face whose donor is outside the domain the geometric flux is replaced by the algebraic C_donor*a (wyFaceFluxBc), because a uniform prescribed ghost band has no usable MYC normal and the inflow colour is boundary DATA. Setting this also makes the property ghosts of that face follow the inflow colour through the registered closures (rho_ghost = rho(C_inflow)) instead of copying the interior — without it the inlet FACE density is the interior's, wrong by up to the density ratio at a liquid inlet into a gas domain. Default (nothing set): the zero-gradient copy, i.e. today's behaviour, bit for bit.
`set_vof_inflow_profile`	set_vof_inflow_profile(self, face: str, profile: ndarray[dtype=float64, order='C']) -> None Per-position inflow colour (Nb,Nc) on the INNER grid of the face's two perpendicular axes, resampled with the same clamp rule set_domain_bc_profile uses for the velocity. This is how a liquid distributor over part of an inlet is expressed: the colour profile says WHERE liquid enters, the velocity profile says how fast.
`set_wall_slip_length`	set_wall_slip_length(self, lambda_cells: float) -> None Rung V6b (WO-V6b) - the VELOCITY half of the dynamic contact line. In the Robust-Scaled cut-cell IBM closure the TANGENTIAL wall datum stops being no-slip and becomes the NAVIER condition u_t(wall) - u_body = lambda * d(u_t)/dn with `lambda_cells` the slip length in the caller's own units (CELLS without a physical domain, where the cell size is 1); the model needs lambda < one cell. The wall-NORMAL component is untouched: the wall stays impermeable, with the moving-body velocity as its datum. 0 (the default) restores the validated no-slip closure BIT-IDENTICALLY - the closure polynomials are the lambda = 0 members of the same family, and the lambda > 0 branch is never entered. It is ONE lambda with set_contact_angle_dynamic: either call overwrites the shared value (last call wins), so the Cox-Voinov inner cut-off and the momentum wall closure cannot disagree. Only this call switches the MOMENTUM half on, so every WO-V6 result taken with the angle half alone is unchanged. Why it exists, measured: with no-slip the contact line advances only at the scheme's own numerical slip, ~1/180 of Lucas-Washburn (WO-V6 finding 6), and WO-V7's pore-scale campaign found every IMBIBITION verdict inverted below the capillary number at which the imposed velocity crosses that numerical slip velocity. Typical values are 0.01-0.5 cells; state lambda with every dynamic-wetting result. Staggered grids only; needs set_solid(..., cutcell_pressure=True). Changing it rebuilds the three cut-cell overlays and the momentum operator in place (the velocity field is NOT reset). APPROXIMATIONS, both stated because they are invisible from the call: the tangential projector is taken DIAGONAL (the cross terms -n_c n_j u_j are dropped; they vanish exactly for an axis-aligned wall and are O(n_c n_j) otherwise), and an axis whose BOTH neighbours are solid (a one-cell fluid gap) keeps the no-slip closure - counted by diagnostics.wall_slip_sandwich_cells().
`size`	size(self) -> int The number of ranks in MPI_COMM_WORLD.
`spacing`	The cell size [dx, dy, dz] = extent/cells as three floats; [1.0, 1.0, 1.0] without a physical domain (the cell-unit default). Three spacings within 1e-12 relative are snapped to one value, so an isotropic box is exactly isotropic.
`step`	step(self) -> None Advance the solver one time step (semi-implicit: diffusion + projection). ATOMIC across the two explicit two-phase stability caps: the Weymouth-Yue boundedness cap and the Brackbill capillary limit are evaluated at the HEAD of the call, so a rejected dt leaves every field bitwise as it was and a retry at a smaller dt is exact. (Before 2026-09-04 the momentum half had already advanced by the rejected dt, so catch-and-halve desynchronised colour and momentum.) The WY pre-check reads the field the call starts from; where the projection itself accelerates the field inside the step, use step_adaptive().
`step_adaptive`	step_adaptive(self, cfl_target: float = 0.4, capillary_cfl: float = 0.4, dt_max: float = inf) -> float One step at a dt re-picked from the CURRENT state; returns the dt used. dt = min(cfl_target * cfl_dt, capillary_cfl * capillary_dt, dt_max) with both limits taken from vof_step_limits() at the moment of the call, then set_dt(dt) + step(). Needs enable_vof(). WHY: both limits are INSTANTANEOUS. Re-picking every ten steps is safe for a free bubble or a falling film and is NOT safe in a packing -- a throat jet can grow max
`surface_tension`	surface_tension(self) -> float The surface-tension coefficient (0 when the CSF is off).
`sync_porous_prev`	sync_porous_prev(self) -> None Reseed eps^n = eps^{n+1} (d(eps)/dt = 0 for the next step). Call it once the first void-fraction field has been written into 'eps' and before the first step, so that step sees no spurious d(eps)/dt source; any later call simply restarts the d(eps)/dt bookkeeping from the current field.
`unit_scales`	The factors between the caller's units and the ones the solver computes in. The solver keeps computing on the UNIT LATTICE and folds the metric into constants at its API boundary, so every array reached through the RAW field registry (diagnostics.field_view / get_field / set_field) is in those INTERNAL units, unlike get_u/get_p, which convert. `identity` is True when no extent was given, and then every factor is exactly 1.0 and internal == physical. Multiply a physical value by *_to_internal, and an internal one by force_to_physical (a total force) to come back.
`update_properties`	update_properties(self) -> None Apply all registered property/force closures now (also done at the top of step()).
`velocity_residual_tolerance`	velocity_residual_tolerance(self) -> float The momentum residual tolerance in force (resolves the follow-the-pressure default).
`vof_bc_volumes`	vof_bc_volumes(self) -> list[float] Signed liquid volume that crossed each of the six domain faces during the LAST colour advection, in cell-volume units, POSITIVE for liquid ENTERING the domain. Local to this rank. Six entries in face order (-x,+x,-y,+y,-z,+z); only faces carrying a domain BC are meaningful.
`vof_bc_volumes_total`	vof_bc_volumes_total(self) -> list[float] The same, accumulated since enable_vof() or the last reset_vof_bc_volumes().
`vof_block_assign`	vof_block_assign(self) -> str The current master-assignment mode ('round-robin'
`vof_block_color`	vof_block_color(self, id: int) -> numpy.ndarray[dtype=float64] One marker's OWN inner color as a Fortran-order (nx,ny,nz) float64 array over its block box (diagnostics.vof_block_stats()['lo'/'hi']); empty on a rank that does not master it. This is the block's ONLY state, so {box, color} per marker is a COMPLETE checkpoint of the container -- and it is the only exact one: re-seeding from the union color field (enable_vof_blocks_from_field) hands each of two TOUCHING markers a slice of the other, because that gather is a copy of the union clipped to the seed extent. enable_vof_blocks_from_colors(boxes, colors) restarts the container from these.
`vof_blocks_enabled`	vof_blocks_enabled(self) -> bool True while the block container exists (enable_vof_blocks with seeds; None drops it).
`vof_cfl_limit`	vof_cfl_limit(self) -> float The Weymouth-Yue boundedness cap currently in force.
`vof_curvature`	vof_curvature(self) -> numpy.ndarray[dtype=float64] The curvature field's inner region as a Fortran-order (nx,ny,nz) float64 array, in 1/h in 1/LENGTH. Call compute_vof_curvature() first; get_field('kappa') is the RAW solver-internal field, which is 1/h.
`vof_curvature_branch`	vof_curvature_branch(self) -> numpy.ndarray[dtype=float64] Which cascade branch produced each cell's curvature, as a Fortran-order (nx,ny,nz) float64 array: 0 not interfacial, 1 height function, 2 height function in a non-preferred direction, 3 mixed height-position fit, 4 PLIC-volumetric paraboloid fit, 5 the same with the rank-deficient 3-parameter model, 6 NO estimate. == get_field('kappa_branch').
`vof_has_geometry`	vof_has_geometry(self) -> bool True when the colour advection is running the CUT-CELL (openness-weighted) kernels, i.e. an immersed solid is present and set_solid ran with cutcell_pressure=True. False means the uncut rung-V1 kernels are running, byte-identically to a solid-free build.
`vof_interface_area`	vof_interface_area(self) -> float Total interfacial area of the colour field, in CELLS squared (h^2), summed over the inner region and globally reduced under MPI. Uses the geometry diagnostics.set_phase_change_area selects, so the number a page quotes and the number the phase change integrates are the same one. Needs enable_vof; phase change need not be on. A sphere of radius R cells reads 4 pi R^2 to about 0.2-0.8 % once its colour field is resolved (WO-P3c).
`vof_max_courant`	vof_max_courant(self) -> float INTERFACE-LOCAL Courant number max
`vof_momentum_enabled`	vof_momentum_enabled(self) -> bool Whether momentum-consistent transport is on.
`vof_step_limits`	vof_step_limits(self) -> dict Both explicit two-phase step limits at the current state, and which one binds: the largest dt the interface-local Weymouth-Yue CFL cap admits ('cfl_dt'), the Brackbill capillary limit ('capillary_dt'), and min(cfl_dt, capillary_cfl*capillary_dt) with a 'capillary_binds' flag. At pore-scale capillary numbers the capillary limit is expected to bind first.
`wall_slip_length`	wall_slip_length(self) -> float The Navier slip length in force, in the caller's units (0 = no-slip).

`SolverColocated`¶

One incompressible Navier-Stokes solver on a block of the grid (staggered MAC, the reference) -- or, for SolverColocated, the same API with cell-centered velocities and the ABC approximate projection. Construct on a (physical) domain; set_rho / set_mu / set_dt and the domain BCs; give it a geometry (set_solid, set_pressure_geometry or a scene); step(); read out with the get_* arrays. Developer instruments and tuning live on .diagnostics.

Method / property	Description
`add_field`	add_field(self, name: str) -> None Register a new zero-initialised cell-centred field on the grid (for transported scalars or material properties). Idempotent.
`add_scalar`	add_scalar(self, name: str, diffusivity: float = 0.0, scheme: str = 'koren', iters: int = 50) -> None Register a transported scalar (temperature/concentration/...): constant diffusivity (grid units), advection scheme 'fou' (first-order upwind), 'koren' (Koren TVD, the default) or 'sou' (second-order upwind), and RB-GS diffusion sweeps. The scalar is a registered field (get_field/set_field/diagnostics.field_view). Requires geometry (set_solid/set_pressure_geometry) for the openness-weighted operators.
`advance_scalars`	advance_scalars(self) -> None Advance all registered scalars one dt with the current velocity (also done by step()).
`advect_vof`	advect_vof(self, dt: float) -> None KINEMATIC colour advection: advance C ONCE with the solver's CURRENT face velocity over `dt` (seconds), with NO Navier-Stokes step. This is the advection-benchmark entry point (Zalesak, LeVeque, and the cut-cell conservation gates): a frozen projected velocity advecting a colour field is a pure statement about the transport scheme, with the momentum and pressure solves out of the picture.\n\nIt THROWS if the current velocity is not discretely divergence-free to 1e-10, measured with max_open_divergence_projected() (the NON-mutating diagnostic). That is not a nicety: Weymouth-Yue's exact conservation is CONDITIONAL on sum_f o_f u_f = 0 per cell, because the dilation term adds H(C-1/2) times that residual to every full cell's volume budget. Run step() to a steady state (or project()) and advect with the solver's own output; never with an analytic sample, which is solenoidal only to O(h^2). It ALSO throws when no cut-cell pressure operator exists (WO-R2): without one the divergence guard measures nothing at all — max_open_divergence() returns 0.0 — and a cell-centre-sampled LeVeque field (true max
`advect_vof_blocks`	advect_vof_blocks(self, dt: float) -> None Advance every block by dt with the CURRENT (projected) face velocity and union the result into 'C' — the block twin of advect_vof(dt), and it carries the same precondition: it RAISES unless the face field is discretely divergence-free to 1e-10, because Weymouth-Yue's exact conservation is conditional on it. Per step: gather the face velocity from the ranks that own each block's cells to its master (plain Isend/Irecv — the block table and the decomposition are both replicated, so every rank knows every message size and an NBX handshake would only rediscover it), run the three sweeps on the dense block, re-centre the box, replicate the table, and scatter the inner colour back with UNPACK_MAX. NOTE on the union: max() from an empty union CLIPS the negative Weymouth-Yue round-off residue to an exact 0 — measured up to 6.2e-17 on the LeVeque gate, and every measured union/global-field difference was exactly that and nothing else.
`apply_phase_change`	apply_phase_change(self, dt: float) -> None KINEMATIC phase-change step: build mdot / the PLIC areas / the normals, deposit the divergence source (census only — nothing is projected) and apply the interface regression over `dt`. No Navier-Stokes step, no advection. This is the P0a / P1 driver: an interface regressing under a prescribed or thermally-computed flux, with the momentum and pressure solves out of the picture.
`balanced_force_projection`	Whether the balanced-force projection runs (see set_balanced_force_projection).
`block_origin`	block_origin(self) -> list[int] This rank's inner-block origin in GLOBAL cells ([0,0,0] single-rank). Shift the coupling deposit origin by this so particles in global coordinates land in the local block.
`capillary_dt`	capillary_dt(self) -> float The Brackbill (1992) capillary time-step limit sqrt((rho_1+rho_2) h^3 / (4 pi sigma)), +inf when surface tension is off. Denner & van Wachem (2015) verified this IS the stability boundary of an explicit CSF — the prefactor 1/(4 pi), the h^{3/2} scaling and the SUM of the phase densities (both phases oscillate). The density sum is the declared phase pair under enable_vof_momentum, otherwise min(rho)+max(rho) over the field.
`cell_centers`	cell_centers(self) -> tuple The physical cell-center coordinates of THIS rank's inner block as three 1-D float64 arrays (x, y, z), i.e. origin + (i + 1/2) * spacing per axis. np.meshgrid(s.cell_centers(), indexing='ij') is the grid an SDF for set_solid is sampled on. Under MPI these are the LOCAL block's centers in GLOBAL coordinates. CELL-UNIT TRAP: without an extent this is i + 1/2, while set_scene puts cell i's centre at i -- the same body given as a scene and as an SDF sampled here lands half a cell apart. With an extent both use origin + (i + 1/2) spacing. Unifying the cell-unit convention changes results; deferred to 2.0.0.
`cells`	This rank's inner block resolution (nx, ny, nz) -- the shape of every field array (get_u, get_p, the SDF for set_solid). Under MPI it is the LOCAL block from mpi_block().
`clear_divergence_sink`	clear_divergence_sink(self) -> None Drop the auto-balanced sink region.
`clear_divergence_source`	clear_divergence_source(self) -> None Drop the prescribed divergence source (the field stays registered).
`compute_vof_curvature`	compute_vof_curvature(self) -> dict VoF rung V3: compute the interface curvature from the CURRENT colour field and store it in the registered cell fields 'kappa' and 'kappa_branch'. Returns THIS RANK's branch census as a dict. 'kappa' is the RAW registered field, kappa = 2H in 1/h (cell units); vof_curvature() returns the same thing in 1/length. POSITIVE for a convex blob of liquid: a sphere of liquid of radius R cells reads +2/R, an infinite cylinder +1/R, a plane 0. The cascade is Popinet's (2009): a standard height function on 7-cell column sums of C over the 3x3 transverse patch in the direction of the largest
`contact_angle`	contact_angle(self) -> float The prescribed static contact angle in degrees (90 if none was set).
`decomposition_levels`	This solver's decomposition mode (0 = aligned ORB, >= 2 = coarse-first with that depth).
`decomposition_max_imbalance`	The load-imbalance budget the coarse-first search stays within.
`diagnostics`	The developer tier of this solver (instruments, tuning, ablations, internals).
`enable_cell_force`	enable_cell_force(self) -> None Allocate + register the per-cell body-force fields force_x/force_y/force_z and route them into the momentum RHS, for an external writer (e.g. CFD-DEM drag feedback) to fill directly via diagnostics.field_view('force_z'). They persist across steps until overwritten.
`enable_drag`	enable_drag(self) -> None Enable implicit (semi-implicit) linear drag for CFD-DEM: allocate the per-cell 'drag_beta' field (added to the momentum diagonal so a -beta(u-u_p) source is treated implicitly -> unconditionally stable for the stiff beta of a dense bed) plus force_x/y/z (which carry betau_p, the RHS target). Fill 'drag_beta' and 'force_*' via diagnostics.field_view each step.
`enable_phase_change`	enable_phase_change(self, rho_gas: float, rho_liquid: float, h_lv: float = 1.0) -> None Enable VoF phase change (VOF_PLAN section 9, rungs P0/P1; the Boyd & Ling 2023 / Malan 2021 pattern). Registers 'mdot' (the interfacial mass flux, kg m^-2 s^-1 in solver units, POSITIVE = evaporation) and 'pc_source' (the deposited divergence source, 1/s). What it does per step, at the HEAD of step(): (1) reconstruct the PLIC plane of every interfacial cell and take its POLYGON AREA A analytically (A =
`enable_vof`	enable_vof(self) -> None Enable the geometric (PLIC + Weymouth-Yue) VoF colour field 'C'. Registers 'C' as an ordinary cell field (get_field/set_field/diagnostics.field_view/closure input) and allocates the colour field's OWN g=3 working block + halo; the solver's G=2 is untouched. C is advected at the end of every step() with the just-projected face velocities. RUNG V2a SCOPE — read before using: there is NO surface tension (rung V4) and NO momentum-consistent transport (rung V2b). Mass and momentum are advected by different fluxes, so a mixed cell carries a spurious interfacial momentum source of order delta-rho; the literature (Rudman 1998, Arrufat 2021) is unambiguous that this breaks down around density ratio 1000 unless the resolution is absurd. V2a is valid only at MODEST density ratios for cases with motion. A high-ratio case AT REST (the hydrostatic acid test) is exact, because there is no momentum to mis-advect. COLLOCATED (SolverColocated) is supported since rung V8 (2026-09-02), ALL-FLUID only: the colour is advected by the PROJECTED face field uf_/vf_/wf_ (the field the ABC approximate projection makes exactly divergence-free, which is what Weymouth-Yue's conservation proof needs), and the CSF enters the implicit momentum predictor as the mass-adjoint face integral (doc/collocated_varrho_forces.md). An immersed solid on that grid THROWS (the cut-cell colour transport of rung V5a is staggered-only), and so does enable_vof_momentum. On the collocated grid the SPURIOUS-CURRENT number to read is the FACE field (get_uf/get_vf/get_wf): the CELL field additionally carries the approximate projection's invisible odd-even mode, which centerToFace annihilates and the projection therefore cannot remove. An IMMERSED SOLID is supported since rung V5a: the geometric fluxes are openness-weighted and the update is done in fluid-volume units, so sum eps_eff*C is conserved exactly against the projection's own openness-weighted divergence; it needs set_solid(..., cutcell_pressure=True) (the staircase operator has no face openness to weight with) and it approximates the solid-clipped flux polygon by (whole-cell PLIC slab) x (open area) — see diagnostics.vof_diagnostics()['clipped_volume'].
`enable_vof_block_csf`	enable_vof_block_csf(self) -> None Form the surface-tension force PER BLOCK (rung W2): each marker runs its own curvature cascade on its own dense box and forms the V4 balanced-force face force there (sigma kappa_f (C(i) - C(i-s))/h with the same selective face-curvature rule addCsfRhs uses), and the three face fields are scattered into the global RHS with UNPACK_SUM — TBFsolver's VOF.f90::computeSurfaceTension structure on the suite's kernels. WHY THE FORCE AND NOT THE CURVATURE IS SCATTERED: kappa is not additive and the union colour is a max, so a face between two OVERLAPPING markers has no single (kappa, dC) pair to build a force from. The force is the additive quantity, and forming it where each marker's own colour still exists is the only place the balanced-force pairing (the same face difference the projection's gradient uses) is available per marker. Requires set_surface_tension(sigma) and the block container; STAGGERED only. Once on, step() drives the whole two-phase stage through the blocks: the union C feeds the closures exactly as before, and the colour is advected by the blocks in the advect_vof slot.
`enable_vof_blocks`	enable_vof_blocks(self, seeds: collections.abc.Sequence[collections.abc.Sequence[float]]
`enable_vof_blocks_from_colors`	enable_vof_blocks_from_colors(self, boxes: collections.abc.Sequence[collections.abc.Sequence[int]], colors: list) -> None Restart the block container from a vof_block_color() checkpoint: one marker per (lo_x, lo_y, lo_z, hi_x, hi_y, hi_z) INNER box (exactly the 'lo'/'hi' of diagnostics.vof_block_stats(), NOT grown by the margin again) with its colour written directly. Exact whatever the markers are doing -- nothing is gathered out of the union field.
`enable_vof_blocks_from_field`	enable_vof_blocks_from_field(self, boxes: collections.abc.Sequence[collections.abc.Sequence[int]]) -> None Seed the block container from the CURRENT colour field: one marker per given global index box (lo_x, lo_y, lo_z, hi_x, hi_y, hi_z; half-open), its colour copied out of the field set_vof installed. The general seeding path — enable_vof_blocks(spheres) is a convenience over it — and the way a marker of any shape enters (a quasi-2-D Hysing cylinder, a scanned bubble, an arbitrary blob found with a host connected-component labelling). The boxes are the BUBBLE extents; the container grows each by its own 3-cell margin. They must not share a cell: the seeding gather is a copy, not a union, so a cell inside two boxes would be handed to both markers.
`enable_vof_momentum`	enable_vof_momentum(self, rho_gas: float, rho_liquid: float) -> None Advect rhou on the HALF-SHIFTED momentum control volumes with the SAME geometric fluxes, the same sweep order and one frozen dilation flag as the colour field, then recover u = (rho^c u)/rho^c. This is what makes the scheme usable above density ratio ~100 (Rudman 1998; Arrufat 2021: a raindrop accurate within 15% at 15 cells/diameter instead of ~200). The two phase densities are REQUIRED, and must be the ones the rho closure produces at C=0 and C=1 — the momentum flux is rho_gas(a-F) + rho_liquidF with F the geometric liquid flux, so the scheme has to know which density each phase carries. Turning this on MOVES the VoF advection from the end of step() to its head (the momentum advection must precede the predictor that consumes it), so the advecting field is u^n — the previous step's projected output. u at step 0 must therefore be discretely divergence-free (rest or a uniform field both are). Requires the variable-density path, staggered layout, explicit advection, no immersed solid and no porous continuity, and is refused together with the per-bubble block container (enable_vof_blocks, in either order): the momentum rides the structured colour's fluxes, not the markers'.
`extent`	The GLOBAL domain side lengths. Equal to the cell counts without a physical domain.
`field_names`	field_names(self) -> list[str] Names of all registered fields (velocity u/v/w, p, sdf, plus any added), sorted.
`get_field`	get_field(self, name: str) -> numpy.ndarray[dtype=float64] Return a registered field's inner region as a Fortran-order (nx,ny,nz) float64 array.
`get_ox`	get_ox(self) -> numpy.ndarray[dtype=float64] The -x face openness field of the cut-cell operator: the fluid area fraction (0..1, dimensionless) of the face between cells (i-1,j,k) and (i,j,k), as a Fortran-order (nx,ny,nz) float64 array over THIS rank's inner cells. Computed by set_solid from the SDF (1 everywhere without a solid); it is what the momentum operator and the geometric cut-cell projection see. get_ox_proj is the openness the PROJECTION actually conserves fluxes through, which differs under diagnostics.set_ghost_projection (binary/coupled faces).
`get_ox_proj`	get_ox_proj(self) -> numpy.ndarray[dtype=float64] -x face openness whose fluxes the projection CONSERVES (binary/COUPLED under diagnostics.set_ghost_projection, geometric cut-cell otherwise). Use for flux bookkeeping (peclet.pnm extract_network_flow).
`get_oy`	get_oy(self) -> numpy.ndarray[dtype=float64] The geometric -y face openness per inner cell, (nx,ny,nz); see get_ox.
`get_oy_proj`	get_oy_proj(self) -> numpy.ndarray[dtype=float64] -y face openness the projection conserves (see get_ox_proj).
`get_oz`	get_oz(self) -> numpy.ndarray[dtype=float64] The geometric -z face openness per inner cell, (nx,ny,nz); see get_ox.
`get_oz_proj`	get_oz_proj(self) -> numpy.ndarray[dtype=float64] -z face openness the projection conserves (see get_ox_proj).
`get_p`	get_p(self) -> numpy.ndarray[dtype=float64] Return the physical pressure as a Fortran-order (nx,ny,nz) float64 array (index [x,y,z]).
`get_u`	get_u(self) -> numpy.ndarray[dtype=float64] Return the x-velocity component as a Fortran-order (nx,ny,nz) float64 array (index [x,y,z]).
`get_uf`	get_uf(self) -> numpy.ndarray[dtype=float64] Return the divergence-free FACE x-velocity (collocated: projected MAC field; staggered: == get_u).
`get_v`	get_v(self) -> numpy.ndarray[dtype=float64] Return the y-velocity component as a Fortran-order (nx,ny,nz) float64 array (index [x,y,z]).
`get_vf`	get_vf(self) -> numpy.ndarray[dtype=float64] Return the divergence-free FACE y-velocity (collocated: projected MAC field; staggered: == get_v).
`get_vof`	get_vof(self) -> numpy.ndarray[dtype=float64] The colour field's inner region as a Fortran-order (nx,ny,nz) float64 array (== get_field('C')).
`get_w`	get_w(self) -> numpy.ndarray[dtype=float64] Return the z-velocity component as a Fortran-order (nx,ny,nz) float64 array (index [x,y,z]).
`get_wf`	get_wf(self) -> numpy.ndarray[dtype=float64] Return the divergence-free FACE z-velocity (collocated: projected MAC field; staggered: == get_w).
`ghost_width`	ghost_width(self) -> int Ghost-layer width g of the velocity block (diagnostics.field_view returns an (n+2g) buffer).
`global_cells`	The GLOBAL grid resolution (gnx, gny, gnz); == `cells` single-rank, the numbers init_mpi gets under MPI (the CFD-DEM co-decomposition weight field is this shape).
`has_cutcell_pressure`	has_cutcell_pressure(self) -> bool True once the cut-cell pressure operator exists (set_solid or set_pressure_geometry was called). The porous continuity requires it; the coupling driver auto-sets an all-fluid geometry when absent.
`has_field`	has_field(self, name: str) -> bool Whether a field of this name is registered.
`has_moving_instance`	has_moving_instance(self) -> bool True when at least one scene instance carries a nonzero velocity.
`has_scalar`	has_scalar(self, name: str) -> bool Whether a transported scalar of this name is registered.
`has_scene`	has_scene(self) -> bool True once set_scene has installed an analytic CSG scene on this solver. False for geometry set any other way (set_solid with a sampled SDF, set_pressure_geometry), and false again if no scene has been given yet. set_solid_from_scene and the scene-derived per-body hydrodynamic force/motion bookkeeping all require this to be True first.
`hydro_force_torque`	hydro_force_torque(self) -> numpy.ndarray[dtype=float64] DEPRECATED (warns) -- the canonical names are hydro_force_torque_reaction() for the force and torque on the bodies, and diagnostics.hydro_force_torque_traction() for what this returns. It is NOT the force on the bodies: it is the reconstructed surface-traction integral, which under-reads the drag by ~30 % at any resolution (the viscous part; traction / reaction 0.685-0.730 over phi 0.008-0.45, N 24-128). Returns exactly what diagnostics.hydro_force_torque_traction() returns (shape (4, n_instances, 3): force, torque, pressure part, viscous part); see that method for the details. Where the reaction refuses (SolverColocated, porous, variable properties, domain BCs, implicit advection) the traction is the only estimate there is -- call it by its diagnostic name and allow for the under-read. Kept working under the suite's rule for a shipped name (docs/NAMING.md section 0); removal needs a major release.
`hydro_force_torque_reaction`	hydro_force_torque_reaction(self) -> numpy.ndarray[dtype=float64] THE hydrodynamic force and torque on each scene instance: the DISCRETE REACTION -- the momentum the fluid actually lost to each body, assembled from the composed step's budget (time term, body force, and the viscous fluxes at u; the pressure telescopes inside each owner region and needs no field). Shape (2, n_instances, 3): [0] force, [1] torque about the instance centre. The force is exactly conservative (sum = f V_fluid at steady state, to the solver residual) and as accurate as the flow solution it sustains: Zick & Homsy's drag within 0.7 % at N = 24 (tests/python/test_hydro_force_units.py). The torque carries the closed-form transposed-stress wall term; measured 2026-09-30 on a steadily spinning sphere it reads 1.019-1.027 of 8 pi mu a^3 Omega / (1 - phi) at R/h 6-12 (over-reads ~2 %, falling with resolution), and on a translating sphere whose true torque is zero it is zero to 1e-15 on a symmetric placement and
`init_mpi`	init_mpi(self, gnx: int, gny: int, gnz: int) -> None Wire the multi-rank step: pass the GLOBAL grid dims (gnx,gny,gnz). The Solver must have been constructed with this rank's LOCAL block dims (from mpi_block). MPI_Init is called if needed. Call it BEFORE the geometry (set_solid / set_pressure_geometry / set_solid_from_scene): the geometry builds the pressure and velocity multigrids for the partition in force when it runs, so init_mpi raises if a geometry already exists (built first, every rank would solve its own block-periodic pressure problem).
`instance_center`	instance_center(self, i: int) -> list[float] The resolved centre of rotation of instance i (world coordinates).
`max_open_divergence`	max_open_divergence(self) -> float Return the max cut-cell velocity-flux divergence max
`max_open_divergence_projected`	max_open_divergence_projected(self) -> float max
`max_porous_residual`	max_porous_residual(self) -> float Residual of the volume-averaged continuity max
`num_scene_instances`	num_scene_instances(self) -> int Number of instances in the scene set by set_scene (0 before it). Instance indices run 0..num_scene_instances()-1 in set_instance_transform / set_instance_motion, scene_owner() and the hydrodynamic force/torque accessors.
`origin`	The physical lower corner of the GLOBAL domain (0,0,0 by default).
`physical_units`	True when the solver was given an extent, i.e. lengths are physical rather than cells.
`pressure_bottom_extent`	The agglomerated-bottom extent threshold (see set_pressure_bottom_extent).
`pressure_telescope`	pressure_telescope(self) -> bool Whether the pressure multigrid's coarse levels telescope (merge ORB siblings onto fewer ranks once a level's block turns odd) rather than stop coarsening in place. Reads back set_pressure_telescope's current setting; ON by default since 2026-09-02 (the FoxBerry ladder), so this is normally True even though it was never called explicitly -- check it rather than assuming the legacy in-place-only hierarchy.
`rank`	rank(self) -> int This rank's index in MPI_COMM_WORLD.
`rebuild_geometry`	rebuild_geometry(self) -> None Re-derive all geometry from the current instance transforms. Velocity and pressure are PRESERVED across the rebuild (set_solid zeroes them by design). Full rebuild: the measured cost is ~65% momentum/IBM stencils, ~35% pressure/MG, scene sampling in the noise. Cells uncovered by the motion inherit zero, not an extrapolated fluid value.
`refresh_wall_velocity`	refresh_wall_velocity(self) -> None Re-derive ONLY the wall-velocity fields and the momentum operator that folds them in, for a driver that changes an instance's VELOCITY every step while its transform stays put -- the linearised oscillating body, a shear cell driven by counter-moving plates. set_instance_motion alone does not reach those fields (they are built in set_solid_from_scene), so the alternative is a full rebuild_geometry() per step to update a boundary condition. Refuses if a transform has changed since the last geometry build, because it does NOT re-sample the SDF, the apertures, the ownership field or the pressure operator. Velocity and pressure are untouched.
`reset_vof_bc_volumes`	reset_vof_bc_volumes(self) -> None Zero the per-face boundary liquid ledger (both the per-step and the running totals).
`set_advection`	set_advection(self, on: bool) -> None Enable/disable explicit high-order momentum advection (default scheme SOU). Off ⇒ Stokes.
`set_advection_scheme`	set_advection_scheme(self, scheme: str) -> None High-order momentum advection scheme: 'sou' (second-order upwind, the default) or 'koren' (Koren TVD).
`set_backflow_stabilization`	set_backflow_stabilization(self, beta: float) -> None Outflow backflow-stabilization coefficient (Bazilevs 2009 / Esmaily-Moghadam 2011): beta in [0,1] scales the dissipative outflow term that prevents backflow divergence when flow reverses at the outlet (e.g. a separated wake / BFS recirculation). Default 0.2; 0 = off. Inert where the outlet is purely outgoing.
`set_balanced_force_projection`	set_balanced_force_projection(self, enabled: bool) -> None Balanced-force projection (flow doc/collocated_varrho_forces.md §4.6), on both grids. Once per step, before the momentum predictor, the gradient part of the forces (set_body_force, the per-cell force fields, the CSF) is moved into the pressure with one extra solve of the projection's own operator, so static balances -- a hydrostatic column, a constant-curvature drop -- are exact from the first step at every viscosity, dt, density ratio and immersed-solid aperture. The solve does not depend on the flow state, so stability and the converged steady state are identical on and off; off leaves a decaying transient balance residue. get_p() stays the total physical pressure; its balanced part is the registered field 'p_balanced', which is RESTART STATE: registered as soon as the option is active, restoring it with 'p' (set_field) keeps the split, and omitting it (or switching the option on after steps with it off) re-splits P -- never a second copy of the balanced pressure. Cost: one extra pressure solve per step while the forces change, solved for the increment of the previous split and skipped when that is already converged -- a static interface costs no iterations. Default: ON on SolverColocated with variable density or surface tension, OFF otherwise; with an inflow/outflow face that default resolves OFF (a one-time stderr notice says so); an explicit call always wins. Refused (named error) with porous continuity, the ghost projection, set_fluid_only_constraint(2), the block CSF, inflow/outflow faces, the non-incremental pressure, and on the constant-density SolverColocated path (nothing to balance there).
`set_body_force`	set_body_force(self, force: collections.abc.Sequence[float]) -> None Set the body force per unit volume (fx, fy, fz) as one 3-sequence — e.g. a mean pressure gradient.
`set_capillary_cfl`	set_capillary_cfl(self, factor: float) -> None Safety factor on the capillary limit: step() raises when dt > factor * capillary_dt(). Default 1.0 — the Brackbill formula was measured to be the boundary itself, so it carries no built-in margin. Set it huge to disable the check, the same escape hatch set_vof_cfl_limit is for the Weymouth-Yue cap.
`set_collocated_scheme`	set_collocated_scheme(self, name: str) -> None Collocated cut-cell projection scheme (no effect on the staggered solver). DEFAULT since 2026-08-25: AUTO = 'ghost' where supported, falling back to 'gauge-exact' with a stderr notice on porous/variable-rho/domain-BC/Chebyshev configurations (ghost v1 limits); any explicit selection here disables AUTO. 'ghost' -- the fluid-only constraint scheme (binary-openness divergence + directional closures + the same gauge-exact gradient). Family-free, unconditionally stable (dt 60..1e20, no stabilizer), protocol-independent (C2), both-bed clean ladder -1.4% (R=8) -> +0.22% asymptote, Z&H anchor -0.018% at N=128. Costs: nonsymmetric BiCGStab pressure solve (~2.3-2.7x), ~1.6 KB/cell overlay (single-GPU size cap), fragmentation guard. MPI-capable (np=1,2,4 ctests); at-scale np>=16 hardening in progress. Equivalent to diagnostics.set_ghost_projection(True, 2, 2); the (1, 2) mixed mode stays quarantined there. 'gauge-exact' (the AUTO fallback) -- the aperture constraint with the directional gauge-exact pressure gradient. Cheapest scheme measured (symmetric MG-PCG, no fragmentation guard). CAVEATS from the 2026-08 campaign (doc/collocated_invisible_subspace.md): possesses an attractor FAMILY of steady states (support-inconsistent gradient/constraint pair) selected by march protocol, and a rotational-update instability whose wall-blend stopgap (diagnostics.set_rotational_wall_weight) has a resolution-dependent margin -- UNSTABLE at (R>=16 cells/radius, dt>=600); use fixed dt<=60 protocols at high resolution. Clean-protocol bias vs the staggered reference: -2.5% (R=8) -> ~+0.2% asymptote. 'embed' -- the complete Basilisk embed.h port (doc/history/collocated_embed_port_plan.md): the finite-volume momentum operator with the true-normal dirichlet_gradient wall drag applied as a defect correction on the IBM matrix, the openness-weighted -grad(P) predictor and cell correction, the wall-aware face map and the solid-cut-cell sliver mask. The live candidate for removing the collocated accuracy ceiling; its two intermediate rungs are diagnostics.set_face_interp(5
`set_contact_angle`	set_contact_angle(self, theta: float) -> None Rung V5b (WO-S). Prescribe a STATIC contact angle, in DEGREES, measured THROUGH THE LIQUID, on every immersed SDF wall. Needs set_solid(..., cutcell_pressure=True) + enable_vof; with no call the wall stays at the neutral 90-degree fill of rung V5a and every V5a number is byte-identical. What it changes is PASS 1 of the solid-band fill and nothing else: the band cells receive the volume fractions of the plane that continues the fluid-side interface into the solid at angle theta (m . n_w = cos theta, with n_w = grad(sdf)/
`set_contact_angle_dynamic`	set_contact_angle_dynamic(self, enabled: bool, theta_e: float = 90.0, slip_length_cells: float = 0.0, mu_liquid: float = 0.0, sigma: float = 0.0) -> None Rung V6 (WO-V6). The DYNAMIC contact angle -- `enabled=False` turns it (and the hysteresis) off so the static angle stands again; `enabled=True` imposes at the grid scale the Cox-Voinov apparent angle of a contact line moving at speed U_cl, with the outer cut-off at the CELL SIZE and an EXPLICIT slip length lambda (Afkhami, Zaleski & Bussmann, JCP 228:5370, 2009): theta_Delta^3 = theta_e^3 + 9 Ca_cl ln(Delta/lambda), Ca_cl = mu_l U_cl / sigma angles in radians internally, arguments in DEGREES; Ca_cl > 0 advancing, < 0 receding; theta_Delta is clamped into [1, 179] degrees (diagnostics.set_contact_angle_clamp). `slip_length_cells` is lambda/Delta and must lie in (0, 1). It is not a tuning knob you may omit from a report: a VoF contact line's NUMERICAL slip is proportional to the cell size, so without an explicit lambda the imposed angle is silently grid-dependent (VOF_PLAN section 6). ALWAYS state lambda alongside a dynamic-wetting result. `mu_liquid` is the liquid dynamic viscosity entering Ca_cl (solver units); `sigma` defaults to whatever set_surface_tension holds. theta_e also becomes the static base angle, so this call replaces set_contact_angle. U_cl is measured as -u.t_hat at the anchor fluid cell of each band cell, with t_hat the IN-WALL direction of the fluid-side PLIC normal (m points into the gas, so the liquid advances along -t_hat), then smoothed with a 3-point mean along t_hat. Needs set_solid(..., cutcell_pressure=True) + enable_vof; with no call every V5b number is byte-identical.
`set_contact_angle_field`	set_contact_angle_field(self, theta: ndarray[dtype=float64, order='F']) -> None Per-cell static contact angle in DEGREES, as a Fortran-order (nx,ny,nz) array. Only the value AT the solid band cell being filled is read, so cells away from a wall are irrelevant. theta is a field so the dynamic-angle rung (V6) changes only what fills it.
`set_contact_angle_hysteresis`	set_contact_angle_hysteresis(self, theta_a: float, theta_r: float) -> None Rung V6 (WO-V6). Advancing / receding contact-angle hysteresis, in DEGREES. While the measured apparent angle lies in [theta_r, theta_a] the contact line is PINNED and the fill imposes the APPARENT angle itself — which reproduces the current interface exactly, because the V5b fill is idempotent (WO-S finding 1), so no Young force appears and the line does not move. Above theta_a the advancing angle is imposed and below theta_r the receding one, each with the Cox-Voinov correction of set_contact_angle_dynamic when that is also configured. Sets the static base to (theta_a+theta_r)/2 if no angle was set yet.
`set_decomposition`	set_decomposition(self, levels: int, max_imbalance: float = 1.05) -> None Choose how this solver's shared MPI decomposition is built. levels=0 (DEFAULT) = the aligned ORB: split positions on the FINE grid are snapped to a power of two (capped at 16, i.e. 5 nested multigrid levels). levels >= 2 = COARSE-FIRST: build the ORB on the grid coarsened levels-1 times and refine the partition upward, so every block is a multiple of the coarsening factor BY CONSTRUCTION and the hierarchy nests for the full requested depth. Each candidate depth is built and MEASURED, and the deepest one whose max-block/min-block ratio stays within max_imbalance (default 1.05) is taken, else the aligned ORB. Call BEFORE init_mpi() (a later call raises), and pass the SAME levels / max_imbalance to flow.mpi_block() -- both derive the same partition and must agree.
`set_divergence_sink`	set_divergence_sink(self, weights: ndarray[dtype=float64, order='F']) -> None WO-P23: an AUTO-BALANCED sink region for the phase-change divergence source. Give a non-negative weight per cell (Fortran-order (nx,ny,nz)); after each deposit the solver subtracts (global sum of the phase-change source) * w / (global sum of w), so a CLOSED domain's Poisson RHS is exactly compatible every step with no user bookkeeping. Put the weights in the LIQUID far from the interface, where the exact solution simply has the liquid leaving. This is what makes a closed-box phase-change run possible without the variable-density outflow operator (whose density-ratio inconsistency is WO-R2's subject). Registered as the field 'div_sink'.
`set_divergence_source`	set_divergence_source(self, array: ndarray[dtype=float64, order='F']) -> None A PRESCRIBED extra divergence source (1/s) on the inner cells, added to the Poisson RHS beside the phase-change deposit, so the projection solves div(open u) = S. Registered as the field 'div_source'. Its use is to make a CLOSED (periodic) domain compatible with a net vapour production by carrying a balancing sink; with an outflow face the outflow does that job and this is unnecessary.
`set_domain_bc`	set_domain_bc(self, face: str, type: str, velocity: collections.abc.Sequence[float] = [0.0, 0.0, 0.0]) -> None Set a per-face domain BC. `face` is one of '-x', '+x', '-y', '+y', '-z', '+z'; `type` is 'periodic' (the default on every face), 'wall' (no-slip), 'inflow' (Dirichlet `velocity` -- an inlet, or a lid with a tangential velocity), 'outflow' (zero-gradient velocity, p = 0) or 'slip' (free-slip / symmetry plane: zero normal velocity, zero normal derivative of the tangential components, pressure Neumann like a wall; `velocity` is ignored). The TYPE must be set BEFORE the geometry (set_solid / set_pressure_geometry / set_solid_from_scene): it is folded into the operators when the geometry is built, and a call that would CHANGE it raises afterwards. A VALUE update on a face whose type is unchanged (ramping an inflow jet, a lid's tangential speed) is allowed at any time and takes effect the next step. A half-domain closed by a 'slip' face reproduces the full symmetric domain pointwise (tests/kokkos/test_freeslip). Both grids (Solver and SolverColocated).
`set_domain_bc_profile`	set_domain_bc_profile(self, face: str, profile: ndarray[dtype=float64, order='C']) -> None Prescribe a per-position inlet velocity profile (Nb,Nc,3) over a face (one of '-x', '+x', '-y', '+y', '-z', '+z'; sets it to 'inflow'). Like set_domain_bc, the TYPE (making a face inflow for the first time) must be set BEFORE the geometry; updating the profile on a face that is ALREADY inflow is allowed after the geometry.
`set_dt`	set_dt(self, dt: float) -> None Set the time step dt (the caller's time unit); the momentum solve is scaled by 1/dt (well-conditioned at large dt). Under a physical domain the FIRST call pins the reference time the internal scales are built on; any later call (step_adaptive makes them) rebuilds the momentum operator at the next step.
`set_field`	set_field(self, name: str, array: ndarray[dtype=float64, order='F']) -> None Write a Fortran-order (nx,ny,nz) float64 array into a registered field's inner region (ghosts refilled on the next diagnostics.exchange_field/step).
`set_implicit_advection`	set_implicit_advection(self, on: bool) -> None Use implicit-FOU advection with deferred-correction TVD.
`set_instance_motion`	set_instance_motion(self, instance: int, lin_vel: collections.abc.Sequence[float] = [0.0, 0.0, 0.0], ang_vel: collections.abc.Sequence[float] = [0.0, 0.0, 0.0], center: collections.abc.Sequence[float]
`set_instance_transform`	set_instance_transform(self, instance: int, translation: collections.abc.Sequence[float], quat: collections.abc.Sequence[float] = [0.0, 0.0, 0.0, 1.0]) -> None Move one scene instance (quat as (x,y,z,w)). Takes effect at the next rebuild_geometry(): the SDF, the cut-cell overlay, the apertures and the pressure operator are all derived from the transforms, so moving one without rebuilding would run the solver on stale geometry.
`set_mass_flux`	set_mass_flux(self, array: ndarray[dtype=float64, order='F']) -> None P0: prescribe a per-cell interfacial mass flux from a Fortran-order (nx,ny,nz) float64 array (== set_field('mdot') plus the ghost fill). Turns the thermal mass flux OFF.
`set_mass_flux_uniform`	set_mass_flux_uniform(self, mdot: float) -> None P0: prescribe a UNIFORM interfacial mass flux (kg m^-2 s^-1, solver units; positive = evaporation, liquid consumed). Turns the thermal mass flux OFF.
`set_mu`	set_mu(self, mu: float) -> None Set the dynamic viscosity mu (the caller's units); a change after the geometry rebuilds the momentum operator at the next step.
`set_phase_change_energy`	set_phase_change_energy(self, enabled: bool, rho_cp_gas: float = 0.0, rho_cp_liquid: float = 0.0) -> None CONSISTENT rhoc_pT transport (VOF_PLAN section 9 item 6) for the scalar named by set_phase_change_thermal; `enabled=False` returns to the constant-diffusivity scalar operator and the Koren TVD advective term (the rung P0/P1 energy path). Two things change together: (1) TRANSPORT: H = (rho c_p) T is advected with the colour advection's OWN geometric fluxes, sweep order and frozen dilation flag (vof/energy_advect.hpp), and T is recovered as H / (rho c_p)(C^{n+1}) -- so a uniform temperature is preserved EXACTLY at any heat capacity ratio. With two different fluxes the heat carried into a mixed cell is divided by a capacity built from another flux: an error of order d(rho c_p), i.e. ~2000x at water/steam, which is the artificial interfacial heating this removes. (2) DIFFUSION: the implicit operator becomes A_C = rho c_p(C)/dt + sum_f k_f open_f with k_f the arithmetic mean of the cells' k(C) (the k_gas/k_liquid of set_phase_change_thermal) EXCEPT at a face touching a per-cell Dirichlet (interfacial) cell, where the pure neighbour's own k is used -- a Dirichlet row is an identity row, so that coefficient's only job is the conductance with which the pure cell reaches a boundary condition that already sits at the interface. Units: rho*c_p in J/(cell^3 K), k in W/(cell K). Requires set_phase_change_thermal.
`set_phase_change_thermal`	set_phase_change_thermal(self, enabled: bool, scalar: str = '', t_sat: float = 0.0, k_gas: float = 0.0, k_liquid: float = 0.0, r_int: float = 0.0) -> None P1: compute mdot each step from a registered scalar (the temperature) instead of prescribing it; `enabled=False` stops that (back to the prescribed mdot field): mdot = ( k_g dT_g/dn - k_l dT_l/dn ) / h_lv , n = m/
`set_porous_continuity`	set_porous_continuity(self, on: bool = True) -> None Enable the volume-averaged (porous) continuity for unresolved CFD-DEM (staggered only): the projection enforces d(eps)/dt + div(eps u) = 0 instead of div(u)=0, so the fluid velocity is NOT solenoidal where the void fraction changes (bubbling/expansion). Binds the 'eps' field (void fraction from the particle deposition, created seeded to 1 if absent; the coupling writes it each step BEFORE step()). eps=1 everywhere reduces exactly to div(u)=0. Pair with max_porous_residual() for the meaningful convergence check.
`set_pressure_bottom`	set_pressure_bottom(self, mode: str) -> None Coarse-level (bottom) solve of the pressure multigrid. A V-cycle converges at a domain-independent rate only if its COARSEST level is effectively solved, and a geometric hierarchy cannot always get small enough: an axis stops coarsening once it turns odd, and under MPI once ANY rank's block turns odd -- so at fixed cells/rank the coarsest GLOBAL grid grows with the rank count and the bottom is progressively under-solved. 'auto' (DEFAULT) agglomerates the coarsest level onto a global operator and solves it exactly whenever it exceeds set_pressure_bottom_extent (4) cells on any axis, and uses the cheap smoothed bottom otherwise (byte-identical to 'smoother' then). 'smoother' = never agglomerate (legacy). 'agglomerated' = always. Measured on one GPU (2048x64x64 channel): a smoothed bottom needs 13.5 pressure iterations/step at 4 levels and 6.0 at 6, against 4.4 at full geometric depth; agglomerated it is 4.0 at BOTH depths, and faster in wall-clock than the deep hierarchy. Works on the cut-cell IBM and ghost-projection paths (per-fluid-component null-space projection, 2026-08-13).
`set_pressure_bottom_extent`	set_pressure_bottom_extent(self, cells: int) -> None The threshold set_pressure_bottom('auto') applies: agglomerate the coarsest level onto a global operator once the coarsest GLOBAL grid exceeds this many cells on any AXIS (not cell count -- what a few smoother sweeps cannot fix is a mode spanning many cells along one axis). DEFAULT 4.
`set_pressure_chebyshev`	set_pressure_chebyshev(self, on: bool, max_iter: int = 120, rtol: float = 1e-09) -> None Use the communication-light Chebyshev pressure accelerator. Mutually exclusive with the two CG drivers (on=True clears an FCG selection); on=False deselects Chebyshev and falls back to FCG if that is selected, otherwise to MG-PCG.
`set_pressure_fcg`	set_pressure_fcg(self, on: bool, max_iter: int = 200, rtol: float = 1e-08) -> None Use the FLEXIBLE MG-CG pressure accelerator: the same V-cycle-preconditioned CG as set_pressure_pcg, the same stopping estimate/mean removal, and the same cap+tolerance (they are shared -- it is the same solve), differing only in the beta recurrence (Polak-Ribiere r^T(z_{k+1}-z_k)/r^T z_k instead of Fletcher-Reeves). Costs one extra level-0 vector and one extra global dot per iteration; buys tolerance of a preconditioner that is not symmetric with respect to the fine operator. On a symmetric preconditioner the two betas are algebraically identical, so FCG reproduces PCG's iteration count. Unlike set_pressure_pcg this flag GENUINELY selects: on=True clears the Chebyshev selection (so it works after diagnostics.set_density_mode/set_porous), on=False returns to MG-PCG, and a later set_pressure_chebyshev(True, ...) wins over it.
`set_pressure_geometry`	set_pressure_geometry(self, sdf: ndarray[dtype=float64, order='F']) -> None Set an all-fluid SDF for the cut-cell pressure operator without an immersed solid (the channel/BFS domain-BC path). For a no-slip immersed BODY in an inflow/outflow domain, call set_solid(sdf, cutcell_pressure=True) instead -- do NOT also call this (a second geometry setter overwrites the SDF and wipes the solid).
`set_pressure_multigrid`	set_pressure_multigrid(self, on: bool, levels: int = 4) -> None Set the pressure multigrid depth (levels=1 => pure RB-GS, no coarse grid).
`set_pressure_pcg`	set_pressure_pcg(self, on: bool, max_iter: int = 200, rtol: float = 1e-08) -> None Use the MG-PCG pressure accelerator (single-GPU default) and set its iteration cap and relative tolerance. on=True GENUINELY selects: it clears both competing selections (Chebyshev and FCG), so it works after diagnostics.set_density_mode / set_porous -- last set wins. on=False RAISES: MG-PCG is the terminal fallback of the driver dispatch, so it cannot be deselected on its own -- select the driver you want instead (set_pressure_chebyshev(True, ...) or set_pressure_fcg(True, ...)). Under diagnostics.set_ghost_projection the nonsymmetric operator is solved by BiCGStab; this call still sets that solve's cap/tolerance (they are shared) and does not change its method.
`set_pressure_solver_params`	set_pressure_solver_params(self, iters: int) -> None Set the pressure smoother iteration count.
`set_pressure_telescope`	set_pressure_telescope(self, on: bool) -> None Coarse-level TELESCOPING of the pressure multigrid (multi-rank). The geometric hierarchy coarsens a level in place, which needs every rank's block even on that axis; once a block turns odd the hierarchy used to STOP, leaving a coarsest grid far too large to solve (384^3 on 1536 ranks: blocks die at 3x6x4 with the coarsest extent still 48, and the pressure iteration count grows 16.6 -> 38.7 across the ladder). With this on, a blocked level instead merges ORB siblings onto fewer ranks (the merged block's origin is its parent's split value, so parity is restored), moves the residual down / correction up within rank groups, and the roots keep coarsening on a sub-communicator to a tiny bottom. Idle ranks below the telescope point cost nothing that matters; the iteration count becomes a property of the problem rather than the rank count. ON by default since 2026-09-02 (the FoxBerry ladder); False restores the in-place-only hierarchy, byte-identical to before telescoping existed. See suite/docs/archive/MG_TELESCOPING_PLAN.md.
`set_property_model`	set_property_model(self, target: str, kind: str, field: str, params: collections.abc.Sequence[float] = [], field2: str = '') -> None Register a device closure writing a property/body-force field from input field(s). target: a registered field (a property 'mu'/'rho'/… or a body-force component 'force_x'/'force_y'/'force_z'). kind: 'linear' (params [p0,p1,p2]: p0+p1field+p2field2), 'boussinesq' (params [rho0,g,beta,T0]: rho0gbeta(field-T0) buoyancy), 'arrhenius' (params [mu_ref,B,Tref]: mu_refexp(B*(1/field-1/Tref))). Applied at the top of step().
`set_property_table`	set_property_table(self, target: str, field: str, x: collections.abc.Sequence[float], y: collections.abc.Sequence[float]) -> None Register a tabulated property: target = piecewise-linear interpolation of (x, y) at the input field value (x ascending, clamped at the ends).
`set_rho`	set_rho(self, rho: float) -> None Set the fluid density rho (the caller's units). Under a physical domain the FIRST call pins the reference density the internal scales are built on; a later change rebuilds the momentum operator at the next step.
`set_scalar_bc`	set_scalar_bc(self, name: str, face: str, type: str, value: float = 0.0) -> None Scalar boundary condition on a domain face ('-x', '+x', '-y', '+y', '-z', '+z'): type 'periodic', 'neumann' (zero flux, adiabatic) or 'dirichlet' (the given value). Single-rank.
`set_scene`	set_scene(self, node_ints: ndarray[dtype=int32, order='C'], node_reals: ndarray[dtype=float64, order='C'], instance_ints: ndarray[dtype=int32, order='C'], instance_reals: ndarray[dtype=float64, order='C'], periodic: bool = False) -> None Install an analytic geometry scene from core's flat encoding (3 ints + 16 reals per node, 2 ints + 18 reals per instance -- reals[17] is the centre-pinned flag; legacy 17-real records are accepted, reading an all-zero centre as 'follows the body'). Coordinates are the caller's own PHYSICAL ones when the solver was given an extent -- so ONE scene serves flow, dem and voro unchanged -- and CELL UNITS on the global inner grid otherwise (cell (i,j,k)'s centre at (i,j,k) -- half a cell from cell_centers(), which returns i + 1/2 in cell units). The scene is replicated on every rank, so scene-derived geometry needs no communication and -- unlike diagnostics.set_exact_crossings -- is NOT single-rank only. periodic=True treats the scene as min-image periodic over the global grid (one instance per body, no images); periodic=False leaves images to the caller.
`set_solid`	set_solid(self, sdf: ndarray[dtype=float64, order='F'], cutcell_pressure: bool = False) -> None Set the solid SDF as a Fortran-order (nx,ny,nz) float64 array sampled at cell_centers() (negative inside the solid, positive in fluid, in the solver's length unit). cutcell_pressure=True enables the open-face-weighted cut-cell pressure operator (proper no-slip); it composes with domain BCs, so this is the single call for a no-slip immersed body in an inflow/outflow domain. The multigrid coarse operators are always the rediscretized per-level cut-cell operators (the former `pressure_coarse=` mode selector is gone).
`set_solid_from_scene`	set_solid_from_scene(self, cutcell_pressure: bool = True) -> None Sample the installed scene onto this rank's inner grid and install it as the solid, entirely on device (no nxnynz float64 host round trip).
`set_state`	set_state(self, u: ndarray[dtype=float64, order='F'], v: ndarray[dtype=float64, order='F'], w: ndarray[dtype=float64, order='F']) -> None Upload an initial velocity field (u,v,w each a Fortran-order (nx,ny,nz) float64 array).
`set_superficial_velocity`	set_superficial_velocity(self, enabled: bool, axis: str = 'x', velocity: float = 0.0) -> None Hold the SUPERFICIAL velocity along a periodic axis at `velocity` (physical units): the volume mean of that velocity component, i.e. the volume flux of the MIXTURE through any cross-section divided by the cross-section's full area. In multiphase flow 'velocity' alone is ambiguous between superficial and interstitial (per-phase); this is the superficial one. At the end of every step() one uniform shift is added to that component on every face; on an axis with no boundary a uniform shift of the normal velocity leaves every cell divergence unchanged, so the field stays discretely divergence-free. `velocity = 0` is a closed (batch) column: the flow can circulate but not pass through (TBFsolver's constant flow rate at 0). The mean is a device reduction, global under MPI; single-rank the step makes no host round trip for it. Staggered Solver, all-fluid domain (set_pressure_geometry), periodic `axis` only -- step() raises otherwise. set_superficial_velocity(False) turns it off. diagnostics.last_superficial_velocity_shift() is the shift the last step applied.
`set_surface_tension`	set_surface_tension(self, sigma: float) -> None VoF rung V4: turn on the BALANCED-FORCE continuum surface force with coefficient sigma (0 = off, the default). Enables VoF and the curvature cascade, which then runs once per step at the head, from the same colour field the density closure sees. The force at the staggered velocity unknown u_c(i) is F = sigma * kappa_f * (C(i) - C(i - s_c)) / h i.e. the colour difference is the projection's OWN face difference — the identical discrete gradient operator the pressure uses (Francois et al. 2006; Popinet 2009). That is the whole content of 'balanced force': with a constant kappa the force is exactly a discrete gradient, so the projection annihilates it and a static drop stays at machine zero. Face-interpolating a cell-centred sigmakappagrad(C) instead — the obvious way to reuse the per-cell body force — is not a discrete gradient of anything and leaves spurious currents of order sigmakappa/mu that no curvature accuracy removes. UNITS: sigma is a force per unit LENGTH in the caller's own units -- the system rho, mu and set_body_force are stated in (so per CELL without a physical domain, where the cell size is 1). SIGN: kappa is positive for a convex blob of liquid and C is the liquid fraction, so the equilibrium pressure is P = sigmakappa*C + const — the Young-Laplace overpressure INSIDE the drop. Surface tension is EXPLICIT, so step() enforces the Brackbill capillary time step (see capillary_dt()).
`set_velocity`	set_velocity(self, c: int, array: ndarray[dtype=float64, order='F']) -> None Write component c's inner velocity from a Fortran-order (nx,ny,nz) float64 array (solid rows re-masked). Initial-condition hook: a uniform stream around a fixed body is the Galilean twin of the same body towed through fluid at rest.
`set_velocity_multigrid`	set_velocity_multigrid(self, on: bool, levels: int = 4, vcycles: int = 8) -> None Enable velocity (momentum) multigrid for the implicit diffusion solve.
`set_velocity_residual_tolerance`	set_velocity_residual_tolerance(self, rtol: float) -> None Residual-based momentum stop: a component's implicit solve ends once max
`set_vof`	set_vof(self, array: ndarray[dtype=float64, order='F']) -> None Set the colour field from a Fortran-order (nx,ny,nz) float64 array: C = the LIQUID volume fraction of the cell, in [0,1]. Enables VoF if it is not on yet. Initialise it SHARP (0/1 with the exact fraction in the interface cells) — geometric VoF keeps a sharp interface and a diffuse initial profile is simply a different problem.
`set_vof_backflow`	set_vof_backflow(self, face: str, value: float = 0.0) -> None inletOutlet colour on OUTFLOW face 'face' ('-x' ... '+z'; default 0 = gas): where the boundary face velocity points back INTO the domain, the colour ghost carries this value instead of the zero-gradient copy. This is OpenFOAM's inletOutletFvPatchField (Rusche 2002 thesis section 4), the standard VoF outlet — an outlet is a place where you know what leaves (whatever is inside) but must STATE what comes back. Zero-gradient on a reversing face re-injects whatever happens to be sitting at the outlet, which for a draining film feeds the film back into the domain. Where the fluid leaves, zero-gradient is kept.
`set_vof_block_assign`	set_vof_block_assign(self, mode: str, every: int = 0) -> None Master assignment of the block container (rung W1). 'round-robin' by block id (rung W0's, blind to block SIZE), 'lpt' = LONGEST-PROCESSING-TIME greedy on the block cell counts, 'orb' = weighted ORB over a 1-D block space (core's BlockDecomposer<1> with the cell counts as weights). 'every' > 0 re-runs the assignment every 'every' block steps and MIGRATES the colour of any block that changed master (nothing else in a block is state, so the migration is exact and the bitwise gates hold across it). All three are pure functions of the REPLICATED block table, so every rank computes the same assignment without an exchange -- which is what lets a re-assignment happen mid-run without breaking a bitwise gate. LPT is the measured winner (see the WO-W12 findings): the ORB's blocks must be CONTIGUOUS in block id, which LPT is free of.
`set_vof_cfl_limit`	set_vof_cfl_limit(self, value: float) -> None Weymouth-Yue boundedness cap on the interface-local Courant number. Default 0.25 — the PROVEN 3D bound 1/(2(N-1)) from Weymouth's thesis eq. A.33; the widely-quoted 0.5 is the 2D value. step() throws when it is exceeded. Conservation is INDEPENDENT of boundedness (the telescoping is algebraic), so an over-CFL run loses 0<=C<=1, not volume.
`set_vof_inflow`	set_vof_inflow(self, face: str, value: float) -> None Colour of the fluid ENTERING through inflow face 'face' ('-x', '+x', '-y', '+y', '-z', '+z'); 1 = liquid, 0 = gas. The face must already be an inflow (set_domain_bc(face, 'inflow', ...)) or this raises. The value may be FRACTIONAL, and it then means 'this fraction of the incoming flux is liquid' — a flux statement, not a sub-cell interface position. That is exactly how it is implemented: on a face whose donor is outside the domain the geometric flux is replaced by the algebraic C_donor*a (wyFaceFluxBc), because a uniform prescribed ghost band has no usable MYC normal and the inflow colour is boundary DATA. Setting this also makes the property ghosts of that face follow the inflow colour through the registered closures (rho_ghost = rho(C_inflow)) instead of copying the interior — without it the inlet FACE density is the interior's, wrong by up to the density ratio at a liquid inlet into a gas domain. Default (nothing set): the zero-gradient copy, i.e. today's behaviour, bit for bit.
`set_vof_inflow_profile`	set_vof_inflow_profile(self, face: str, profile: ndarray[dtype=float64, order='C']) -> None Per-position inflow colour (Nb,Nc) on the INNER grid of the face's two perpendicular axes, resampled with the same clamp rule set_domain_bc_profile uses for the velocity. This is how a liquid distributor over part of an inlet is expressed: the colour profile says WHERE liquid enters, the velocity profile says how fast.
`set_wall_slip_length`	set_wall_slip_length(self, lambda_cells: float) -> None Rung V6b (WO-V6b) - the VELOCITY half of the dynamic contact line. In the Robust-Scaled cut-cell IBM closure the TANGENTIAL wall datum stops being no-slip and becomes the NAVIER condition u_t(wall) - u_body = lambda * d(u_t)/dn with `lambda_cells` the slip length in the caller's own units (CELLS without a physical domain, where the cell size is 1); the model needs lambda < one cell. The wall-NORMAL component is untouched: the wall stays impermeable, with the moving-body velocity as its datum. 0 (the default) restores the validated no-slip closure BIT-IDENTICALLY - the closure polynomials are the lambda = 0 members of the same family, and the lambda > 0 branch is never entered. It is ONE lambda with set_contact_angle_dynamic: either call overwrites the shared value (last call wins), so the Cox-Voinov inner cut-off and the momentum wall closure cannot disagree. Only this call switches the MOMENTUM half on, so every WO-V6 result taken with the angle half alone is unchanged. Why it exists, measured: with no-slip the contact line advances only at the scheme's own numerical slip, ~1/180 of Lucas-Washburn (WO-V6 finding 6), and WO-V7's pore-scale campaign found every IMBIBITION verdict inverted below the capillary number at which the imposed velocity crosses that numerical slip velocity. Typical values are 0.01-0.5 cells; state lambda with every dynamic-wetting result. Staggered grids only; needs set_solid(..., cutcell_pressure=True). Changing it rebuilds the three cut-cell overlays and the momentum operator in place (the velocity field is NOT reset). APPROXIMATIONS, both stated because they are invisible from the call: the tangential projector is taken DIAGONAL (the cross terms -n_c n_j u_j are dropped; they vanish exactly for an axis-aligned wall and are O(n_c n_j) otherwise), and an axis whose BOTH neighbours are solid (a one-cell fluid gap) keeps the no-slip closure - counted by diagnostics.wall_slip_sandwich_cells().
`size`	size(self) -> int The number of ranks in MPI_COMM_WORLD.
`spacing`	The cell size [dx, dy, dz] = extent/cells as three floats; [1.0, 1.0, 1.0] without a physical domain (the cell-unit default). Three spacings within 1e-12 relative are snapped to one value, so an isotropic box is exactly isotropic.
`step`	step(self) -> None Advance the solver one time step (semi-implicit: diffusion + projection). ATOMIC across the two explicit two-phase stability caps: the Weymouth-Yue boundedness cap and the Brackbill capillary limit are evaluated at the HEAD of the call, so a rejected dt leaves every field bitwise as it was and a retry at a smaller dt is exact. (Before 2026-09-04 the momentum half had already advanced by the rejected dt, so catch-and-halve desynchronised colour and momentum.) The WY pre-check reads the field the call starts from; where the projection itself accelerates the field inside the step, use step_adaptive().
`step_adaptive`	step_adaptive(self, cfl_target: float = 0.4, capillary_cfl: float = 0.4, dt_max: float = inf) -> float One step at a dt re-picked from the CURRENT state; returns the dt used. dt = min(cfl_target * cfl_dt, capillary_cfl * capillary_dt, dt_max) with both limits taken from vof_step_limits() at the moment of the call, then set_dt(dt) + step(). Needs enable_vof(). WHY: both limits are INSTANTANEOUS. Re-picking every ten steps is safe for a free bubble or a falling film and is NOT safe in a packing -- a throat jet can grow max
`surface_tension`	surface_tension(self) -> float The surface-tension coefficient (0 when the CSF is off).
`sync_porous_prev`	sync_porous_prev(self) -> None Reseed eps^n = eps^{n+1} (d(eps)/dt = 0 for the next step). Call it once the first void-fraction field has been written into 'eps' and before the first step, so that step sees no spurious d(eps)/dt source; any later call simply restarts the d(eps)/dt bookkeeping from the current field.
`unit_scales`	The factors between the caller's units and the ones the solver computes in. The solver keeps computing on the UNIT LATTICE and folds the metric into constants at its API boundary, so every array reached through the RAW field registry (diagnostics.field_view / get_field / set_field) is in those INTERNAL units, unlike get_u/get_p, which convert. `identity` is True when no extent was given, and then every factor is exactly 1.0 and internal == physical. Multiply a physical value by *_to_internal, and an internal one by force_to_physical (a total force) to come back.
`update_properties`	update_properties(self) -> None Apply all registered property/force closures now (also done at the top of step()).
`velocity_residual_tolerance`	velocity_residual_tolerance(self) -> float The momentum residual tolerance in force (resolves the follow-the-pressure default).
`vof_bc_volumes`	vof_bc_volumes(self) -> list[float] Signed liquid volume that crossed each of the six domain faces during the LAST colour advection, in cell-volume units, POSITIVE for liquid ENTERING the domain. Local to this rank. Six entries in face order (-x,+x,-y,+y,-z,+z); only faces carrying a domain BC are meaningful.
`vof_bc_volumes_total`	vof_bc_volumes_total(self) -> list[float] The same, accumulated since enable_vof() or the last reset_vof_bc_volumes().
`vof_block_assign`	vof_block_assign(self) -> str The current master-assignment mode ('round-robin'
`vof_block_color`	vof_block_color(self, id: int) -> numpy.ndarray[dtype=float64] One marker's OWN inner color as a Fortran-order (nx,ny,nz) float64 array over its block box (diagnostics.vof_block_stats()['lo'/'hi']); empty on a rank that does not master it. This is the block's ONLY state, so {box, color} per marker is a COMPLETE checkpoint of the container -- and it is the only exact one: re-seeding from the union color field (enable_vof_blocks_from_field) hands each of two TOUCHING markers a slice of the other, because that gather is a copy of the union clipped to the seed extent. enable_vof_blocks_from_colors(boxes, colors) restarts the container from these.
`vof_blocks_enabled`	vof_blocks_enabled(self) -> bool True while the block container exists (enable_vof_blocks with seeds; None drops it).
`vof_cfl_limit`	vof_cfl_limit(self) -> float The Weymouth-Yue boundedness cap currently in force.
`vof_curvature`	vof_curvature(self) -> numpy.ndarray[dtype=float64] The curvature field's inner region as a Fortran-order (nx,ny,nz) float64 array, in 1/h in 1/LENGTH. Call compute_vof_curvature() first; get_field('kappa') is the RAW solver-internal field, which is 1/h.
`vof_curvature_branch`	vof_curvature_branch(self) -> numpy.ndarray[dtype=float64] Which cascade branch produced each cell's curvature, as a Fortran-order (nx,ny,nz) float64 array: 0 not interfacial, 1 height function, 2 height function in a non-preferred direction, 3 mixed height-position fit, 4 PLIC-volumetric paraboloid fit, 5 the same with the rank-deficient 3-parameter model, 6 NO estimate. == get_field('kappa_branch').
`vof_has_geometry`	vof_has_geometry(self) -> bool True when the colour advection is running the CUT-CELL (openness-weighted) kernels, i.e. an immersed solid is present and set_solid ran with cutcell_pressure=True. False means the uncut rung-V1 kernels are running, byte-identically to a solid-free build.
`vof_interface_area`	vof_interface_area(self) -> float Total interfacial area of the colour field, in CELLS squared (h^2), summed over the inner region and globally reduced under MPI. Uses the geometry diagnostics.set_phase_change_area selects, so the number a page quotes and the number the phase change integrates are the same one. Needs enable_vof; phase change need not be on. A sphere of radius R cells reads 4 pi R^2 to about 0.2-0.8 % once its colour field is resolved (WO-P3c).
`vof_max_courant`	vof_max_courant(self) -> float INTERFACE-LOCAL Courant number max
`vof_momentum_enabled`	vof_momentum_enabled(self) -> bool Whether momentum-consistent transport is on.
`vof_step_limits`	vof_step_limits(self) -> dict Both explicit two-phase step limits at the current state, and which one binds: the largest dt the interface-local Weymouth-Yue CFL cap admits ('cfl_dt'), the Brackbill capillary limit ('capillary_dt'), and min(cfl_dt, capillary_cfl*capillary_dt) with a 'capillary_binds' flag. At pore-scale capillary numbers the capillary limit is expected to bind first.
`wall_slip_length`	wall_slip_length(self) -> float The Navier slip length in force, in the caller's units (0 = no-slip).

`SolverDiagnostics`¶

Developer tier of a Solver, reached as s.diagnostics: instruments (censuses, timers, budgets, ledgers), solver TUNING beyond the public driver selection, ablation switches and the raw zero-copy / MPI internals. Nothing here is needed to set up, run or read out a simulation. It holds a reference to the solver and no state of its own; the setters that CHANGE RESULTS say so in their docstring.

Method / property	Description
`advection_wall_velocity`	Whether the wall-velocity advection inputs are in force.
`anderson_accelerator`	anderson_accelerator(self, window: int = 5, mixing: float = 1.0) -> peclet.flow._flow._AndersonAccelerator A new Anderson accelerator of this solver's steady march (doc/steady_acceleration.md; the public driver is peclet.flow.march_to_steady). window in [1, 8] (default 5), mixing in (0, 1] (default 1.0). Raises on an unsupported configuration (a collocated scheme other than 'ghost', VoF, phase change, scalars, porous, variable properties, drag, cell forces, moving scenes, superficial velocity, pressure warm start, the balanced-force projection) and when the history, (2m+3) * n_fields * 8 * n_padded bytes, does not fit. Collective under MPI.
`balanced_force_failures`	balanced_force_failures(self) -> int How many balanced-force solves have failed over the solver's lifetime (see last_balanced_force_failed).
`bcast_from_root`	bcast_from_root(self, value: object) -> object Broadcast a value from rank 0 (identity in the single-rank module; mirrors the MPI API).
`comm_avoiding`	The communication-avoiding smoothing mode in force (see set_comm_avoiding).
`contact_angle_diagnostics`	contact_angle_diagnostics(self) -> dict The solid-band census of the CURRENT colour field on this rank: how many band cells each branch of pass 1 wrote ('contact_cells' the theta plane of the cell's own anchor, 'neighbour_cells' the mean of the anchor's MIXED neighbours' planes where the anchor itself is pure phase, 'pure_cells' the pure-phase continuation C_s = C_f, 'parallel_cells' an interface parallel to the wall where no rotation is defined, 'neutral_cells' the WO-Q neutral-mean fallback, 'unfilled_cells' left untouched), plus 'mean_apparent_angle' — the mean angle the fluid-only normal reported at the contact cells BEFORE the rotation, in degrees. That last number is the direct read-out of how far the fluid-side interface still is from the prescribed angle, measured on the fill's own data rather than on a post-processed shape. Rung V6 (WO-V6) adds, all zero unless set_contact_angle_dynamic / set_contact_angle_hysteresis is configured: 'dynamic_cells' (band cells the V6 pass produced an angle for), 'pinned_cells' / 'advancing_cells' / 'receding_cells' (the hysteresis branch census), 'mean_imposed_theta' and 'mean_apparent_theta' in degrees, 'max_Ca_cl' = max
`csf_diagnostics`	csf_diagnostics(self) -> dict Census of the CSF face force on THIS RANK: per component the max
`exchange_field`	exchange_field(self, name: str) -> None Fill a registered field's ghost cells (cross-rank + periodic under MPI; periodic single-rank).
`exchange_field_add`	exchange_field_add(self, name: str) -> None Add-reduce halo: fold ghost-layer deposits back onto their owner (cross-rank + periodic). The particle->grid deposition primitive for MPI CFD-DEM; single-rank non-periodic no-op.
`field_view`	field_view(self, name: str) -> numpy.ndarray[dtype=float64] Zero-copy view of a registered field's full padded buffer as a Fortran-order (nx+2g, ny+2g, nz+2g) array (g = ghost_width); host → NumPy, device → DLPack (CuPy).
`fluid_momentum_cells`	fluid_momentum_cells(self) -> list[int] Unmasked (fluid) staggered momentum cells per component -- the exact discrete datum of the reaction identity: at steady state, sum over bodies of F_c = f_c * N_c.
`get_cut_owner`	get_cut_owner(self) -> numpy.ndarray[dtype=int32] Per-inner-cell owning scene instance (the argmin behind the sampled SDF), (nx,ny,nz) int32; -1 before set_solid_from_scene. Moving geometry reads its wall velocity off the owner and resolved CFD-DEM posts the hydrodynamic force back to it.
`hydro_force_torque_traction`	hydro_force_torque_traction(self) -> numpy.ndarray[dtype=float64] DIAGNOSTIC -- not the force on the bodies; that is Solver.hydro_force_torque_reaction(). The reconstructed surface-TRACTION loads, shape (4, n_instances, 3): [0] force, [1] torque about the instance centre, [2] the pressure part, [3] the viscous part ([0] == [2] + [3]). The cut-cell surface integral of (-p I + mu(grad u + grad u^T)) against the exact aperture wall-area vector, with a central-difference velocity gradient that spans the wall. It UNDER-READS THE VISCOUS PART: traction / reaction drag = 0.685-0.730 over phi 0.008-0.45 and N 24-128 (measured 2026-09-30), i.e. ~30 % low at ANY resolution -- refinement does not fix it. The torque is low by the same mechanism: a steadily spinning sphere reads 0.59-0.72 of 8 pi mu a^3 Omega / (1 - phi) at R/h 6-12 (the reaction: 1.019-1.027). What it offers that the reaction does not: the pressure/viscous split, and it runs where the reaction refuses (the collocated solver, porous, variable properties, domain BCs, implicit advection). A consistent traction (from the momentum operator's small-cell-robust wall reconstruction) is planned. Was Solver.hydro_force_torque() (deprecated). Atomics: tolerance-reproducible, not bitwise.
`instance_center_pinned`	instance_center_pinned(self, i: int) -> bool True if instance i's centre of rotation is PINNED (an explicit finite centre in the encoding, or set_instance_motion(center=...)); False if it follows the body's translation (NaN in the encoding, the builder's default; legacy all-zero raw arrays).
`last_balanced_force_failed`	last_balanced_force_failed(self) -> bool Whether the last step's balanced-force solve FAILED (driver breakdown or a non-finite result). A failed solve keeps P and p_balanced unchanged for that step; the first one warns on stderr.
`last_balanced_force_iterations`	last_balanced_force_iterations(self) -> int Pressure-driver iterations of the last step's balanced-force solve (set_balanced_force_projection); 0 when it is off, the force divergence was zero, or the previous split already met the stop.
`last_momentum_residual`	last_momentum_residual(self) -> float max over components of max
`last_outer_iterations`	last_outer_iterations(self) -> int Return the outer-iteration count from the last step().
`last_pressure_iterations`	last_pressure_iterations(self) -> int Return the pressure-solver iteration count from the last step(). A solve that BROKE DOWN (non-finite preconditioner output) reports the iteration CAP, so the usual rule-3b 'a capped pressure solve invalidates the run' check sees it; pressure_solve_failed() distinguishes the two.
`last_step_timers`	last_step_timers(self) -> dict Per-phase wall times (seconds, THIS rank, device-fenced) of the last step(): 'predictor' (ghost fills + RHS/advection/stencil builds), 'momentum' (implicit-diffusion solves), 'projection' (cut-cell pressure projection), 'step' (whole step). 'pressure_allreduce'/'pressure_allreduce_count' = time in / number of global reductions (MPI_Allreduce) inside the pressure solve — the latency-bound term of the distributed solve (0 on a single rank).
`last_superficial_velocity_shift`	last_superficial_velocity_shift(self) -> float The uniform velocity shift set_superficial_velocity applied at the end of the last step() (physical velocity units; 0 when it is off). Reading it synchronises with the device.
`moving_instance_cut_cells`	moving_instance_cut_cells(self) -> numpy.ndarray[dtype=float64] Per instance: inner cells on this rank it owns that touch a fractional face aperture, recounted whenever any instance moves. ZERO for a moving instance means its surface sits on grid planes (or is sub-cell) and its wall velocity is silently inert (set_solid_from_scene warns). Empty when no instance moves.
`moving_instance_degenerate_points`	moving_instance_degenerate_points(self) -> numpy.ndarray[dtype=float64] Per instance: staggered velocity points where the sampled SDF is EXACTLY zero -- a face on a lattice plane. Those points are fluid to the mask and not ghosts to the cut-cell fold, so a moving body's datum never enters there (set_solid_from_scene warns). Empty when no instance moves.
`multigrid_aspect_threshold`	The anisotropic-coarsening threshold in force (see set_multigrid_aspect_threshold).
`outflow_backflow`	outflow_backflow(self) -> dict Outflow REVERSAL census over the outflow faces (collective under MPI): 'max_reverse' = the largest reversed normal velocity max(0, -u.n), 'fraction' = reversed / all outlet faces, 'energy_influx' = sum over the reversed faces of rho
`outflow_rho_correction`	outflow_rho_correction(self) -> bool Whether the 1/rho_f factor is applied to the outflow face correction (WO-R item 4, reversed by WO-R2's operator fix; default True).
`periodic_image_overlap_cells`	periodic_image_overlap_cells(self) -> int Cells on this rank whose solid/fluid sign was decided by a periodic IMAGE of an instance wider than the box (set_solid_from_scene warns when nonzero): the scene evaluates the UNION of images, so a slab wider than the box refills any cavity carved from it. 0 when no instance is that wide or the images agree.
`phase_change_area`	phase_change_area(self) -> str The set_phase_change_area mode in force, as its name.
`phase_change_budget`	phase_change_budget(self) -> dict WO-P3f: the ENERGY BUDGET of the LAST energy solve (units: J for the enthalpies, W for `q_gfm`; cell volume = 1). h_open / h_open_new sum rho c_p (T - T_sat) over the UNMASKED cells, before / after the solve h_liquid, h_masked the same over the pure-liquid and the interfacial cells d_overwrite sum rho c_p (dval - T) over masked cells: what the Dirichlet rows inject when they overwrite the transported temperature (negative = destroyed) d_overwrite_new the part of it on cells that were NOT masked last step e_enter / e_leave the enthalpy carried OUT of / INTO the solved set by the cells that changed class this step q_gfm the heat the plane-anchored rows deliver INTO the unmasked set (negative while a bubble grows); the energy equation's own interfacial flux, to be compared with the regression's `mdot h_lv A_Gamma` n_enter_/n_leave_ the class-change census The discrete balance the entries close over one solve is `sum rho c_p (T^{n+1} - T*)/dt = q_gfm + (domain boundary flux) + (solve residual)`.
`phase_change_carry_ledger`	phase_change_carry_ledger(self) -> dict WO-P3f: the last `set_phase_change_carry_conserve` pass — the enthalpy actually handed back ('deposited') and the enthalpy of the interfacial cells that had NO neighbour left in the solve ('lost', which stays destroyed). Both 0 when the option is off.
`phase_change_diagnostics`	phase_change_diagnostics(self) -> dict Per-rank phase-change census of the LAST step: mdot extrema/mean and the interfacial cell count, the total PLIC interface area, the liquid volume the regression removed, the clip-and-redistribute ledger (
`pressure_coarse_ghost`	Whether the Neumann coarse ghost is in force (see set_pressure_coarse_ghost).
`pressure_exact_residual`	pressure_exact_residual(self) -> bool Whether the exact level-0 operator apply is in force (see set_pressure_exact_residual).
`pressure_mg_level_ratios`	pressure_mg_level_ratios(self) -> list[list[int]] The pressure multigrid's per-level coarsening ratio, one (rx, ry, rz) per level (read-only). On an isotropic domain every coarsenable axis halves, as it always has. On an ANISOTROPIC one (extent giving different spacings per axis) the hierarchy coarsens the FINEST coarsenable axis and defers one that is already at least set_multigrid_aspect_threshold (default 2) times coarser, so the coarse operators stay close to isotropic and the point smoother keeps its rate -- e.g. cells (N, 2N, N/2) on a cube gives (1,2,1), (2,2,1), (2,2,2), ... Empty until set_solid/set_pressure_geometry has built the operator. See flow/doc/anisotropic_metric.md §5.
`pressure_solve_failed`	pressure_solve_failed(self) -> bool Did the last pressure solve break down on a non-finite recurrence scalar (a preconditioner or operator that produced NaN/Inf)? The driver then zeroes the correction and continues -- so the projection was handed NOTHING for that step and the run is invalid. Until 2026-09-04 the only trace was one line on stdout ('preconditioner produced non-finite z') while last_pressure_iterations() reported 0, i.e. a perfectly healthy-looking solve (measured: examples/pore-scale-imbibition, a 3.1-cell throat). Set set_pressure_strict(True) to raise instead.
`pressure_strict`	Whether a non-finite pressure preconditioner output raises (see set_pressure_strict).
`reaction_budget_terms`	reaction_budget_terms(self) -> numpy.ndarray[dtype=float64] Diagnostic decomposition of the reaction identity, shape (2, 3): [0] the unsteady sum sum_i (rho/dt)(u_i - u^n_i) and [1] the advective sum sum_i A_i, both over the FLUID momentum cells of each component. The full discrete identity is sum_bodies F_c = f_c*N_c + sum fb + [1]_c - [0]_c; the Stokes form drops both because they vanish at steady state with advection off. sum_i A_i is the advection operator's net momentum flux through the cut walls -- an O(h) property of that operator, not of the budget.
`rebalance_by_weights`	rebalance_by_weights(self, weights: collections.abc.Sequence[float]) -> int Dynamic load balancing: redistribute the solver's state onto the weighted ORB of per-cell weights (global x-fastest, gnxgnygnz), its split planes aligned to 2^a cells with a the largest alignment whose weight imbalance stays within 1.05, so the pressure multigrid keeps coarsening in place for a levels. Returns that alignment (1 = the plain weighted ORB). Pass fluid work + gamma*particle_count; a coupled dem migrates onto the SAME partition with migrate_to_weights(weights, align=). State-preserving (bit-exact at np=1, reduction floor at np>1).
`reset_vof_timing`	reset_vof_timing(self) -> None Zero the VoF stage timers and the step counter without disarming them.
`set_advection_wall_velocity`	set_advection_wall_velocity(self, on: bool) -> None A0: fill the advection inputs' masked (solid) rows with the instantaneous WALL velocity instead of zeros, so a MOVING body's advective term sees the body's own motion. DEFAULT True; False is the pre-A0 ablation (measurement only). Inert on a static scene and on the collocated grid.
`set_aperture_order`	set_aperture_order(self, order: int) -> None Face-aperture estimator (DEFAULT 2 since 2026-08-26): 1 = legacy one-sample linear model, 2 = marching-squares (5 samples/face, O(h^2); removes the convexity bias measured at +0.59/+0.27% bed permeability at R=8/12 -- see doc/collocated_paper_plan.md row 51). For analytic geometry, exact apertures via set_openness_override are better still. Call before set_solid (a later call raises).
`set_comm_avoiding`	set_comm_avoiding(self, mode: str) -> None Communication-avoiding red-black smoothing (multi-rank): exchange a 2-deep ghost layer once per red-black PAIR instead of a 1-deep layer before every colour, redundantly re-smoothing the first colour's ghost ring so the second colour reads exactly what a per-colour exchange would have delivered. BIT-IDENTICAL, at half the halo events (measured: np=32 weak efficiency 35% -> 62%). 'both' (DEFAULT)
`set_contact_angle_clamp`	set_contact_angle_clamp(self, lo: float, hi: float) -> None The clamp on the Cox-Voinov cube, in degrees (default 1 / 179). The cubic has no solution beyond the maximum receding capillary number (film entrainment) and the fill's plane construction degenerates at 0/180.
`set_contact_angle_pivot`	set_contact_angle_pivot(self, mode: str) -> None ABLATION -- how the theta-plane is anchored in the fluid cell. 'volume' (DEFAULT) match the anchor cell's liquid volume with plicAlpha; 'centroid' pass through the PLIC centroid p_f (the Afkhami-Bussmann / Basilisk contact.h rule); 'projected-centroid' the work order's c = p_f - sdf(p_f) n_w; 'contact-line' the contact line on the wall. 'volume', 'centroid' and 'contact-line' are exactly idempotent (1e-15 on gate G0a); 'projected-centroid' is NOT -- projecting the centroid along n_w shifts the plane by -sdf(p_f) cos(theta), measured 0.26 in cell fraction at theta = 60, with the wrong sign (it removes liquid from the band for a wetting angle).
`set_contact_angle_smoothing`	set_contact_angle_smoothing(self, on: bool) -> None ABLATION - the 3-point in-wall mean of U_cl (default ON). A MAC velocity next to a wall is noisy cell to cell and the cube root of the Cox-Voinov relation puts that noise straight into the imposed angle; off measures what the smoothing is worth.
`set_csf_mode`	set_csf_mode(self, mode: str) -> None ABLATION. 'face' (default, the only production mode) evaluates the surface-tension force as sigmakappa_f(C(i)-C(i-s))/h at the face -- the projection's own gradient operator. 'cell' evaluates a CELL-CENTRED sigmakappagrad(C) and face-interpolates it with the arithmetic mean, exactly as the per-cell body-force machinery carries a rho*g field: consistent, convergent, and wrong for surface tension, because the result is not in the range of the operator the projection inverts. Shipped so the difference is a measured number (see the vof_surface_tension ctest), not an argument.
`set_deferred_correction`	set_deferred_correction(self, on: bool) -> None Deferred-correction advection: True (default) = 2nd order (implicit FOU + explicit high-order correction, the high-order scheme being SOU by default or Koren TVD via set_advection_scheme); False = pure implicit FOU (1st-order upwind, more dissipative but unconditionally stable at sharp shear layers).
`set_density_mode`	set_density_mode(self, mode: str = 'variable') -> None Enable variable density: binds the 'rho' field (get/set_field('rho'), created seeded with set_rho's value if absent) into the momentum time term, the advection weight, the per-cell body force (face-interpolated), and the pressure projection (face coefficient opennessrho0/rho_f with the matching 1/rho_f velocity correction; rho0 = set_rho's value, so a uniform field reduces exactly to the constant solver). A closure targeting 'rho' (e.g. a linear mixture of a transported phase fraction) enables this automatically. For gravity, register a closure force_z = linear(rho, params=[0, -g]). COLLOCATED (SolverColocated, rung V8): supported since 2026-09-02, ALL-FLUID only (set_pressure_geometry; an immersed solid, the ghost projection and set_rho_face_harmonic throw). The face coefficient and the face correction are the same as on the staggered grid, applied to the ABC projection's MAC face field, which is the MOMENTUM-weighted map (rho_L u_L + rho_R u_R)/(rho_L + rho_R) of the cell velocities; the CELL correction is the average of the two face corrections of each axis (never a cell-centred grad(phi)/rho_c). The pressure, the CSF and set_body_force (a mean pressure gradient) enter the implicit momentum predictor as the matching face integral rho_c avg_faces((f - grad_f P)/rho_f); a per-cell force enters at the cell value. The pair is adjoint, so the step is stable at every dt and density ratio and its steady state is dt-independent (doc/collocated_varrho_forces.md). Rated to density ratio ~10 for cases WITH MOTION (a translating drop exceeds the VoF CFL cap at ratio 100 and 1000; momentum consistency needs Favre face states and is not in this rung); exact at rest at ratio 1000.
`set_exact_crossings`	set_exact_crossings(self, t: collections.abc.Sequence[float]) -> None Analytic-SDF capability: exact wall-crossing fractions overriding the linear-interpolated theta in the momentum cut-cell overlay AND the ghost-projection closures. Flat array of 9nxnynz values, blocks [(c3+k)]: component c's staggered point at inner cell i toward its +k neighbour; NaN = no crossing. Call BEFORE set_solid (a later call raises); empty list clears. Single-rank, staggered momentum placement.
`set_exact_crossings_from_scene`	set_exact_crossings_from_scene(self) -> None Compute EXACT wall crossings from the installed scene, on device, on every rank -- the in-solver replacement for set_exact_crossings + scripts/exact_apertures_spheres.py. Bisection, not Newton: only SIGN correctness is guaranteed for bound-only leaves.
`set_face_interp`	set_face_interp(self, mode: int) -> None The two INTERMEDIATE rungs of the Basilisk embed.h port on the collocated grid (doc/history/collocated_embed_port_plan.md), reachable for the ladder that measures what each piece of the complete port ('embed' = mode 7) is worth: 5 = the FV momentum operator with the true-normal dirichlet_gradient wall drag as a defect correction, on the wall-aware face map with the transpose-gradient pressure force; 6 = plus the openness-weighted -grad(P) predictor and cell correction, on the plain face map (7 adds the wall-aware face map and the solid-cut-cell sliver mask). CHANGES RESULTS; no effect on the staggered solver. Modes 1-4 and 10-13 (ablations) were deleted with their kernels at 1.0.0.
`set_fluid_only_constraint`	set_fluid_only_constraint(self, mode: str) -> None Fluid-only pressure constraint (route 2b, a collocated mechanism instrument; call before set_solid): 'filter' = Design A (openness filter), 'star' = Design B (SPD Kron star elimination), 'off'. Any setting other than 'off' disables the AUTO scheme.
`set_fresh_cell_seed`	set_fresh_cell_seed(self, on: bool) -> None Fresh-cell policy for MOVING geometry. A point a body has just uncovered inherits whatever the solid held there -- zero, or a stale masked value -- which is one of the two textbook mechanisms behind spurious force oscillations in a moving-boundary IBM. True seeds it with the LOCAL WALL VELOCITY instead, so it starts moving with the surface that released it. Bounded (no extrapolation), no new field, and exactly the old behaviour when nothing moves. DEFAULT TRUE since 2026-08-30: measured on an oscillating sphere translating through the grid, it removes a resolution-INDEPENDENT +2.6..2.9% drag bias, cuts the spurious force oscillation 20-50x to within 17% of the non-moving floor, and improves the resolved CFD-DEM loop's total-momentum conservation 95x. Pass False for the pre-2026-08-30 behaviour.
`set_ghost_projection`	set_ghost_projection(self, on: bool, matrix_order: int = 2, rhs_order: int = 2) -> None QUARANTINED 2026-08-18 (verification only, unsupported): superseded by the gauge-exact collocated scheme, which matches its accuracy at 5-6x lower cost. Kept as the independent second discretization behind the cross-IBM physics gate. Enabling it on the collocated grid silently selects the plain face map, since it owns the operators the gauge-exact scheme replaces. EXPERIMENTAL directional ghost-cell projection (second staggered IBM): point-based FD divergence with wall-anchored directional closures instead of the openness-weighted cut-cell projection; solved by MG-preconditioned BiCGStab. Call BEFORE set_solid (a later call raises). Closure orders (1=linear, 2=quadratic): (matrix_order, rhs_order) = (2,2) full quadratic 13-point matrix; (1,1) linear 7-point; (1,2) mixed/deferred — 2nd-order steady constraint on a 7-point matrix. Collocated: the same closures/matrix on the face-averaged field, plus a directional (one-sided 2nd-order) cell gradient for the -grad(P) predictor and cell correction; requires face_interp 0. v1: single-rank, periodic + IBM, stationary walls; incompatible with porous/variable-rho/domain-BC/Chebyshev.
`set_incremental_pressure`	set_incremental_pressure(self, on: bool) -> None Toggle the rotational incremental-pressure projection.
`set_multigrid_aspect_threshold`	set_multigrid_aspect_threshold(self, theta: float) -> None Anisotropic-coarsening aspect threshold theta shared by the pressure and velocity multigrids (doc/anisotropic_metric.md §5.1): an axis is coarsened only while its cell size stays within theta of the smallest one, so a stretched grid semi-coarsens instead of building levels whose operator is dominated by one direction. DEFAULT 2.0; must be > 1; 1e9 restores full coarsening on every axis (the §8.5 item (c) ablation). Only ever consulted on an ANISOTROPIC metric, so it is inert on cubic cells.
`set_openness_override`	set_openness_override(self, ox: collections.abc.Sequence[float], oy: collections.abc.Sequence[float], oz: collections.abc.Sequence[float]) -> None Analytic-SDF capability: exact face-aperture (openness) fields for the cut-cell projection, overriding the sampled-SDF openness. Inner nxnynz arrays, x-fastest; ox[i] = fluid fraction of the -x face of cell i. Call BEFORE set_solid (a later call raises); empty ox clears. Single-rank.
`set_outer_iterations`	set_outer_iterations(self, n: int) -> None Set the number of Picard/outer iterations per step.
`set_outer_tolerance`	set_outer_tolerance(self, tol: float) -> None Set the outer (Picard) convergence tolerance.
`set_outflow_rho_correction`	set_outflow_rho_correction(self, on: bool = True) -> None The 1/rho_f mobility factor on the HIGH-side outflow face correction. DEFAULT TRUE since WO-R2; pass False for the ablation. A projection correction cancels the discrete divergence only if it uses the SAME face coefficient the operator row used. Until WO-R2 the multigrid re-imposed the literal openness 1.0 at every Dirichlet domain face, overwriting buildRhoCoeff's open_f*rho0/rho_f, so the plain phi difference was the consistent partner and WO-R measured this factor making things seven orders WORSE. WO-R2 fixed the operator (CutcellMG::setOutflowCoefficient) and the verdict inverted. Stratified duct, ratio 10, 5 steps, max
`set_phase_change_area`	set_phase_change_area(self, mode: str) -> None WO-P3c: WHICH GEOMETRY the interfacial area A_Gamma comes from. A_Gamma sets the plane shift dV = mdot A dt / rho_l and the divergence source S = mdot A (1/rho_g - 1/rho_l), so a bubble grows as int mdot dA and a biased area is a biased growth rate. 'plic' (DEFAULT, rungs P0/P1): plicArea =
`set_phase_change_budget`	set_phase_change_budget(self, on: bool = True) -> None WO-P3f INSTRUMENT: turn on the energy budget of the phase-change energy solve. Costs one extra cell field and two reductions per energy solve, and is skipped entirely when off (the solve is then bit-identical). Read it with `phase_change_budget()`. What it exists for: interfacial cells are Dirichlet rows, i.e. they are OUTSIDE the energy solve, so the set over which enthalpy is conserved changes membership every step as the interface sweeps. A liquid cell that becomes interfacial LEAVES that set carrying its superheat rho c_p (T - T_sat) and an interfacial cell that becomes pure RE-ENTERS it carrying `pcCarriedValue`; neither transfer appears in the latent-heat book-keeping.
`set_phase_change_carry_conserve`	set_phase_change_carry_conserve(self, on: bool = True) -> None WO-P3f OPTION (default OFF): make the interfacial cells' per-cell Dirichlet OVERWRITE enthalpy-conserving. An interfacial cell's row is the identity `T = dval`, so whatever the geometric energy transport left there is discarded every step. Over a cell's interfacial lifetime that telescopes to rho c_p (T_entry - dval_exit): a liquid cell the interface sweeps enters with its superheat and leaves as vapour at T_sat, and the difference goes nowhere. Measured on Scriven 128^3 by `phase_change_budget()['d_overwrite']`: -0.7 % of mdot h_lv A_Gamma at Ja = 0.5 and -4.3 % at Ja = 2, one-signed. With this on, that enthalpy is handed to the interfacial cell's face neighbours that are still in the solve, weighted by n_d^2 (the clip-and-redistribute allocation), as a fixed-order gather so it is decomposition-independent. Which side receives is decided per axis by which pure neighbour deviates from T_Gamma in the same direction as the interfacial cell itself, i.e. the phase the enthalpy came from -- the superheated liquid on an evaporating bubble, the superheated vapour on the Stefan problem.
`set_phase_change_curvature_distance`	set_phase_change_curvature_distance(self, on: bool = True) -> None WO-P3g item 3 (default OFF): measure the GFM row's theta and the one-sided fits' sample distances to the CURVED interface, with the mean curvature taken PER CELL from the V3 curvature cascade (the same kappa surface tension uses; positive for a convex blob of liquid, i.e. -2/R for a gas bubble). Supersedes `set_phase_change_fit_curvature`, which prescribes one curvature for the whole field; where both are set the cascade wins.
`set_phase_change_deposit_fallback`	set_phase_change_deposit_fallback(self, on: bool = True) -> None WO-P3f open item 6 / WO-P3g (OFF by default): give an interfacial cell whose two along-the-normal deposit candidates are BOTH still interfacial a target from the 5^3 box instead of leaving the divergence source in place. A cell that keeps its source carries div(open u) = S on its own faces, so Weymouth-Yue advects the colour with a field that is not the liquid velocity -- read it out with `phase_change_diagnostics()['band_div']` and the 'fallback_cells' count.
`set_phase_change_energy_muscl`	set_phase_change_energy_muscl(self, on: bool) -> None MinMod-limited donor reconstruction of the face temperature in the CONSISTENT energy flux, instead of the default plain donor-cell (first-order upwind) value. OFF by default. The consistency is what the geometric flux buys; its first-order upwind numerical diffusion
`set_phase_change_energy_order`	set_phase_change_energy_order(self, order: int) -> None WO-P3g: the ORDER of the interfacial energy operator. 1 (the shipped WO-P23...P3f scheme) or 2. `order = 2` turns on, TOGETHER, the four pieces WO-P3f's instruments indicated: 1. `set_phase_change_mdot_operator(True)` -- mdot is the energy operator's OWN interfacial flux q/(h_lv A) instead of a separate least-squares fit, so the heat the energy equation loses and the mass the regression produces are one discrete quantity (and the interfacial AREA cancels out of the mass balance entirely); 2. `set_phase_change_gfm_order(2)` -- the Gibou-Fedkiw three-point ghost-fluid row (2/((1+theta) theta), 2/(1+theta)) instead of the two-point (1/theta, 1), which is exact on a quadratic profile at every theta; 3. `set_phase_change_curvature_distance(True)` -- the row's theta and the one-sided fits' sample distances are measured to the CURVED interface, with kappa taken per cell from the V3 curvature cascade; 4. `set_phase_change_carry_conserve(True)` -- WO-P3f's enthalpy-conserving Dirichlet overwrite. Why together and not one at a time: WO-P3f measured the shipped scheme's 1 % Scriven error to be the residue of a CANCELLATION between the fit's +6 % curvature bias, the two-point row's -5 % flux deficit and the overwrite's -0.7...-4.3 % enthalpy destruction, so repairing any ONE alone makes the gate worse.
`set_phase_change_fit_curvature`	set_phase_change_fit_curvature(self, kappa: float) -> None WO-P3f INSTRUMENT (default 0, bitwise inert): prescribe the interface curvature `kappa = div(n)` that the one-sided temperature fits use to measure a sample's distance to the CURVED interface instead of to the interfacial cell's tangent plane (`vof::pcCurvedDistance`). For a spherical GAS bubble of radius R, whose PLIC normal points inward, kappa = -2/R. Why it exists: the tangent-plane distance makes the fit FIRST ORDER in h/R on a curved interface, and every off-axis sample is hotter than the plane model expects, so the fitted dT/dn and with it mdot come out HIGH. Measured a priori on an exact sphere with an exactly linear profile: +19.2 / +12.1 / +8.8 / +6.2 % at R = 6 / 10 / 14 / 20, observed order 0.91-0.98 in h/R. This entry point takes a PRESCRIBED kappa so that bias can be measured against a known geometry; it is not a curvature estimator.
`set_phase_change_gfm_order`	set_phase_change_gfm_order(self, order: int) -> None WO-P3g item 2 (default 1): the order of the one-sided (ghost-fluid) Dirichlet row. 1 = the shipped two-point form k o (T_Gamma - T_i)/theta; 2 = Gibou-Fedkiw's three-point form through (T_behind, T_i, T_Gamma), which reproduces a quadratic temperature profile exactly at every theta.
`set_phase_change_mdot_operator`	set_phase_change_mdot_operator(self, on: bool = True) -> None WO-P3g item 1 (default OFF): take mdot from the energy operator's own interfacial flux -- the sum of the ghost-fluid rows' Dirichlet couplings evaluated with the converged T -- instead of the one-sided least-squares fit, which stays as `phase_change_diagnostics()['mdot_fit']`.
`set_phase_change_plane_dirichlet`	set_phase_change_plane_dirichlet(self, on: bool) -> None The PLANE-ANCHORED (ghost-fluid) interfacial Dirichlet condition. ON by default. Rungs P0/P1 pinned the whole interfacial CELL at T_G, so the numerical thermal boundary sat at the cell CENTRE while the mass-flux gradient is fitted from the PLIC PLANE — a mismatch of up to half a cell that CHANGES SIGN as the interface sweeps through a cell, and the first-order component of the P1 Stefan error (WO-P01 finding 6). With this on, the condition becomes PER FACE: for a pure cell i whose neighbour j is interfacial, that face's row carries kopen(T_i - T_G)/theta with theta the distance in cells from i's centre to j's PLIC plane along the face's own axis, and the interfacial cell's own value is never read by any neighbour. A per-CELL value (the literal reading of the work order) does NOT work and the measurement is in the findings: giving the interfacial cell the value the one-sided profile takes at its centre is right for the side the fit came from and wrong for the other side — on the P1 Stefan ladder it heats the saturated liquid through the interfacial cell's liquid-side face, which gives the liquid a spurious gradient that feeds straight back into mdot: +6.20/+5.62/+5.42 % at N = 64/128/256, order 0.10, against +1.31/+0.59/+0.20 % and order 1.37 for the cell-centre pinning it was meant to improve on. The per-face form has no such asymmetry (the liquid-side face reads T_G at its own theta, i.e. exactly zero flux for a saturated liquid). False restores the rung P0/P1 behaviour bit-for-bit.
`set_phase_change_quadratic_fit`	set_phase_change_quadratic_fit(self, on: bool) -> None Fit the one-sided interfacial temperature gradients QUADRATICALLY through the interface value (T - T_G = G phi + Q phi^2) instead of linearly. This is VOF_PLAN section 9 item 1's Aslam quadratic extrapolation in least-squares form — the same 5^3 pure-cell samples, the same Malan collinearity weights, one more basis function, no PDE sweeps. Once the plane-anchored Dirichlet has removed the cell-centre mismatch, the linear fit's O(T'' h) curvature bias is the leading error of the rung: its samples start about one cell from the plane and reach two and a half, so a curved profile tilts the straight line through the origin. ON by default (WO-P23): with the plane-anchored Dirichlet it takes the P1 Stefan interface position from +0.195 % to +0.003 % at N = 256, and the mdot kernel itself from order 1.1 to order 2.0. set(False) is the ablation.
`set_porous_conservative`	set_porous_conservative(self, on: bool) -> None eps-conservative porous momentum + projection pair (default True): time term (eps_f rho/dt) u, eps rho-weighted advective form, projection coefficients open*(eps rho idt)/(eps rho idt+beta) with matching correction. False = the legacy plain-u pair (A/B only; it lets the projection drag gas with the moving porosity at zero inertia cost — a spurious late-time energy source in clustering flows).
`set_porous_deps_dt`	set_porous_deps_dt(self, on: bool) -> None Include (default True) or drop the d(eps)/dt source in the porous projection RHS. Drop it to enforce div(eps u)=0 when the per-cell eps deposit's time-derivative is too jagged and destabilizes the eps-weighted pressure solve.
`set_pressure_bottom_solver`	set_pressure_bottom_solver(self, engine: str) -> None Engine of the AGGLOMERATED pressure bottom (set_pressure_bottom 'auto' / 'agglomerated'). 'auto' (DEFAULT): on a GPU backend, single rank, the singular (periodic / wall) operator, a bottom of at most 8192 cells with at most 64 fluid components (solids allowed) and an axis whose planes hold at most 192 cells, 'direct': one device launch runs flexible CG in FP64 (inner tolerance 1e-5, cap 100) preconditioned by a block-tridiagonal FP32 direct factor of the bottom (planes along one axis, explicit Schur-complement inverses), refactored once per operator change -- no host transfer; everywhere else the host GraphAMG solve. 'direct' forces the device engine and raises at the bottom solve, naming the failed condition, where it is not eligible; 'algebraic' forces GraphAMG (the A/B instrument). Host backends always run GraphAMG under 'auto'. doc/vof_step_performance_design.md §13.
`set_pressure_coarse_ghost`	set_pressure_coarse_ghost(self, on: bool) -> None Zero-gradient (Neumann) coarse ghost on wall/inflow faces before the pressure multigrid's prolongation -- the WO-H symmetry repair. DEFAULT True. False restores the pre-2026-08-30 periodic-wrap ghost, which reinstates the asymmetry that stalls MG-PCG on every 3-D wall-bounded grid: a MEASUREMENT ABLATION, never a production setting. Inert on periodic / IBM problems, where the two are the same.
`set_pressure_exact_residual`	set_pressure_exact_residual(self, on: bool = True) -> None Apply the level-0 pressure operator EXACTLY (matrix-free, double, flux form) in the residual and the Krylov matvec instead of reading the float band storage. P1 of the suite defect-correction campaign (suite/docs/archive/DEFECT_CORRECTION_PLAN.md). Per solver, OFF at construction. enable_vof() turns it ON, because a two-phase coefficient contrast is exactly what amplifies the float operator's broken row-sum identity A*1 = 0. Measured on Hysing case 2 (64x128x4, adaptive dt, nvidia-cuda): max
`set_pressure_graph_amg`	set_pressure_graph_amg(self, on: bool) -> None Solve the pressure MG's coarsest level with an agglomerated mesh-agnostic algebraic multigrid (core GraphAMG), decomposition-agnostic: with levels=1 this gives a mesh-independent pressure solve that works under a WEIGHTED ORB (where the geometric coarse levels can't cleanly coarsen). Applied at the next set_solid.
`set_pressure_mean_removal`	set_pressure_mean_removal(self, scope: str) -> None Nullspace (mean) removal scope in the pressure solve: 'fine' (DEFAULT — only the projections the Krylov iteration needs: rhs/residual, the fine-level V-cycle exit, the final iterate; ~3x fewer global-reduction latency hits per iteration, the measured winner of the multi-GPU ablation) or 'all' (legacy — every V-cycle level + after every matvec). Iteration counts are identical (A preserves mean-freeness); results equal within solver tolerance, not bit-identical.
`set_pressure_strict`	set_pressure_strict(self, on: bool) -> None Raise instead of report when the pressure preconditioner returns a non-finite correction. DEFAULT False: the shipped behaviour caps the solve, prints one line, and raises the pressure_solve_failed() flag, so a check can catch it without changing control flow. True turns the same event into an exception.
`set_pressure_underrelax`	set_pressure_underrelax(self, omega: float) -> None Pressure under-relaxation factor omega_p in (0,1] for the incremental accumulation (MFIX §10.1); 1.0 = off (default). <1 damps the incremental predictor overshoot on stiff porous+drag.
`set_pressure_warmstart`	set_pressure_warmstart(self, on: bool) -> None Seed each pressure solve from the previous step's phi (default off).
`set_property_mode`	set_property_mode(self, mode: str = 'variable', harmonic: bool = False) -> None Enable variable-coefficient momentum (variable viscosity): mode 'variable' binds the 'mu' field (get/set_field('mu')) into the diffusion operator; 'constant' reverts. harmonic = harmonic face-viscosity mean (continuous shear stress across a jump) vs arithmetic. A closure targeting 'mu' enables this automatically. The incremental-rotational pressure scheme (large-dt / steady-Stokes) stays active — see set_variable_rotational.
`set_rho_face_harmonic`	set_rho_face_harmonic(self, on: bool = True) -> None Use the HARMONIC instead of the arithmetic face mean of rho in the pressure projection (both the operator coefficient and the velocity correction, so the projection stays exact). DEFAULT OFF, and it should stay off: the arithmetic mean of rho IS the harmonic mean of the mobility 1/rho — the series-correct choice for a normal flux — and it is what makes the discrete hydrostatic balance EXACT, because the momentum time term and the face body force interpolate rho arithmetically and are not switched by this flag. Shipped as a measured knob for the coefficient-coarsening question (VOF_PLAN S3), not as an alternative scheme.
`set_rotational_filter`	set_rotational_filter(self, on: bool, eps: float = 0.05) -> None Experimental filtered rotational update: smooth div(u) (mask-aware axis-wise 1-2-1, one-sided into the fluid at walls) before accumulating -mudiv into P. Keeps the O(1) pressure-relaxation gain, removes the checkerboard feedback channel.
`set_rotational_pressure`	set_rotational_pressure(self, on: bool) -> None PM I ablation (Guy-Fogelson 2005): False drops the rotational -mudiv(u) term from the incremental pressure accumulation (constant-mu path only). Default True = shipped rotational (Timmermans) update.
`set_rotational_wall_weight`	set_rotational_wall_weight(self, w0: float) -> None Wall-banded rotational blend: at fluid cells with a solid axis-neighbour use P += (rho/dt + w0mu/dx^2)phi - (1-w0)mudiv(u*); bulk keeps the full rotational update. Stabilizes the boundary rows without slowing bulk pressure relaxation, and keeps them relaxing at dt->infinity. 0 (default) = off.
`set_rotational_weight`	set_rotational_weight(self, w: float) -> None Under-relax the rotational term: P += ctphi - wmudiv(u). 1 = shipped, 0 = PM I; small w raises the boundary-mode stability threshold ~1/w at ~1/w slower smooth-mode pressure relaxation. phi=0 stays the unique fixed point for any w>0 at every dt.
`set_uf_advection`	set_uf_advection(self, on: bool) -> None Advect the COLLOCATED momentum with the projected, divergence-free MAC face field uf/vf/wf of the last projection instead of the un-projected cell->face average 0.5*(u_i+u_j). DEFAULT True -- the shipped scheme: doc/flow_colocated_plan.md step 3, the Almgren-Bell-Colella prescription (the field the projection just made solenoidal IS the conservative advective flux), and what the FOU operator's conservative row-sum identity needs. False is the ablation that restores the phase-2 average (measurement only; see doc/uf_advection.md). Inert on the staggered grid (the stored velocity already IS the projected face velocity) and before the first projection or set_velocity has built a face field.
`set_variable_rotational`	set_variable_rotational(self, mode: str = 'min', chi: float = 1.0) -> None Rotational-pressure term under variable viscosity (the constant-mu Timmermans term -mudiv(u) is only valid for homogeneous viscosity — Deteix & Yakoubi 2018). 'min' (default): constant coefficient chimu_min — provably stable at any contrast, exact fallback to the constant-mu scheme for uniform mu. 'full': pointwise chimu(i) — better pressure consistency at MILD contrast only. 'off': plain incremental (no rotational term). All modes keep the incremental predictor (large-dt / steady-Stokes capability).
`set_velocity_mg_smoother`	set_velocity_mg_smoother(self, name: str, degree: int = 0, eig_ratio: float = 10.0) -> None Which smoother the velocity multigrid uses on every level -- orthogonal to set_velocity_solver, and read only when that selects 'multigrid'. 'gauss_seidel' (DEFAULT) or 'chebyshev'. ABLATION, AND CHEBYSHEV IS MEASURED SLOWER: at degree 4 / eig_ratio 6 it IS the stronger smoother per cycle (8.0 V-cycles against red-black's 9.3) and still loses 1.21x on wall time, because matching red-black's smoothing strength costs the same four halo exchanges two sweeps cost, while Chebyshev needs a separate residual and update pass where a Gauss-Seidel sweep fuses the two. Kept because it is the apparatus for a question that gets asked repeatedly. degree=0 follows the V-cycle's own pre/post/bottom counts.
`set_velocity_multigrid_auto`	set_velocity_multigrid_auto(self, cells_per_rank: int, min_global_cells: int = -1, min_block_extent: int = -1) -> None Narrow the rule that selects the momentum solver. Since 2026-09-15 the DEFAULT is the velocity V-cycle on every validated operator mode and at every rank count: measured on the 1.0.0 scaling benchmark it is faster at every rung of both the CPU and the GPU ladder (2.23x at 147k cells/rank, 2.19x on one H100), and the red-black smoother it replaced was running at its sweep cap -- not converging -- on every rung of both. cells_per_rank sets an upper bound in cells per rank (0 = never use the V-cycle); set_velocity_multigrid_auto(65536, 1 << 23) restores the 1.0.0 behaviour exactly. min_block_extent (default 16) is the shortest per-rank inner extent worth building a hierarchy on. An explicit set_velocity_multigrid() overrides all of it.
`set_velocity_solver`	set_velocity_solver(self, name: str) -> None Which solver runs the implicit momentum (diffusion) equation. 'auto' (DEFAULT) picks by the operator's condition number kappa = 1 + 4 dt mu (sum w)/rho -- the multigrid at kappa >= 13, the smoother below -- and names no geometry, so an immersed body and a domain-BC box at the same dt get the same solver. 'gauss_seidel' is the red-black smoother; 'multigrid' the 3-level velocity V-cycle; 'chebyshev' Chebyshev semi-iteration over a Gershgorin interval. ALL FOUR SOLVE THE SAME EQUATION TO THE SAME TOLERANCE, so this changes cost and not the answer -- except that 'gauss_seidel' runs out of sweeps above kappa ~ 50 and then returns a solution that never met it. MEASURED (384^3 cut-cell bed, Snellius): multigrid beats gauss_seidel 1.66x-2.23x across the CPU ladder and 2.19x on one H100; chebyshev beats gauss_seidel by the sqrt(kappa) factor theory predicts but loses to multigrid except at the smallest per-rank blocks, where the evidence is suggestive and not conclusive.
`set_velocity_solver_params`	set_velocity_solver_params(self, iters: int, rtol: float = 0.0, min_iters: int = 2) -> None Momentum smoother control: `iters` RB-GS sweeps per component, or with rtol > 0 a TOLERANCE STOP — end the loop once the sweep's max velocity increment has contracted to rtol of the first sweep's (iters becomes the cap, min_iters the floor). Easy regimes (small nu*dt/dx^2) exit after ~3-5 sweeps; stiff regimes run to the cap unchanged. rtol = 0 (default) is the legacy fixed count, byte-identical.
`set_vof_block_debris`	set_vof_block_debris(self, enabled: bool) -> None Per-step marker DEBRIS removal of the block container (doc/vof_overlap_design.md 5.3, 11). Default: ON under enable_vof_block_csf, off otherwise. A debris cell is an interfacial cell of a marker (1e-8 < C < 1-1e-8) with every cell of that marker within two cells below C = 1/2: colour with no resolvable interface near it (what two overlapping markers tear off each other). Its colour is zeroed and the volume is returned to the marker's attached interfacial cells weighted C(1-C), in a fixed summation order, so the marker volume is exact and the state decomposition-independent. Ledger in vof_block_stats(): debris_cells / debris_volume (this step), debris_returned, debris_lost, debris_unresolved (cumulative, migrate with the block). enabled=False keeps the census (what WOULD be removed) and moves nothing.
`set_vof_block_device_staging`	set_vof_block_device_staging(self, on: bool) -> None Pack/unpack the block gather and scatter in the block's own MEMORY SPACE (rung W1 item b, default True), with a host staging copy only per MPI MESSAGE — and none at all for the master's own cells, which are a device-to-device copy. False selects rung W0's host-staged path (a full mirror of the local patch per step), which is what the device-vs-host measurement compares against. Every step is a copy of a double, so the two paths are BITWISE identical; the ctest gates that rather than asserting it.
`set_vof_block_kappa_clip`	set_vof_block_kappa_clip(self, enabled: bool, kappa_max: float
`set_vof_block_overlap_density`	set_vof_block_overlap_density(self, enabled: bool) -> None OPT-IN MODEL, default False (the MAX union verbatim). With it the union colour that feeds rho(C) and mu(C) -- and the registered 'C' -- becomes C_eff = max_k C_k where S = sum_k C_k <= 1 + 1e-8 and max(0, 2 - S) where markers overlap: the overlap volume is treated as the sub-grid FILM liquid (the 'tent'), so the phantom interface of one marker inside another is an ordinary liquid-gas interface under the base capillary bound. Its cost: the lens holds liquid the real state does not (buoyancy and mixture weight off by 2 V_overlap during a contact), and it changes collision dynamics. See doc/vof_overlap_design.md 5.6.
`set_vof_block_pool`	set_vof_block_pool(self, on: bool) -> None Recycle the advectors a re-centring retires, keyed by the exact box extent (rung W1 item c, default True). A translating bubble keeps its box SIZE and only moves its origin, so the hit rate is ~100 % and the ten Views of the new box are not allocated. A recycled advector is handed back in the state a freshly initialised one is in (colour and the three face-velocity fields zeroed), so the pool is bitwise inert.
`set_vof_curvature_mixed_height_fit`	set_vof_curvature_mixed_height_fit(self, on: bool = True) -> None Enable cascade branch 3, the mixed height-position fit (a paraboloid through the interface positions of whichever columns closed). DEFAULT OFF and it should stay off: measured on an exact-fraction sphere at 16/32/64 it takes over the 19.5-59.6% of cells the height function cannot serve and DESTROYS the convergence of the max curvature error (order 0.00 vs 1.86 with the PLIC-volumetric fallback instead), because its data set is the columns that closed - a slope-selected, asymmetric subset whose lever-arm bias is scale invariant. Shipped as a re-measurable instrument, not a configuration.
`set_vof_curvature_weight_width`	set_vof_curvature_weight_width(self, width: float) -> None Wendland support width of the PLIC-volumetric fallback fit, in CELL units. Default 2.5, which is Han et al.'s recommended pairing with the 5^3 stencil. Their translating droplet recovers first-order convergence of the spurious currents at 3.5 and loses convergence entirely at 4.5 (over-smoothing), so this is a real trade-off between the locality of the estimate and its robustness to transport error — not a free knob.
`set_vof_curvature_worklist`	set_vof_curvature_worklist(self, on: bool) -> None WO-V9: run the V3 curvature cascade (the PLIC plane pass, the height functions and the PV fallback) over a COMPACTED list of the interfacial cells instead of over the whole inner region. Default ON. The cascade is the most divergent kernel in the pipeline -- on a resolved droplet the interfacial cells are well under 1 % of the block, so a dense parallel_for puts one or two active lanes in a warp and the other thirty run the guard and then idle for the whole height function. Compaction is a pure re-ordering (each cell reads the same neighbours and writes the same value), so the two paths are BIT-IDENTICAL; tests/kokkos/test_vof_timing.cpp gates that and WO-V9 records the gain.
`set_vof_cutcell_flux_clamp`	set_vof_cutcell_flux_clamp(self, on: bool) -> None Ablation: Weymouth's admissible-interval clamp on the openness-weighted cut-cell flux (ON by default). The clamp bounds
`set_vof_flux_clamp`	set_vof_flux_clamp(self, on: bool = True) -> None ABLATION: drop the Weymouth flux clamp max(0,
`set_vof_interface_eps`	set_vof_interface_eps(self, eps: float) -> None Wisp threshold on the curvature cascade's interfacial predicate once surface tension is on: a cell carries an interface only while eps < C < 1-eps. Default 1e-8. This is NOT optional and it is not cosmetic. Weymouth-Yue leaves round-off colour residue (measured down to -3e-35) in every cell its sweeps touch; those cells satisfy 0 < C < 1, so the cascade builds a PLIC polygon of area ~0 for them and returns
`set_vof_kappa_constant`	set_vof_kappa_constant(self, kappa: float) -> None INSTRUMENT: set kappa to a constant everywhere (inner + ghosts), mark every cell as carrying a valid estimate, and freeze it. The CSF force is then EXACTLY the discrete gradient of sigmakappaC, so the projection must annihilate it to round-off from any colour field whatsoever — the exactness gate of rung V4, independent of curvature accuracy and of resolution.
`set_vof_kappa_frozen`	set_vof_kappa_frozen(self, on: bool = True) -> None INSTRUMENT: stop recomputing the curvature at the head of each step and use whatever is in the 'kappa'/'kappa_branch' fields. With set_vof_kappa_constant this isolates the balanced-force identity from the curvature estimator.
`set_vof_momentum_cell_flag`	set_vof_momentum_cell_flag(self, on: bool = True) -> None ABLATION: use the PRESSURE-cell frozen dilation flag H(C^n-1/2) on the shifted control volume instead of its structural analogue H(C^c,n-1/2). The flag that must be shared is the one of the pair that telescopes, and both members of that pair live on the shifted volume — this switch is the literal reading of the work order, kept as a measurement.
`set_vof_momentum_muscl`	set_vof_momentum_muscl(self, on: bool = True) -> None MinMod-limited linear reconstruction of the donor velocity in the momentum flux, instead of the DEFAULT plain donor-cell upwind. Both preserve the uniform-velocity identity exactly (a uniform field has a zero slope bit for bit), but on a control volume a sweep EMPTIES the slope's deviation from the volume's own velocity is amplified by drho*F/rho^c, which is unbounded in the density ratio. Measured at ratio 1e4, 50 steps: with the slope the uniform-velocity residual grows to 2.2e-10, without it it is flat at 6.7e-16; at ratio 1e3 the slope is harmless. Turn it on deliberately and re-run the ratio sweep if you do.
`set_vof_phantom_capillary_bound`	set_vof_phantom_capillary_bound(self, enabled: bool = True) -> None The phantom-aware capillary limit of the block container (doc/vof_overlap_design.md 5.6, 11.3), DEFAULT True: while two block markers overlap (max S > 1 + 1e-2) the capillary limit of vof_step_limits() and step() uses 2 rho_min in place of rho_l + rho_g -- the part of one marker's surface inside another is an interface with gas on both sides. False restores the pre-5.6 bound (the ablation; a fixed-dt driver at a density ratio above ~31 otherwise throws at first contact).
`set_vof_rho_floor`	set_vof_rho_floor(self, fraction: float) -> None Floor on rho^c in the recovery divide u = (rho^c u)/rho^c, as a FRACTION of min(rho_gas, rho_liquid). Default 1e-6. rho^c leaves [rho_gas, rho_liquid] only through a wisp in the half-shifted colour and reaching zero would need C^c ~ -1/(ratio-1), so this is a guard, not a model — vof_momentum_diagnostics()['floored'] reports how many control volumes it actually touched.
`set_vof_step_parity`	set_vof_step_parity(self, n: int) -> None Set the sweep-permutation counter of the NEXT colour advection: the Weymouth-Yue sweep order is kWySweepPerm[n % 6], cycled so no axis is systematically favoured. Exposed so a benchmark can hold the permutation fixed (n constant) or resume a run across a restart. Default: it increments once per advection from 0.
`set_vof_timing`	set_vof_timing(self, on: bool) -> None WO-V9: arm (or disarm) the VoF pipeline's per-stage timers, and reset them. OFF by default. When armed, every stage boundary calls Kokkos::fence() before reading the clock -- the same rule the step's three coarse phase timers already follow, because on a device backend queued work would otherwise be billed to whichever stage next reads the clock. A fence changes WHEN work happens, never WHAT it computes: a run with the timers armed is bit-identical to the same run without them (gated in tests/kokkos/test_vof_timing.cpp). When disarmed the cost is one branch per stage and no fence at all.
`set_vof_wisp_eps`	set_vof_wisp_eps(self, eps: float) -> None Wisp tolerance on the Weymouth-Yue mixed-cell predicate: a cell counts as carrying an interface only while eps < C < 1 - eps, and one outside that band is fluxed ALGEBRAICALLY as C_donor * a — its ACTUAL colour, so the exact telescoping conservation is untouched. The same threshold gates the interface-local Courant band (
`set_vof_worklist`	set_vof_worklist(self, on: bool) -> None WyAdvector::useWorklist -- compact the PLIC reconstruction pass onto the mixed cells (a parallel_scan over the reconstruction region, then a dense parallel_for over the compacted list) instead of a guarded parallel_for over the whole region. Default ON. Pure optimization: the flux never reads a non-mixed cell's plane, so the two paths produce the same field BIT FOR BIT (gated in tests/kokkos/test_vof_advect.cpp and measured in WO-V9).
`set_wall_flux_divergence`	set_wall_flux_divergence(self, on: bool) -> None Rung 3 on/off (default on = correct physics). Off leaves the moving wall's no-slip datum in the momentum operator but drops the wall's own volume flux from the cut-cell projection -- the configuration the Galilean gate uses to exhibit the failure the term fixes.
`uf_advection`	Whether the projected face field is the collocated advecting velocity.
`velocity_mg_smoother`	velocity_mg_smoother(self) -> str The velocity multigrid's smoother.
`velocity_multigrid_active`	velocity_multigrid_active(self) -> bool Whether the velocity (momentum) solve uses the multigrid V-cycle rather than plain RB-GS. Reads back the OUTCOME, not the setting: an explicit set_velocity_multigrid(on=...) call always wins, but absent one the AUTO rule (set_velocity_multigrid_auto) decides it when the geometry is built (set_solid / set_pressure_geometry / set_solid_from_scene) -- on for a distributed run (np > 1) in an eligible operator mode with a global problem of at least min_global_cells (default 8M) and fewer than cells_per_rank (default 65536) cells per rank, off otherwise. So the value can differ from what the caller last set explicitly, and is only meaningful after geometry has been installed (False on a fresh solver).
`velocity_solver`	velocity_solver(self) -> str The momentum solver actually in force, resolving 'auto'. Only meaningful after the first step(), which is where the choice is made (dt and mu are final only then).
`vof_advected_velocity`	vof_advected_velocity(self, component: int) -> numpy.ndarray[dtype=float64] The recovered advected velocity (rho^c u_c)/rho^c of component c on the inner cells — the momentum RHS's time base. Exposed so the uniform-velocity consistency identity can be gated on the ADVECTION ALONE, with the projection and the momentum solve out of the picture.
`vof_block_census`	vof_block_census(self) -> dict Per-rank load census of the block container: 'masters'[r] = blocks rank r masters, 'cells'[r] = the inner cells those blocks carry (the actual VoF work).
`vof_block_curvature_stats`	vof_block_curvature_stats(self) -> dict Branch census of the last block-CSF curvature pass, SUMMED over this rank's blocks (local to the rank). 'no_estimate' must be 0 on any gated case. 'clipped' = cells whose
`vof_block_force`	vof_block_force(self, component: int) -> numpy.ndarray[dtype=float64] The block CSF face force of component c (0,1,2) as scattered (summed over markers) into the solver's LOW-face slot of each cell, internal units; the field the momentum RHS adds. A diagnostic.
`vof_block_imbalance`	vof_block_imbalance(self) -> float max/mean of the per-rank block-cell load under the CURRENT master assignment (round robin by block id at rung W0). 1.0 is perfect; the weighted-ORB assignment of rung W1 is what this number is here to grade.
`vof_block_imbalance_of`	vof_block_imbalance_of(self, mode: int) -> float The max/mean per-rank block-cell load the given assignment mode WOULD give on the current blocks, without applying it — so a study can put the three modes side by side on one swarm without perturbing the run.
`vof_block_kappa`	vof_block_kappa(self, id: int) -> numpy.ndarray[dtype=float64] One marker's OWN curvature (internal units, 1/cell) over its block box, the layout of vof_block_color; empty unless this rank masters it and the block CSF is on. A diagnostic.
`vof_block_pool_stats`	vof_block_pool_stats(self) -> dict Block-pool census: 'hits' = advectors recycled, 'misses' = advectors allocated.
`vof_block_stats`	vof_block_stats(self) -> list Per-bubble Lagrangian census, one dict per block, in block-id order. 'lo'/'hi' (the global index box) and 'master' are replicated on every rank; the MEASURED entries — 'volume' (sum of C over the box, cell volumes), 'centroid', 'velocity' (d(centroid)/dt of the last advection), 'moments' (the central second moments xx, yy, zz, xy, xz, yz divided by the volume, i.e. the deformation) — are filled only on the block's MASTER and are zero elsewhere. 'discarded' is the colour a re-centring dropped, cumulatively: Weymouth-Yue leaves round-off residue in every cell its sweeps touch and the box tracks the BUBBLE, not that wake, so the residue falling outside the new box is discarded. Never physical liquid (measured -9.5e-17 over a 20-cell translation, against a bubble volume of 524), but it is reported rather than hidden — a container that silently loses mass is not acceptable. 'debris_cells' / 'debris_volume' are this step's marker DEBRIS removal (block CSF on only; master only; see set_vof_block_debris): interfacial cells with every cell of the marker within two cells below C = 1/2, and their summed colour. 'debris_returned', 'debris_lost' and 'debris_unresolved' are the cumulative ledger (doc/vof_overlap_design.md 5.3): volume + discarded + debris_lost is conserved. 'full_axis' counts re-centrings whose box had to span a whole periodic axis (it is then [0, L) and copied with the wrap, losslessly; 12.3). 'residue_returned' is the cumulative SIGNED sub-wisp-threshold residue cleared after each advection and returned to the interface with the debris (13): volume + discarded + debris_lost is exact. 'overlap_*' are the GLOBAL marker-overlap census of the last block advection (vof_overlap_design 5.6; block CSF or the overlap density on), the same in every dict: 'overlap_max_sum' = max over cells of S = sum_k C_k, 'overlap_excess' = sum of (S - 1)^+ (cell volumes), 'overlap_cells' = cells with S > 1 + 1e-8, and 'overlap_bound_active' = the phantom-aware capillary limit (2 rho_min in place of rho_l + rho_g) is in force for the next step (triggered by max S > 1 + 1e-2).
`vof_curvature_worklist`	vof_curvature_worklist(self) -> bool Whether the curvature-cascade compaction is on.
`vof_diagnostics`	vof_diagnostics(self) -> dict Colour census over THIS RANK's inner cells: sum (cell-volume units), min, max, the number of mixed cells (0<C<1) and of wisps (C within 1e-8 of 0 or 1). No clipping is applied at this rung, so min/max may leave [0,1] if the CFL cap is raised. With an immersed solid (rung V5a) the dict also carries the cut-cell quantities, all zero without one: 'volume' = sum eps_effC over fluid cells, the functional the openness-weighted scheme conserves EXACTLY (eps_eff = max(eps, 1/64), the 4^3 subsampling resolution — see src/vof/cutcell.hpp rule 1); 'raw_volume' = sum epsC, which differs only on eps==0 cells that still own an open face and is therefore NOT the conserved functional; 'solid_sum' = sum of C over solid cells (0 by construction — the neutral band fill lives on the working block, not on 'C'); 'min_fluid'/'max_fluid' over UNCUT fluid cells (eps==1), where Weymouth's boundedness applies verbatim; 'clipped_volume' / 'clipped_signed' = the liquid volume the cut-cell clip moved during the last advection (a TRIPWIRE on the flux approximation, not a mechanism: if it is not negligible the fix is the solid-clipped flux polygon); 'cut_cells' / 'solid_cells'. 'inflow_volume' / 'outflow_volume' are the liquid volume that ENTERED / LEFT through the domain faces during the LAST colour advection, in the same cell-volume units as 'sum' (0 unless a VoF boundary colour is set — rung V-BC). They are the advector's OWN face fluxes, so sum(C) - sum(C_0) = integral(inflow - outflow) holds to round-off whatever the interface does; vof_bc_volumes() breaks them out per face.
`vof_dynamic_field`	vof_dynamic_field(self, which: int) -> numpy.ndarray[dtype=float64] The V6 per-cell dynamic-wetting state on the inner region, as (nx,ny,nz): 0 the IMPOSED angle in degrees, 1 the measured APPARENT angle in degrees, 2 the smoothed U_cl (positive = the liquid ADVANCES), 3 Ca_cl, 4 the state (0 not a contact cell, 1 Cox-Voinov on the static base, 2 PINNED, 3 advancing, 4 receding). Regenerates the fill, so it is also the direct gate on the pass's decomposition independence.
`vof_filled_color`	vof_filled_color(self) -> numpy.ndarray[dtype=float64] The colour field INCLUDING the neutral solid-band fill — what the MYC and height-function stencils actually read — as a Fortran-order (nx,ny,nz) array. get_vof()/'C' is the canonical field and carries EXACTLY 0 in solid cells; the fill is a stencil device regenerated at every ghost exchange (three passes with a shrinking depth budget, src/vof/cutcell.hpp), and it is what makes the wall look 90-degree neutral instead of perfectly non-wetting.
`vof_geometry`	vof_geometry(self, which: int) -> numpy.ndarray[dtype=float64] The cut-cell geometry the colour block runs on, as a Fortran-order (nx,ny,nz) array: 0 = the cell fluid fraction eps (4^3-subsampled, a multiple of 1/64), 1/2/3 = the openness of the +x/+y/+z face of each cell (the ADVECTOR's high-face convention, one cell shifted from the solver's ox/oy/oz), 4 = the classification (1 = SOLID, i.e. eps == 0 AND all six faces closed). On an ALL-FLUID solver (set_pressure_geometry, or no cut-cell pressure operator) there is no geometry block and this returns the TRIVIAL geometry the V1 transport kernels actually execute: eps = 1, openness = 1, classification = 0. It used to raise, which meant one diagnostic could not serve a packed scene and its all-fluid control.
`vof_interface_eps`	vof_interface_eps(self) -> float The wisp threshold set by set_vof_interface_eps.
`vof_last_courant`	vof_last_courant(self) -> float The interface-local Courant number of the step just taken (0 before the first step).
`vof_momentum_diagnostics`	vof_momentum_diagnostics(self) -> dict Census over THIS RANK's momentum control volumes: per-component min/max of the half-shifted colour C^c, the summed momentum rho^c u_c (the conservation census), the minimum rho^c before the floor, and the number of control volumes the floor touched.
`vof_rho_floor`	vof_rho_floor(self) -> float The absolute rho^c floor used by the last recovery.
`vof_step_parity`	vof_step_parity(self) -> int The sweep-permutation counter of the next colour advection.
`vof_timing`	vof_timing(self) -> dict WO-V9: cumulative VoF stage times in SECONDS since the last set_vof_timing/reset_vof_timing, on this rank. `steps` is the number of step() calls in the window, so every entry divides down to a per-step cost. `step`/`predictor`/`momentum_solve`/`projection` are the step's three coarse phases summed over the SAME window (the remainder of `step` is BC re-imposition, Picard bookkeeping and the scalars). `vof_advect` is the whole colour stage and `vof_bridge` the part of it spent in the G=2 <-> g=3 bridges and the C ghost policy; `vof_momentum_advect` is its V2b twin. `k_*` are the advector's own kernels, shared by the colour, momentum-consistent and consistent-energy drivers, so they are the per-KERNEL breakdown of whichever of the two stages is running: `k_reconstruct` the MYC + plicAlpha pass, `k_fluxes` the geometric face fluxes, `k_sweep` the update, `k_clip` the V5a cut-cell clip, `k_exchange` the g = 3 ghost exchange, `k_sweeps` the number of sweeps (3 per advect call). Returns zeros unless the timers are armed.
`vof_wisp_eps`	vof_wisp_eps(self) -> float The wisp tolerance currently in force (see set_vof_wisp_eps).
`vof_worklist`	vof_worklist(self) -> bool Whether the reconstruction-pass compaction is on.
`wall_area_probe`	wall_area_probe(self) -> list[float] Diagnostic: [sum_c x_cA_wall_x, sum_c y_cA_wall_y, sum_c z_c*A_wall_z] over all cut cells. Must equal -V_solid componentwise if the aperture wall-area vectors are right, so it separates a force-integral error in the GEOMETRY from one in the traction.
`wall_flux_imbalance`	wall_flux_imbalance(self) -> float Sum over cells of u_wall . A_wall -- the compatibility datum of the singular pressure problem. Exactly zero for a translating body in a periodic box; small but nonzero for rotation. Reported, not corrected.
`wall_slip_sandwich_cells`	wall_slip_sandwich_cells(self) -> tuple Per velocity component, the number of cut-cell AXES at which a one-cell fluid gap made the Navier closure inapplicable and the no-slip one was kept. Nonzero means part of the wall is silently no-slip; report it.

`SolverColocatedDiagnostics`¶

Developer tier of a Solver, reached as s.diagnostics: instruments (censuses, timers, budgets, ledgers), solver TUNING beyond the public driver selection, ablation switches and the raw zero-copy / MPI internals. Nothing here is needed to set up, run or read out a simulation. It holds a reference to the solver and no state of its own; the setters that CHANGE RESULTS say so in their docstring.

Method / property	Description
`advection_wall_velocity`	Whether the wall-velocity advection inputs are in force.
`anderson_accelerator`	anderson_accelerator(self, window: int = 5, mixing: float = 1.0) -> peclet.flow._flow._AndersonAcceleratorColocated A new Anderson accelerator of this solver's steady march (doc/steady_acceleration.md; the public driver is peclet.flow.march_to_steady). window in [1, 8] (default 5), mixing in (0, 1] (default 1.0). Raises on an unsupported configuration (a collocated scheme other than 'ghost', VoF, phase change, scalars, porous, variable properties, drag, cell forces, moving scenes, superficial velocity, pressure warm start, the balanced-force projection) and when the history, (2m+3) * n_fields * 8 * n_padded bytes, does not fit. Collective under MPI.
`balanced_force_failures`	balanced_force_failures(self) -> int How many balanced-force solves have failed over the solver's lifetime (see last_balanced_force_failed).
`bcast_from_root`	bcast_from_root(self, value: object) -> object Broadcast a value from rank 0 (identity in the single-rank module; mirrors the MPI API).
`comm_avoiding`	The communication-avoiding smoothing mode in force (see set_comm_avoiding).
`contact_angle_diagnostics`	contact_angle_diagnostics(self) -> dict The solid-band census of the CURRENT colour field on this rank: how many band cells each branch of pass 1 wrote ('contact_cells' the theta plane of the cell's own anchor, 'neighbour_cells' the mean of the anchor's MIXED neighbours' planes where the anchor itself is pure phase, 'pure_cells' the pure-phase continuation C_s = C_f, 'parallel_cells' an interface parallel to the wall where no rotation is defined, 'neutral_cells' the WO-Q neutral-mean fallback, 'unfilled_cells' left untouched), plus 'mean_apparent_angle' — the mean angle the fluid-only normal reported at the contact cells BEFORE the rotation, in degrees. That last number is the direct read-out of how far the fluid-side interface still is from the prescribed angle, measured on the fill's own data rather than on a post-processed shape. Rung V6 (WO-V6) adds, all zero unless set_contact_angle_dynamic / set_contact_angle_hysteresis is configured: 'dynamic_cells' (band cells the V6 pass produced an angle for), 'pinned_cells' / 'advancing_cells' / 'receding_cells' (the hysteresis branch census), 'mean_imposed_theta' and 'mean_apparent_theta' in degrees, 'max_Ca_cl' = max
`csf_diagnostics`	csf_diagnostics(self) -> dict Census of the CSF face force on THIS RANK: per component the max
`exchange_field`	exchange_field(self, name: str) -> None Fill a registered field's ghost cells (cross-rank + periodic under MPI; periodic single-rank).
`exchange_field_add`	exchange_field_add(self, name: str) -> None Add-reduce halo: fold ghost-layer deposits back onto their owner (cross-rank + periodic). The particle->grid deposition primitive for MPI CFD-DEM; single-rank non-periodic no-op.
`field_view`	field_view(self, name: str) -> numpy.ndarray[dtype=float64] Zero-copy view of a registered field's full padded buffer as a Fortran-order (nx+2g, ny+2g, nz+2g) array (g = ghost_width); host → NumPy, device → DLPack (CuPy).
`fluid_momentum_cells`	fluid_momentum_cells(self) -> list[int] Unmasked (fluid) staggered momentum cells per component -- the exact discrete datum of the reaction identity: at steady state, sum over bodies of F_c = f_c * N_c.
`get_cut_owner`	get_cut_owner(self) -> numpy.ndarray[dtype=int32] Per-inner-cell owning scene instance (the argmin behind the sampled SDF), (nx,ny,nz) int32; -1 before set_solid_from_scene. Moving geometry reads its wall velocity off the owner and resolved CFD-DEM posts the hydrodynamic force back to it.
`hydro_force_torque_traction`	hydro_force_torque_traction(self) -> numpy.ndarray[dtype=float64] DIAGNOSTIC -- not the force on the bodies; that is Solver.hydro_force_torque_reaction(). The reconstructed surface-TRACTION loads, shape (4, n_instances, 3): [0] force, [1] torque about the instance centre, [2] the pressure part, [3] the viscous part ([0] == [2] + [3]). The cut-cell surface integral of (-p I + mu(grad u + grad u^T)) against the exact aperture wall-area vector, with a central-difference velocity gradient that spans the wall. It UNDER-READS THE VISCOUS PART: traction / reaction drag = 0.685-0.730 over phi 0.008-0.45 and N 24-128 (measured 2026-09-30), i.e. ~30 % low at ANY resolution -- refinement does not fix it. The torque is low by the same mechanism: a steadily spinning sphere reads 0.59-0.72 of 8 pi mu a^3 Omega / (1 - phi) at R/h 6-12 (the reaction: 1.019-1.027). What it offers that the reaction does not: the pressure/viscous split, and it runs where the reaction refuses (the collocated solver, porous, variable properties, domain BCs, implicit advection). A consistent traction (from the momentum operator's small-cell-robust wall reconstruction) is planned. Was Solver.hydro_force_torque() (deprecated). Atomics: tolerance-reproducible, not bitwise.
`instance_center_pinned`	instance_center_pinned(self, i: int) -> bool True if instance i's centre of rotation is PINNED (an explicit finite centre in the encoding, or set_instance_motion(center=...)); False if it follows the body's translation (NaN in the encoding, the builder's default; legacy all-zero raw arrays).
`last_balanced_force_failed`	last_balanced_force_failed(self) -> bool Whether the last step's balanced-force solve FAILED (driver breakdown or a non-finite result). A failed solve keeps P and p_balanced unchanged for that step; the first one warns on stderr.
`last_balanced_force_iterations`	last_balanced_force_iterations(self) -> int Pressure-driver iterations of the last step's balanced-force solve (set_balanced_force_projection); 0 when it is off, the force divergence was zero, or the previous split already met the stop.
`last_momentum_residual`	last_momentum_residual(self) -> float max over components of max
`last_outer_iterations`	last_outer_iterations(self) -> int Return the outer-iteration count from the last step().
`last_pressure_iterations`	last_pressure_iterations(self) -> int Return the pressure-solver iteration count from the last step(). A solve that BROKE DOWN (non-finite preconditioner output) reports the iteration CAP, so the usual rule-3b 'a capped pressure solve invalidates the run' check sees it; pressure_solve_failed() distinguishes the two.
`last_step_timers`	last_step_timers(self) -> dict Per-phase wall times (seconds, THIS rank, device-fenced) of the last step(): 'predictor' (ghost fills + RHS/advection/stencil builds), 'momentum' (implicit-diffusion solves), 'projection' (cut-cell pressure projection), 'step' (whole step). 'pressure_allreduce'/'pressure_allreduce_count' = time in / number of global reductions (MPI_Allreduce) inside the pressure solve — the latency-bound term of the distributed solve (0 on a single rank).
`last_superficial_velocity_shift`	last_superficial_velocity_shift(self) -> float The uniform velocity shift set_superficial_velocity applied at the end of the last step() (physical velocity units; 0 when it is off). Reading it synchronises with the device.
`moving_instance_cut_cells`	moving_instance_cut_cells(self) -> numpy.ndarray[dtype=float64] Per instance: inner cells on this rank it owns that touch a fractional face aperture, recounted whenever any instance moves. ZERO for a moving instance means its surface sits on grid planes (or is sub-cell) and its wall velocity is silently inert (set_solid_from_scene warns). Empty when no instance moves.
`moving_instance_degenerate_points`	moving_instance_degenerate_points(self) -> numpy.ndarray[dtype=float64] Per instance: staggered velocity points where the sampled SDF is EXACTLY zero -- a face on a lattice plane. Those points are fluid to the mask and not ghosts to the cut-cell fold, so a moving body's datum never enters there (set_solid_from_scene warns). Empty when no instance moves.
`multigrid_aspect_threshold`	The anisotropic-coarsening threshold in force (see set_multigrid_aspect_threshold).
`outflow_backflow`	outflow_backflow(self) -> dict Outflow REVERSAL census over the outflow faces (collective under MPI): 'max_reverse' = the largest reversed normal velocity max(0, -u.n), 'fraction' = reversed / all outlet faces, 'energy_influx' = sum over the reversed faces of rho
`outflow_rho_correction`	outflow_rho_correction(self) -> bool Whether the 1/rho_f factor is applied to the outflow face correction (WO-R item 4, reversed by WO-R2's operator fix; default True).
`periodic_image_overlap_cells`	periodic_image_overlap_cells(self) -> int Cells on this rank whose solid/fluid sign was decided by a periodic IMAGE of an instance wider than the box (set_solid_from_scene warns when nonzero): the scene evaluates the UNION of images, so a slab wider than the box refills any cavity carved from it. 0 when no instance is that wide or the images agree.
`phase_change_area`	phase_change_area(self) -> str The set_phase_change_area mode in force, as its name.
`phase_change_budget`	phase_change_budget(self) -> dict WO-P3f: the ENERGY BUDGET of the LAST energy solve (units: J for the enthalpies, W for `q_gfm`; cell volume = 1). h_open / h_open_new sum rho c_p (T - T_sat) over the UNMASKED cells, before / after the solve h_liquid, h_masked the same over the pure-liquid and the interfacial cells d_overwrite sum rho c_p (dval - T) over masked cells: what the Dirichlet rows inject when they overwrite the transported temperature (negative = destroyed) d_overwrite_new the part of it on cells that were NOT masked last step e_enter / e_leave the enthalpy carried OUT of / INTO the solved set by the cells that changed class this step q_gfm the heat the plane-anchored rows deliver INTO the unmasked set (negative while a bubble grows); the energy equation's own interfacial flux, to be compared with the regression's `mdot h_lv A_Gamma` n_enter_/n_leave_ the class-change census The discrete balance the entries close over one solve is `sum rho c_p (T^{n+1} - T*)/dt = q_gfm + (domain boundary flux) + (solve residual)`.
`phase_change_carry_ledger`	phase_change_carry_ledger(self) -> dict WO-P3f: the last `set_phase_change_carry_conserve` pass — the enthalpy actually handed back ('deposited') and the enthalpy of the interfacial cells that had NO neighbour left in the solve ('lost', which stays destroyed). Both 0 when the option is off.
`phase_change_diagnostics`	phase_change_diagnostics(self) -> dict Per-rank phase-change census of the LAST step: mdot extrema/mean and the interfacial cell count, the total PLIC interface area, the liquid volume the regression removed, the clip-and-redistribute ledger (
`pressure_coarse_ghost`	Whether the Neumann coarse ghost is in force (see set_pressure_coarse_ghost).
`pressure_exact_residual`	pressure_exact_residual(self) -> bool Whether the exact level-0 operator apply is in force (see set_pressure_exact_residual).
`pressure_mg_level_ratios`	pressure_mg_level_ratios(self) -> list[list[int]] The pressure multigrid's per-level coarsening ratio, one (rx, ry, rz) per level (read-only). On an isotropic domain every coarsenable axis halves, as it always has. On an ANISOTROPIC one (extent giving different spacings per axis) the hierarchy coarsens the FINEST coarsenable axis and defers one that is already at least set_multigrid_aspect_threshold (default 2) times coarser, so the coarse operators stay close to isotropic and the point smoother keeps its rate -- e.g. cells (N, 2N, N/2) on a cube gives (1,2,1), (2,2,1), (2,2,2), ... Empty until set_solid/set_pressure_geometry has built the operator. See flow/doc/anisotropic_metric.md §5.
`pressure_solve_failed`	pressure_solve_failed(self) -> bool Did the last pressure solve break down on a non-finite recurrence scalar (a preconditioner or operator that produced NaN/Inf)? The driver then zeroes the correction and continues -- so the projection was handed NOTHING for that step and the run is invalid. Until 2026-09-04 the only trace was one line on stdout ('preconditioner produced non-finite z') while last_pressure_iterations() reported 0, i.e. a perfectly healthy-looking solve (measured: examples/pore-scale-imbibition, a 3.1-cell throat). Set set_pressure_strict(True) to raise instead.
`pressure_strict`	Whether a non-finite pressure preconditioner output raises (see set_pressure_strict).
`reaction_budget_terms`	reaction_budget_terms(self) -> numpy.ndarray[dtype=float64] Diagnostic decomposition of the reaction identity, shape (2, 3): [0] the unsteady sum sum_i (rho/dt)(u_i - u^n_i) and [1] the advective sum sum_i A_i, both over the FLUID momentum cells of each component. The full discrete identity is sum_bodies F_c = f_c*N_c + sum fb + [1]_c - [0]_c; the Stokes form drops both because they vanish at steady state with advection off. sum_i A_i is the advection operator's net momentum flux through the cut walls -- an O(h) property of that operator, not of the budget.
`rebalance_by_weights`	rebalance_by_weights(self, weights: collections.abc.Sequence[float]) -> int Dynamic load balancing: redistribute the solver's state onto the weighted ORB of per-cell weights (global x-fastest, gnxgnygnz), its split planes aligned to 2^a cells with a the largest alignment whose weight imbalance stays within 1.05, so the pressure multigrid keeps coarsening in place for a levels. Returns that alignment (1 = the plain weighted ORB). Pass fluid work + gamma*particle_count; a coupled dem migrates onto the SAME partition with migrate_to_weights(weights, align=). State-preserving (bit-exact at np=1, reduction floor at np>1).
`reset_vof_timing`	reset_vof_timing(self) -> None Zero the VoF stage timers and the step counter without disarming them.
`set_advection_wall_velocity`	set_advection_wall_velocity(self, on: bool) -> None A0: fill the advection inputs' masked (solid) rows with the instantaneous WALL velocity instead of zeros, so a MOVING body's advective term sees the body's own motion. DEFAULT True; False is the pre-A0 ablation (measurement only). Inert on a static scene and on the collocated grid.
`set_aperture_order`	set_aperture_order(self, order: int) -> None Face-aperture estimator (DEFAULT 2 since 2026-08-26): 1 = legacy one-sample linear model, 2 = marching-squares (5 samples/face, O(h^2); removes the convexity bias measured at +0.59/+0.27% bed permeability at R=8/12 -- see doc/collocated_paper_plan.md row 51). For analytic geometry, exact apertures via set_openness_override are better still. Call before set_solid (a later call raises).
`set_comm_avoiding`	set_comm_avoiding(self, mode: str) -> None Communication-avoiding red-black smoothing (multi-rank): exchange a 2-deep ghost layer once per red-black PAIR instead of a 1-deep layer before every colour, redundantly re-smoothing the first colour's ghost ring so the second colour reads exactly what a per-colour exchange would have delivered. BIT-IDENTICAL, at half the halo events (measured: np=32 weak efficiency 35% -> 62%). 'both' (DEFAULT)
`set_contact_angle_clamp`	set_contact_angle_clamp(self, lo: float, hi: float) -> None The clamp on the Cox-Voinov cube, in degrees (default 1 / 179). The cubic has no solution beyond the maximum receding capillary number (film entrainment) and the fill's plane construction degenerates at 0/180.
`set_contact_angle_pivot`	set_contact_angle_pivot(self, mode: str) -> None ABLATION -- how the theta-plane is anchored in the fluid cell. 'volume' (DEFAULT) match the anchor cell's liquid volume with plicAlpha; 'centroid' pass through the PLIC centroid p_f (the Afkhami-Bussmann / Basilisk contact.h rule); 'projected-centroid' the work order's c = p_f - sdf(p_f) n_w; 'contact-line' the contact line on the wall. 'volume', 'centroid' and 'contact-line' are exactly idempotent (1e-15 on gate G0a); 'projected-centroid' is NOT -- projecting the centroid along n_w shifts the plane by -sdf(p_f) cos(theta), measured 0.26 in cell fraction at theta = 60, with the wrong sign (it removes liquid from the band for a wetting angle).
`set_contact_angle_smoothing`	set_contact_angle_smoothing(self, on: bool) -> None ABLATION - the 3-point in-wall mean of U_cl (default ON). A MAC velocity next to a wall is noisy cell to cell and the cube root of the Cox-Voinov relation puts that noise straight into the imposed angle; off measures what the smoothing is worth.
`set_csf_mode`	set_csf_mode(self, mode: str) -> None ABLATION. 'face' (default, the only production mode) evaluates the surface-tension force as sigmakappa_f(C(i)-C(i-s))/h at the face -- the projection's own gradient operator. 'cell' evaluates a CELL-CENTRED sigmakappagrad(C) and face-interpolates it with the arithmetic mean, exactly as the per-cell body-force machinery carries a rho*g field: consistent, convergent, and wrong for surface tension, because the result is not in the range of the operator the projection inverts. Shipped so the difference is a measured number (see the vof_surface_tension ctest), not an argument.
`set_deferred_correction`	set_deferred_correction(self, on: bool) -> None Deferred-correction advection: True (default) = 2nd order (implicit FOU + explicit high-order correction, the high-order scheme being SOU by default or Koren TVD via set_advection_scheme); False = pure implicit FOU (1st-order upwind, more dissipative but unconditionally stable at sharp shear layers).
`set_density_mode`	set_density_mode(self, mode: str = 'variable') -> None Enable variable density: binds the 'rho' field (get/set_field('rho'), created seeded with set_rho's value if absent) into the momentum time term, the advection weight, the per-cell body force (face-interpolated), and the pressure projection (face coefficient opennessrho0/rho_f with the matching 1/rho_f velocity correction; rho0 = set_rho's value, so a uniform field reduces exactly to the constant solver). A closure targeting 'rho' (e.g. a linear mixture of a transported phase fraction) enables this automatically. For gravity, register a closure force_z = linear(rho, params=[0, -g]). COLLOCATED (SolverColocated, rung V8): supported since 2026-09-02, ALL-FLUID only (set_pressure_geometry; an immersed solid, the ghost projection and set_rho_face_harmonic throw). The face coefficient and the face correction are the same as on the staggered grid, applied to the ABC projection's MAC face field, which is the MOMENTUM-weighted map (rho_L u_L + rho_R u_R)/(rho_L + rho_R) of the cell velocities; the CELL correction is the average of the two face corrections of each axis (never a cell-centred grad(phi)/rho_c). The pressure, the CSF and set_body_force (a mean pressure gradient) enter the implicit momentum predictor as the matching face integral rho_c avg_faces((f - grad_f P)/rho_f); a per-cell force enters at the cell value. The pair is adjoint, so the step is stable at every dt and density ratio and its steady state is dt-independent (doc/collocated_varrho_forces.md). Rated to density ratio ~10 for cases WITH MOTION (a translating drop exceeds the VoF CFL cap at ratio 100 and 1000; momentum consistency needs Favre face states and is not in this rung); exact at rest at ratio 1000.
`set_exact_crossings`	set_exact_crossings(self, t: collections.abc.Sequence[float]) -> None Analytic-SDF capability: exact wall-crossing fractions overriding the linear-interpolated theta in the momentum cut-cell overlay AND the ghost-projection closures. Flat array of 9nxnynz values, blocks [(c3+k)]: component c's staggered point at inner cell i toward its +k neighbour; NaN = no crossing. Call BEFORE set_solid (a later call raises); empty list clears. Single-rank, staggered momentum placement.
`set_exact_crossings_from_scene`	set_exact_crossings_from_scene(self) -> None Compute EXACT wall crossings from the installed scene, on device, on every rank -- the in-solver replacement for set_exact_crossings + scripts/exact_apertures_spheres.py. Bisection, not Newton: only SIGN correctness is guaranteed for bound-only leaves.
`set_face_interp`	set_face_interp(self, mode: int) -> None The two INTERMEDIATE rungs of the Basilisk embed.h port on the collocated grid (doc/history/collocated_embed_port_plan.md), reachable for the ladder that measures what each piece of the complete port ('embed' = mode 7) is worth: 5 = the FV momentum operator with the true-normal dirichlet_gradient wall drag as a defect correction, on the wall-aware face map with the transpose-gradient pressure force; 6 = plus the openness-weighted -grad(P) predictor and cell correction, on the plain face map (7 adds the wall-aware face map and the solid-cut-cell sliver mask). CHANGES RESULTS; no effect on the staggered solver. Modes 1-4 and 10-13 (ablations) were deleted with their kernels at 1.0.0.
`set_fluid_only_constraint`	set_fluid_only_constraint(self, mode: str) -> None Fluid-only pressure constraint (route 2b, a collocated mechanism instrument; call before set_solid): 'filter' = Design A (openness filter), 'star' = Design B (SPD Kron star elimination), 'off'. Any setting other than 'off' disables the AUTO scheme.
`set_fresh_cell_seed`	set_fresh_cell_seed(self, on: bool) -> None Fresh-cell policy for MOVING geometry. A point a body has just uncovered inherits whatever the solid held there -- zero, or a stale masked value -- which is one of the two textbook mechanisms behind spurious force oscillations in a moving-boundary IBM. True seeds it with the LOCAL WALL VELOCITY instead, so it starts moving with the surface that released it. Bounded (no extrapolation), no new field, and exactly the old behaviour when nothing moves. DEFAULT TRUE since 2026-08-30: measured on an oscillating sphere translating through the grid, it removes a resolution-INDEPENDENT +2.6..2.9% drag bias, cuts the spurious force oscillation 20-50x to within 17% of the non-moving floor, and improves the resolved CFD-DEM loop's total-momentum conservation 95x. Pass False for the pre-2026-08-30 behaviour.
`set_ghost_projection`	set_ghost_projection(self, on: bool, matrix_order: int = 2, rhs_order: int = 2) -> None QUARANTINED 2026-08-18 (verification only, unsupported): superseded by the gauge-exact collocated scheme, which matches its accuracy at 5-6x lower cost. Kept as the independent second discretization behind the cross-IBM physics gate. Enabling it on the collocated grid silently selects the plain face map, since it owns the operators the gauge-exact scheme replaces. EXPERIMENTAL directional ghost-cell projection (second staggered IBM): point-based FD divergence with wall-anchored directional closures instead of the openness-weighted cut-cell projection; solved by MG-preconditioned BiCGStab. Call BEFORE set_solid (a later call raises). Closure orders (1=linear, 2=quadratic): (matrix_order, rhs_order) = (2,2) full quadratic 13-point matrix; (1,1) linear 7-point; (1,2) mixed/deferred — 2nd-order steady constraint on a 7-point matrix. Collocated: the same closures/matrix on the face-averaged field, plus a directional (one-sided 2nd-order) cell gradient for the -grad(P) predictor and cell correction; requires face_interp 0. v1: single-rank, periodic + IBM, stationary walls; incompatible with porous/variable-rho/domain-BC/Chebyshev.
`set_incremental_pressure`	set_incremental_pressure(self, on: bool) -> None Toggle the rotational incremental-pressure projection.
`set_multigrid_aspect_threshold`	set_multigrid_aspect_threshold(self, theta: float) -> None Anisotropic-coarsening aspect threshold theta shared by the pressure and velocity multigrids (doc/anisotropic_metric.md §5.1): an axis is coarsened only while its cell size stays within theta of the smallest one, so a stretched grid semi-coarsens instead of building levels whose operator is dominated by one direction. DEFAULT 2.0; must be > 1; 1e9 restores full coarsening on every axis (the §8.5 item (c) ablation). Only ever consulted on an ANISOTROPIC metric, so it is inert on cubic cells.
`set_openness_override`	set_openness_override(self, ox: collections.abc.Sequence[float], oy: collections.abc.Sequence[float], oz: collections.abc.Sequence[float]) -> None Analytic-SDF capability: exact face-aperture (openness) fields for the cut-cell projection, overriding the sampled-SDF openness. Inner nxnynz arrays, x-fastest; ox[i] = fluid fraction of the -x face of cell i. Call BEFORE set_solid (a later call raises); empty ox clears. Single-rank.
`set_outer_iterations`	set_outer_iterations(self, n: int) -> None Set the number of Picard/outer iterations per step.
`set_outer_tolerance`	set_outer_tolerance(self, tol: float) -> None Set the outer (Picard) convergence tolerance.
`set_outflow_rho_correction`	set_outflow_rho_correction(self, on: bool = True) -> None The 1/rho_f mobility factor on the HIGH-side outflow face correction. DEFAULT TRUE since WO-R2; pass False for the ablation. A projection correction cancels the discrete divergence only if it uses the SAME face coefficient the operator row used. Until WO-R2 the multigrid re-imposed the literal openness 1.0 at every Dirichlet domain face, overwriting buildRhoCoeff's open_f*rho0/rho_f, so the plain phi difference was the consistent partner and WO-R measured this factor making things seven orders WORSE. WO-R2 fixed the operator (CutcellMG::setOutflowCoefficient) and the verdict inverted. Stratified duct, ratio 10, 5 steps, max
`set_phase_change_area`	set_phase_change_area(self, mode: str) -> None WO-P3c: WHICH GEOMETRY the interfacial area A_Gamma comes from. A_Gamma sets the plane shift dV = mdot A dt / rho_l and the divergence source S = mdot A (1/rho_g - 1/rho_l), so a bubble grows as int mdot dA and a biased area is a biased growth rate. 'plic' (DEFAULT, rungs P0/P1): plicArea =
`set_phase_change_budget`	set_phase_change_budget(self, on: bool = True) -> None WO-P3f INSTRUMENT: turn on the energy budget of the phase-change energy solve. Costs one extra cell field and two reductions per energy solve, and is skipped entirely when off (the solve is then bit-identical). Read it with `phase_change_budget()`. What it exists for: interfacial cells are Dirichlet rows, i.e. they are OUTSIDE the energy solve, so the set over which enthalpy is conserved changes membership every step as the interface sweeps. A liquid cell that becomes interfacial LEAVES that set carrying its superheat rho c_p (T - T_sat) and an interfacial cell that becomes pure RE-ENTERS it carrying `pcCarriedValue`; neither transfer appears in the latent-heat book-keeping.
`set_phase_change_carry_conserve`	set_phase_change_carry_conserve(self, on: bool = True) -> None WO-P3f OPTION (default OFF): make the interfacial cells' per-cell Dirichlet OVERWRITE enthalpy-conserving. An interfacial cell's row is the identity `T = dval`, so whatever the geometric energy transport left there is discarded every step. Over a cell's interfacial lifetime that telescopes to rho c_p (T_entry - dval_exit): a liquid cell the interface sweeps enters with its superheat and leaves as vapour at T_sat, and the difference goes nowhere. Measured on Scriven 128^3 by `phase_change_budget()['d_overwrite']`: -0.7 % of mdot h_lv A_Gamma at Ja = 0.5 and -4.3 % at Ja = 2, one-signed. With this on, that enthalpy is handed to the interfacial cell's face neighbours that are still in the solve, weighted by n_d^2 (the clip-and-redistribute allocation), as a fixed-order gather so it is decomposition-independent. Which side receives is decided per axis by which pure neighbour deviates from T_Gamma in the same direction as the interfacial cell itself, i.e. the phase the enthalpy came from -- the superheated liquid on an evaporating bubble, the superheated vapour on the Stefan problem.
`set_phase_change_curvature_distance`	set_phase_change_curvature_distance(self, on: bool = True) -> None WO-P3g item 3 (default OFF): measure the GFM row's theta and the one-sided fits' sample distances to the CURVED interface, with the mean curvature taken PER CELL from the V3 curvature cascade (the same kappa surface tension uses; positive for a convex blob of liquid, i.e. -2/R for a gas bubble). Supersedes `set_phase_change_fit_curvature`, which prescribes one curvature for the whole field; where both are set the cascade wins.
`set_phase_change_deposit_fallback`	set_phase_change_deposit_fallback(self, on: bool = True) -> None WO-P3f open item 6 / WO-P3g (OFF by default): give an interfacial cell whose two along-the-normal deposit candidates are BOTH still interfacial a target from the 5^3 box instead of leaving the divergence source in place. A cell that keeps its source carries div(open u) = S on its own faces, so Weymouth-Yue advects the colour with a field that is not the liquid velocity -- read it out with `phase_change_diagnostics()['band_div']` and the 'fallback_cells' count.
`set_phase_change_energy_muscl`	set_phase_change_energy_muscl(self, on: bool) -> None MinMod-limited donor reconstruction of the face temperature in the CONSISTENT energy flux, instead of the default plain donor-cell (first-order upwind) value. OFF by default. The consistency is what the geometric flux buys; its first-order upwind numerical diffusion
`set_phase_change_energy_order`	set_phase_change_energy_order(self, order: int) -> None WO-P3g: the ORDER of the interfacial energy operator. 1 (the shipped WO-P23...P3f scheme) or 2. `order = 2` turns on, TOGETHER, the four pieces WO-P3f's instruments indicated: 1. `set_phase_change_mdot_operator(True)` -- mdot is the energy operator's OWN interfacial flux q/(h_lv A) instead of a separate least-squares fit, so the heat the energy equation loses and the mass the regression produces are one discrete quantity (and the interfacial AREA cancels out of the mass balance entirely); 2. `set_phase_change_gfm_order(2)` -- the Gibou-Fedkiw three-point ghost-fluid row (2/((1+theta) theta), 2/(1+theta)) instead of the two-point (1/theta, 1), which is exact on a quadratic profile at every theta; 3. `set_phase_change_curvature_distance(True)` -- the row's theta and the one-sided fits' sample distances are measured to the CURVED interface, with kappa taken per cell from the V3 curvature cascade; 4. `set_phase_change_carry_conserve(True)` -- WO-P3f's enthalpy-conserving Dirichlet overwrite. Why together and not one at a time: WO-P3f measured the shipped scheme's 1 % Scriven error to be the residue of a CANCELLATION between the fit's +6 % curvature bias, the two-point row's -5 % flux deficit and the overwrite's -0.7...-4.3 % enthalpy destruction, so repairing any ONE alone makes the gate worse.
`set_phase_change_fit_curvature`	set_phase_change_fit_curvature(self, kappa: float) -> None WO-P3f INSTRUMENT (default 0, bitwise inert): prescribe the interface curvature `kappa = div(n)` that the one-sided temperature fits use to measure a sample's distance to the CURVED interface instead of to the interfacial cell's tangent plane (`vof::pcCurvedDistance`). For a spherical GAS bubble of radius R, whose PLIC normal points inward, kappa = -2/R. Why it exists: the tangent-plane distance makes the fit FIRST ORDER in h/R on a curved interface, and every off-axis sample is hotter than the plane model expects, so the fitted dT/dn and with it mdot come out HIGH. Measured a priori on an exact sphere with an exactly linear profile: +19.2 / +12.1 / +8.8 / +6.2 % at R = 6 / 10 / 14 / 20, observed order 0.91-0.98 in h/R. This entry point takes a PRESCRIBED kappa so that bias can be measured against a known geometry; it is not a curvature estimator.
`set_phase_change_gfm_order`	set_phase_change_gfm_order(self, order: int) -> None WO-P3g item 2 (default 1): the order of the one-sided (ghost-fluid) Dirichlet row. 1 = the shipped two-point form k o (T_Gamma - T_i)/theta; 2 = Gibou-Fedkiw's three-point form through (T_behind, T_i, T_Gamma), which reproduces a quadratic temperature profile exactly at every theta.
`set_phase_change_mdot_operator`	set_phase_change_mdot_operator(self, on: bool = True) -> None WO-P3g item 1 (default OFF): take mdot from the energy operator's own interfacial flux -- the sum of the ghost-fluid rows' Dirichlet couplings evaluated with the converged T -- instead of the one-sided least-squares fit, which stays as `phase_change_diagnostics()['mdot_fit']`.
`set_phase_change_plane_dirichlet`	set_phase_change_plane_dirichlet(self, on: bool) -> None The PLANE-ANCHORED (ghost-fluid) interfacial Dirichlet condition. ON by default. Rungs P0/P1 pinned the whole interfacial CELL at T_G, so the numerical thermal boundary sat at the cell CENTRE while the mass-flux gradient is fitted from the PLIC PLANE — a mismatch of up to half a cell that CHANGES SIGN as the interface sweeps through a cell, and the first-order component of the P1 Stefan error (WO-P01 finding 6). With this on, the condition becomes PER FACE: for a pure cell i whose neighbour j is interfacial, that face's row carries kopen(T_i - T_G)/theta with theta the distance in cells from i's centre to j's PLIC plane along the face's own axis, and the interfacial cell's own value is never read by any neighbour. A per-CELL value (the literal reading of the work order) does NOT work and the measurement is in the findings: giving the interfacial cell the value the one-sided profile takes at its centre is right for the side the fit came from and wrong for the other side — on the P1 Stefan ladder it heats the saturated liquid through the interfacial cell's liquid-side face, which gives the liquid a spurious gradient that feeds straight back into mdot: +6.20/+5.62/+5.42 % at N = 64/128/256, order 0.10, against +1.31/+0.59/+0.20 % and order 1.37 for the cell-centre pinning it was meant to improve on. The per-face form has no such asymmetry (the liquid-side face reads T_G at its own theta, i.e. exactly zero flux for a saturated liquid). False restores the rung P0/P1 behaviour bit-for-bit.
`set_phase_change_quadratic_fit`	set_phase_change_quadratic_fit(self, on: bool) -> None Fit the one-sided interfacial temperature gradients QUADRATICALLY through the interface value (T - T_G = G phi + Q phi^2) instead of linearly. This is VOF_PLAN section 9 item 1's Aslam quadratic extrapolation in least-squares form — the same 5^3 pure-cell samples, the same Malan collinearity weights, one more basis function, no PDE sweeps. Once the plane-anchored Dirichlet has removed the cell-centre mismatch, the linear fit's O(T'' h) curvature bias is the leading error of the rung: its samples start about one cell from the plane and reach two and a half, so a curved profile tilts the straight line through the origin. ON by default (WO-P23): with the plane-anchored Dirichlet it takes the P1 Stefan interface position from +0.195 % to +0.003 % at N = 256, and the mdot kernel itself from order 1.1 to order 2.0. set(False) is the ablation.
`set_porous_conservative`	set_porous_conservative(self, on: bool) -> None eps-conservative porous momentum + projection pair (default True): time term (eps_f rho/dt) u, eps rho-weighted advective form, projection coefficients open*(eps rho idt)/(eps rho idt+beta) with matching correction. False = the legacy plain-u pair (A/B only; it lets the projection drag gas with the moving porosity at zero inertia cost — a spurious late-time energy source in clustering flows).
`set_porous_deps_dt`	set_porous_deps_dt(self, on: bool) -> None Include (default True) or drop the d(eps)/dt source in the porous projection RHS. Drop it to enforce div(eps u)=0 when the per-cell eps deposit's time-derivative is too jagged and destabilizes the eps-weighted pressure solve.
`set_pressure_bottom_solver`	set_pressure_bottom_solver(self, engine: str) -> None Engine of the AGGLOMERATED pressure bottom (set_pressure_bottom 'auto' / 'agglomerated'). 'auto' (DEFAULT): on a GPU backend, single rank, the singular (periodic / wall) operator, a bottom of at most 8192 cells with at most 64 fluid components (solids allowed) and an axis whose planes hold at most 192 cells, 'direct': one device launch runs flexible CG in FP64 (inner tolerance 1e-5, cap 100) preconditioned by a block-tridiagonal FP32 direct factor of the bottom (planes along one axis, explicit Schur-complement inverses), refactored once per operator change -- no host transfer; everywhere else the host GraphAMG solve. 'direct' forces the device engine and raises at the bottom solve, naming the failed condition, where it is not eligible; 'algebraic' forces GraphAMG (the A/B instrument). Host backends always run GraphAMG under 'auto'. doc/vof_step_performance_design.md §13.
`set_pressure_coarse_ghost`	set_pressure_coarse_ghost(self, on: bool) -> None Zero-gradient (Neumann) coarse ghost on wall/inflow faces before the pressure multigrid's prolongation -- the WO-H symmetry repair. DEFAULT True. False restores the pre-2026-08-30 periodic-wrap ghost, which reinstates the asymmetry that stalls MG-PCG on every 3-D wall-bounded grid: a MEASUREMENT ABLATION, never a production setting. Inert on periodic / IBM problems, where the two are the same.
`set_pressure_exact_residual`	set_pressure_exact_residual(self, on: bool = True) -> None Apply the level-0 pressure operator EXACTLY (matrix-free, double, flux form) in the residual and the Krylov matvec instead of reading the float band storage. P1 of the suite defect-correction campaign (suite/docs/archive/DEFECT_CORRECTION_PLAN.md). Per solver, OFF at construction. enable_vof() turns it ON, because a two-phase coefficient contrast is exactly what amplifies the float operator's broken row-sum identity A*1 = 0. Measured on Hysing case 2 (64x128x4, adaptive dt, nvidia-cuda): max
`set_pressure_graph_amg`	set_pressure_graph_amg(self, on: bool) -> None Solve the pressure MG's coarsest level with an agglomerated mesh-agnostic algebraic multigrid (core GraphAMG), decomposition-agnostic: with levels=1 this gives a mesh-independent pressure solve that works under a WEIGHTED ORB (where the geometric coarse levels can't cleanly coarsen). Applied at the next set_solid.
`set_pressure_mean_removal`	set_pressure_mean_removal(self, scope: str) -> None Nullspace (mean) removal scope in the pressure solve: 'fine' (DEFAULT — only the projections the Krylov iteration needs: rhs/residual, the fine-level V-cycle exit, the final iterate; ~3x fewer global-reduction latency hits per iteration, the measured winner of the multi-GPU ablation) or 'all' (legacy — every V-cycle level + after every matvec). Iteration counts are identical (A preserves mean-freeness); results equal within solver tolerance, not bit-identical.
`set_pressure_strict`	set_pressure_strict(self, on: bool) -> None Raise instead of report when the pressure preconditioner returns a non-finite correction. DEFAULT False: the shipped behaviour caps the solve, prints one line, and raises the pressure_solve_failed() flag, so a check can catch it without changing control flow. True turns the same event into an exception.
`set_pressure_underrelax`	set_pressure_underrelax(self, omega: float) -> None Pressure under-relaxation factor omega_p in (0,1] for the incremental accumulation (MFIX §10.1); 1.0 = off (default). <1 damps the incremental predictor overshoot on stiff porous+drag.
`set_pressure_warmstart`	set_pressure_warmstart(self, on: bool) -> None Seed each pressure solve from the previous step's phi (default off).
`set_property_mode`	set_property_mode(self, mode: str = 'variable', harmonic: bool = False) -> None Enable variable-coefficient momentum (variable viscosity): mode 'variable' binds the 'mu' field (get/set_field('mu')) into the diffusion operator; 'constant' reverts. harmonic = harmonic face-viscosity mean (continuous shear stress across a jump) vs arithmetic. A closure targeting 'mu' enables this automatically. The incremental-rotational pressure scheme (large-dt / steady-Stokes) stays active — see set_variable_rotational.
`set_rho_face_harmonic`	set_rho_face_harmonic(self, on: bool = True) -> None Use the HARMONIC instead of the arithmetic face mean of rho in the pressure projection (both the operator coefficient and the velocity correction, so the projection stays exact). DEFAULT OFF, and it should stay off: the arithmetic mean of rho IS the harmonic mean of the mobility 1/rho — the series-correct choice for a normal flux — and it is what makes the discrete hydrostatic balance EXACT, because the momentum time term and the face body force interpolate rho arithmetically and are not switched by this flag. Shipped as a measured knob for the coefficient-coarsening question (VOF_PLAN S3), not as an alternative scheme.
`set_rotational_filter`	set_rotational_filter(self, on: bool, eps: float = 0.05) -> None Experimental filtered rotational update: smooth div(u) (mask-aware axis-wise 1-2-1, one-sided into the fluid at walls) before accumulating -mudiv into P. Keeps the O(1) pressure-relaxation gain, removes the checkerboard feedback channel.
`set_rotational_pressure`	set_rotational_pressure(self, on: bool) -> None PM I ablation (Guy-Fogelson 2005): False drops the rotational -mudiv(u) term from the incremental pressure accumulation (constant-mu path only). Default True = shipped rotational (Timmermans) update.
`set_rotational_wall_weight`	set_rotational_wall_weight(self, w0: float) -> None Wall-banded rotational blend: at fluid cells with a solid axis-neighbour use P += (rho/dt + w0mu/dx^2)phi - (1-w0)mudiv(u*); bulk keeps the full rotational update. Stabilizes the boundary rows without slowing bulk pressure relaxation, and keeps them relaxing at dt->infinity. 0 (default) = off.
`set_rotational_weight`	set_rotational_weight(self, w: float) -> None Under-relax the rotational term: P += ctphi - wmudiv(u). 1 = shipped, 0 = PM I; small w raises the boundary-mode stability threshold ~1/w at ~1/w slower smooth-mode pressure relaxation. phi=0 stays the unique fixed point for any w>0 at every dt.
`set_uf_advection`	set_uf_advection(self, on: bool) -> None Advect the COLLOCATED momentum with the projected, divergence-free MAC face field uf/vf/wf of the last projection instead of the un-projected cell->face average 0.5*(u_i+u_j). DEFAULT True -- the shipped scheme: doc/flow_colocated_plan.md step 3, the Almgren-Bell-Colella prescription (the field the projection just made solenoidal IS the conservative advective flux), and what the FOU operator's conservative row-sum identity needs. False is the ablation that restores the phase-2 average (measurement only; see doc/uf_advection.md). Inert on the staggered grid (the stored velocity already IS the projected face velocity) and before the first projection or set_velocity has built a face field.
`set_variable_rotational`	set_variable_rotational(self, mode: str = 'min', chi: float = 1.0) -> None Rotational-pressure term under variable viscosity (the constant-mu Timmermans term -mudiv(u) is only valid for homogeneous viscosity — Deteix & Yakoubi 2018). 'min' (default): constant coefficient chimu_min — provably stable at any contrast, exact fallback to the constant-mu scheme for uniform mu. 'full': pointwise chimu(i) — better pressure consistency at MILD contrast only. 'off': plain incremental (no rotational term). All modes keep the incremental predictor (large-dt / steady-Stokes capability).
`set_velocity_mg_smoother`	set_velocity_mg_smoother(self, name: str, degree: int = 0, eig_ratio: float = 10.0) -> None Which smoother the velocity multigrid uses on every level -- orthogonal to set_velocity_solver, and read only when that selects 'multigrid'. 'gauss_seidel' (DEFAULT) or 'chebyshev'. ABLATION, AND CHEBYSHEV IS MEASURED SLOWER: at degree 4 / eig_ratio 6 it IS the stronger smoother per cycle (8.0 V-cycles against red-black's 9.3) and still loses 1.21x on wall time, because matching red-black's smoothing strength costs the same four halo exchanges two sweeps cost, while Chebyshev needs a separate residual and update pass where a Gauss-Seidel sweep fuses the two. Kept because it is the apparatus for a question that gets asked repeatedly. degree=0 follows the V-cycle's own pre/post/bottom counts.
`set_velocity_multigrid_auto`	set_velocity_multigrid_auto(self, cells_per_rank: int, min_global_cells: int = -1, min_block_extent: int = -1) -> None Narrow the rule that selects the momentum solver. Since 2026-09-15 the DEFAULT is the velocity V-cycle on every validated operator mode and at every rank count: measured on the 1.0.0 scaling benchmark it is faster at every rung of both the CPU and the GPU ladder (2.23x at 147k cells/rank, 2.19x on one H100), and the red-black smoother it replaced was running at its sweep cap -- not converging -- on every rung of both. cells_per_rank sets an upper bound in cells per rank (0 = never use the V-cycle); set_velocity_multigrid_auto(65536, 1 << 23) restores the 1.0.0 behaviour exactly. min_block_extent (default 16) is the shortest per-rank inner extent worth building a hierarchy on. An explicit set_velocity_multigrid() overrides all of it.
`set_velocity_solver`	set_velocity_solver(self, name: str) -> None Which solver runs the implicit momentum (diffusion) equation. 'auto' (DEFAULT) picks by the operator's condition number kappa = 1 + 4 dt mu (sum w)/rho -- the multigrid at kappa >= 13, the smoother below -- and names no geometry, so an immersed body and a domain-BC box at the same dt get the same solver. 'gauss_seidel' is the red-black smoother; 'multigrid' the 3-level velocity V-cycle; 'chebyshev' Chebyshev semi-iteration over a Gershgorin interval. ALL FOUR SOLVE THE SAME EQUATION TO THE SAME TOLERANCE, so this changes cost and not the answer -- except that 'gauss_seidel' runs out of sweeps above kappa ~ 50 and then returns a solution that never met it. MEASURED (384^3 cut-cell bed, Snellius): multigrid beats gauss_seidel 1.66x-2.23x across the CPU ladder and 2.19x on one H100; chebyshev beats gauss_seidel by the sqrt(kappa) factor theory predicts but loses to multigrid except at the smallest per-rank blocks, where the evidence is suggestive and not conclusive.
`set_velocity_solver_params`	set_velocity_solver_params(self, iters: int, rtol: float = 0.0, min_iters: int = 2) -> None Momentum smoother control: `iters` RB-GS sweeps per component, or with rtol > 0 a TOLERANCE STOP — end the loop once the sweep's max velocity increment has contracted to rtol of the first sweep's (iters becomes the cap, min_iters the floor). Easy regimes (small nu*dt/dx^2) exit after ~3-5 sweeps; stiff regimes run to the cap unchanged. rtol = 0 (default) is the legacy fixed count, byte-identical.
`set_vof_block_debris`	set_vof_block_debris(self, enabled: bool) -> None Per-step marker DEBRIS removal of the block container (doc/vof_overlap_design.md 5.3, 11). Default: ON under enable_vof_block_csf, off otherwise. A debris cell is an interfacial cell of a marker (1e-8 < C < 1-1e-8) with every cell of that marker within two cells below C = 1/2: colour with no resolvable interface near it (what two overlapping markers tear off each other). Its colour is zeroed and the volume is returned to the marker's attached interfacial cells weighted C(1-C), in a fixed summation order, so the marker volume is exact and the state decomposition-independent. Ledger in vof_block_stats(): debris_cells / debris_volume (this step), debris_returned, debris_lost, debris_unresolved (cumulative, migrate with the block). enabled=False keeps the census (what WOULD be removed) and moves nothing.
`set_vof_block_device_staging`	set_vof_block_device_staging(self, on: bool) -> None Pack/unpack the block gather and scatter in the block's own MEMORY SPACE (rung W1 item b, default True), with a host staging copy only per MPI MESSAGE — and none at all for the master's own cells, which are a device-to-device copy. False selects rung W0's host-staged path (a full mirror of the local patch per step), which is what the device-vs-host measurement compares against. Every step is a copy of a double, so the two paths are BITWISE identical; the ctest gates that rather than asserting it.
`set_vof_block_kappa_clip`	set_vof_block_kappa_clip(self, enabled: bool, kappa_max: float
`set_vof_block_overlap_density`	set_vof_block_overlap_density(self, enabled: bool) -> None OPT-IN MODEL, default False (the MAX union verbatim). With it the union colour that feeds rho(C) and mu(C) -- and the registered 'C' -- becomes C_eff = max_k C_k where S = sum_k C_k <= 1 + 1e-8 and max(0, 2 - S) where markers overlap: the overlap volume is treated as the sub-grid FILM liquid (the 'tent'), so the phantom interface of one marker inside another is an ordinary liquid-gas interface under the base capillary bound. Its cost: the lens holds liquid the real state does not (buoyancy and mixture weight off by 2 V_overlap during a contact), and it changes collision dynamics. See doc/vof_overlap_design.md 5.6.
`set_vof_block_pool`	set_vof_block_pool(self, on: bool) -> None Recycle the advectors a re-centring retires, keyed by the exact box extent (rung W1 item c, default True). A translating bubble keeps its box SIZE and only moves its origin, so the hit rate is ~100 % and the ten Views of the new box are not allocated. A recycled advector is handed back in the state a freshly initialised one is in (colour and the three face-velocity fields zeroed), so the pool is bitwise inert.
`set_vof_curvature_mixed_height_fit`	set_vof_curvature_mixed_height_fit(self, on: bool = True) -> None Enable cascade branch 3, the mixed height-position fit (a paraboloid through the interface positions of whichever columns closed). DEFAULT OFF and it should stay off: measured on an exact-fraction sphere at 16/32/64 it takes over the 19.5-59.6% of cells the height function cannot serve and DESTROYS the convergence of the max curvature error (order 0.00 vs 1.86 with the PLIC-volumetric fallback instead), because its data set is the columns that closed - a slope-selected, asymmetric subset whose lever-arm bias is scale invariant. Shipped as a re-measurable instrument, not a configuration.
`set_vof_curvature_weight_width`	set_vof_curvature_weight_width(self, width: float) -> None Wendland support width of the PLIC-volumetric fallback fit, in CELL units. Default 2.5, which is Han et al.'s recommended pairing with the 5^3 stencil. Their translating droplet recovers first-order convergence of the spurious currents at 3.5 and loses convergence entirely at 4.5 (over-smoothing), so this is a real trade-off between the locality of the estimate and its robustness to transport error — not a free knob.
`set_vof_curvature_worklist`	set_vof_curvature_worklist(self, on: bool) -> None WO-V9: run the V3 curvature cascade (the PLIC plane pass, the height functions and the PV fallback) over a COMPACTED list of the interfacial cells instead of over the whole inner region. Default ON. The cascade is the most divergent kernel in the pipeline -- on a resolved droplet the interfacial cells are well under 1 % of the block, so a dense parallel_for puts one or two active lanes in a warp and the other thirty run the guard and then idle for the whole height function. Compaction is a pure re-ordering (each cell reads the same neighbours and writes the same value), so the two paths are BIT-IDENTICAL; tests/kokkos/test_vof_timing.cpp gates that and WO-V9 records the gain.
`set_vof_cutcell_flux_clamp`	set_vof_cutcell_flux_clamp(self, on: bool) -> None Ablation: Weymouth's admissible-interval clamp on the openness-weighted cut-cell flux (ON by default). The clamp bounds
`set_vof_flux_clamp`	set_vof_flux_clamp(self, on: bool = True) -> None ABLATION: drop the Weymouth flux clamp max(0,
`set_vof_interface_eps`	set_vof_interface_eps(self, eps: float) -> None Wisp threshold on the curvature cascade's interfacial predicate once surface tension is on: a cell carries an interface only while eps < C < 1-eps. Default 1e-8. This is NOT optional and it is not cosmetic. Weymouth-Yue leaves round-off colour residue (measured down to -3e-35) in every cell its sweeps touch; those cells satisfy 0 < C < 1, so the cascade builds a PLIC polygon of area ~0 for them and returns
`set_vof_kappa_constant`	set_vof_kappa_constant(self, kappa: float) -> None INSTRUMENT: set kappa to a constant everywhere (inner + ghosts), mark every cell as carrying a valid estimate, and freeze it. The CSF force is then EXACTLY the discrete gradient of sigmakappaC, so the projection must annihilate it to round-off from any colour field whatsoever — the exactness gate of rung V4, independent of curvature accuracy and of resolution.
`set_vof_kappa_frozen`	set_vof_kappa_frozen(self, on: bool = True) -> None INSTRUMENT: stop recomputing the curvature at the head of each step and use whatever is in the 'kappa'/'kappa_branch' fields. With set_vof_kappa_constant this isolates the balanced-force identity from the curvature estimator.
`set_vof_momentum_cell_flag`	set_vof_momentum_cell_flag(self, on: bool = True) -> None ABLATION: use the PRESSURE-cell frozen dilation flag H(C^n-1/2) on the shifted control volume instead of its structural analogue H(C^c,n-1/2). The flag that must be shared is the one of the pair that telescopes, and both members of that pair live on the shifted volume — this switch is the literal reading of the work order, kept as a measurement.
`set_vof_momentum_muscl`	set_vof_momentum_muscl(self, on: bool = True) -> None MinMod-limited linear reconstruction of the donor velocity in the momentum flux, instead of the DEFAULT plain donor-cell upwind. Both preserve the uniform-velocity identity exactly (a uniform field has a zero slope bit for bit), but on a control volume a sweep EMPTIES the slope's deviation from the volume's own velocity is amplified by drho*F/rho^c, which is unbounded in the density ratio. Measured at ratio 1e4, 50 steps: with the slope the uniform-velocity residual grows to 2.2e-10, without it it is flat at 6.7e-16; at ratio 1e3 the slope is harmless. Turn it on deliberately and re-run the ratio sweep if you do.
`set_vof_phantom_capillary_bound`	set_vof_phantom_capillary_bound(self, enabled: bool = True) -> None The phantom-aware capillary limit of the block container (doc/vof_overlap_design.md 5.6, 11.3), DEFAULT True: while two block markers overlap (max S > 1 + 1e-2) the capillary limit of vof_step_limits() and step() uses 2 rho_min in place of rho_l + rho_g -- the part of one marker's surface inside another is an interface with gas on both sides. False restores the pre-5.6 bound (the ablation; a fixed-dt driver at a density ratio above ~31 otherwise throws at first contact).
`set_vof_rho_floor`	set_vof_rho_floor(self, fraction: float) -> None Floor on rho^c in the recovery divide u = (rho^c u)/rho^c, as a FRACTION of min(rho_gas, rho_liquid). Default 1e-6. rho^c leaves [rho_gas, rho_liquid] only through a wisp in the half-shifted colour and reaching zero would need C^c ~ -1/(ratio-1), so this is a guard, not a model — vof_momentum_diagnostics()['floored'] reports how many control volumes it actually touched.
`set_vof_step_parity`	set_vof_step_parity(self, n: int) -> None Set the sweep-permutation counter of the NEXT colour advection: the Weymouth-Yue sweep order is kWySweepPerm[n % 6], cycled so no axis is systematically favoured. Exposed so a benchmark can hold the permutation fixed (n constant) or resume a run across a restart. Default: it increments once per advection from 0.
`set_vof_timing`	set_vof_timing(self, on: bool) -> None WO-V9: arm (or disarm) the VoF pipeline's per-stage timers, and reset them. OFF by default. When armed, every stage boundary calls Kokkos::fence() before reading the clock -- the same rule the step's three coarse phase timers already follow, because on a device backend queued work would otherwise be billed to whichever stage next reads the clock. A fence changes WHEN work happens, never WHAT it computes: a run with the timers armed is bit-identical to the same run without them (gated in tests/kokkos/test_vof_timing.cpp). When disarmed the cost is one branch per stage and no fence at all.
`set_vof_wisp_eps`	set_vof_wisp_eps(self, eps: float) -> None Wisp tolerance on the Weymouth-Yue mixed-cell predicate: a cell counts as carrying an interface only while eps < C < 1 - eps, and one outside that band is fluxed ALGEBRAICALLY as C_donor * a — its ACTUAL colour, so the exact telescoping conservation is untouched. The same threshold gates the interface-local Courant band (
`set_vof_worklist`	set_vof_worklist(self, on: bool) -> None WyAdvector::useWorklist -- compact the PLIC reconstruction pass onto the mixed cells (a parallel_scan over the reconstruction region, then a dense parallel_for over the compacted list) instead of a guarded parallel_for over the whole region. Default ON. Pure optimization: the flux never reads a non-mixed cell's plane, so the two paths produce the same field BIT FOR BIT (gated in tests/kokkos/test_vof_advect.cpp and measured in WO-V9).
`set_wall_flux_divergence`	set_wall_flux_divergence(self, on: bool) -> None Rung 3 on/off (default on = correct physics). Off leaves the moving wall's no-slip datum in the momentum operator but drops the wall's own volume flux from the cut-cell projection -- the configuration the Galilean gate uses to exhibit the failure the term fixes.
`uf_advection`	Whether the projected face field is the collocated advecting velocity.
`velocity_mg_smoother`	velocity_mg_smoother(self) -> str The velocity multigrid's smoother.
`velocity_multigrid_active`	velocity_multigrid_active(self) -> bool Whether the velocity (momentum) solve uses the multigrid V-cycle rather than plain RB-GS. Reads back the OUTCOME, not the setting: an explicit set_velocity_multigrid(on=...) call always wins, but absent one the AUTO rule (set_velocity_multigrid_auto) decides it when the geometry is built (set_solid / set_pressure_geometry / set_solid_from_scene) -- on for a distributed run (np > 1) in an eligible operator mode with a global problem of at least min_global_cells (default 8M) and fewer than cells_per_rank (default 65536) cells per rank, off otherwise. So the value can differ from what the caller last set explicitly, and is only meaningful after geometry has been installed (False on a fresh solver).
`velocity_solver`	velocity_solver(self) -> str The momentum solver actually in force, resolving 'auto'. Only meaningful after the first step(), which is where the choice is made (dt and mu are final only then).
`vof_advected_velocity`	vof_advected_velocity(self, component: int) -> numpy.ndarray[dtype=float64] The recovered advected velocity (rho^c u_c)/rho^c of component c on the inner cells — the momentum RHS's time base. Exposed so the uniform-velocity consistency identity can be gated on the ADVECTION ALONE, with the projection and the momentum solve out of the picture.
`vof_block_census`	vof_block_census(self) -> dict Per-rank load census of the block container: 'masters'[r] = blocks rank r masters, 'cells'[r] = the inner cells those blocks carry (the actual VoF work).
`vof_block_curvature_stats`	vof_block_curvature_stats(self) -> dict Branch census of the last block-CSF curvature pass, SUMMED over this rank's blocks (local to the rank). 'no_estimate' must be 0 on any gated case. 'clipped' = cells whose
`vof_block_force`	vof_block_force(self, component: int) -> numpy.ndarray[dtype=float64] The block CSF face force of component c (0,1,2) as scattered (summed over markers) into the solver's LOW-face slot of each cell, internal units; the field the momentum RHS adds. A diagnostic.
`vof_block_imbalance`	vof_block_imbalance(self) -> float max/mean of the per-rank block-cell load under the CURRENT master assignment (round robin by block id at rung W0). 1.0 is perfect; the weighted-ORB assignment of rung W1 is what this number is here to grade.
`vof_block_imbalance_of`	vof_block_imbalance_of(self, mode: int) -> float The max/mean per-rank block-cell load the given assignment mode WOULD give on the current blocks, without applying it — so a study can put the three modes side by side on one swarm without perturbing the run.
`vof_block_kappa`	vof_block_kappa(self, id: int) -> numpy.ndarray[dtype=float64] One marker's OWN curvature (internal units, 1/cell) over its block box, the layout of vof_block_color; empty unless this rank masters it and the block CSF is on. A diagnostic.
`vof_block_pool_stats`	vof_block_pool_stats(self) -> dict Block-pool census: 'hits' = advectors recycled, 'misses' = advectors allocated.
`vof_block_stats`	vof_block_stats(self) -> list Per-bubble Lagrangian census, one dict per block, in block-id order. 'lo'/'hi' (the global index box) and 'master' are replicated on every rank; the MEASURED entries — 'volume' (sum of C over the box, cell volumes), 'centroid', 'velocity' (d(centroid)/dt of the last advection), 'moments' (the central second moments xx, yy, zz, xy, xz, yz divided by the volume, i.e. the deformation) — are filled only on the block's MASTER and are zero elsewhere. 'discarded' is the colour a re-centring dropped, cumulatively: Weymouth-Yue leaves round-off residue in every cell its sweeps touch and the box tracks the BUBBLE, not that wake, so the residue falling outside the new box is discarded. Never physical liquid (measured -9.5e-17 over a 20-cell translation, against a bubble volume of 524), but it is reported rather than hidden — a container that silently loses mass is not acceptable. 'debris_cells' / 'debris_volume' are this step's marker DEBRIS removal (block CSF on only; master only; see set_vof_block_debris): interfacial cells with every cell of the marker within two cells below C = 1/2, and their summed colour. 'debris_returned', 'debris_lost' and 'debris_unresolved' are the cumulative ledger (doc/vof_overlap_design.md 5.3): volume + discarded + debris_lost is conserved. 'full_axis' counts re-centrings whose box had to span a whole periodic axis (it is then [0, L) and copied with the wrap, losslessly; 12.3). 'residue_returned' is the cumulative SIGNED sub-wisp-threshold residue cleared after each advection and returned to the interface with the debris (13): volume + discarded + debris_lost is exact. 'overlap_*' are the GLOBAL marker-overlap census of the last block advection (vof_overlap_design 5.6; block CSF or the overlap density on), the same in every dict: 'overlap_max_sum' = max over cells of S = sum_k C_k, 'overlap_excess' = sum of (S - 1)^+ (cell volumes), 'overlap_cells' = cells with S > 1 + 1e-8, and 'overlap_bound_active' = the phantom-aware capillary limit (2 rho_min in place of rho_l + rho_g) is in force for the next step (triggered by max S > 1 + 1e-2).
`vof_curvature_worklist`	vof_curvature_worklist(self) -> bool Whether the curvature-cascade compaction is on.
`vof_diagnostics`	vof_diagnostics(self) -> dict Colour census over THIS RANK's inner cells: sum (cell-volume units), min, max, the number of mixed cells (0<C<1) and of wisps (C within 1e-8 of 0 or 1). No clipping is applied at this rung, so min/max may leave [0,1] if the CFL cap is raised. With an immersed solid (rung V5a) the dict also carries the cut-cell quantities, all zero without one: 'volume' = sum eps_effC over fluid cells, the functional the openness-weighted scheme conserves EXACTLY (eps_eff = max(eps, 1/64), the 4^3 subsampling resolution — see src/vof/cutcell.hpp rule 1); 'raw_volume' = sum epsC, which differs only on eps==0 cells that still own an open face and is therefore NOT the conserved functional; 'solid_sum' = sum of C over solid cells (0 by construction — the neutral band fill lives on the working block, not on 'C'); 'min_fluid'/'max_fluid' over UNCUT fluid cells (eps==1), where Weymouth's boundedness applies verbatim; 'clipped_volume' / 'clipped_signed' = the liquid volume the cut-cell clip moved during the last advection (a TRIPWIRE on the flux approximation, not a mechanism: if it is not negligible the fix is the solid-clipped flux polygon); 'cut_cells' / 'solid_cells'. 'inflow_volume' / 'outflow_volume' are the liquid volume that ENTERED / LEFT through the domain faces during the LAST colour advection, in the same cell-volume units as 'sum' (0 unless a VoF boundary colour is set — rung V-BC). They are the advector's OWN face fluxes, so sum(C) - sum(C_0) = integral(inflow - outflow) holds to round-off whatever the interface does; vof_bc_volumes() breaks them out per face.
`vof_dynamic_field`	vof_dynamic_field(self, which: int) -> numpy.ndarray[dtype=float64] The V6 per-cell dynamic-wetting state on the inner region, as (nx,ny,nz): 0 the IMPOSED angle in degrees, 1 the measured APPARENT angle in degrees, 2 the smoothed U_cl (positive = the liquid ADVANCES), 3 Ca_cl, 4 the state (0 not a contact cell, 1 Cox-Voinov on the static base, 2 PINNED, 3 advancing, 4 receding). Regenerates the fill, so it is also the direct gate on the pass's decomposition independence.
`vof_filled_color`	vof_filled_color(self) -> numpy.ndarray[dtype=float64] The colour field INCLUDING the neutral solid-band fill — what the MYC and height-function stencils actually read — as a Fortran-order (nx,ny,nz) array. get_vof()/'C' is the canonical field and carries EXACTLY 0 in solid cells; the fill is a stencil device regenerated at every ghost exchange (three passes with a shrinking depth budget, src/vof/cutcell.hpp), and it is what makes the wall look 90-degree neutral instead of perfectly non-wetting.
`vof_geometry`	vof_geometry(self, which: int) -> numpy.ndarray[dtype=float64] The cut-cell geometry the colour block runs on, as a Fortran-order (nx,ny,nz) array: 0 = the cell fluid fraction eps (4^3-subsampled, a multiple of 1/64), 1/2/3 = the openness of the +x/+y/+z face of each cell (the ADVECTOR's high-face convention, one cell shifted from the solver's ox/oy/oz), 4 = the classification (1 = SOLID, i.e. eps == 0 AND all six faces closed). On an ALL-FLUID solver (set_pressure_geometry, or no cut-cell pressure operator) there is no geometry block and this returns the TRIVIAL geometry the V1 transport kernels actually execute: eps = 1, openness = 1, classification = 0. It used to raise, which meant one diagnostic could not serve a packed scene and its all-fluid control.
`vof_interface_eps`	vof_interface_eps(self) -> float The wisp threshold set by set_vof_interface_eps.
`vof_last_courant`	vof_last_courant(self) -> float The interface-local Courant number of the step just taken (0 before the first step).
`vof_momentum_diagnostics`	vof_momentum_diagnostics(self) -> dict Census over THIS RANK's momentum control volumes: per-component min/max of the half-shifted colour C^c, the summed momentum rho^c u_c (the conservation census), the minimum rho^c before the floor, and the number of control volumes the floor touched.
`vof_rho_floor`	vof_rho_floor(self) -> float The absolute rho^c floor used by the last recovery.
`vof_step_parity`	vof_step_parity(self) -> int The sweep-permutation counter of the next colour advection.
`vof_timing`	vof_timing(self) -> dict WO-V9: cumulative VoF stage times in SECONDS since the last set_vof_timing/reset_vof_timing, on this rank. `steps` is the number of step() calls in the window, so every entry divides down to a per-step cost. `step`/`predictor`/`momentum_solve`/`projection` are the step's three coarse phases summed over the SAME window (the remainder of `step` is BC re-imposition, Picard bookkeeping and the scalars). `vof_advect` is the whole colour stage and `vof_bridge` the part of it spent in the G=2 <-> g=3 bridges and the C ghost policy; `vof_momentum_advect` is its V2b twin. `k_*` are the advector's own kernels, shared by the colour, momentum-consistent and consistent-energy drivers, so they are the per-KERNEL breakdown of whichever of the two stages is running: `k_reconstruct` the MYC + plicAlpha pass, `k_fluxes` the geometric face fluxes, `k_sweep` the update, `k_clip` the V5a cut-cell clip, `k_exchange` the g = 3 ghost exchange, `k_sweeps` the number of sweeps (3 per advect call). Returns zeros unless the timers are armed.
`vof_wisp_eps`	vof_wisp_eps(self) -> float The wisp tolerance currently in force (see set_vof_wisp_eps).
`vof_worklist`	vof_worklist(self) -> bool Whether the reconstruction-pass compaction is on.
`wall_area_probe`	wall_area_probe(self) -> list[float] Diagnostic: [sum_c x_cA_wall_x, sum_c y_cA_wall_y, sum_c z_c*A_wall_z] over all cut cells. Must equal -V_solid componentwise if the aperture wall-area vectors are right, so it separates a force-integral error in the GEOMETRY from one in the traction.
`wall_flux_imbalance`	wall_flux_imbalance(self) -> float Sum over cells of u_wall . A_wall -- the compatibility datum of the singular pressure problem. Exactly zero for a translating body in a periodic box; small but nonzero for rotation. Reported, not corrected.
`wall_slip_sandwich_cells`	wall_slip_sandwich_cells(self) -> tuple Per velocity component, the number of cut-cell AXES at which a one-cell fluid gap made the Navier closure inapplicable and the no-slip one was kept. Nonzero means part of the wall is silently no-slip; report it.

Module attributes¶

Attribute	Value in this build
`execution_space`	`'OpenMP'`
`has_mpi`	`True`

`finalize`¶

finalize() -> None

Release every live object and zero-copy array of this module, then Kokkos::finalize() (deterministic teardown; also run automatically at interpreter exit). Idempotent. After it, the module's solver objects and *_view arrays must not be used.

`march_to_steady`¶

March ``solver`` to its steady state; return a :class:`MarchResult`.

    Parameters
    ----------
    solver : peclet.flow.Solver or peclet.flow.SolverColocated
        A configured solver (geometry, properties, dt, forcing). Its current state is the start.
    monitor : callable
        ``monitor() -> float``: the scalar the stop instrument watches, e.g. the mean velocity along
        the forcing. Under MPI it must return the same value on every rank (a global reduction).
    rtol : float
        Relative bound on the change of ``monitor()`` still to come.
    max_steps : int
        Cap on the total number of solver steps.
    accelerate : bool
        Anderson acceleration (``True``) or the plain march (``False``). Raises on a configuration
        the accelerator refuses, with the reason and the hint to pass ``accelerate=False``.
    window : int
        Anderson window, 1 to 8.
    check_every, num_passes, slow_rate, roundoff
        The stop instrument: block length in steps, consecutive passing blocks required, the
        slowest per-step contraction assumed, and the relative block change accepted as round-off.
    callback : callable, optional
        ``callback(steps, phase)`` after every step, ``phase`` in ``{"accelerate", "certify",
        "plain"}`` (logging, checkpoints).

`mpi_block`¶

mpi_block(gnx: int, gny: int, gnz: int, levels: int = 0, max_imbalance: float = 1.05) -> tuple[list[int], list[int]]

Return this MPI rank's ORB block of the global (gnx,gny,gnz) grid as (origin, size), each a length-3 list [x,y,z]. Use it to slice the global SDF into this rank's local block for a distributed Solver, in this order: Solver(*size), init_mpi(gnx, gny, gnz), then the geometry (init_mpi raises after set_solid / set_pressure_geometry). MPI_Init is called if needed.

levels / max_imbalance choose how the partition is built and MUST match the values given to Solver.set_decomposition() -- both derive the same partition and have to agree. levels=0 (DEFAULT) = the aligned ORB: split positions on the FINE grid are snapped to a power of two (capped at 16, i.e. 5 nested multigrid levels). levels >= 2 = COARSE-FIRST: build the ORB on the grid coarsened levels-1 times and refine the partition upward, so every block is a multiple of the coarsening factor BY CONSTRUCTION and the hierarchy nests for the full requested depth. Coarse-first also balances better: the aligned ORB picks a split and then rounds it (a balanced 96|96 can round to 128|64), whereas on the coarse grid one cell IS the quantum, so the rank count no longer has to be a power of two -- what matters is that the COARSE grid divides among the ranks. Each candidate depth is BUILT and measured, and the deepest one whose max-block/min-block ratio stays within max_imbalance (default 1.05) is taken, else the aligned ORB.

`predict_hierarchy`¶

predict_hierarchy(gnx: int, gny: int, gnz: int, np: int, levels: int, telescope: bool = False, min_extent: int = 4, decomposition_levels: int = 0, max_imbalance: float = 1.05, weights: collections.abc.Sequence[float] | None = None) -> object

The pressure-multigrid hierarchy init_mpi WOULD build for this grid / rank count / level request, for the level-0 decomposition that decomposition_levels / max_imbalance describe (pass what you will pass to Solver.set_decomposition) and with or without coarse-level telescoping -- a pure function, no MPI, no allocation, any rank count. Returns one row per level: (global dims, ranks holding the level, block-0 dims, ratio to the next level, telescopes-out). Pre-flight any job with it.

weights: per-cell weights, exactly what diagnostics.rebalance_by_weights takes (global x-fastest, gnx*gny*gnz values). Given, the forecast is the hierarchy AFTER that rebalance: level 0 is its aligned weighted ORB (split planes on multiples of 2^a, a the largest alignment whose weight imbalance stays within 1.05; decomposition_levels / max_imbalance are then ignored, as the rebalance ignores set_decomposition), and a stage may be a repartition onto fewer ranks (reported as telescopes-out). The return is then the pair (rows, align), align = 2^a -- the value rebalance_by_weights returns and a co-decomposing dem takes as migrate_to_weights(weights, align=...). None (default): the rows alone, unchanged.

peclet.flow — Eulerian Navier–Stokes solver¶

peclet.flow¶

Solver¶

SolverColocated¶

SolverDiagnostics¶

SolverColocatedDiagnostics¶