Discrete Hydrostatic Consistency in horiz_press_grad

Creation date: 2026/07/28

Contributors: Xylar Asay-Davis, Claude

Summary

The Omega design document for the higher-order horizontal pressure gradient (components/omega/doc/design/PGradHighOrder.md, §2.3 and §3.7) requires discrete hydrostatic consistency: for horizontally uniform conservative temperature and absolute salinity, with coordinate surfaces tilted arbitrarily between adjacent columns, the scheme must return exactly zero horizontal pressure-gradient acceleration (HPGA), to machine precision, with no background reference state. Its §5.2 makes that a gating test and names the Polaris two-column task as one of the places it should be exercised.

This design adds two variants to the ocean/horiz_press_grad task family that measure that property, together with the reference-free analysis they need. The new capability is:

• Two exact resting states in which the true HPGA is identically zero, so the model’s HPGA is the error and no reference solution is required: hydrostatic_consistency (tilted coordinate, flat sea floor) and bathymetry_step (flat coordinate, stepped sea floor).

• A tilt scan: at fixed resolution, the tilt is swept over roughly three decades and the exponent \(q\) in \(\lVert\text{HPGA}\rVert \sim (\text{tilt})^{q}\) is measured. The exponent, not a single pass/fail, is the discriminating output.

• A reference-free analysis step, RestingAnalysis, which compares against zero rather than against ReferenceColumn, retains the Omega-vs-Polaris consistency check, and guards the resting-state premise itself.

The primary challenge is not the numerics but the severity of the configuration. The existing ztilde_gradient variant is already structurally a resting state with a tilted coordinate, and the centered scheme passes it with an RMS HPGA of \(\sim 7\times10^{-12}\) m s\(^{-2}\) — six orders of magnitude below the \(\sim 3\times10^{-6}\) m s\(^{-2}\) bottom-layer error seen in realistic global Omega runs. A new variant that inherits that configuration would pass trivially and measure nothing. The design is successful if the centered scheme is driven to an RMS HPGA of at least \(10^{-6}\) m s\(^{-2}\) somewhere in the sweep (the sensitivity gate), if the measured tilt exponent is stable and interpretable, and if the resulting diagnostic is ready to be pointed at the FiniteVolume scheme when it exists.

Two stages, and a note on the word “phase”

This design was written and first implemented before Omega’s FiniteVolume scheme existed, and extended once it did. The two Polaris stages are named rather than numbered, because PGradHighOrder.md uses Phase 1 and Phase 2 for the scheme’s order — Phase 1 being HorzOrder: 2 with a linear vertical reconstruction, Phase 2 order 4 with PPM — and two independent pairs of ordinals in one document is a reliable source of confusion.

stage	what it delivered
the calibration stage	the two resting-state variants, the reference-free analysis, and the Centered baseline measured in Omega: tilt exponent \(q = 1.000\), sensitivity gate passing by 21–23×
the rollout stage	scheme selection and the §5.1 gates set from a Centered-versus-FiniteVolume comparison

Everything here concerns Omega Phase 1. “Phase” below always means the Omega sense.

The calibration stage’s output is what makes the rollout readable: with the centered exponent measured at 1.000 to five digits, a FiniteVolume measurement is read against a known-good value rather than against an expectation.

What the rollout stage adds

• A scheme axis. Each variant runs both Centered and FiniteVolume over one shared initial condition, so the two are compared at an identical state.

• A second Python kernel to compare against. init.py writes HPGAFiniteVolume alongside HPGA, and each analysis step compares Omega against whichever matches the scheme that produced its output.

The two variants deliberately isolate two different mechanisms, which the sizing work below showed are not the same mechanism:

	hydrostatic_consistency	bathymetry_step
what tilts	the p-star reference grid, via z_tilde_bot_grad	nothing; only the sea floor moves
where the error lives	every layer, growing with depth	the deepest layer valid in both columns, only
behaviour vs. the swept parameter	a clean power law	a staircase, set by where maxLevelCell changes
what it represents	thin, steeply sloped layers in a general ALE coordinate	the bottom partial cell at a bathymetry step

Requirements

Requirement: The test configuration is an exact resting state

Date last modified: 2026/07/28

Contributors: Xylar Asay-Davis, Claude

Each new variant shall be a configuration whose true, continuous HPGA is identically zero at every point in the water column, so that the model’s computed HPGA is entirely error and no reference solution is needed. The resting-state property shall be a consequence of the configuration (horizontally uniform conservative temperature and absolute salinity specified as functions of pseudo-height, with no surface-pressure gradient), not an approximation.

Requirement: The geometry under test is verified, not assumed

Date last modified: 2026/07/28

Contributors: Xylar Asay-Davis, Claude

The analysis shall confirm at run time that the geometry actually built is the geometry the configuration asked for, and shall fail loudly if it is not. A sweep whose swept parameter has been silently altered by the vertical-grid construction reports an error level for a configuration that was never tested, in either direction — too large or too small.

Requirement: The test reaches the severity at which the centered scheme fails

Date last modified: 2026/07/28

Contributors: Xylar Asay-Davis, Claude

At its most severe configuration, the sweep shall drive the existing centered scheme to an RMS HPGA of at least \(10^{-6}\) m s\(^{-2}\) — the order of the bottom-layer HPGA error observed in realistic global Omega runs with a uniform temperature and salinity state. A sweep that does not reach this level is not exercising the failure mode the higher-order scheme is meant to fix, and any conclusion drawn from it about the new scheme would be unsupported.

Requirement: The measurement discriminates between mechanisms, not just pass and fail

Date last modified: 2026/07/28

Contributors: Xylar Asay-Davis, Claude

The analysis shall report the exponent \(q\) in \(\lVert\text{HPGA}\rVert \sim (\text{tilt})^{q}\), restricted to the portion of the sweep where the result is well above the round-off floor, in addition to absolute error levels. A single pass/fail number cannot distinguish “the scheme represents specific volume as piecewise constant in pressure” from “the scheme has a higher-order residual” from “the scheme is exact”, and those are the outcomes the test exists to separate.

Requirement: The deepest layer valid in both columns is measured, not discarded

Date last modified: 2026/07/28

Contributors: Xylar Asay-Davis, Claude

The error metric shall include the deepest layer that is valid in both columns bounding the edge. That layer is the bottom partial cell, and it is where the pressure-gradient error concentrates; excluding it removes the signal these tests exist to measure.

This applies to the reference-based Analysis step as well as to the new resting-state analysis. Analysis currently drops that layer, which is a leftover from the previous five-column finite-difference reference whose stencil could not be formed where a neighbouring column’s collocation point fell below its own bathymetry. The reference is now a single analytic column at the edge, documented as valid to the seafloor, so the exclusion no longer has a basis. Fixing it is in scope here because it is the same decision, and because leaving the two analysis steps disagreeing about which layers count would be incoherent.

Requirement: The existing variants and their pass criteria are untouched

Date last modified: 2026/08/02

Contributors: Xylar Asay-Davis, Claude

The four existing variants (salinity_gradient, surface_pressure_gradient, temperature_gradient, ztilde_gradient), the ReferenceColumn evaluator, and the convergence-rate and reference thresholds in horiz_press_grad.cfg shall behave exactly as they do today. The new capability shall not change their step names, work directories, or results.

Amended for the rollout stage. This requirement held for the calibration stage and is deliberately relaxed by the rollout, in two respects:

• the scheme axis puts a scheme label in every forward step’s directory, so the step names do change and the recorded centered baselines have to be regenerated. This is the cost of running both schemes over one shared Init rather than duplicating the initial condition per scheme;

• the convergence bands become per-variant and per-scheme. A single band cannot serve four variants whose measured centered rates are 1.6, 1.8, 2.0 and 1.1.

What the rollout does not change is the four variants’ results under Centered, or the meaning of HPGA. Any change there is a regression.

Requirement: The cross-implementation check compares the same algorithm

Date last modified: 2026/08/02

Contributors: Xylar Asay-Davis, Claude

omega_vs_polaris_rms_threshold is the only check in the family that compares Omega’s arithmetic against an independent implementation of the same scheme, and it is worth exactly as much as the claim that the two sides are the same algorithm. Where the Python kernel makes a discretization choice that Omega makes differently — the quadrature rule and its order, which end of the column the scan is anchored at, how the reconstruction is closed at the ends of a column, how the edge-shared expansion point is formed — the check stops measuring implementation disagreement and starts measuring the gap between two different algorithms, which is not distinguishable from a bug in either.

Every such choice shall therefore be reconciled against PGradHighOrder.md before the check is run in anger, and each resolution recorded with its provenance. Two constraints on how:

1. The design is the source, not the C++. The Python kernel was written from PGradHighOrder.md precisely so that this comparison has content. Reading the Omega implementation to settle a question of intent the design leaves open is legitimate; transcribing it is not, because a Python kernel derived from the C++ would compare an implementation against itself.

2. Reconcile up front rather than by bisection. The check produces one scalar per resolution. It can report that the two sides disagree; it cannot report which of several choices disagrees, and several of them first bite in the same configurations.

Where a reconciliation finds that the Python follows the design and Omega does not, it stops being a Polaris change and is handed back rather than papered over by adjusting the Python to match.

Algorithm Design

Algorithm Design: Why these configurations are exact resting states

Date last modified: 2026/07/28

Contributors: Xylar Asay-Davis, Claude

Omega’s vertical coordinate is pseudo-height \(\tilde z = -p/(\rho_0 g)\), so pressure is an exactly linear function of \(\tilde z\) and a surface of constant \(\tilde z\) is a surface of constant \(p\). Init interpolates conservative temperature \(\Theta\) and absolute salinity \(S_A\) against \(\tilde z\) from piecewise profiles whose node values and node pseudo-heights are set by *_mid and *_grad configuration options. Setting every *_grad to zero makes \(\Theta\) and \(S_A\) the same functions of \(\tilde z\), and therefore of \(p\), in both columns. Specific volume is then a function of pressure alone,

\[ \alpha = \alpha\bigl(\Theta(p), S_A(p), p\bigr) \equiv \alpha(p), \]

isobars are level surfaces, and the true HPGA is identically zero everywhere in the fluid — for any sea-floor shape and any tilt of the coordinate surfaces. This is what makes both variants valid resting states, and it is the discrete analogue of the property §3.7 of the Omega design asserts.

Two configuration choices are therefore forced:

• temperature_grad, salinity_grad and z_tilde_grad are all zero (inherited from the base horiz_press_grad.cfg), and surface_pressure_grad is zero, so the sea surface is flat.

• The deepest profile node lies below the deepest reference-grid interface in either column, so the PCHIP interpolant is never asked to extrapolate when the tilt makes one column’s reference grid deeper than the other’s.

Algorithm Design: What actually tilts the coordinate in a p-star column

Date last modified: 2026/07/28

Contributors: Xylar Asay-Davis, Claude

This is the design’s central finding, and it determined the shape of both variants.

init_pstar_vertical_coord builds each column by clipping a reference 1-D grid at that column’s BottomPressure. In Init, the reference grid for column \(i\) is uniform over \([0, -\tilde z^{\text{bot}}_i]\) with a level count fixed by the mid-column value, so z_tilde_bot_grad stretches one column’s reference grid relative to the other’s. The cross-edge offset of interface \(k\) then grows linearly with depth,

\[ \Delta \tilde z_k \;=\; \frac{k}{N}\, \bigl(\tilde z^{\text{bot}}_1 - \tilde z^{\text{bot}}_0\bigr), \]

which is the tilted-coordinate configuration §3.7 describes. This is the only mechanism in the two-column task that tilts interior coordinate surfaces.

A sea-floor step does not tilt them. With z_tilde_bot_grad = 0 both columns share one reference grid, and moving the sea floor changes only where that grid is clipped. Every interior interface sits at an identical pseudo-height — and hence an identical pressure — in both columns, so the cross-edge pressure difference is exactly zero above the bottom. Offline replication of the Polaris Init arithmetic confirms this: with a sea-floor gradient of 10–200 m km\(^{-1}\) and a flat coordinate, the RMS HPGA over the interior layers is \(\sim 10^{-15}\) m s\(^{-2}\) (round-off), while the deepest layer valid in both columns carries \(2\times10^{-6}\) to \(8\times10^{-5}\) m s\(^{-2}\).

Two consequences follow:

1. bathymetry_step is a test of the bottom partial cell, not of interior layer tilt, and its whole signal is in the one layer the reference-based analysis discards. Hence the requirement to include that layer.

2. Its error is not a power law in the sea-floor gradient. The signal is set by how the two columns’ maxLevelCell and clipped bottom thicknesses relate, which changes in steps. Fitting an exponent to it would be meaningless, so the fit is configurable and is disabled for that variant.

Algorithm Design: The tilt scan and the exponent

Date last modified: 2026/07/28

Contributors: Xylar Asay-Davis, Claude

Because the exact answer is zero there is nothing to converge to, so the discriminating measurement is not a convergence rate against a reference but a scan of the tilt at fixed resolution. For hydrostatic_consistency the swept parameter is z_tilde_bot_grad; the analysis fits

\[ \log_{10} \lVert \text{HPGA} \rVert_{\text{RMS}} \;=\; q \, \log_{10} (\text{tilt}) + \text{const} \]

over the points below tilt_fit_max, and reports \(q\). The readings are:

measured \(q\)	reading
\(q \approx 1\)	\(\alpha\) is effectively piecewise constant in pressure — no better than the centered scheme
\(q \approx 2\)	consistent with a residual set by the per-cell Taylor expansion point of PGradHighOrder.md §3.3
round-off, independent of tilt	consistent with §3.7 holding as written

An exponent that is neither, or one that drifts across the sweep, is a real possible outcome and would say the mechanism is not one of the three above.

Restricting the fit to tilt_fit_max matters for two reasons: at small tilt the error approaches the round-off floor and flattens the fit, and at large tilt the two columns’ maxLevelCell values begin to differ, which adds a staircase component. This mirrors the existing omega_vs_reference_convergence_fit_max_resolution option.

The scan is meaningful immediately, against the existing Centered scheme and before any new Omega code exists. Offline replication of the Polaris arithmetic gives \(q = 1.00\) for the centered scheme, which is the expected reading for a scheme that holds \(\alpha\) constant within a layer, and which calibrates the diagnostic.

Algorithm Design: Sizing the sweep so that it is severe enough

Date last modified: 2026/07/28

Contributors: Xylar Asay-Davis, Claude

The existing variants use a 500 m water column resolved at 0.5–4 m, which is between one and three orders of magnitude finer than the ~250 m layers in the deep ocean of a realistic global run. Offline replication of the Polaris Init arithmetic (which reproduces the recorded \(7\times10^{-12}\) m s\(^{-2}\) for the existing ztilde_gradient to within 4%) gives, at a fixed coordinate tilt of 5 m km\(^{-1}\):

layer thickness	RMS HPGA (centered)
256 m	\(1.6\times10^{-6}\)
128 m	\(5.0\times10^{-7}\)
64 m	\(9.6\times10^{-8}\)
16 m	\(6.4\times10^{-9}\)

i.e. roughly second order in layer thickness, and only the coarsest rows are near the target severity. The new variants therefore use a deep, coarsely resolved column: a sea floor at \(-3500\) m, a reference pseudo-bottom at \(-4096\) m, and 16, 32 or 64 levels (256, 128 or 64 m). The level counts are multiples of 16, which Omega prefers.

The gap between the sea floor and the reference pseudo-bottom is deliberate. The p-star iteration converges to a pseudo-bottom depth slightly greater than the geometric water-column thickness (because in-situ density exceeds \(\rho_0\) at depth), and the reference grid must extend below it in both columns. The 576 m of head room here permits tilts up to a few hundred m km\(^{-1}\); without it the shallower column’s reference grid runs out, the p-star iteration diverges, and the reported HPGA jumps to \(O(1)\) m s\(^{-2}\) — a failure mode the Init guard below is designed to catch.

A further measured result shapes the sweep: at a fixed tilt gradient, the RMS HPGA is independent of horizontal resolution to five digits (1.5570e-06 at 8, 4, 2 and 1 km). The offset and the cell spacing scale together, so the difference quotient is unchanged. The horizontal sweep is therefore reduced to two resolutions — enough to confirm the same independence holds in Omega, which has not been checked — rather than the seven the existing variants use.

Algorithm Design: The bathymetry guard

Date last modified: 2026/07/28

Contributors: Xylar Asay-Davis, Claude

The swept tilt has to be the geometry actually under test. The analysis therefore checks

\[ \max_i \lvert \texttt{bottomDepth}_i - \texttt{bottomDepthRequested}_i \rvert \;<\; \texttt{resting\_state\_max\_bathy\_error} \]

at every point in the sweep and fails otherwise, with Init writing the requested bathymetry alongside the achieved one.

The thing this guards against is partial-cell snapping. min_pc_fraction forbids bottom cells thinner than a set fraction of a layer, so the achievable pseudo-bottom depths are quantized and a requested bathymetry can fall in a gap between them; _snap_partial_cells then moves the sea floor to the nearest representable depth.

That does not break the resting state. Since run_pstar_init anchors the geometric column at the prescribed sea surface, ssh is exact by construction, the sea surface stays level, and the state remains exactly at rest. What breaks is the meaning of the sweep. This matters most for bathymetry_step, where the swept parameter is the sea floor itself: offline replication shows that at a nominal gradient of 100 m km\(^{-1}\) over 4 km, snapping moves both columns to the same representable depth, so a nominal 400 m step becomes no step, and the measured RMS HPGA falls to \(4\times10^{-16}\) m s\(^{-2}\) — which would read as a spectacular pass rather than a test that had stopped testing anything.

A guard on ssh would be useless here, and was: an earlier version of this design checked \(\lvert \eta_1 - \eta_0 \rvert\), which is identically zero once the column is surface-anchored and therefore can never fire.

Both new variants set partial_cell_type = None, which leaves the requested bathymetry untouched. This does not cost partial-cell coverage: in a p-star column the reference grid is clipped at the pseudo-bottom depth regardless, so the deepest layer is a partial cell either way. The option controls only whether the topography is snapped. The guard remains in place so that anyone who later enables snapping finds out immediately.

Historical note. Before run_pstar_init anchored the column at the prescribed sea surface, the same snapping surfaced as a spurious sea-surface tilt — up to 12 m between adjacent columns in this configuration, producing an apparent “error” of \(3\times10^{-2}\) m s\(^{-2}\) that was entirely real barotropic physics. That is fixed in the framework (see the p-star partial-cell work), and this branch is rebased onto that fix.

Implementation

Implementation: shared metric helpers

Date last modified: 2026/07/28

Contributors: Xylar Asay-Davis, Claude

analysis.py holds several private helpers that the new analysis needs verbatim: locating the internal edge, RMS over finite values, the log-log power law fit, and the netCDF writer. Rather than importing private names across step modules, move them to a new dependency-light leaf module polaris/tasks/ocean/horiz_press_grad/metrics.py with public names (get_internal_edge, rms, power_law_fit, write_metric_dataset, format_value_list, format_value_error_pairs) and have analysis.py import them. This commit is a pure refactor with no behaviour change.

Implementation: including the bottom layer in the reference-based analysis

Date last modified: 2026/07/28

Contributors: Xylar Asay-Davis, Claude

In analysis.py, take one more interface and one more layer:

# was: [: max_level_index + 1]  and  hpga_forward[:max_level_index]
z_tilde_inter_for_ref = z_tilde_inter_edge[: max_level_index + 2]
hpga_ref_diff = hpga_forward[: max_level_index + 1] - ref_layer_mean

This is safe with respect to interface masking: ZTildeInterface is valid for indices 0 .. maxLevelCell in each column and max_level_index is min(maxLevelCell) - 1, so interfaces 0 .. max_level_index + 1 are finite in both columns by construction. The omega_vs_python comparison in the same function already used [: max_level_index + 1], so after this change the two comparisons agree about which layers count.

Measured effect. Running the real Init arithmetic and the real ReferenceColumn against a synthetic two-cell mesh (only mesh generation bypassed), with the Polaris HPGA standing in for Omega’s tendency:

variant	fitted slope	RMS at 0.5 km
salinity_gradient	1.787 → 1.799	2.600e-07 → 2.601e-07
surface_pressure_gradient	2.014 → 2.012	1.137e-11 → 1.136e-11
temperature_gradient	1.625 → 1.625	2.770e-07 → 2.770e-07
ztilde_gradient	2.003 → 1.119	6.779e-12 → 4.954e-11

The first three are unaffected. ztilde_gradient is the only variant with a tilted coordinate, so its two columns clip the bottom layer at different thicknesses; that cell carries -1.55e-09 m s\(^{-2}\) at 0.5 km against an interior RMS of 6.8e-12, converges at roughly first order, and pulls the fitted slope down. That is a real first-order error in the bottom partial cell which the exclusion was hiding — the same behaviour hydrostatic_consistency measures directly, and what the higher-order scheme exists to fix.

ztilde_gradient.cfg therefore takes its own convergence band, provisionally [0.9, 1.4], marked in the file as requiring re-tuning from the first real Omega run. The shared omega_vs_reference_high_res_rms_threshold of 1e-6 needs no change.

Implementation: the reference-free analysis step

Date last modified: 2026/07/28

Contributors: Xylar Asay-Davis, Claude

Add polaris/tasks/ocean/horiz_press_grad/resting_analysis.py defining RestingAnalysis, a sibling of Analysis rather than a mode of it, so the existing step is untouched. Per sweep point it:

1. locates the internal edge with get_internal_edge;

2. takes maxLevelCell for the two bounding cells and forms the valid range 0 .. min(maxLevelCell) - 1 inclusive of the deepest common layer;

3. reads Omega’s NormalVelocityTend at that edge and takes its RMS over the valid layers — the truth is zero, so this is the error;

4. does the same for the Polaris-side HPGA from init.nc;

5. RMS-differences the two and applies omega_vs_polaris_rms_threshold, which is retained: it confirms Omega implements the discretization Polaris thinks it does, independently of the resting-state property;

6. applies the bathymetry guard, comparing bottomDepth from vert_coord.nc against bottomDepthRequested from init.nc.

It then groups the sweep by resolution pair, fits power_law_fit over the points with tilt \(\le\) tilt_fit_max when tilt_fit is true, writes resting_state.nc and resting_state.png (log-log, one series per resolution pair), logs the exponents, and applies the gates:

• sensitivity gate — the largest RMS over the whole sweep must be at least resting_state_sensitivity_min_rms; this is a hard prerequisite, and failing it means the sweep must be redesigned before anything is concluded from it;

• consistency gate — RMS below resting_state_max_rms, deliberately disabled (none) until the rollout stage measures a value to set it from;

• Omega-vs-Polaris and resting-state checks as above.

Implementation: the tilt axis and the resting-state task

Date last modified: 2026/07/28

Contributors: Xylar Asay-Davis, Claude

Init and Forward gain an optional subdir_suffix, and Init gains optional tilt_option / tilt. When a tilt is given, Init.run sets that option in its own config before building the columns, exactly as it already sets vertical_grid:vert_levels. Init.run also gains a guard that raises when a column’s reference pseudo-bottom does not extend below its converged pseudo-bottom depth — the divergence mode described above — with an error naming the offending column and the head room available.

Task wiring goes in a new module resting_state_task.py defining HorizPressGradRestingStateTask, rather than branching HorizPressGradTask. The existing task pairs horiz_resolutions[i] with vert_resolutions[i] and keys its step dictionaries by horizontal resolution alone, so it cannot express the repeated horizontal resolutions the new sweep needs; the new task keys by the (horiz_res, vert_res, tilt) triple. Init and Forward are reused unchanged apart from the optional arguments above.

This is a deliberate departure from the working plan, which leaned toward adding the tilt axis inside HorizPressGradTask._setup_steps. Doing so would have required renaming the existing variants’ steps to disambiguate the vertical resolution, changing their work directories for no benefit to them.

Implementation: configuration

Date last modified: 2026/07/28

Contributors: Xylar Asay-Davis, Claude

horiz_press_grad.cfg gains the resting-state options with inert defaults, so the four existing variants are unaffected. Which task class a variant uses is decided by name in __init__.py, so no separate “this is a resting state” flag is needed. The new variant .cfg files set the deep column, the coarse vertical resolutions, and the sweep:

# hydrostatic_consistency.cfg (abridged)
geom_z_bot_mid = -3500.0
z_tilde_bot_mid = -4096.0
horiz_resolutions = [4.0, 4.0, 4.0, 1.0]
vert_resolutions = [256.0, 128.0, 64.0, 256.0]
tilt_option = z_tilde_bot_grad
tilt_values = [0.05, 0.1, 0.5, 1.0, 5.0, 10.0, 50.0]
tilt_fit = True
tilt_fit_max = 10.0

# bathymetry_step.cfg (abridged)
geom_z_bot_mid = -3500.0
z_tilde_bot_mid = -4096.0
horiz_resolutions = [4.0, 4.0]
vert_resolutions = [256.0, 128.0]
tilt_option = geom_z_bot_grad
tilt_values = [1.0, 5.0, 10.0, 25.0, 50.0, 100.0, 200.0]
tilt_fit = False

Both set partial_cell_type = None in [vertical_grid], for the reason given under the bathymetry guard above.

forward.yaml was unchanged for the calibration stage: every variant ran the existing Centered scheme. The rollout gives it a PressureGrad block, filled in per step — see below.

Implementation: the scheme axis

Date last modified: 2026/08/02

Contributors: Xylar Asay-Davis, Claude

Forward gains a scheme argument, and forward.yaml a PressureGrad block filled in through add_yaml_file’s template_replacements:

PressureGrad:
   PressureGradType: {{ pressure_grad_type }}
   QuadraturePoints: {{ quadrature_points }}

HorzOrder and VerticalReconstruction are left at their Omega defaults of 2 and linear, which are the only Phase 1 values; the Phase 2 values exist as keys but are rejected with an error, so a configuration written for Phase 2 cannot quietly run as Phase 1. An unrecognized PressureGradType is likewise fatal rather than falling back to Centered, so a typo aborts the run instead of producing centered answers that read as a pass.

Both task classes key their step dictionaries by (resolution, scheme) — or (horiz_res, vert_res, tilt, scheme) for the resting-state sweep — and build one Forward per scheme over one shared Init. Init is scheme-independent: it writes both HPGA and HPGAFiniteVolume from the same state, so duplicating it per scheme would cost run time and, worse, would leave the two schemes’ comparison open to the objection that they were run on different initial conditions.

This is the departure from the calibration stage’s structure noted under the amended requirement above: it changes the existing variants’ step directories. The alternative — a task per (variant, scheme) pair — preserves them, but runs Init twice per state and leaves no analysis step that sees both schemes, which is the comparison the rollout exists to make.

quadrature_points is a single config option that drives both this template and the Python kernel’s Gauss order. It is deliberately not two options: the two sides using different quadrature would make them different algorithms, which is what the cross-implementation requirement above forbids. Sweeping it for accuracy remains available and moves both sides together.

Testing

Testing and Validation: the sensitivity gate is the prerequisite

Date last modified: 2026/07/28

Contributors: Xylar Asay-Davis, Claude

Run both variants with Centered and confirm the sensitivity gate passes. Offline replication of the Polaris arithmetic predicts the centered RMS HPGA listed below; the first job of the real run is to confirm Omega agrees, via the retained omega_vs_polaris_rms_threshold check.

variant	configuration	predicted centered RMS
hydrostatic_consistency	4 km, 256 m, tilt 50 m km\(^{-1}\)	\(2.1\times10^{-5}\)
hydrostatic_consistency	4 km, 64 m, tilt 50 m km\(^{-1}\)	\(9.8\times10^{-7}\)
bathymetry_step	4 km, 256 m, 100 m km\(^{-1}\)	\(2.3\times10^{-5}\)
bathymetry_step	4 km, 128 m, 25 m km\(^{-1}\)	\(3.8\times10^{-6}\)

If the gate fails, the sweep is not exercising the failure mode and must be redesigned before the variants are used to judge any scheme. This is the single most consequential check in the calibration stage.

Testing and Validation: the tilt exponent

Date last modified: 2026/07/28

Contributors: Xylar Asay-Davis, Claude

Confirm that hydrostatic_consistency returns \(q \approx 1\) for Centered at every resolution pair, and that the value is stable when tilt_fit_max is varied. Offline replication gives \(q = 0.95\)–\(1.00\), the droop below 1 coming from the largest tilts where maxLevelCell starts to differ between columns — which is what tilt_fit_max is there to exclude. A centered exponent far from 1 should be understood before the scan is used on the new scheme.

Confirm that bathymetry_step produces the expected staircase and that its interior layers sit at round-off, isolating the signal to the bottom cell.

Testing and Validation: resting-state and consistency checks

Date last modified: 2026/07/28

Contributors: Xylar Asay-Davis, Claude

• Bathymetry guard. Confirm it passes for both variants as configured, and confirm it fires for a deliberately broken configuration — set partial_cell_type = partial on bathymetry_step with a sea-floor gradient of 32 m km\(^{-1}\), which offline replication shows moves the sea floor 12.33 m from the request.

• Zero-tilt sanity. With the swept parameter set to zero the RMS HPGA must be exactly zero (offline replication gives 0.0, not merely round-off), since the two columns become identical.

• Omega-vs-Polaris. The retained check at omega_vs_polaris_rms_threshold must pass at every sweep point. Its present value of \(10^{-10}\) was tuned for errors of order \(10^{-7}\); it may need re-tuning for a signal of order \(10^{-5}\) and should be set from the observed differences.

• Regression. Run the four existing variants and confirm their step names and work directories are unchanged. Their results change, because the bottom layer is now included: three are unaffected in practice, while ztilde_gradient’s convergence slope drops from ~2.0 to ~1.1 and its provisional band must be re-tuned from the observed value.

Testing and Validation: the rollout stage

Date last modified: 2026/08/02

Contributors: Xylar Asay-Davis, Claude

Two corrections to how this section read before Omega’s FiniteVolume landed, both from omega_finite_volume_handoff.md. The keys are HorzOrder: 2 and VerticalReconstruction: linear, not ReconstructionOrder and constant — constant truncates the equation-of-state expansion to \(\alpha_0\) and would run the centered scheme under a new name. And there is no “reduced configuration” or centered limit of the new scheme: the two are separate implementations and no setting reduces one to the other. Phase 1 offers exactly two configurations, Centered and FiniteVolume with its defaults.

The Python counterpart also already exists — it was written during the calibration stage, from the design rather than from the C++ — so what this stage tests is not whether it can be built but whether it and Omega agree.

Reconciliation, before anything is run. Work through every discretization choice named in the cross-implementation requirement above, record its resolution and provenance, and give each confirmed mismatch its own fix and its own test. Two are already confirmed: the Python kernel defaults to a 4-point Gauss rule where Omega defaults to 2, and it anchors its column scan at the sea surface where Omega anchors at the sea floor. Neither shows up on the exact set, where the integrand is zero at every quadrature point; both bite off it, which is where the §5.1 gates live.

Gates, set from measurement and not before. Three points, of which the first two would otherwise have been got wrong:

• bands are per variant, not per scheme alone. The measured Centered rates are 1.6, 1.8, 2.0 and 1.1 across the four gradient variants, so a single “≈ 2 for centered” band fails three of them;

• the accuracy gate is a ratio at the coarse end, set from the measured advantage — 6.5× at 256 m layers, 2.4× at 64 m, ~400× on stepped bathymetry;

• do not gate on relative convergence order. The two schemes converge at similar rates on smooth profiles; the advantage is a constant factor that narrows under refinement, so a gate requiring the new scheme’s slope to beat the centered one would fail where nothing is wrong.

resting_state_max_rms is set from the measured FiniteVolume level at a stated tilt, as an explicit multiple. The improvement gate — FiniteVolume RMS below Centered RMS at every sweep point — needs no prior measurement and ships with the analysis change.

Report absolutes. \(\max\lvert\mathcal{S}\rvert\) behaved unexpectedly under vertical refinement on the linear variant during the calibration stage, so absolute residuals and tendencies in m s\(^{-2}\) are reported alongside any ratio, not in place of them.

The variant to expect surprises from is bathymetry_step. Every Omega-side measurement of the new scheme so far used equal layer counts in both columns, so this is where the rule for evaluating a column’s reconstruction below its own floor is first genuinely exercised.

Finally, the result — whichever way it falls — is written back into PGradHighOrder.md as a measured property with its configuration and tilt range attached. That is an Omega-side commit in a separate repository and is not folded into the Polaris work.