Cut-cell normal-mode stability machinery + E4_1 alignment with Brady-Livescu 2021 - #30
Merged
Conversation
…ection, round-trip, RK4 containment)
New stencil_gen/cutcell_eval.py: make_cutcell_evaluator(scheme, bc) derives the cut-cell closure once per scheme (lru_cache; mathematica path for zeros-schemes, matching psi_sweep_group_velocity's dispatch), lambdifies every Dirichlet-block entry over (psi, *alpha), and exposes rows() (wall column dropped, offsets 0..T-2, row_indices 0..R-2), dims, alpha_symbols, denominators (irreducible non-constant factors for feasibility pre-checks), and interior(). E4_1 tests are @pytest.mark.slow (derivation ~5.5 min); warm-cache 20-psi evaluation measured ~0.2 ms.
… gks_kreiss substrate defects Four false-unstable mechanisms diagnosed by the blocked 1.2 attempt: sigma_min/det divergence at large |s| for augmented m>p matrices, unbounded NM witness runaway, near-origin anomalies routed to the zero-discrimination classifier, and the uncovered Dirichlet tail-constant kappa=1 structural zero. Dense eigenvalues (N=20..160) confirm E2_1 alpha=0 stable at all five psi.
Fixes the four false-positive mechanisms diagnosed by the blocked 1.2
attempt (E2_1 alpha=0 Dirichlet rows, dense-eigenvalue-confirmed stable):
(a) column-equilibrate M(s) before the SVD on the augmented m>p path
(det-zero locations invariant; sigma_min stays O(1) at large |s|
instead of decaying like 1/(4|s|)) and bound the witness search to
the provable eigensolution region |s| <= S_max (sup-norm/Gershgorin
bound, new eigensolution_s_bound): grid points beyond are masked
unevaluated, Nelder-Mead iterates beyond are rejected (stops the
|s| ~ 1e155 runaway/underflow);
(b) reflect NM probes at Re(s) = eps_imag (no more sub-sonic-axis
count_mismatch anomalies from probes at Re ~ 1e-13);
(c) route near-origin anomalies (|s| < 0.1) through the structural
benign-zero discrimination, never _classify_imag_axis;
(d) generalize the structural test: benign iff the sigma_min right
singular vector aligns with x_const (conservative full-constant
zero, row sums ~ 0) OR x_tail = (0_nodal, e_{kappa->1}) (Dirichlet
tail-constant zero, tail row sums ~ 0), row-sum identity verified;
vanishing-order fit/deflated margin cover the double zero, so the
psi->1 row-coincidence degeneracy no longer trips sigma_tol.
Tests: the 1.2 diagnostic as regression (five psi stable, dense-eig
cross-check at N=40), equilibrated sigma O(1) at |s|=1e3, the 0.2
engineered det-zero still caught, witness |s| <= S_max, NM runaway and
axis-probe regressions, structural triage units; all 0.1-0.4 verdicts
unchanged.
New stencil_gen/reflection.py: mode_split(interior_weights, interior_offsets, s) -> ModeSplit partitions the interior kappa roots into rightgoing (|kappa| < 1 - tol), leftgoing (|kappa| > 1 + tol), and borderline (unit-modulus, resolved by the Trefethen perturbation test); enforces the exact p/p split for centered 2p+1 interiors at Re(s) > 0 and raises DefectiveKappaError on coalescing borderline (sonic) roots. The per-root perturbation logic is factored out of gks_kreiss._classify_imag_axis into the shared helper gks_kreiss.perturbed_root_directions (behavior-preserving refactor) rather than duplicated, per the item text.
New stencil_gen/fully_discrete.py: cached RK4 stability-region boundary
(Newton continuation over theta in [0, 8 pi) — arg R winds 4x around dS
by the argument principle, so the plan's naive 2 pi sweep would trace a
quarter loop and overestimate the margin; closure asserted), containment
margin min_{z in dS} sigma_min(zI - dt L) with the 1-Lipschitz-certified
margin_certified = margin - max_arc_spacing/2, eigenvalue-in-S check
(closed region), and rk4_containment_over_psi reusing the item-3.1
interval assembler at fixed N=40 plus evaluation at the worst
(psi_L, psi_R, N) triple(s) from the item-3.2 scan. Resolvent
quantities only — no matrix powers, no expm, no stepping.
Tests cover the plan's (a) L = -I polyline-exact margin, (b) normal
advection circulant at CFL 0.8/3.0 (with the honest zero margin of the
raw circulant's boundary constant mode), (c) non-normal perturbation
shrinking the margin, boundary anchors/winding/coverage, caching, and
the E2_1 cutcell sweep wiring (E4_1 variant slow-marked).
Post-review corrections to the 4.5 findings entry and calibration script: - Requote all three probe rows' realized psi from the _realized_psi_set replication at the hardcoded (rounded) coordinates, binary-verified against logs/geometry.csv: P_large 0.4988/0.9298, P_mid 0.0695/0.755, P_small 8.84e-4/0.868 (the earlier draft quoted pre-rounding search coordinates; C++ rates/deltas unaffected). - Extend the part-0 scan to all 8 single-alpha +/-0.01 directions and correct the quoted range to 12-20/20 (zero-width conclusion unchanged). - Fix the script docstring/comments to match what is actually quoted. Full script rerun (exit 0) reproduces every reported figure.
Replace the mis-ported one-node-downwind E4_1 cut-cell construction with
a direct symbolic implementation of the published Brady & Livescu 2021
(JCP 426:109794) Sec 3.2/3.3 rules: new derive_cut_cell_paper builds the
5x7 block with row 0 AT the wall, rows centered (wall, f0..f3), a
transition row with no wall support, the paper blend rule
alpha_ij = psi*alphau_ij + (1-psi)*alphau_{i-1,j-1}, and the t+1-equation
cut conservation (fluid-column sums + wall-column equation
w0*B[0,0] = wu^T Bu[:,0] = -1) solved linearly with all 4 uniform alphas
surviving symbolically.
Deleted (no shims, no dual dispatch): derive_cut_cell_mathematica,
construct_cut_cell_with_free_gammas, blend_free_gammas, _safe_get,
zeroed_col_default, zeroed_col_aza, and the zeros path in
derive_cut_cell_scheme (which now forwards zeros schemes to
derive_cut_cell_paper, ending the old 2-param alpha space).
cutcell_eval and group_velocity dispatch switched to the new function.
Acceptance (new slow TestPaperE41Construction in test_e4_cut_cell.py):
blocks/dirichlet slice/weights match tests/fixtures/
e4_1_paper_reference.json at all 5 psi to <=1e-10 (observed ~5e-15);
psi=1 identity B(1) = [conserved uniform closure; interior] is an exact
symbolic zero; wall-column identity -1 exact; Taylor moments hold at the
paper row centers. E2_1/E2_2 verified bit-identical (srepr diff) and the
E2 test subsets pass unchanged; test_temo.py 127/127.
Legacy E4 slow tests anchored to the mis-port now fail as expected;
A7 overhauls them. C++ regen (E4_1.cpp) is a later item.
cutcell_eval.py: document the corrected E4_1 contract on the paper
construction — canonical alpha order alpha_0..alpha_3 = (alphau_04,
alphau_14, alphau_24, alphau_25), the C++ alpha[] contract; Dirichlet rows
at grid nodes 0..3 with integer offsets after the wall-column drop;
denominator structure (psi+1, psi+2, psi+3; psi-independent
uniform-conservation factor 66*alpha_0 - 73; one mixed cut-conservation
polynomial, degree 8 in psi) replacing the dead D6/197-288 lore. Stale
'>170 s' derivation timing updated (~7 s with derive_cut_cell_paper).
gks_kreiss/reflection/interval_operator needed no changes — they consume
the evaluator generically (dims/alpha_symbols) and carry no 2-alpha or
old-layout assumptions.
test_cutcell_eval.py money tests (slow, fixture id=1 alphas from
tests/fixtures/e4_1_paper_reference.json):
(a) kreiss_cutcell_check stable at psi in {0.05, 0.3, 0.7}; at 0.97 the
interior branch-point artifact (witness s=1.17600..., defective kappa,
sigma=inf) is adjudicated via the investigation's deflated winding
(W=0, swap-invariant G=(det M/(k1-k2))^2 on an RHP rectangle) and a
dense fallback;
(b) dense eigenvalues of the assembled N=60 operator: max Re < 0 at all
fixture psi (measured -1.5e-4..-1.6e-5);
(c) default 20-psi sweep runs; recorded verdict: any_unstable=True driven
ENTIRELY by the branch-point artifact at the 6 near-one grid psi, all
adjudicating clean; worst finite margin 2.10e-2 at psi=0.0720; no
genuine bad-psi pocket.
Plus structure pins: denominator factor set asserted exactly, no
denominator zero on the C++ clamp range at alpha=0 or id=1; alpha-order
pin reproducing the fixture Dirichlet blocks at all fixture psi.
Focused run: tests/test_cutcell_eval.py --run-slow -k 'e4 or E4' →
19 passed; full file --run-slow 56 passed; consumer suites
(gks_kreiss/reflection/interval_operator) 192 passed, 5 skipped.
…7/288)
Regenerate src/stencils/E4_1.cpp and E4_1.t.cpp (and their
scripts/stencil_gen/output/ sources) from the paper construction
(derive_cut_cell_paper), replacing the mis-ported 2-alpha scheme.
- std::array<real, 4> alpha in the A2 canonical order
alpha_0..3 = (alphau_04, alphau_14, alphau_24, alphau_25); the plain
copy_zero_padded constructor now takes all four from the Lua alpha table.
- REMOVED the 197/288 constructor guard, its comment block, the
<stdexcept> include, and the psi std::clamp in both nbs methods. A2's
denominator structure ((psi+1)(psi+2)(psi+3), psi-independent
66*alpha_0-73, one degree-8 mixed polynomial) has no psi-pole in (0,1)
at feasible alpha, so no guard/clamp is emitted. A dense scan of the
floating block over (0,1] at the id=1 alphas stays finite
(max|coeff| = 3.09, no pole), confirming the endpoints are safe without
the old clamp.
- Floating nbs = paper 5x7 (wall row + wall column); Dirichlet nbs =
floating rows 1..4 (wall row dropped, wall column kept, 4x7), matching
the operator's split of column 0 into the boundary-coupling B matrix.
- .t.cpp default alpha = fixture id=1 values; asserts the C++ nbs output
equals the committed e4_1_paper_reference.json floating rows at psi in
{0.3, 0.7} plus the Dirichlet block at psi=0.7, margin 1e-8.
Derivation reproduces the fixture to 2.9e-15. t-E4_1 passes; the full
stencils label (11 tests, incl. unchanged E2_1) passes.
Update the Python->C++ L8 bridge for the paper-aligned 4-alpha E4_1 cut
closure and remove the bridge handling built for the mis-port's 197/288
std::invalid_argument guard (the throw no longer exists in src/; grep of
src/ confirms no invalid_argument anywhere).
cpp_bridge.py: drop BridgeResult.infeasible/.reason and the _detect_infeasible
std::invalid_argument classifier; the nonzero-exit path returns a plain
failure result. brady2d_stability.py: cut-template verdict is now
growth_rate < rate_tol (no infeasible gate); layer8 result dict and the
cfl_values aggregate no longer carry infeasible/reason. _L8_CUT_SCHEME_TYPE
is unchanged (E4->make_E4_1); its comment documents the 4 free boundary
alphas (alphau_04/14/24/25) passed via params["alpha"] and that no alpha
vector is infeasible. Tests: remove TestInfeasibleDetection and the
propagation/below-floor cases; E4 cut-dispatch test now exercises 4 alphas;
new slow smoke test runs E4_1 at the paper id=1 alphas end-to-end.
Consumers of the removed fields (the two findings scripts that read
BridgeResult.infeasible/.reason) updated in lockstep so both run again:
- calibration_plan48_4_4.py: fmt()/as_dict() drop the infeasible/reason
keys; part (d) is retargeted from the deleted 197/288 guard probe to the
corrected 4-alpha E4_1 at the fixture id=1 alphas on the cut template
(N in {31, 61}, CFL 0.8), with an E4_1 uncut-template uniform baseline
and the cut-vs-uniform delta. No 197/288 / infeasible content remains.
- calibration_plan48_4_5.py: bridge_row() drops the infeasible key
(nothing downstream read it). Verified by grep that the only remaining
.infeasible accesses in stencil_gen are on CutCellInfeasible objects, a
different type, not BridgeResult.
Cut rung at fixture id=1 alphas, CFL 0.8, t_final=10 (exit=0, no throw):
N=31: cut growth_rate=+0.006461 final_linf=7.47e-3;
uniform baseline growth_rate=-0.001489 -> delta=+0.00795
N=61: cut growth_rate=+0.148840 final_linf=1.36e-3;
uniform baseline growth_rate=+0.001843 -> delta=+0.14700
final_linf falls with N (converging); per plan-48 the reportable cut-cell
signal is the cut-vs-uniform delta, the template's floating-wall growth
having been divided out (template-BC issue not chased here).
Remove everything specific to the retired one-node-downwind E4_1 mis-port from the Python optimizer/sweep stack and reference docs, per the e4-align directive (no shims, no dual dispatch): - stencil_gen/optimizer.py: delete _E4_CLASSICAL_ALPHA1_CPP_FLOOR and _record_cpp_cutcell_diagnostic (the alpha[1] >= 197/288 diagnostic, which referenced a C++ guard that e4-align A5 already removed from src/), and both call sites in run_staged_optimize. Replace the optimizer.py:61-74 stale comment block with one sentence pointing to the paper construction. - sweeps/optimize.py, sweeps/pareto.py, sweeps/bo.py: drop the cpp_cutcell_violates_197_288 plumbing (imports, persistence, docstrings) that consumed the removed diagnostic. - tests/test_optimizer.py, tests/test_sweep_bo.py, tests/test_sweep_pareto.py: delete tests that exercised the removed flag/function directly; lightly correct TestClassicalAlphaBounds docstrings that stated the removed C++ guard as current fact. TestAlphaBasinSurvey is untouched (self-contained via monkeypatched run_staged_optimize, unaffected by the removal). - scripts/stencil_gen/docs/optimization_reference.md, cutcell_normal_mode_reference.md: remove 197/288 and D6 lore, correct the DEFAULT_BOUNDS/extras-schema/basin-survey descriptions, and mark the legacy-closure GKS-instability finding as historical (superseded by e4-align, detail preserved in cutcell_trefethen_theory.md / e4_1_determinant_analysis.md). - docs/reference/py-derivation.md: rewrite the E4_1 cut-cell construction description to the paper convention — derive_cut_cell_ mathematica and its supporting dead code no longer exist; derive_cut_cell_scheme now dispatches internally to derive_cut_cell_paper for zeros-schemes. - group-velocity-analysis skill: checked, does not mention 197/288 or old E4_1 specifics — no change needed. Focused tests: tests/test_optimizer.py tests/test_sweep_bo.py tests/test_sweep_pareto.py -x -q -> 175 passed, 13 skipped. Full stencil_gen suite (non-slow): 1228 passed, 249 skipped, 1 xfailed. Slow TestStagedClassicalAlpha::test_staged_classical_e4_single_seed verified separately with --run-slow -> 1 passed.
Overhaul the E4_1 cut-cell test module to the published Brady & Livescu 2021
(JCP 426:109794) construction, deleting everything that pinned the removed
mis-port.
Deleted (mis-port specifics, no back-compat): the zeros={3,4} uniform hack
(TestE4UniformBoundaryWithZeros), the whole one-node-downwind cut construction
via construct_cut_cell_stencil on E4_1 (TestE4ZeroConstrainedCutCell,
TestE4CutCellConservationSolution, TestE4TEMOConstruction), the 2-parameter
alpha space + {0.1,0.7} default codegen tests (TestE4CodeGeneration,
TestE4TestFileGeneration), the infeasibility / approach-B psi(psi-1)-pole
narratives (test_e4_1_conservation_fails/with_zeros, _conservative_*,
_constant_weights_infeasible_r5, _psi_dependent_conservation_infeasible,
TestCutCellConservationAfterUniform, TestPolynomialStructure,
TestPolynomialBoundaryRows, TestFractionFreeConservation,
TestApproachAInfeasibility, TestPolynomialFullStencil), and the E4_1-specific
2-alpha assertions in TestDeriveCutCellScheme / TestE4CutCellSchemeWithZeros /
TestBuildCutCellConservationSystem.
Kept and re-anchored (construction-agnostic): the generic uniform-boundary
derivation and its conservation (TestE4UniformBoundary,
TestE4UniformConservation — shared with E2, uniform stage matches the paper);
the paper uniform stage (TestMathematicaUniformConservation); the shared
cut-cell conservation machinery via E2_1 (TestBuildCutCellConservationSystem);
and E2_1/E2_2 reproduction plus a new E4_1-forwards-to-paper dispatch check
(TestDeriveCutCellScheme).
Consolidated paper anchors into TestPaperE41Construction against
tests/fixtures/e4_1_paper_reference.json (single module-scoped
derive_cut_cell_paper fixture, folding the former duplicate TestMathematicaWorkflow):
5x7/4x7 dims, wall row 0, transition row 4 (no wall support), 4 surviving
alphas, block/dirichlet/weight match at all 5 psi, psi=1 uniform identity,
wall-column -1 identity, pole-free denominators, and q+1 Taylor moments about
the paper row centers (wall, f0..f3).
2892 -> 744 lines. Focused run: tests/test_e4_cut_cell.py --run-slow -> 38
passed. Full fast suite tests/ -> 1216 passed, 163 skipped.
Review verdict: the A0-A7 alignment is correct (details in the item report). Construction re-derived independently from the paper PDF (pp. 8-11): temo.derive_cut_cell_paper satisfies all paper constraints exactly (5x7 shape, per-row Taylor accuracy, Bd,0 zero pattern, blend rule, psi=0/psi=1 limits, t+1-equation conservation with wall column = -1) and reproduces the FULL printed Table A.4 (35 entries) to 2.8e-14 when fed the corrupted E2_1.db alphas. Fixture id=1 alphas match the independently extracted Appendix B text digit-for-digit. Production C++ E4_1.cpp matches Python at off-fixture psi=0.42 to 1e-15. Fast suite 1215 passed / 163 skipped; slow paper anchors 9/9. Fixes for defects found (residual mis-port references that survived A3/A6, per the review's zero-references gate): - docs/reference/stencils.md: stale-as-current 2-alpha E4_1 lore (alpha[2], the 197/288 constructor throw, psi clamp, "re-add guards" workflow, 361-line guard-section test description) rewritten to the paper-aligned 4-alpha reality (src file is verbatim codegen output). - stencil_gen/benchmarks/alpha_basin_survey.py: remove the dead cpp_cutcell_violates_197_288 plumbing (A6 removed its optimizer producer, so the flag could never be set); 4-column survey table. - tests/test_optimizer.py: drop the flag-propagation test and flag fixtures accordingly; neutralize a 197/288 literal used as an arbitrary params-roundtrip value. - scripts/stencil_gen/calibration_plan48_1_4.py: DELETE - the plan-48 1.4 legacy-closure findings script asserts the retired 2-parameter zeros-path parameterization (derive_cut_cell_scheme now returns 4 alphas, so its assert fires); wholly mis-port-specific; its numbers remain in the plan-48 findings log. - docs/reference/py-brady2d.md, py-derivation.md, ONBOARDING.md: stale claims that the flag is "fully-computed" in optimizer.py, that output/E4_1.cpp diverges from src (now byte-identical), and that the E4_1 conservation xfail still exists (retired in A7).
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
1 participant
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
Summary
Two bodies of work, sequential on this branch:
Plan 48 — cut-cell normal-mode stability machinery (Trefethen program). Adds a theorem-backed stability analysis of ψ-parameterized cut-cell closures to the Python pipeline:
gks_kreiss.py: generalized (augmented) Kreiss determinant for closures with more evolution rows than decaying modes; defective/anomalous κ made observable (was silently "stable"); structural benign-zero discrimination (conservative full-constant and Dirichlet tail-constant classes); column-equilibrated, provably bounded witness search. Fixes a pre-existing bug where L2 fed only p of r closure rows.cutcell_eval.py(cached lambdified closure evaluator),reflection.py(reflection matrices R(s;ψ), per-end A_L/A_R via the mirror convention, Trefethen-1985 round-trip operator),interval_operator.py(two-cut-end interval assembly + resolvent scans),fully_discrete.py(RK4 stability-region containment, resolvent-only).E4 paper alignment. The repo's E4_1 cut construction was a one-node-downwind mis-port of the published construction (no wall row; rows centered at fluid nodes instead of (wall, f0..f3)). This branch replaces it with the Brady & Livescu 2021 (JCP 426:109794) construction exactly:
derive_cut_cell_paperintemo.py: wall row, paper blend rule, t+1-equation cut conservation with the wall-column identity, Dirichlet = drop wall row. All mis-ported paths deleted; E2_1/E2_2 verified bit-identical.E4_1.cppregenerated:std::array<real,4>alphas (αu04, αu14, αu24, αu25), 197/288 guard removed (it was a property of the mis-ported family's denominator; the paper family has no ψ-pole on (0,1] at feasible alphas).tests/fixtures/e4_1_paper_reference.json(id=1 parameters from the paper's Appendix B text; reconstruction verified to ~1e-15 by two independent implementations).Test plan
--run-slow).stencilslabel 11/11;operators3/4 with only the pre-existing, CI-excludedt-laplacianfailure (docs/CLEANUP_PLAN.md §0a).