Functions
	sdf_sphere (x, y, z)

	stokes_u (x, y, z)

	periodic_u (x, y, z)

	phi_man (x, y, z)

	poly_D (t)

	poly_Nc (t)

	poly_Nnb (t)

	build_geo (N, sdf=sdf_sphere, mode="center")

	overlay_cells (geo)

	crossing_points (geo, a, side, cells_flat)

	closure_weights (stf, thf)

	as_order (stf, order)

	row_rescale (geo, order=2)

	gather_face (u, a, m, cells_flat)

	divergence (geo, u3, ubc_fn, rho=None, u_explicit=None)

	assemble (geo, rho, order=2)

	binary_openness_op (geo)

	face_fields (geo, ufun)

	test_extrapolation (Ns)

	test_divergence (Ns)

	solve_system (A, b, n_active)

	test_solve (Ns, mode="center")

	test_probes (N)

	sdf_bed (npz_path, jitter=0.0)

	sdf_pair (gap, rad=0.22)

	solve_set (geo, M)

	probe_split (geo, verbose=True, power_iters=200, tol=1e-6, order=2)

	compat_bias (pr, geo, exact_leftnull=True)

	test_split (args)

	test_degenerate ()

Variables
float	R = 0.3102

	C0 = np.array([0.013, -0.007, 0.004])

int	THETA_MIN = 1e-4

	COUPLED

	QUAD

	LIN

	BC_ONLY

	EXPLICIT

dict	STATE_NAMES

	ap = argparse.ArgumentParser()

	action

	help

	type

	int

	default

	float

	nargs

	choices

	args = ap.parse_args()

list	Ns_eval = [16, 32, 64] if args.quick else [16, 32, 64, 128]

list	Ns_solve = [16, 24, 32] if args.quick else [16, 24, 32, 48]

	s1 = test_extrapolation(Ns_eval)

	s2 = test_divergence(Ns_eval)

	s3 = test_solve(Ns_solve, mode=args.classify)

	s4 = test_probes(args.probe_n)

	s5 = test_degenerate()

list	gates

int	npass = 0

Detailed Description

A-priori validation of the directional ghost-cell IBM projection (plan Phase 0).

The proposed scheme (staggered, point-based, NO openness): the divergence of a fluid-centered
pressure cell uses plain face differences; a face whose staggered velocity point is solid gets its
velocity from the momentum IBM's 1-D wall-anchored quadratic along the face's own axis
(poly_D/poly_Nc/poly_N_nb/poly_Nbc of src/cut_cell_ibm.hpp, reused verbatim):

    poly_D(th) * u_ghost = 2*u_bc + poly_Nc(th)*u_near + poly_N_nb(th)*u_far
    th = sdf_near/(sdf_near - sdf_ghost)  in (0,1], clamped [1e-4, 1]
    (wall at distance th below the near point; ghost point at distance 1; unit grid spacing)

Substituting corrected velocities u = u* - grad(phi) makes the closure implicit in phi: the
Poisson row gains couplings to phi(+/-1), phi(+/-2) along the axis (13-point, nonsymmetric).
Per-row conditioning rescale rho = min(1, min_f D_f) — the D_rescale analog.

Face-state cascade per (cell, axis, side)  [sdf >= 0 fluid]:
  COUPLED   face point fluid AND neighbor center fluid          -> standard +/- (phi_i - phi_nb)
  (sandwich) both face points of the axis solid, center fluid   -> BC_ONLY both sides
  GHOST_QUAD face point solid, near+far sources exist            -> quadratic closure, th in (0,1]
  GHOST_LIN  face point solid, only near source                  -> th*u_g = u_bc + (th-1)*u_near
  SLIVER     face point fluid but neighbor center solid: same quadratic with EXTENDED
             th = 1 + sdf_g/(sdf_g - sdf_beyond) in (1,2)  (evaluation INSIDE the data hull,
             D = th(1+th) > 2: well conditioned); falls back like QUAD->LIN
  BC_ONLY    no usable fluid source (or sandwich)                -> u_face = u_bc, no phi coupling
  EXPLICIT   sliver with no crossing on the u-line               -> face flux = u* (no phi term)

Tests (gates in main):
  1. extrapolation accuracy at ghost faces vs the analytic Stokes-sphere field, anchored at the
     linearized-SDF crossing (the scheme; expect O(h^2)) AND at the exact sphere crossing
     (pure polynomial truncation; expect O(h^3))
  2. closed divergence of the exact solenoidal field on near-IB cells (physical units):
     localized boundary truncation, expect O(h) — the SAME structure as the momentum IBM,
     global 2nd order then comes from elliptic damping and is measured by test 3
  3. assembled sparse solve: u* = u_periodic_exact + discrete_grad(phi_man); solve
     A phi = -div*(u*) with inhomogeneous u_bc = exact wall values; expect phi and the corrected
     velocity 2nd order; verify diagnostic div(u_corr) == residual identically (round-off)
  4. solver probes (dense, small N): A@1 = 0 on active rows; left-null compatibility gap;
     spectrum of the binary-openness-MG-surrogate-preconditioned operator (BiCGStab health);
     deferred-correction rate max|1-lambda|
  5. degenerate geometries: slab channel, sandwich slit, one-cell gap, wall through a face point

This file is the reference implementation for src/ghost_projection.hpp: the C++ gpFillEntry must
reproduce these closure coefficients to float tolerance.

Function Documentation

◆ sdf_sphere()

ghost_projection_apriori.sdf_sphere	(	x,
		y,
		z
	)

Definition at line 55 of file ghost_projection_apriori.py.

◆ stokes_u()

ghost_projection_apriori.stokes_u	(	x,
		y,
		z
	)

Exact Stokes flow past the sphere, U_inf = (1,0,0). Solenoidal, no-slip at r=R.

Definition at line 58 of file ghost_projection_apriori.py.

Referenced by test_extrapolation().

◆ periodic_u()

ghost_projection_apriori.periodic_u	(	x,
		y,
		z
	)

Manufactured periodic solenoidal field (period-1 box), nonzero at the sphere wall.
Each component varies along its own axis (non-degenerate discrete divergence).

Definition at line 70 of file ghost_projection_apriori.py.

◆ phi_man()

ghost_projection_apriori.phi_man	(	x,
		y,
		z
	)

Definition at line 78 of file ghost_projection_apriori.py.

Referenced by compat_bias(), test_probes(), and test_solve().

◆ poly_D()

ghost_projection_apriori.poly_D ( t )

Definition at line 84 of file ghost_projection_apriori.py.

Referenced by closure_weights().

◆ poly_Nc()

ghost_projection_apriori.poly_Nc ( t )

Definition at line 87 of file ghost_projection_apriori.py.

Referenced by closure_weights().

◆ poly_Nnb()

ghost_projection_apriori.poly_Nnb ( t )

Definition at line 90 of file ghost_projection_apriori.py.

Referenced by closure_weights().

◆ build_geo()

ghost_projection_apriori.build_geo	(	N,
		sdf = `sdf_sphere`,
		mode = `"center"`
	)

Classify every (cell, axis, side) face. Returns dict with sdf samples, states, thetas.

mode="center" (the shipped scheme): a cell carries a pressure unknown iff its CENTRE is fluid,
and a face is COUPLED iff the face point is fluid AND both centres are fluid. A fluid face
whose neighbour centre is solid therefore has nowhere to couple, which is what the SLIVER /
EXTENDED-theta / EXPLICIT branches exist to paper over.

mode="face" (the proposal, and the rule in doc/Robust_Scaled_IBM_Solver.tex: "pressure values
are required in cells where any face contains a fluid velocity"): a cell carries a pressure
unknown iff at least one of its six faces has a FLUID velocity point, and a face is COUPLED iff
its own velocity point is fluid. Every fluid face then has a live unknown on both sides, so the
sliver branches vanish; the wall-anchored closure is used only where it is really needed, at a
SOLID face point. Solid-centred cells that own fluid faces (the throats threading between
spheres) get their continuity equation back.

Definition at line 99 of file ghost_projection_apriori.py.

Referenced by test_degenerate(), test_divergence(), test_extrapolation(), test_probes(), test_solve(), and test_split().

◆ overlay_cells()

ghost_projection_apriori.overlay_cells ( geo )

Active cells with at least one non-COUPLED face.

Definition at line 213 of file ghost_projection_apriori.py.

Referenced by test_divergence(), and test_solve().

◆ crossing_points()

ghost_projection_apriori.crossing_points	(	geo,
		a,
		side,
		cells_flat
	)

Definition at line 221 of file ghost_projection_apriori.py.

Referenced by divergence(), and test_extrapolation().

◆ closure_weights()

ghost_projection_apriori.closure_weights	(	stf,
		thf
	)

(w_bc, w_n1, w_n2, D) for QUAD/LIN/BC_ONLY flat state/theta arrays.

Definition at line 241 of file ghost_projection_apriori.py.

References poly_D(), poly_Nc(), and poly_Nnb().

Referenced by assemble(), divergence(), row_rescale(), and test_extrapolation().

◆ as_order()

ghost_projection_apriori.as_order	(	stf,
		order
	)

Closure ORDER selection, mirroring gpOrderWeights: order 1 evaluates a QUAD face with the
linear closure (that is what matrix_order=1 does to the implicit phi couplings).

Definition at line 261 of file ghost_projection_apriori.py.

Referenced by assemble(), and row_rescale().

◆ row_rescale()

ghost_projection_apriori.row_rescale	(	geo,
		order = `2`
	)

rho = min(1, min over ghost faces of D_f) per cell (flat array). D comes from the MATRIX
weights, so `order` is matrix_order.

Definition at line 270 of file ghost_projection_apriori.py.

References as_order(), and closure_weights().

Referenced by probe_split(), test_degenerate(), test_probes(), and test_solve().

◆ gather_face()

ghost_projection_apriori.gather_face	(	u,
		a,
		m,
		cells_flat
	)

Definition at line 283 of file ghost_projection_apriori.py.

Referenced by divergence(), and test_extrapolation().

◆ divergence()

ghost_projection_apriori.divergence	(	geo,
		u3,
		ubc_fn,
		rho = `None`,
		u_explicit = `None`
	)

Closed point divergence (grid units: sum of face differences) on active cells,
row-rescaled by rho. u3 = 3 face fields; ubc_fn(x,y,z)->(3 components); u_explicit
supplies the field read at EXPLICIT faces (defaults to u3).

Definition at line 286 of file ghost_projection_apriori.py.

References closure_weights(), crossing_points(), float, and gather_face().

Referenced by compat_bias(), test_divergence(), test_probes(), and test_solve().

◆ assemble()

ghost_projection_apriori.assemble	(	geo,
		rho,
		order = `2`
	)

Sparse A (N^3 x N^3): binary-openness base + closure deltas, overlay rows scaled by rho.
Inactive rows = identity. Convention: A phi = -div(u*) (positive diagonal).

Definition at line 330 of file ghost_projection_apriori.py.

References as_order(), closure_weights(), and float.

Referenced by probe_split(), test_degenerate(), test_probes(), and test_solve().

◆ binary_openness_op()

ghost_projection_apriori.binary_openness_op ( geo )

The symmetric MG surrogate: 7-point op with o=1 on COUPLED faces, 0 otherwise.

Definition at line 380 of file ghost_projection_apriori.py.

Referenced by probe_split(), and test_probes().

◆ face_fields()

ghost_projection_apriori.face_fields	(	geo,
		ufun
	)

Definition at line 408 of file ghost_projection_apriori.py.

Referenced by compat_bias(), test_divergence(), test_extrapolation(), test_probes(), and test_solve().

◆ test_extrapolation()

ghost_projection_apriori.test_extrapolation ( Ns )

QUAD ghost-face closures vs the smooth continuation of the analytic Stokes field.
Variant 'scheme':     u_bc = 0 (what the solver knows)      -> O(h^2) (wall-anchoring error)
Variant 'consistent': u_bc = field at the linearized anchor -> O(h^3) (pure poly truncation)

Definition at line 411 of file ghost_projection_apriori.py.

References build_geo(), closure_weights(), crossing_points(), face_fields(), float, gather_face(), int, and stokes_u().

◆ test_divergence()

ghost_projection_apriori.test_divergence ( Ns )

Definition at line 458 of file ghost_projection_apriori.py.

References build_geo(), divergence(), face_fields(), float, and overlay_cells().

◆ solve_system()

ghost_projection_apriori.solve_system	(	A,
		b,
		n_active
	)

Definition at line 484 of file ghost_projection_apriori.py.

◆ test_solve()

ghost_projection_apriori.test_solve	(	Ns,
		mode = `"center"`
	)

Definition at line 495 of file ghost_projection_apriori.py.

References assemble(), build_geo(), divergence(), face_fields(), float, int, overlay_cells(), phi_man(), and row_rescale().

◆ test_probes()

ghost_projection_apriori.test_probes ( N )

Definition at line 584 of file ghost_projection_apriori.py.

References assemble(), binary_openness_op(), build_geo(), divergence(), face_fields(), float, int, phi_man(), and row_rescale().

◆ sdf_bed()

ghost_projection_apriori.sdf_bed	(	npz_path,
		jitter = `0.0`
	)

Periodic union-of-spheres SDF from a pack_bed.py packing, mapped onto the unit box.

Requires a CUBIC packing box (s100/s101/s108); the harness grid is the unit cube, so the
sphere radius in cells is N/box. Periodic images are included, so the geometry the harness
classifies is the same one the solver sees on that bed at that resolution.

Definition at line 631 of file ghost_projection_apriori.py.

References float.

Referenced by test_split().

◆ sdf_pair()

ghost_projection_apriori.sdf_pair	(	gap,
		rad = `0.22`
	)

Two spheres on the x-axis separated by `gap` (in units of the box side), periodic images
included. gap -> 0 is the near-tangent pathological configuration.

Definition at line 662 of file ghost_projection_apriori.py.

Referenced by test_split().

◆ solve_set()

ghost_projection_apriori.solve_set	(	geo,
		M
	)

The cells the solver actually solves on: active, phi-coupled, and in the LARGEST connected
component of the COUPLED graph (flow's fragmentation guard, flow_ibm.hpp set_solid).

Definition at line 675 of file ghost_projection_apriori.py.

Referenced by probe_split().

◆ probe_split()

ghost_projection_apriori.probe_split	(	geo,
		verbose = `True`,
		power_iters = `200`,
		tol = `1e-6`,
		order = `2`
	)

A4: rho(S^-1 N) by power iteration, S = binary-openness 7-point op, N = gp overlay delta.
A5: the compatibility bias of the nonsymmetric A under the solver's mean removal.
Sparse throughout, so it runs on beds the dense probe cannot touch.

Definition at line 689 of file ghost_projection_apriori.py.

References assemble(), binary_openness_op(), float, int, row_rescale(), and solve_set().

Referenced by test_split().

◆ compat_bias()

ghost_projection_apriori.compat_bias	(	pr,
		geo,
		exact_leftnull = `True`
	)

A5 (second half): the per-solve bias the mean removal leaves behind. Build the physical RHS
b = -div(u*) the solver would see, project it the way the solver does (remove the mean), and
measure the component that the TRUE left null vector still sees -- that part is unreachable by
any Krylov iteration and is what the incremental-rotational pressure accumulates.

Definition at line 809 of file ghost_projection_apriori.py.

References divergence(), face_fields(), float, and phi_man().

Referenced by test_split().

◆ test_split()

ghost_projection_apriori.test_split ( args )

Phase-A4/A5 driver: run the sparse split + compatibility probes on the requested geometry.

Definition at line 847 of file ghost_projection_apriori.py.

References build_geo(), compat_bias(), probe_split(), sdf_bed(), and sdf_pair().

◆ test_degenerate()

ghost_projection_apriori.test_degenerate ( )

Definition at line 873 of file ghost_projection_apriori.py.

References assemble(), build_geo(), float, int, and row_rescale().

Variable Documentation

◆ R

float ghost_projection_apriori.R = 0.3102

Definition at line 51 of file ghost_projection_apriori.py.

◆ C0

ghost_projection_apriori.C0 = np.array([0.013, -0.007, 0.004])

Definition at line 52 of file ghost_projection_apriori.py.

◆ THETA_MIN

int ghost_projection_apriori.THETA_MIN = 1e-4

Definition at line 53 of file ghost_projection_apriori.py.

◆ COUPLED

ghost_projection_apriori.COUPLED

Definition at line 94 of file ghost_projection_apriori.py.

◆ QUAD

ghost_projection_apriori.QUAD

Definition at line 94 of file ghost_projection_apriori.py.

◆ LIN

ghost_projection_apriori.LIN

Definition at line 94 of file ghost_projection_apriori.py.

◆ BC_ONLY

ghost_projection_apriori.BC_ONLY

Definition at line 94 of file ghost_projection_apriori.py.

◆ EXPLICIT

ghost_projection_apriori.EXPLICIT

Definition at line 94 of file ghost_projection_apriori.py.

◆ STATE_NAMES

dict ghost_projection_apriori.STATE_NAMES

Initial value:

1= {COUPLED: "COUPLED", QUAD: "QUAD", LIN: "LIN",

2 BC_ONLY: "BC_ONLY", EXPLICIT: "EXPLICIT"}

Definition at line 95 of file ghost_projection_apriori.py.

◆ ap

ghost_projection_apriori.ap = argparse.ArgumentParser()

Definition at line 905 of file ghost_projection_apriori.py.

◆ action

ghost_projection_apriori.action

Definition at line 906 of file ghost_projection_apriori.py.

◆ help

ghost_projection_apriori.help

Definition at line 906 of file ghost_projection_apriori.py.

◆ type

ghost_projection_apriori.type

Definition at line 907 of file ghost_projection_apriori.py.

◆ int

ghost_projection_apriori.int

Definition at line 907 of file ghost_projection_apriori.py.

Referenced by probe_split(), test_degenerate(), test_extrapolation(), test_probes(), and test_solve().

◆ default

ghost_projection_apriori.default

Definition at line 907 of file ghost_projection_apriori.py.

◆ float

ghost_projection_apriori.float

Definition at line 912 of file ghost_projection_apriori.py.

Referenced by assemble(), compat_bias(), divergence(), probe_split(), sdf_bed(), test_degenerate(), test_divergence(), test_extrapolation(), test_probes(), and test_solve().

◆ nargs

ghost_projection_apriori.nargs

Definition at line 914 of file ghost_projection_apriori.py.

◆ choices

ghost_projection_apriori.choices

Definition at line 916 of file ghost_projection_apriori.py.

◆ args

ghost_projection_apriori.args = ap.parse_args()

Definition at line 922 of file ghost_projection_apriori.py.

◆ Ns_eval

list ghost_projection_apriori.Ns_eval = [16, 32, 64] if args.quick else [16, 32, 64, 128]

Definition at line 928 of file ghost_projection_apriori.py.

◆ Ns_solve

list ghost_projection_apriori.Ns_solve = [16, 24, 32] if args.quick else [16, 24, 32, 48]

Definition at line 929 of file ghost_projection_apriori.py.

◆ s1

ghost_projection_apriori.s1 = test_extrapolation(Ns_eval)

Definition at line 931 of file ghost_projection_apriori.py.

◆ s2

ghost_projection_apriori.s2 = test_divergence(Ns_eval)

Definition at line 932 of file ghost_projection_apriori.py.

◆ s3

ghost_projection_apriori.s3 = test_solve(Ns_solve, mode=args.classify)

Definition at line 933 of file ghost_projection_apriori.py.

◆ s4

ghost_projection_apriori.s4 = test_probes(args.probe_n)

Definition at line 934 of file ghost_projection_apriori.py.

◆ s5

ghost_projection_apriori.s5 = test_degenerate()

Definition at line 935 of file ghost_projection_apriori.py.

◆ gates

list ghost_projection_apriori.gates

Initial value:

=  [
        ("extrapolation order (consistent)  >= 2.5", s1["consistent"] >= 2.5),
        ("extrapolation order (scheme)      >= 1.7", s1["scheme"] >= 1.7),
        ("near-IB divergence order          >= 0.8", s2 >= 0.8),
        ("solve: phi order                  >= 1.5", s3["o_phi"] >= 1.5),
        ("solve: corrected-velocity order   >= 1.7", s3["o_u"] >= 1.7),
        ("diagnostic == residual identity", s3["ident"]),
        ("preconditioned spectrum Re > 0", s4["re_min"] > 0.0),
        ("compatibility gap < 1e-3", s4["gap"] < 1e-3),
        ("degenerate geometries sane", s5),
    ]

Definition at line 938 of file ghost_projection_apriori.py.

◆ npass

int ghost_projection_apriori.npass = 0

Definition at line 949 of file ghost_projection_apriori.py.

Functions

Variables

Detailed Description

Function Documentation

◆ sdf_sphere()

◆ stokes_u()

◆ periodic_u()

◆ phi_man()

◆ poly_D()

◆ poly_Nc()

◆ poly_Nnb()

◆ build_geo()

◆ overlay_cells()

◆ crossing_points()

◆ closure_weights()

◆ as_order()

◆ row_rescale()

◆ gather_face()

◆ divergence()

◆ assemble()

◆ binary_openness_op()

◆ face_fields()

◆ test_extrapolation()

◆ test_divergence()

◆ solve_system()

◆ test_solve()

◆ test_probes()

◆ sdf_bed()

◆ sdf_pair()

◆ solve_set()

◆ probe_split()

◆ compat_bias()

◆ test_split()

◆ test_degenerate()

Variable Documentation

◆ R

◆ C0

◆ THETA_MIN

◆ COUPLED

◆ QUAD

◆ LIN

◆ BC_ONLY

◆ EXPLICIT

◆ STATE_NAMES

◆ ap

◆ action

◆ help

◆ type

◆ int

◆ default

◆ float

◆ nargs

◆ choices

◆ args

◆ Ns_eval

◆ Ns_solve

◆ s1

◆ s2

◆ s3

◆ s4

◆ s5

◆ gates

◆ npass