Axis-agnostic dynamics + terrain-following coordinates (Gal-Chen/SLEVE, Schär)

.coveragerc

@@ -1,17 +1,20 @@
 [run]
 source = src
+parallel = True
+concurrency = multiprocessing
+sigterm = True
 omit = 
     src/pybella/data_assimilation/*
     src/pybella/inputs/*
     src/pybella/utils/debug_helpers.py
     src/pybella/utils/slices.py
-    */tests/*
-    */test_*
+    src/pybella/interfaces/postprocessing/strip_target_file.py
+    src/pybella/utils/operators/laplacian/lap3D.py
+    src/pybella/utils/operators/laplacian/lap2D_numba.py
     setup.py
 
 [report]
 exclude_lines =
-    pragma: no cover
     def __repr__
     raise AssertionError
     raise NotImplementedError

.github/workflows/deploy.yml

@@ -1,6 +1,10 @@
 name: deploy
 
-on: [push, pull_request, workflow_dispatch]
+on:
+  push:
+    branches: [develop]
+  pull_request:
+  workflow_dispatch:
 
 permissions:
   contents: write
@@ -36,6 +40,22 @@ jobs:
           # pytest --cov=src --cov-report=xml --cov-report=html --cov-report=term
           pytest ./test_scripts/test_blending.py
           pytest ./test_scripts/test_flow_solver.py
+          pytest ./test_scripts/test_3d_elliptic_oracle.py
+          pytest ./test_scripts/test_igw_analytic.py
+          pytest ./test_scripts/test_3d_coriolis_oracle.py
+          pytest ./test_scripts/test_axes.py
+          pytest ./test_scripts/test_permutation_oracle.py
+          pytest ./test_scripts/test_zvert_smoke.py
+          pytest ./test_scripts/test_terrain_transform.py
+          pytest ./test_scripts/test_terrain_identity_oracle.py
+          pytest ./test_scripts/test_terrain_elliptic_oracle.py
+          pytest ./test_scripts/test_terrain_2d_oracle.py
+          pytest ./test_scripts/test_terrain_resting_atmosphere.py
+          pytest ./test_scripts/test_agnesi_smoke.py
+          pytest ./test_scripts/test_agnesi_analytic.py
+          pytest ./test_scripts/test_agnesi_2d_equivalence.py
+          pytest ./test_scripts/test_schaer_smoke.py
+          pytest ./test_scripts/test_schaer_analytic.py
 
       # - name: Upload coverage HTML report
       #   if: ${{ !env.ACT }}

changelog.d/agnesi_2d_equivalence.added.md

@@ -0,0 +1,6 @@
+Native-2D Agnesi equivalence proof: the 2D terrain path reproduces the
+quasi-2D 3D golden-master run to the pre-existing lap2D/lap3D wall-convention
+floor (~10% on the small wave fields, solver-tolerance independent, same gap
+as a flat wall-bounded impulse) and — the strong gate — passes the Smith
+(1980) analytic oracle with the same thresholds and near-identical
+calibration (w 0.397, u' 0.432, drag 0.981, flux 3.9%) at ~5x less runtime.

changelog.d/axial_phase0_axes_module.added.md

@@ -0,0 +1,5 @@
+Added the axis-geometry module `pybella/utils/axes.py` (role-space convention:
+cyclic permutation mapping (h1, v, h2) roles onto array axes, vertical-axis
+accessors, slab/profile/permutation helpers) with unit tests, plus the
+bit-for-bit H5 run comparator `test_scripts/compare_h5_runs.py` used to gate
+the pure-refactor phases of the axial-agnosticity work. No behaviour change.

changelog.d/axial_phase1_gravity_config.changed.md

@@ -0,0 +1,5 @@
+Axial-agnosticity Phase 1: `gravity_direction` is now a configuration input
+(default 1 = y-vertical, validated; 2D runs require 1), `gravity_strength` and
+`coriolis_strength` are computed role-based from it, and all seven hardcoded
+`ud.gravity_strength[1]` reads route through `axes.vertical_axis(ud)`.
+Bit-identical for all existing cases (verified at tol=0 on full run outputs).

changelog.d/axial_phase2_role_canonical.changed.md

@@ -0,0 +1,9 @@
+Axial-agnosticity Phase 2: Coriolis and explicit dynamics are now
+role-canonical — `multiply_inverse_terms` binds (wh1, wv, wh2) and the
+momentum components in (h1, vertical, h2) role order via the cyclic axis
+permutation (njit kernels unchanged); new `compute_inverse_coefficients`
+exposes the cached H^-1 fields for the upcoming tensor elliptic operator;
+buoyancy and the rhoX stratification coupling act on the configured vertical
+momentum; the explicit momentum rows are written in role symbols with
+expression trees preserved; the advection pwchi special case keys on the
+vertical axis. Bit-identical for all existing cases (verified at tol=0).

changelog.d/axial_phase3_boundaries_hydrostates.changed.md

@@ -0,0 +1,10 @@
+Axial-agnosticity Phase 3: hydrostatic state, vertical profiles, and
+boundary handling are axis-generic — `States` profiles carry their vertical
+axis (`expand_profile`), `integrated_state`/`analytical_state` integrate
+along the configured vertical, the nodal-divergence wall zeroing loops over
+all WALL/RAYLEIGH axes (fixing the previously broken x-WALL elliptic path),
+the gravity ghost-cell handler threads the physical vertical axis (also
+fixing a latent 3D bug where it read `gravity_strength[2] = 0` during
+advection sweeps), `_set_boundary` mirrors the wall-normal momentum of the
+actual wall axis, and the quasi-2D nodal broadcast generalises to any
+degenerate axis. Bit-identical for all existing cases (verified at tol=0).

changelog.d/axial_phase4_lap3d_full_tensor.changed.md

@@ -0,0 +1,13 @@
+Axial-agnosticity Phase 4: the 3D elliptic operator now carries the full
+H^-1 tensor coefficients (C_ij = (Gamma^-1 P Theta) h[role(i),role(j)], the
+same H^-1 the momentum correction applies), replacing the legacy ad-hoc x-z
+`corrf` cross terms; the 3D preconditioner uses the C_ii diagonals. With
+identity H^-1 the operator reduces bit-exactly to the previous one (Oracle
+A), and with full Coriolis a y-uniform 3D solve now matches the trusted 2D
+path under the pseudovector axis mapping to ~1e-6 (new
+test_scripts/test_3d_coriolis_oracle.py). That oracle also exposed and
+fixed the implicit-side sibling of the 2D out-of-plane Coriolis defect (the
+pressure-correction w-row was ndim==3-only). Regenerated golden masters:
+travelling_vortex_3d_coriolis (operator upgrade), igw_baldauf_brdar,
+internal_long_wave, unstable_lamb (w-row fix; the igw analytic-oracle
+out-of-plane error improves 0.060 -> 0.045). All other cases bit-identical.

changelog.d/axial_phase5_permutation_oracle.added.md

@@ -0,0 +1,11 @@
+Axial-agnosticity Phase 5: the permutation oracle
+(`test_scripts/test_permutation_oracle.py`, CI-wired) proves the endgame —
+the 2D internal-long-wave reference (gravity, stratification, walls, full
+Coriolis) embedded as z-vertical (gravity_direction=2) and x-vertical
+(gravity_direction=0) quasi-2D 3D twins reproduces the sigma-mapped 2D
+fields through full solver steps to <=1e-6, with exact uniformity along
+the degenerate axis and the Coriolis pseudovector mapping produced
+automatically by the role-based configuration. Also fixes a latent bug the
+oracle exposed: SpaceDiscr stored `dxyz/ig/ic/stride` as class-level shared
+arrays, so two grids coexisting in one process corrupted each other; they
+are now per-instance. Bit-identical for all existing cases.

changelog.d/axial_phase6_straka_walls_smoke.changed.md

@@ -0,0 +1,8 @@
+Axial-agnosticity Phase 6 wrap-up: Straka now runs its faithful free-slip
+wall configuration on all boundaries, exercising the repaired x-WALL
+elliptic path (target regenerated for the wall config); a z-vertical
+plumbing smoke case (`smoke_zvert` + `test_scripts/test_zvert_smoke.py`)
+pushes `gravity_direction = 2` through the production `-ic` entry path; the
+axis conventions, proof obligations, defects fixed, and known limitations
+are documented in `dev_notes/axial_agnosticity.md`. A permanent z-vertical
+golden-master case is deferred to the terrain-following work.

changelog.d/cross_platform_tolerances.fixed.md

@@ -0,0 +1,6 @@
+Golden-master tolerances recalibrated for cross-platform CI: the first GitHub
+runner pass deviated from locally generated targets by 2.3e-6 (igw rhou) to
+7.0e-5 (Agnesi rhou) — different CPU/BLAS/numba reorder the bicgstab
+reductions, ~100x the same-machine scatter the old gates were tuned to. igw
+returns to the 1e-5 default; the two terrain cases (long elliptic iteration
+chains) gate at 5e-4. Physics remains guarded by the analytic oracles.

changelog.d/fix_2d_coriolis_w_row.fixed.md

@@ -0,0 +1,12 @@
+Fixed the 2D out-of-plane Coriolis defect: the `rhow` momentum row in
+`explicit_euler.do_forward_step` was guarded by `if ndim == 3`, so 2D runs
+applied only the implicit half of the out-of-plane Coriolis rotation. Found
+by the Baldauf-Brdar analytic oracle (out-of-plane velocity error pinned at
+~0.44 rel-L2 independent of dt, sim/ref amplitude ratio ~0.6, with O(f t)
+feedback into u); after the fix the error drops to 0.06 (dt-convergent to
+0.03) and the amplitude ratio to 1.03. Affects only 2D runs with Coriolis
+components in the x/y slots: the `igw_baldauf_brdar` and
+`internal_long_wave` golden-master targets were deliberately regenerated;
+all other cases are bit-identical (Lamb cases use only `strength[2]`, which
+does not enter the w-row). Oracle gates tightened accordingly
+(`test_igw_analytic.GATES`: vo 0.60 -> 0.10).

changelog.d/fix_3d_elliptic_path.fixed.md

@@ -0,0 +1,12 @@
+Fixed the full-3D (`inz > 1`) implicit/elliptic solver path, broken since the
+package restructure: re-wired `lap3D` to the interior-sized `npf` arrays
+(missing imports, coefficient slicing, 3D preconditioner via
+`preconditioner.prepare_diag`), fixed the flat-vector memory layout
+(C-order `[x, y, z]`; the old reshape silently transposed x and z on
+non-cubic grids), fixed the 3D `rhs` shape mismatch and a sign error on the
+x-component of the 3D nodal divergence (legacy, dating to Oct 2021), and made
+the pressure diagnostic kernel dimension-agnostic. The 3D path is validated
+against the 2D solver on a y-uniform quasi-2D problem to ~1e-10
+(`test_scripts/test_3d_elliptic_oracle.py`), and the
+`test_travelling_vortex_3d_coriolis` regression case now runs with a
+committed golden-master target.

changelog.d/igw_analytic_oracle.added.md

@@ -0,0 +1,14 @@
+Added a physics oracle for the Baldauf & Brdar (2013) internal-gravity-wave
+case: `pybella/tests/baldauf_brdar_analytic.py` builds a numerically-exact
+solution of the linearised compressible Euler equations about the isothermal
+background (Bretherton-transformed constant-coefficient system, staggered
+vertical collocation, exact-in-time eigenpropagation per x-Fourier mode;
+energy drift ~1e-13) and evolves the simulation's own initial condition.
+`test_scripts/test_igw_analytic.py` gates the regression configuration
+against it (catching sign/dispersion/amplitude *wrongness*, complementing
+the golden masters which catch *change*) and writes ref/sim/diff PNGs.
+Refinement studies (dt 500->125 s, dx 20->10 km, f on/off) decompose the
+measured sim-vs-linear residual and surface a known solver limitation: in 2D
+runs the out-of-plane momentum receives only the implicit half of the
+Coriolis rotation (`explicit_euler.do_forward_step` skips the `rhow` row for
+`ndim == 2`), pinning that component's error at ~0.44 rel-L2.

changelog.d/lap2d_rayleigh_wall.fixed.md

@@ -0,0 +1,6 @@
+`lap2D_manual` now treats RAYLEIGH boundaries as walls, matching `lap3D`,
+`lap2D_numba` and the divergence slab-zeroing. Previously a sponged top in a
+native-2D run got a periodic-in-y elliptic stencil. No-op for the existing 2D
+golden masters (their vertical handling goes through the atmospheric-extension
+branch — verified bit-identical, and the stable lamb wave still propagates at
+0.99 Cs with stable amplitude at twice the regression horizon).

changelog.d/schaer_ridge_case.added.md

@@ -0,0 +1,8 @@
+Schär (2002) ridge golden-master case (`test_schaer_ridge`): two-scale
+orography (h0 250 m, a 5 km, lambda 4 km) under SLEVE coordinates, native 2D,
+with analytic gradients and the exact cos^2 smooth/residual split. New linear
+FFT mountain-wave oracle (`schaer_linear_analytic.py`, self-tested against
+Smith's closed-form drag) gates the wave field, drag and flux constancy, and
+the Gal-Chen-vs-SLEVE discriminator: spurious small-scale w aloft collapses
+~17x under SLEVE (E_ss 0.004 vs 0.070) where the true lambda-scale response
+is evanescent-dead. Smoke + analytic tests wired into CI.

changelog.d/sleve_transform.added.md

@@ -0,0 +1,8 @@
+SLEVE vertical transform (Schär et al. 2002; Leuenberger et al. 2010
+exponent n): two-scale orography split `ud.orography_smooth` (large-scale)
++ residual against the total `ud.orography`, per-component sinh decay,
+analytic Jacobian — the first eta-dependent J through the operators.
+Selected via `ud.vertical_transform = SLEVETransform(s1, s2, n)`; the
+activation contract and the Gal-Chen path are untouched. Proven by
+transform identities, scale-separation property, and SLEVE-parametrized
+elliptic-composition + resting-atmosphere oracles (3D and native 2D).

changelog.d/straka_density_current.added.md

@@ -0,0 +1,11 @@
+Added the Straka density current (Straka et al. 1993) as regression case
+`test_straka` — the suite's first nonlinear, advection-dominated gravity case
+(256x32 at 200 m, dt = 4 s to t = 900 s; front position, peak winds and
+symmetry verified against the published benchmark). Includes a new explicit
+constant-coefficient diffusion module (`flow_solver/numerics/diffusion.py`,
+enabled per-case via `ud.diffusion` / `ud.diffusion_coeff`, off by default)
+implementing the benchmark's fixed K = 75 m^2/s on velocity and potential
+temperature. The case runs periodic in x: the x-WALL elliptic path is
+currently broken/untested (wall momentum zeroing in the nodal divergence
+handles the vertical axis only) — known limitation, to be fixed with the
+axial-agnosticity refactor.

changelog.d/terrain_2d_native.changed.md

@@ -0,0 +1,6 @@
+Native 2D terrain-following runs: the 2D metric divergence (contravariant/
+J-weighted fluxes), the 2x2 elliptic tensor `J A^T H^-1 A` folded into the
+lap2D cross-term coefficient slots, and a geometry-aware 2D preconditioner
+diagonal. Proven by a 2D div∘correction composition oracle (flat/terrain,
+with/without Coriolis), forced-flat bit-identity, and the native-2D
+resting-atmosphere / conservation gates.

changelog.d/terrain_phase0_transform_scaffolding.added.md

@@ -0,0 +1,6 @@
+Terrain-following coordinates, phase 0: `VerticalTransform`/`GalChenTransform` +
+`MetricFields` scaffolding (`flow_solver/discretisation/terrain.py`), built in
+`grid_init` and attached as `elem.metric`/`node.metric` (`None` without
+`ud.orography` — uniform-Cartesian path untouched). Transform unit tests, h≡0
+identity oracle, and a quasi-2D mountain-wave smoke case (`smoke_agnesi`)
+through the full-tensor 3D elliptic path.

changelog.d/terrain_phase1_metric_divergence.changed.md

@@ -0,0 +1,5 @@
+Terrain phase 1: the nodal divergence computes `J∇·F` when terrain is active —
+J-weighted horizontal fluxes and the contravariant vertical flux
+`θ(mom_v − G1·mom_h1 − G2·mom_h2)` in role space, with unchanged differencing
+stencils. Ghost-slab wall zeroing covers the contravariant fluxes automatically.
+Forced-flat oracle (≤1e-13 vs plain path) and role-wiring checks added.

changelog.d/terrain_phase2_metric_gradients.changed.md

@@ -0,0 +1,5 @@
+Terrain phase 2: pressure gradients map to physical space via the terrain
+gradient matrix A (slope correction of the horizontal rows, 1/J on the
+vertical) in both the explicit forward step and the implicit pressure
+correction; the explicit π update divides the J-weighted divergence by the
+node Jacobian. Bypassed entirely without terrain.

changelog.d/terrain_phase3_elliptic_tensor.changed.md

@@ -0,0 +1,5 @@
+Terrain phase 3: the 3D elliptic operator folds the metric into its tensor —
+`C_ij = wplus ⊙ (J Aᵀ H⁻¹ A)` in role space (lap3D cross terms self-activate)
+and the Helmholtz centre term is J-weighted at nodes. New elliptic oracle
+proves the operator equals the discrete divergence∘momentum-correction
+composition, flat and over an Agnesi hill, to ~1e-12.

changelog.d/terrain_phase4_hydrostates_fields.changed.md

@@ -0,0 +1,6 @@
+Terrain phase 4: hydrostates at physical height. `States` gains a field mode
+(full per-column fields when terrain is active); `analytical_state` evaluates
+its closed form at z(ξ,η) with local dz = J·dη, `integrated_state` gains a
+fine-grid quadrature branch. Resting-atmosphere oracle: balanced atmosphere
+over a 400 m Agnesi hill stays at rest to ~1e-10 m/s (the solve floor),
+identical to flat. `column`/`initial_pressure` assert no terrain.

changelog.d/terrain_phase5_bottom_bc.changed.md

@@ -0,0 +1,6 @@
+Terrain phase 5: gravity-axis ghost cells respect the terrain — the
+CONTRAVARIANT momentum (mom_v − G·mom_h) is reflected oddly and the
+Cartesian vertical momentum rebuilt with the ghost cell's slope terms;
+ghost hydrostatics use the physical image height and local dz = J·dη.
+Advection sweeps flip the metric alongside the solution arrays. Gates:
+uniform flow over a forced-flat metric matches the plain path to 1e-12.

changelog.d/terrain_phase6_advection_cfl.changed.md

@@ -0,0 +1,9 @@
+Terrain phase 6: advection uses contravariant vertical / J-weighted horizontal
+mass fluxes (HLL upwinding follows automatically), the cell update divides by
+J, recovery's Courant velocity divides J back out, and the CFL gains the
+metric vertical signal speed. Orography is evaluated on periodic-wrapped
+coordinates — a non-periodic hill on a periodic axis otherwise makes the
+elliptic system inconsistent at the duplicated nodes (found via dense
+null-space analysis; bicgstab diverged). smoke_agnesi now runs the 400 m
+witch-of-Agnesi hill end-to-end: J-weighted mass/P conserved to 1e-10,
+max |w| ≈ 0.44 m/s vs the ~0.5 m/s linear estimate.

changelog.d/terrain_phase7_sponge_generalisation.changed.md

@@ -0,0 +1,6 @@
+Terrain phase 7: Rayleigh sponge generalised — profile builders read the
+configured vertical axis via `axes.coords_along` (η-based taper, correct in
+terrain-following coordinates), damping broadcasts axis-aware for 3D fields,
+the background Y uses field-mode hydrostates under terrain, and the third
+velocity component is damped in 3D. Lifts the documented vertical=1
+limitation; the 2D path (lamb golden masters) is bit-identical.

changelog.d/terrain_phase8_agnesi_case_oracle.added.md

@@ -0,0 +1,7 @@
+Agnesi hydrostatic mountain-wave case (`test_agnesi_hydrostatic`): N=0.01 1/s,
+U=10 m/s over a 100 m witch-of-Agnesi hill with a=10 km in Gal-Chen
+terrain-following coordinates, RAYLEIGH sponge above 12 km. Golden-master
+regression run (15 steps) plus a Smith-1980 analytic oracle
+(`test_agnesi_analytic`): quasi-steady wave field vs the linear hydrostatic
+solution — momentum flux matches the analytic wave drag to 2%, constant with
+height to 4% below the sponge.

changelog.d/terrain_phase9_ci_wiring.infrastructure.md

@@ -0,0 +1,4 @@
+CI runs the terrain test battery: transform units, h≡0 identity oracle,
+elliptic composition oracle, resting-atmosphere/uniform-flow/mountain-wave
+gates, the Agnesi smoke case, the Agnesi-vs-Smith analytic oracle, and the
+`test_agnesi_hydrostatic` golden master in the flow-solver suite.

src/pybella/flow_solver/discretisation/grid.py

@@ -1,6 +1,7 @@
 import numpy as np
 
 from ...utils import options as opts
+from . import terrain
 
 
 def grid_init(ud):
@@ -35,6 +36,11 @@ def grid_init(ud):
     elem = ElemSpaceDiscr(grid, ud)
     node = NodeSpaceDiscr(grid, ud)
 
+    # terrain metric fields (None when ud has no orography: the uniform-
+    # Cartesian path must stay bit-identical and pay no overhead)
+    elem.metric = terrain.build_metric_fields(elem, ud)
+    node.metric = terrain.build_metric_fields(node, ud)
+
     return elem, node
 
 
@@ -138,9 +144,21 @@ def __init__(self, g):
         assert g.iny >= 1
         assert g.inz >= 1
 
+        # per-instance arrays: the class-level definitions above are shared
+        # buffers, and two grids coexisting in one process (e.g. a 2D
+        # reference and its 3D permutation twin) would corrupt each other
+        self.ig = np.zeros((3))
+        self.ic = np.zeros((3))
+        self.stride = np.zeros((3))
+        self.dxyz = np.zeros((3))
+
         self.ndim = g.ndim
         self.normal = big
 
+        # terrain metric fields (discretisation.terrain.MetricFields);
+        # None == uniform Cartesian — every consumer branches on this
+        self.metric: terrain.MetricFields | None = None
+
         self.igx = self.ig[0] = 2
         self.igy = self.ig[1] = 2 if g.iny > 1 else 0
         self.igz = self.ig[2] = 2 if g.inz > 1 else 0