Axis-agnostic dynamics + terrain-following coordinates (Gal-Chen/SLEVE, Schär)#65
Merged
Conversation
…-row - add Straka density current (test_straka, 256x32 @ 200 m) with new explicit diffusion operator (numerics/diffusion.py, K = 75 m^2/s, off by default) - tolerance audit (3 reps/case): 1e-5 default confirmed; igw momenta/rhoX tightened to 1e-8; utility in test_scripts/tolerance_audit.py - add Baldauf-Brdar linear-reference oracle (tests/baldauf_brdar_analytic.py + test_igw_analytic.py): exact-in-time eigenpropagation of the transformed linear system, evolves the sim's own IC; catches wrongness, not just change - fix 2D out-of-plane Coriolis: rhow row in explicit_euler.do_forward_step was ndim==3-only, halving the rotation; vo error 0.44 -> 0.06 vs oracle; igw + internal_long_wave targets deliberately regenerated, others unchanged - plot 3D comparison fields as transverse mid-slice; CI runs all new tests - known limitation: x-WALL elliptic path broken (Straka runs periodic-x)
- add utils/axes.py: role-space (h1,v,h2) convention, cyclic perms only - add test_axes.py (11 tests) and compare_h5_runs.py bit comparator - no behaviour change
…axes - gravity_direction input (default 1), validated; 2D requires y-vertical - gravity/coriolis strength computed role-based in user_data - all 7 gravity_strength[1] reads -> axes.vertical_axis(ud) - bit-identical: all 10 cases verified at tol=0 vs baselines
- multiply_inverse_terms binds (wh1,wv,wh2) + components in (h1,v,h2) role order via the cyclic axis permutation; njit C11 kernels unchanged - new coriolis.compute_inverse_coefficients exposes the cached H^-1 fields for the upcoming full-tensor elliptic operator - buoyancy + rhoX stratification coupling act on the configured vertical momentum; explicit momentum rows rewritten in role symbols (expression trees preserved); advection pwchi keys on the vertical axis - bit-identical: all 10 cases verified at tol=0 vs baselines
- States profiles carry their vertical axis; integrated/analytical_state integrate along the configured vertical; column/initial_pressure are explicitly 2D-only IC helpers - divergence wall zeroing loops over all WALL/RAYLEIGH axes via wall_slabs, fixing the previously broken x-WALL elliptic path - gravity ghost cells thread the physical vertical axis (fixes latent 3D bug reading gravity_strength[2]=0 during sweeps); _set_boundary mirrors the wall-normal momentum of the actual wall axis - quasi-2D nodal broadcast generalised to any degenerate axis; rayleigh + time_update key on the vertical axis - bit-identical: all 10 cases verified at tol=0 vs baselines
- lap3D carries all nine C_ij = (Gamma^-1 P Theta) h[role(i),role(j)] blocks (same H^-1 as the momentum correction), replacing the ad-hoc x-z corrf cross terms; 3D preconditioner uses the C_ii diagonals - bit-exact reduction to the legacy operator at identity H^-1; new 3D-vs-2D full-Coriolis oracle (pseudovector axis mapping) passes <=1e-6 and is CI-wired - fix implicit-side 2D Coriolis defect the oracle exposed: the pressure correction's w-row was ndim==3-only (sibling of the explicit-step fix) - declare float64 to LinearOperator (scipy's int8 dtype probe crashed numba typing on the new in-place cross-term accumulation) - targets regenerated: tv3d (operator), igw/ilw/unstable_lamb (w-row fix); all other cases bit-identical at tol=0; suite 23/23
- new test_permutation_oracle.py (CI-wired): the 2D internal-long-wave reference run as z-vertical and x-vertical quasi-2D 3D twins reproduces the sigma-mapped fields through full solver steps to <=1e-6, with exact transverse uniformity - Coriolis pseudovector mapping emerges automatically from the role-based config (asserted in-test); cyclic-permutations-only rule validated - fix latent SpaceDiscr bug the oracle exposed: dxyz/ig/ic/stride were class-level shared arrays, corrupting any two grids in one process - all 10 cases bit-identical at tol=0; suite 25/25
- Straka runs free-slip walls on all boundaries (Straka et al. 1993 config), proving the repaired x-WALL elliptic path; target regenerated - smoke_zvert case + test pushes gravity_direction=2 through the production -ic path (registry, prepare, time loop); no golden master - lap2D x-y convention note; CI runs all axial oracles + smoke - suite 26/26; all 10 cases bit-identical at tol=0
- VerticalTransform/GalChenTransform (SLEVE-ready) + MetricFields in discretisation/terrain.py; built in grid_init, role-oriented via axes.py - elem.metric/node.metric attach point; None without ud.orography so the uniform-Cartesian path stays bit-identical with zero overhead - transform unit tests (identities, role orientation, flips, J>0 guard) - h==0 identity oracle pins golden-master cases to the bypass path - smoke_agnesi: quasi-2D x-y-vertical slice through the lap3D path (flat for now; gains the witch-of-Agnesi hill with the metric operators)
- contravariant flux kernel: f_v = theta*(mom_v - G1 mom_h1 - G2 mom_h2),
J-weighted horizontal fluxes; same stencils, role-mapped via elem.metric
- ghost-slab wall zeroing unchanged: fluxes are formed from the momenta,
so zeroed ghosts give zero contravariant flux exactly as the flat path
- native-2D metric divergence explicitly NotImplemented (quasi-2D-via-lap3D)
- oracle: forced-flat metric == plain path to <=1e-13; constant-J/G role
wiring checks for gravity_direction in {0,1,2}
- identity gate green (9 golden masters + IGW analytic + permutation)
- terrain.apply_gradient_map: dp <- A dp in role space (horizontal rows slope-corrected before the vertical row is scaled by 1/J) - wired into explicit_euler.do_forward_step and _correction_nodes — the same A the elliptic C_ij will encode, keeping div/projection consistent - explicit pi update uses rhs/J_n (rhs now carries J*div F under terrain) - oracle: flat metric is a bit-exact identity; constant-J/G algebra checks for all gravity_direction values; identity gate green
- terrain.elliptic_tensor folds J/G into the role-indexed H^-1 before the cij assembly in _prepare_3d_system; with H^-1=I it is the classic TFC tensor, with h=0 it reduces exactly to today's coefficients - wcenter *= J_n (node-exact): every term of the solved equation carries exactly one J, matching the J-weighted rhs divergence - preconditioner flows through unchanged via the existing cii hook - test_terrain_elliptic_oracle: operator == div(correction(p)) composition to ~1e-12 flat AND over an Agnesi hill; forced-flat == plain to 1e-13 - identity gate green (golden masters + analytic + permutation oracles)
- States(field_mode=True) when terrain active: full per-column fields; get_S0c passthrough, get_dSdy differences node S0 over physical node heights then averages to cells (field generalisation of the 1D convolve) - analytical_state: closed form at z(xi,eta), dz = J*deta — reduces bit-identically to the legacy 1D profiles without terrain - integrated_state: fine-grid quadrature branch evaluated at z fields - smoke_agnesi init fixed to the perturbation-pressure convention (p2_nodes = 0, rho = rhoY*S0c) — balanced exactly, flat and terrain - resting-atmosphere oracle: ~1e-10 m/s over a 400 m hill after 5 steps (same as flat floor); field-vs-profile equivalence checks at h=0 - identity gate green
- gravity ghost fill reflects the contravariant momentum mom_v - G.mom_h oddly (free slip through the terrain surface), rebuilds the Cartesian vertical momentum with the ghost cell's own slope terms - ghost hydrostatics: stratification at the physical image height, pressure increment over local dz = J*deta - compute_advection flips elem.metric alongside Sol in every sweep; the BC asserts metric orientation matches the sweep - gate: uniform flow with forced-flat metric == plain run to 1e-12 over full steps; resting atmosphere over the hill stays at the solve floor - identity gate green
- advective mass fluxes: contravariant vertical rhoY*(w - G.u_h), J-weighted horizontal; HLL scales every quantity by flux.rhoY so the metric weighting and upwind sign propagate automatically; cell update divides by J - recovery Courant velocity divides J back out (true coordinate velocity) - CFL: vertical eta_dot = (w - G.u_h)/J and signal speed c*sqrt(1+G^2)/J - orography evaluated on periodic-wrapped coordinates: non-periodic h on a periodic axis made the elliptic system inconsistent at duplicated nodes (dense null-space analysis; bicgstab -> NaN). FD slopes differentiate the wrapped orography so the seam stays consistent - smoke_agnesi carries the 400 m Agnesi hill; gates: J-weighted mass and P conserved to 1e-10 over full steps, max|w| 0.44 m/s vs 0.5 m/s linear - identity gate green (golden masters + analytic + permutation oracles)
- get_tau_y/get_bottom_tau_y read coords_along(., vertical_axis) — the eta-based taper is correct in TFC (surfaces flatten aloft) - rayleigh_damping broadcasts 1D profiles axis-aware (2D x-y unchanged, bit-identical), uses field-mode Y0 under terrain, damps the third velocity component in 3D (mountain-wave sponges must absorb w) - lifts the vertical=1 sponge limitation from axial_agnosticity.md - gates green: lamb + unstable_lamb golden masters, full identity suite
- test_agnesi_hydrostatic: N=0.01/s, U=10 m/s, a=10 km, h0=100 m (Na/U=10, Nh0/U=0.1), Gal-Chen TFC, quasi-2D 128x64, RAYLEIGH sponge above 12 km; golden-master target (15-step spin-up, max-abs ~1e-7 vs 1e-5 tolerance) - agnesi_smith_analytic: Smith steady linear solution (FD-verified) + analytic wave drag; anelastic sqrt(rho0) descaling, x-demeaned (periodic mean-flow response is not part of the infinite-domain solution) - test_agnesi_analytic oracle (96x48, t*U/a=12, window 1-4.5 km): calibrated w 0.40 / u' 0.43 / drag ratio 0.98 / flux constancy 4%; gates 0.55 / 0.60 / [0.85,1.15] / 10% - full-res check at t*U/a=20: drag 1.04, flux constant to 2.7% below 5 km
- test_agnesi_hydrostatic joins the flow-solver golden-master suite - deploy.yml runs the terrain tests: transform units, h==0 identity oracle, elliptic composition oracle, resting-atmosphere/uniform-flow/ mountain-wave gates, Agnesi smoke, Agnesi-vs-Smith analytic oracle - full gate green locally: 57 passed in 15:04 (the Smith oracle is the most expensive single test, ~7.6 min on 16 cores)
- lap2D_manual treated RAYLEIGH as periodic in the stencil wrap; lap3D, lap2D_numba and the divergence slab-zeroing all treat it as WALL - prerequisite for native-2D terrain runs with a sponged top - no-op for existing 2D golden masters (lamb cases use the atmospheric-extension branch; verified bit-identical errors vs target) - extended stable-lamb check: wave still travels at 0.994 Cs with stable amplitude at step 630 (2x the regression horizon)
- 2D contravariant/J-weighted divergence kernel (3D minus the h2 leg) replaces the NotImplementedError quasi-2D guard - elliptic_tensor_2d: G2-free 2x2 restriction of J A^T H^-1 A, folded into lap2D's pre-existing cxy/cyx cross-term coefficient slots - 2D preconditioner diag folds the geometric factors only (legacy 2D keeps H^-1 out of the diagonal; forced-flat preconditions bit-identically) - new test_terrain_2d_oracle.py: div∘correction composition == operator (flat/terrain x with/without Coriolis) <= 1e-12, forced-flat <= 1e-13 - resting-atmosphere battery parametrized over native-2D/quasi-2D-3D; 2D divergence wiring + flat-metric identity tests
- agnesi_smith_analytic comparator made ndim-aware (2D or quasi-2D 3D) - path gate: 15 golden steps 2D vs 3D, gates sit on the measured pre-existing lap2D/lap3D wall-convention floor (~10% of wave fields; flat impulse shows 13% — documented, tolerance-independent) - physics gate: native-2D run passes the Smith oracle with the 3D thresholds (w 0.397 / u' 0.432 / drag 0.981 / flux 3.9%), ~5x faster - CI: test_terrain_2d_oracle + test_agnesi_2d_equivalence wired in
- SLEVETransform: z = eta + h1 b1 + h2 b2, sinh decays per scale, analytic eta-derivative Jacobian, linear below-ground extension (keeps fractional Leuenberger n well-defined at ghost rows) - VerticalTransform.n_components: two-component transforms receive (h1, h2)/(dh1, dh2) tuples; builder splits ud.orography_smooth vs total with the SAME periodic wrap (seam consistency), residual h2 - missing ud.orography_smooth raises (no silent degeneration) - 12 new transform tests incl. scale separation, J > 0 rejection, n > 1 uniform surface Jacobian; elliptic + resting oracles parametrized over the transform (first eta-dependent J end to end)
…ator - test_schaer_ridge: native-2D SLEVE golden master (256x64, 300 s gate), exact cos^2 smooth/residual split with analytic gradients - schaer_linear_analytic: steady linear Boussinesq FFT reference for any periodic terrain (radiating + evanescent branches), drag from the analytic surface flux — reproduces Smith's Agnesi drag to 1.5% - discriminator (480 steps, 128x48, both transforms): E_ss(SLEVE) 0.004 vs E_ss(GalChen) 0.070 at 4-9 km where the true small-scale response is evanescent-dead; SLEVE gates w 0.312 / u' 0.307 / drag 1.10 / flux 4.3%; relative assertions + comparison PNG artifact - registered in IC_MODULES, test_flow_solver, CI smoke + analytic
- on: [push, pull_request] double-ran the suite for PR branches - push now triggers on develop only; PRs and manual dispatch unchanged
- CI runner deviates from locally generated targets ~100x more than same-machine scatter (BLAS/numba reduction ordering in bicgstab) - igw: drop the 1e-8 same-machine-audit gates back to the 1e-5 default (CI showed 2.3e-6 on rhou) - agnesi + schaer: 5e-4 on all fields (CI showed 7.0e-5 / 2.4e-5 on rhou; terrain elliptic chains amplify platform rounding) - determinism gates stay ~20x below real-defect scale; physics is guarded by the B&B / Smith / Schaer analytic oracles
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.
Configurable vertical axis and terrain-following vertical coordinates for the flow solver — 26 commits (axial pts 0–6, terrain pts 0–14), bit-identical on the uniform-Cartesian path.
Axial agnosticity (pts 0–6)
ud.gravity_direction ∈ {0,1,2}; dynamics written in role space (h1, v, h2) viautils/axes.py, cyclic permutations onlyTerrain-following coordinates (pts 0–14)
VerticalTransform; active iffud.orographyis set — without itelem.metric is Noneand the solver is bit-identical to beforeJ AᵀH⁻¹A, field-mode hydrostates, contravariant bottom BC, advection/CFL, sponge — in 3D and native 2D (lap2D cross-term folding; ~5× cheaper than the quasi-2D 3D path)test_agnesi_hydrostatic(+ Smith 1980 oracle, drag ratio 0.98) andtest_schaer_ridge(+ linear FFT oracle; Gal-Chen-vs-SLEVE discriminator: spurious small-scale w aloft collapses ~17× under SLEVE)Validation
~7 min added to CI (Agnesi/Schär oracles); demote to nightly if runners choke.