Functions
	_grid (N)

	_minimg (d, N)

	sdf_box (N)

	sdf_cyl (N)

	_hollow_cyl_sdf (X, Y, Z, c, axis, r_out, r_in, H, N)

	sdf_rings (N)

	sdf_pack (N)

	level_set (shape, N)

	rho_field (shape, edge, ratio, N)

	run_one (geom, shape, edge, ratio, driver, case, N, quiet=True, bottom=None, levels=None)

	markdown (records)

	summarize (records)

	cheb_overhead (geom, N, ratio=1e3)

	trace_pcg (geom, shape, edge, ratio, case, N)

	main ()

Variables
int	N_DEFAULT = 32

int	N_PACK = 64

int	STEPS = 20

int	MAXIT = 200

int	RTOL = 1e-8

float	MU = 0.01

float	DT = 1.0

float	GRAV = 0.1

float	ULID = 1.0

int	FBODY = 1e-3

int	VEL_SWEEPS = 40

tuple	RATIOS = (1e2, 1e3, 1e4)

tuple	SHAPES = ("slab", "blob", "tilt")

tuple	GEOM_ORDER = ("box", "cyl", "rings", "pack")

tuple	SHAPE_ORDER = ("const", "slab", "blob", "tilt")

tuple	EDGES = ("smooth", "sharp")

tuple	DRIVERS = ("cheb", "pcg", "fcg")

tuple	DRIVERS_DEFAULT = ("cheb", "pcg")

tuple	CASES = ("hydro", "lid", "per")

dict	GEOMS = {"box": sdf_box, "cyl": sdf_cyl, "rings": sdf_rings, "pack": sdf_pack}

Detailed Description

WO-B (VoF rung S0) -- pressure-driver measurement battery on STATIC two-phase coefficients.

Purpose (see `doc/vof_workorders.md` WO-B and `suite/docs/VOF_PLAN.md` section 5, the S-ladder):
decide, from measurement rather than from the single anecdote in
`doc/variable_density_projection.md` section 2, how far up the S-ladder the variable-density
pressure solver work has to climb.  NO VoF is needed: the density field is manufactured and
static, so this probes the *linear solver* on the exact coefficient operator a VoF interface
would present.

What is swept
-------------
  geometry  x  rho-shape  x  edge  x  density ratio  x  driver  x  forcing case

  geometry   `box`   open box, no immersed solid (`set_pressure_geometry`), walls +-z
             `cyl`   immersed cylinder along y, R = 0.25 N (cut-cell IBM)
             `rings` 3 Raschig rings (the regression suite's `hollow_rings` bed) -- a genuinely
                     tortuous cut-cell pore geometry
             `pack`  `data/packing_ring.vti` strided onto the run grid (opt-in, N = 64)
  rho-shape  `slab`  flat film normal to z, heavy below  (the only true hydrostatic rest state)
             `blob`  sphere of heavy fluid at the domain centre
             `tilt`  film with normal (1,0,1)/sqrt2 -- deliberately NOT aligned with the MG
                     coarsening axes
             `const` THE CONTROL: constant density, `set_density_mode` never called, gravity as
                     a uniform body force.  Same geometry, same BCs, same driver, unit
                     coefficient -- so `pack/const` vs `pack/slab@1e4` is exactly the
                     "geometry alone" vs "geometry x coefficient" separation.  Recorded with
                     ratio 1.
  edge       `smooth` tanh over ~2 cells       `sharp` a one-cell step (what a real VoF
                     interface presents to the operator)
  ratio      1e2, 1e3, 1e4
  driver     `cheb` (`set_pressure_chebyshev`)   `pcg` (`set_pressure_pcg`)
             `fcg` (`set_pressure_fcg`, added by WO-C / rung S1) -- the SAME V-cycle-preconditioned
                     CG as `pcg` with the same cap/tolerance/stopping estimate, differing only in
                     the beta recurrence (Polak-Ribiere instead of Fletcher-Reeves).  `fcg` is NOT
                     in the default `--drivers` set, so the WO-B battery reproduces unchanged; ask
                     for it with `--drivers pcg,fcg` (or `cheb,pcg,fcg`).

**Driver selection (historical trap, REPAIRED 2026-08-30 by WO-H).**  `set_pressure_pcg(on, ...)`
now genuinely selects MG-PCG (it clears both `useChebyshev_` and `useFcg_`), so the plain
`set_pressure_pcg(True, cap, rtol)` after `set_density_mode` is correct.  Until WO-H its `on` flag
was **discarded** -- `IbmSolver::setPressurePcg(bool /*on*/, ...)` only stored
`pcgMaxit_`/`pcgRtol_` -- so a "PCG" run after `set_density_mode("variable")` silently measured
**Chebyshev** (the tell: it capped at exactly **120**, `chebMaxit_`'s default, not the cap passed),
and the only working spelling was `set_pressure_chebyshev(False, ...)`.  `run_one` still issues that
`set_pressure_chebyshev(False, ...)` before `set_pressure_pcg` -- redundant now, kept so this script
reproduces WO-B's and WO-C's recorded batteries byte-for-byte on an older build.
  case       `hydro`  quiescent + gravity closure force_z = -g*rho (the hydrostatic setup);
                      periodic x/y, walls +-z
             `lid`    lid-driven flow: walls -x/+x/-z, the +z face translating at U in x,
                      periodic y.
             `per`    THE PRODUCTION CONFIGURATION: fully periodic + a uniform body force in x
                      -- pore-scale flow through the packing, no domain BC anywhere.  This case
                      is where the two drivers swap places, so it is not optional.  It is
                      skipped for (box, const) and (box, slab), where a z-only stratification in
                      a periodic box leaves u* already x-uniform, div(u*) = 0 identically, and
                      the pressure solve returns after 0 iterations -- a degenerate cell, not a
                      measurement.

What is recorded, per configuration
-----------------------------------
  * pressure iterations per step (`last_pressure_iterations()`) -- median / max / total;
  * "final residual" = `max_open_divergence()` = max|div(open*u)| after the step.  This IS the
    residual of the pressure system in the divergence normalisation: the projection solves
    div(open*(rho0/rho_f) grad phi) = div(open*u*) and then subtracts the gradient, so whatever
    divergence survives is exactly what the linear solve failed to remove;
  * wall time: whole step and the `projection` phase alone (device-fenced per-phase timers);
  * converged: whether the driver stayed strictly below its iteration cap on every step.

Both drivers get the SAME rtol (1e-8) and the SAME cap (200), so "hit the cap" means the same
thing for both and a stalling configuration terminates instead of burning the device.

Run:
    OMP_NUM_THREADS=8 OMP_PROC_BIND=false PYTHONPATH=<build> \
        python tests/study/vardensity_solver_probe.py --out tests/study/vardensity_solver_probe.json
    ... --pack                 add the `data/packing_ring.vti` sub-sweep (N=64)
    ... --trace                print the PCG residual history of one stalling configuration
    ... --geoms box,cyl        restrict the sweep

Function Documentation

◆ _grid()

vardensity_solver_probe._grid ( N )

protected

Definition at line 113 of file vardensity_solver_probe.py.

Referenced by level_set(), sdf_cyl(), and sdf_rings().

◆ _minimg()

vardensity_solver_probe._minimg	(	d,
		N
	)

protected

Definition at line 118 of file vardensity_solver_probe.py.

Referenced by _hollow_cyl_sdf().

◆ sdf_box()

vardensity_solver_probe.sdf_box ( N )

No immersed solid: an all-fluid field for `set_pressure_geometry`.

Definition at line 122 of file vardensity_solver_probe.py.

◆ sdf_cyl()

vardensity_solver_probe.sdf_cyl ( N )

Immersed cylinder along y, centred, R = 0.25 N.

Definition at line 127 of file vardensity_solver_probe.py.

References _grid().

◆ _hollow_cyl_sdf()

vardensity_solver_probe._hollow_cyl_sdf	(	X,
		Y,
		Z,
		c,
		axis,
		r_out,
		r_in,
		H,
		N
	)

protected

One Raschig ring (hollow cylinder): annulus [r_in,r_out] intersect the slab |axial| <= H/2.

Definition at line 134 of file vardensity_solver_probe.py.

References _minimg().

Referenced by sdf_rings().

◆ sdf_rings()

vardensity_solver_probe.sdf_rings ( N )

The regression suite's `hollow_rings` bed: 3 Raschig rings, fixed positions/orientations.

Definition at line 147 of file vardensity_solver_probe.py.

References _grid(), and _hollow_cyl_sdf().

◆ sdf_pack()

vardensity_solver_probe.sdf_pack ( N )

`data/packing_ring.vti` (256^3, physical units) strided onto an N^3 grid, in CELL units.

Definition at line 160 of file vardensity_solver_probe.py.

◆ level_set()

vardensity_solver_probe.level_set	(	shape,
		N
	)

psi > 0 in the HEAVY fluid, in units of cells (so `psi/w` is a cell-width argument).

Definition at line 178 of file vardensity_solver_probe.py.

References _grid().

Referenced by rho_field().

◆ rho_field()

vardensity_solver_probe.rho_field	(	shape,
		edge,
		ratio,
		N
	)

Definition at line 192 of file vardensity_solver_probe.py.

References level_set().

Referenced by cheb_overhead(), and run_one().

◆ run_one()

vardensity_solver_probe.run_one	(	geom,
		shape,
		edge,
		ratio,
		driver,
		case,
		N,
		quiet = `True`,
		bottom = `None`,
		levels = `None`
	)

Definition at line 204 of file vardensity_solver_probe.py.

References rho_field().

Referenced by main(), and trace_pcg().

◆ markdown()

vardensity_solver_probe.markdown ( records )

Definition at line 309 of file vardensity_solver_probe.py.

Referenced by main().

◆ summarize()

vardensity_solver_probe.summarize ( records )

The two questions the S-ladder decision turns on.

Definition at line 347 of file vardensity_solver_probe.py.

Referenced by main().

◆ cheb_overhead()

vardensity_solver_probe.cheb_overhead	(	geom,
		N,
		ratio = `1e3`
	)

Quantify what rung S2 (bound amortization) would actually buy.

Under varRho the coefficient operator is rebuilt every step and `chebBoundsSet_` is
invalidated with it (`flow_ibm.hpp` project()), so `estimateEigenvalues(..., iters=15, ...)`
runs EVERY step -- two 15-step power loops, each applying M^{-1}A once, i.e. **30 extra
V-cycles per step** on top of the ~15-25 the solve itself needs.  Constant density estimates
the bounds once (step 1) and reuses them, so the A/B below isolates the cost.

Definition at line 404 of file vardensity_solver_probe.py.

References rho_field(), and run().

Referenced by main().

◆ trace_pcg()

vardensity_solver_probe.trace_pcg	(	geom,
		shape,
		edge,
		ratio,
		case,
		N
	)

Re-run one configuration with the MG residual trace on (PCG only -- Chebyshev has none).

Definition at line 474 of file vardensity_solver_probe.py.

References run_one().

Referenced by main().

◆ main()

vardensity_solver_probe.main ( )

Definition at line 485 of file vardensity_solver_probe.py.

References cheb_overhead(), main(), markdown(), run_one(), summarize(), and trace_pcg().

Referenced by main().

Variable Documentation

◆ N_DEFAULT

int vardensity_solver_probe.N_DEFAULT = 32

Definition at line 90 of file vardensity_solver_probe.py.

◆ N_PACK

int vardensity_solver_probe.N_PACK = 64

Definition at line 91 of file vardensity_solver_probe.py.

◆ STEPS

int vardensity_solver_probe.STEPS = 20

Definition at line 92 of file vardensity_solver_probe.py.

◆ MAXIT

int vardensity_solver_probe.MAXIT = 200

Definition at line 93 of file vardensity_solver_probe.py.

◆ RTOL

int vardensity_solver_probe.RTOL = 1e-8

Definition at line 94 of file vardensity_solver_probe.py.

◆ MU

float vardensity_solver_probe.MU = 0.01

Definition at line 95 of file vardensity_solver_probe.py.

◆ DT

float vardensity_solver_probe.DT = 1.0

Definition at line 96 of file vardensity_solver_probe.py.

◆ GRAV

float vardensity_solver_probe.GRAV = 0.1

Definition at line 97 of file vardensity_solver_probe.py.

◆ ULID

float vardensity_solver_probe.ULID = 1.0

Definition at line 98 of file vardensity_solver_probe.py.

◆ FBODY

int vardensity_solver_probe.FBODY = 1e-3

Definition at line 99 of file vardensity_solver_probe.py.

◆ VEL_SWEEPS

int vardensity_solver_probe.VEL_SWEEPS = 40

Definition at line 100 of file vardensity_solver_probe.py.

◆ RATIOS

tuple vardensity_solver_probe.RATIOS = (1e2, 1e3, 1e4)

Definition at line 102 of file vardensity_solver_probe.py.

◆ SHAPES

tuple vardensity_solver_probe.SHAPES = ("slab", "blob", "tilt")

Definition at line 103 of file vardensity_solver_probe.py.

◆ GEOM_ORDER

tuple vardensity_solver_probe.GEOM_ORDER = ("box", "cyl", "rings", "pack")

Definition at line 104 of file vardensity_solver_probe.py.

◆ SHAPE_ORDER

tuple vardensity_solver_probe.SHAPE_ORDER = ("const", "slab", "blob", "tilt")

Definition at line 105 of file vardensity_solver_probe.py.

◆ EDGES

tuple vardensity_solver_probe.EDGES = ("smooth", "sharp")

Definition at line 106 of file vardensity_solver_probe.py.

◆ DRIVERS

tuple vardensity_solver_probe.DRIVERS = ("cheb", "pcg", "fcg")

Definition at line 107 of file vardensity_solver_probe.py.

◆ DRIVERS_DEFAULT

tuple vardensity_solver_probe.DRIVERS_DEFAULT = ("cheb", "pcg")

Definition at line 108 of file vardensity_solver_probe.py.

◆ CASES

tuple vardensity_solver_probe.CASES = ("hydro", "lid", "per")

Definition at line 109 of file vardensity_solver_probe.py.

◆ GEOMS

dict vardensity_solver_probe.GEOMS = {"box": sdf_box, "cyl": sdf_cyl, "rings": sdf_rings, "pack": sdf_pack}

Definition at line 174 of file vardensity_solver_probe.py.

Functions

Variables

Detailed Description

Function Documentation

◆ _grid()

◆ _minimg()

◆ sdf_box()

◆ sdf_cyl()

◆ _hollow_cyl_sdf()

◆ sdf_rings()

◆ sdf_pack()

◆ level_set()

◆ rho_field()

◆ run_one()

◆ markdown()

◆ summarize()

◆ cheb_overhead()

◆ trace_pcg()

◆ main()

Variable Documentation

◆ N_DEFAULT

◆ N_PACK

◆ STEPS

◆ MAXIT

◆ RTOL

◆ MU

◆ DT

◆ GRAV

◆ ULID

◆ FBODY

◆ VEL_SWEEPS

◆ RATIOS

◆ SHAPES

◆ GEOM_ORDER

◆ SHAPE_ORDER

◆ EDGES

◆ DRIVERS

◆ DRIVERS_DEFAULT

◆ CASES

◆ GEOMS