Functions
	gpu_mem_mb ()

	rcp_sdf (Ng)

	sc_sphere_sdf (N, phi)

	case_hydro (ratios=(1e1, 1e2, 1e3, 1e4, 1e5, 1e6), N=8, NZ=24, g=0.1, steps=100)

	case_vof (ratios=(1e1, 1e2, 1e3, 1e4), n=32, steps=200)

	case_porous (N=16, beta=4.0, f_drive=0.2, eps=0.6, dt=0.5, steps=200)

	case_contrast (Ngs=(48, 64, 96), steps=4, maxit=300, rtol=1e-8, driver="pcg")

	case_zh (Ns=(32, 48, 64), phi=0.125, mu=0.1, F=1e-3, dt=None, max_steps=600, tol=1e-6)

	case_perm (Ngs=(44, 56), mu=0.1, F=1e-3, dt=80.0, max_steps=3000, tol=1e-6)

	case_cost (Ngs=(64, 96, 128), steps=10, warmup=3)

	case_floor (Ng=96, steps=2, maxit=500, driver="pcg", levels=None)

	case_trace (Ng=96, steps=2, maxit=300, rtol=1e-8, driver="pcg", levels=None)

	main ()

Variables
	HERE = os.path.dirname(os.path.abspath(__file__))

	flow = ensure_flow()

list	ZH_PHI = [0.000125, 0.001, 0.008, 0.027, 0.064, 0.125, 0.216, 0.343, 0.45, 0.5236]

list	ZH_K = [1.096, 1.212, 1.525, 2.008, 2.810, 4.292, 7.442, 15.4, 28.1, 42.1]

	CASES

Detailed Description

WO-M step 2 — what does the solver's fp32 operator storage actually COST us?

The solver stores its multigrid operator (`mac_cutcell_mg.hpp` MReal) and its momentum stencil
(`flow_ibm.hpp` FV) in float. `-DPECLET_FLOW_OPERATOR_DOUBLE` switches both to double. This script is
the A/B: run it against a default build and against a `-DPECLET_FLOW_OPERATOR_DOUBLE` build of the same
commit, and diff the JSON. It measures ACCURACY and COST on cases we actually ship, not on synthetic
ones, because a precision change that costs 2x operator bandwidth for an accuracy nobody needs is
not worth shipping.

Cases
hydro varRho hydrostatic acid test, ratio 1e1..1e6 (walls +-z): steady max|u|, dP/dz error.
vof V2b uniform-velocity identity, ratio 1e1..1e4 (WO-K's decisive gate).
porous porous continuity with a PRESCRIBED uniform drag beta: the steady balance
f = beta*u is exactly solvable, so this measures the porous path against an exact
reference instead of an empirical correlation.
contrast the case that started the campaign: a random close packing (phi 0.63), whose order-2
marching-squares apertures span ~3 decades. Records the FULL residual trace of every
pressure solve (PECLET_FLOW_MG_DEBUG=2) so a floor or a rebound is visible, not just the
iteration count. Rule 3b: a capped solve is INVALID, and this reports the cap directly.
zh Zick & Homsy sphere-array drag vs the published reference (external ground truth).
perm RCP permeability — an actual physical functional, float vs double.
cost ms/step and GPU memory on production-sized beds.

Usage:
PYTHONPATH=<build> python tests/study/precision_ab.py --cases all --json out.json

Function Documentation

◆ gpu_mem_mb()

precision_ab.gpu_mem_mb ( )

Resident GPU memory of THIS process, MiB (0 if no CUDA device / not visible).

Definition at line 47 of file precision_ab.py.

Referenced by case_cost().

◆ rcp_sdf()

precision_ab.rcp_sdf ( Ng )

Definition at line 63 of file precision_ab.py.

Referenced by case_contrast(), case_cost(), case_perm(), and case_trace().

◆ sc_sphere_sdf()

precision_ab.sc_sphere_sdf	(	N,
		phi
	)

Simple-cubic sphere array at solid fraction phi — the Zick & Homsy geometry, identical to
scripts/validate_zick_homsy_sdflow.sc_sdf_xyz so the K it feeds is the recorded one.

Definition at line 79 of file precision_ab.py.

Referenced by case_zh().

◆ case_hydro()

precision_ab.case_hydro	(	ratios = `(1e1, 1e2, 1e3, 1e4, 1e5, 1e6)`,
		N = `8`,
		NZ = `24`,
		g = `0.1`,
		steps = `100`
	)

The hydrostatic acid test: a stratified column at rest must STAY at rest, with the discrete
pressure gradient exactly rho_face*g. Any float defect in the momentum diagonal or in the
projection coefficient shows up here as a spurious velocity.

Definition at line 98 of file precision_ab.py.

◆ case_vof()

precision_ab.case_vof	(	ratios = `(1e1, 1e2, 1e3, 1e4)`,
		n = `32`,
		steps = `200`
	)

WO-K's decisive gate: an arbitrary sharp colour field carried by a UNIFORM velocity must
leave the velocity unchanged. WO-K measured this floored at ~1.2e-7 by the FLOAT momentum
operator (rho_f/dt formed in double in the RHS, stored in float in the stencil) and ~1.2e-15 in
a double build. This re-measures it here so the claim is in one table with its cost.

Definition at line 137 of file precision_ab.py.

◆ case_porous()

precision_ab.case_porous	(	N = `16`,
		beta = `4.0`,
		f_drive = `0.2`,
		eps = `0.6`,
		dt = `0.5`,
		steps = `200`
	)

Porous (volume-averaged) continuity with a PRESCRIBED uniform drag coefficient. At steady
state the momentum balance of a uniform periodic bed is exactly f = beta*u_i (interstitial), so
the reference is analytic — no closure, no empirical correlation, no dem. This is the porous
code path the Ergun fixed-bed benchmark exercises (buildPorousCoeffCons + addDragDiagonal +
projectCorrectPorousCons), measured against an exact answer.

Definition at line 151 of file precision_ab.py.

References mean().

◆ case_contrast()

precision_ab.case_contrast	(	Ngs = `(48, 64, 96)`,
		steps = `4`,
		maxit = `300`,
		rtol = `1e-8`,
		driver = `"pcg"`
	)

THE campaign case: a random close packing whose order-2 (marching-squares) apertures span
about three decades, which is where the collocated session measured float PCG floor at 8e-7 and
REBOUND to 3.8e-5 while a double build converged monotonically to 1e-8.

Reports, per solve: iterations, whether the cap was touched (rule 3b: a capped solve is INVALID,
not degraded), the minimum r/r0 reached and the final r/r0 — a floor shows up as
min == final >> rtol and a rebound as final > min.

Definition at line 185 of file precision_ab.py.

References mean(), and rcp_sdf().

◆ case_zh()

precision_ab.case_zh	(	Ns = `(32, 48, 64)`,
		phi = `0.125`,
		mu = `0.1`,
		F = `1e-3`,
		dt = `None`,
		max_steps = `600`,
		tol = `1e-6`
	)

Zick & Homsy simple-cubic sphere-array drag — the solver's external ground truth.

Definition at line 227 of file precision_ab.py.

References mean(), and sc_sphere_sdf().

◆ case_perm()

precision_ab.case_perm	(	Ngs = `(44, 56)`,
		mu = `0.1`,
		F = `1e-3`,
		dt = `80.0`,
		max_steps = `3000`,
		tol = `1e-6`
	)

The RCP permeability — a physical functional on the high-contrast geometry.

Definition at line 265 of file precision_ab.py.

References mean(), and rcp_sdf().

◆ case_cost()

precision_ab.case_cost	(	Ngs = `(64, 96, 128)`,
		steps = `10`,
		warmup = `3`
	)

ms/step and GPU memory on production-sized high-contrast beds. Operator storage is what the
precision switch changes, so this is where the +28 B/cell of a full fp64 hierarchy shows up.

Definition at line 300 of file precision_ab.py.

References gpu_mem_mb(), and rcp_sdf().

◆ case_floor()

precision_ab.case_floor	(	Ng = `96`,
		steps = `2`,
		maxit = `500`,
		driver = `"pcg"`,
		levels = `None`
	)

ATTAINABLE-ACCURACY FLOOR. Drive the solve with an unreachable tolerance (rtol 1e-14) and let
it run to `maxit`: the residual stops descending where round-off in the CG recurrence takes over,
and that plateau — not the iteration count — is what a precision policy can and cannot move.

Read with the classical bound: CG's attainable relative residual is O(eps * kappa(A)). So the
ratio of the float floor to the double floor should be ~ eps_f32/eps_f64 ~ 5e8 IF storage
precision is the binding term, and the double floor itself measures eps_f64 * kappa, i.e. the
conditioning of the operator — a limit no storage policy can cross. Depth-independence
discriminates the two: a storage/coarsening defect worsens with more levels, a conditioning
limit does not (pass `levels` to check).

Run under PECLET_FLOW_MG_DEBUG=2 and read the [mg] trace; this only sets the configuration up.

Definition at line 341 of file precision_ab.py.

References case_trace().

◆ case_trace()

precision_ab.case_trace	(	Ng = `96`,
		steps = `2`,
		maxit = `300`,
		rtol = `1e-8`,
		driver = `"pcg"`,
		levels = `None`
	)

One high-contrast pressure solve with the operator's own convergence trace. Run under
PECLET_FLOW_MG_DEBUG=2 PECLET_FLOW_MG_DEBUG_SOLVES=<n> and read the `[mg] it ... r/r0` lines:
a FLOOR is min == final >> rtol, a REBOUND is final > min. That distinction is the whole
difference between 'the preconditioner is weak' and 'the operator's A*1=0 defect has stopped
the residual from going anywhere', and an iteration count cannot show it.

Definition at line 358 of file precision_ab.py.

References rcp_sdf().

Referenced by case_floor().

◆ main()

precision_ab.main ( )

Definition at line 391 of file precision_ab.py.

References main().

Referenced by main().

Variable Documentation

◆ HERE

precision_ab.HERE = os.path.dirname(os.path.abspath(__file__))

Definition at line 37 of file precision_ab.py.

◆ flow

precision_ab.flow = ensure_flow()

Definition at line 40 of file precision_ab.py.

◆ ZH_PHI

list precision_ab.ZH_PHI = [0.000125, 0.001, 0.008, 0.027, 0.064, 0.125, 0.216, 0.343, 0.45, 0.5236]

Definition at line 90 of file precision_ab.py.

◆ ZH_K

list precision_ab.ZH_K = [1.096, 1.212, 1.525, 2.008, 2.810, 4.292, 7.442, 15.4, 28.1, 42.1]

Definition at line 91 of file precision_ab.py.

◆ CASES

precision_ab.CASES

Initial value:

1= dict(hydro=case_hydro, trace=case_trace, floor=case_floor, vof=case_vof, porous=case_porous, contrast=case_contrast,

2 zh=case_zh, perm=case_perm, cost=case_cost)

Definition at line 387 of file precision_ab.py.

Functions

Variables