Reworked HPG Reference Solution

Creation date: 2026/06/14

Contributors: Xylar Asay-Davis, Claude

Background: the horizontal pressure gradient and its reference solution

Omega is a non-Boussinesq, hydrostatic ocean model whose prognostic vertical coordinate is pseudo-height, \(\tilde{z} = -p / (\rho_0 g)\), where \(p\) is sea gauge pressure, \(\rho_0\) is a reference seawater density (RhoSw), and \(g\) is gravitational acceleration. The ocean/horiz_press_grad task family in Polaris exists to verify the convergence of Omega’s horizontal pressure-gradient (HPG) operator against a high-fidelity reference solution, as called for in the Omega pressure gradient design document (PGrad.md, “Test: Spatial convergence to exact solution”).

The continuous horizontal pressure-gradient force (HPGF) per unit mass that Omega’s momentum tendency must reproduce is, from the Omega V1 governing equations (OmegaV1GoverningEqns.md),

\[ \mathbf{a}_\text{PGF} = -\left(\alpha\,\nabla p + \nabla \Phi\right), \]

where \(\alpha = 1/\rho\) is specific volume, \(\Phi = g z\) is the geopotential (neglecting tidal and self-attraction-and-loading terms), and \(\nabla\) is the along-layer horizontal gradient. Using the hydrostatic relation, Omega rewrites this through the Montgomery potential \(M = \alpha p + g z\) as

\[ \mathbf{a}_\text{PGF} = -\nabla M + p\,\nabla \alpha, \]

and discretizes it with a centered scheme in which \(\alpha\) is held constant within each layer, the Montgomery gradient is evaluated at layer interfaces and averaged to the layer midpoint, and the along-layer gradient is a two-cell finite difference across the edge.

The two-column task and its comparisons

Each horiz_press_grad task variant (salinity_gradient, temperature_gradient, ztilde_gradient) builds a two-column mesh at a set of paired horizontal and vertical resolutions. Per resolution:

An Init step (a PStarInitStep, see the p-star initialization design) builds the p-star coordinate and the initial conservative temperature (CT) and absolute salinity (SA) from piecewise pseudo-height profiles, and computes a Polaris HPGA from a two-cell finite difference of the Montgomery potential.
A Forward step runs Omega for a single time step and writes Omega’s NormalVelocityTend, which is just the HPGA.
An Analysis step performs two comparisons:
1. omega_vs_reference — Omega’s HPGA against the high-fidelity reference solution, fit to a power law to measure the convergence rate.
2. omega_vs_python — Omega’s HPGA against the Polaris/Python HPGA from Init, a code-correctness check that the two implementations of the same discrete scheme agree to a tight tolerance.

Because each forward run is a single step with no coordinate motion, the forward layer pseudo-heights and geometric heights are identical to those in Init; the omega_vs_python comparison therefore needs no special handling.

The current reference solution and why it must be reworked

The current Reference step (reference.py) builds five columns at \(x = \Delta x_\text{ref}\,[-1.5,-0.5,0,0.5,1.5]\), integrates the hydrostatic relation upward from each seafloor to obtain geometric height \(z(\tilde z)\), forms \(M\), and takes a 4th-order finite difference across columns to obtain \(\partial M/\partial x\) and \(\partial \alpha/\partial x\) at the center column, giving HPGA \(= -\partial M/\partial x + p\,\partial\alpha/\partial x\). It writes the result on a vertical grid deliberately refined (via a greatest-common- divisor construction) so that every test grid’s layer midpoints land exactly on reference grid points. Analysis then samples the reference by exact z̃ matching with no interpolation.

This design has two structural weaknesses:

Grid-alignment fragility. Exact-match sampling only works when the reference and test pseudo-heights coincide. As soon as a nonzero SurfacePressure is applied, the p-star coordinate stretches and squashes the column in proportion to the pressure thickness, and partial bottom cells produce interfaces that no longer align across vertical resolutions. The reference must instead be queryable at arbitrary pseudo-height with no loss of accuracy.
Breakdown near bathymetry. The cross-column finite-difference stencil requires every neighboring column to be valid at a given \(\tilde z\). Near the seafloor, collocation points in some columns fall below their own bathymetry and the stencil cannot be formed, so those layers are masked out. The accuracy is also only as good as the stencil and the matched grid, which is insufficient for assessing the higher-than-second-order HPG operators planned for Omega.

Summary

This design replaces the multi-column, finite-difference reference solution with a single-column, analytic reference computed at the edge between the two columns (\(x = 0\)). The new capability:

Evaluates the exact continuous HPGF directly via a chain-rule / Leibniz expansion of the along-pseudo-height gradient, using TEOS-10 specific-volume derivatives. Because the continuous HPGF is coordinate-invariant, this single central-column expression is exact in the horizontal — there is no finite-difference truncation and no dependence on neighboring columns, so the near-bathymetry breakdown disappears and the reference is fidelity-limited only by a smooth one-dimensional quadrature.
Is a callable function of \(\tilde z\), so it can be sampled at exactly the pseudo-heights a test produces, eliminating the grid-alignment requirement.
Supports fully general horizontal variation of every prescribed input.
Is consumed directly by Analysis through a shared evaluator module; the separate Reference step is removed.

The accuracy assessment compares Omega’s along-layer HPGA to the reference using a layer-averaged target: for each model layer, the reference is averaged over that layer’s pseudo-height extent at the edge. This compares like with like (a layer mean against a layer mean) so that legitimate layer tilt and finite layer thickness are accounted for in the target rather than charged as discretization error.

The primary challenges are (a) formulating the continuous reference correctly so that the along-layer computation is judged on its own merits, (b) differentiating the prescribed profiles with respect to the horizontal coordinate at fixed pseudo-height — including the case where the profile node pseudo-heights themselves tilt — and (c) integrating the resulting expression to high accuracy in the vertical. The design is successful if the reference can be evaluated at any pseudo-height without loss of accuracy, remains valid throughout the water column including immediately above the bathymetry, is accurate enough to serve a beyond-second-order HPG operator, supports general horizontal input variation, and reproduces (with re-tuned thresholds) the convergence behavior the existing regression tests check.

Requirements

Requirement: The reference HPGA can be evaluated at an arbitrary pseudo-height without loss of accuracy

Date last modified: 2026/06/14

Contributors: Xylar Asay-Davis, Claude

The reference solution shall be evaluable at any pseudo-height within the water column produced by a test, to high accuracy, without requiring that the reference and test vertical grids coincide. This must hold when a nonzero sea-surface pressure stretches or squashes the p-star column and when partial bottom cells misalign layer interfaces across resolutions.

Requirement: The reference is valid throughout the water column, including immediately above the bathymetry

Date last modified: 2026/06/14

Contributors: Xylar Asay-Davis, Claude

The reference shall remain valid and accurate for every pseudo-height from the free surface down to the seafloor of the column at the edge between the two test columns. It shall not break down or require masking in the layers nearest the bathymetry.

Requirement: The reference is accurate enough to assess HPG operators beyond second order

Date last modified: 2026/06/14

Contributors: Xylar Asay-Davis, Claude

The reference shall be of sufficient fidelity that, when an HPG operator more accurate than the current second-order centered scheme is introduced in Omega, the reference error does not limit the measured convergence. In particular, the reference shall not introduce a fixed-order horizontal truncation error of its own.

Requirement: The reference supports general horizontal variation of the prescribed inputs

Date last modified: 2026/06/14

Contributors: Xylar Asay-Davis, Claude

The reference shall correctly represent horizontal variation of every prescribed input that the task permits: the geometric sea-surface height, the geometric seafloor depth, the seafloor pseudo-height, and the temperature and salinity profile node values and node pseudo-heights. No input variation supported by the task configuration may be assumed to be horizontally uniform.

Requirement: The accuracy assessment evaluates the along-layer scheme on its own merits

Date last modified: 2026/06/14

Contributors: Xylar Asay-Davis, Claude

The comparison between the model’s along-layer HPGA and the reference shall measure only the genuine discretization error of the along-layer scheme. It shall not attribute legitimate, physically correct effects of tilted layers or of finite layer thickness to discretization error. The reference shall represent the exact continuous pressure-gradient force that the along-layer scheme converges to as both horizontal and vertical resolution are refined.

Requirement: The capability integrates with the existing two-column task without a separate precompute step

Date last modified: 2026/06/14

Contributors: Xylar Asay-Davis, Claude

The reference shall be available to the Analysis step on demand at the pseudo-heights each test produces, without a separate precomputed reference artifact and without changing the structure of the Init/Forward/Analysis workflow or the two comparisons (omega_vs_reference and omega_vs_python).

Algorithm Design

Algorithm Design: The continuous, coordinate-invariant reference and the cancellation of layer-slant terms

Date last modified: 2026/06/14

Contributors: Xylar Asay-Davis, Claude

The reference is the exact continuous pressure-gradient acceleration that Omega’s discrete operator must converge to. The key fact that makes a single central-column evaluation correct is that this acceleration is independent of the vertical coordinate used to evaluate the horizontal gradients.

Write the horizontal gradients along a general layer label \(\sigma\) (constant along a model coordinate surface). Pseudo-height is, by definition, a rescaled pressure, \(p = -\rho_0 g\,\tilde z\), so along constant \(\tilde z\) the pressure is constant. Using the hydrostatic relation \(\partial z/\partial \tilde z = \rho_0 \alpha\), the along-\(\sigma\) gradients of pressure and geometric height are

\[\begin{split} \nabla_\sigma p &= -\rho_0 g\,\left.\partial_x \tilde z\right|_\sigma, \\ \nabla_\sigma z &= \left.\partial_x z\right|_{\tilde z} + \rho_0 \alpha\,\left.\partial_x \tilde z\right|_\sigma . \end{split}\]

Substituting into \(\mathbf{a}_\text{PGF} = -(\alpha\,\nabla_\sigma p + g\,\nabla_\sigma z)\), the two terms proportional to the coordinate slope \(\left.\partial_x\tilde z\right|_\sigma\) carry the same \(\alpha\) at the same point and cancel exactly:

\[ \mathbf{a}_\text{PGF} = -g\,\left.\frac{\partial z}{\partial x}\right|_{\tilde z}. \]

This is the standard statement that the pressure-gradient force does not depend on the coordinate used to compute it. Consequences for the design:

The reference is \(-g\,\partial_x z|_{\tilde z}\), evaluated at the edge (\(x=0\)). It can be computed entirely at constant \(\tilde z\), which is the convenient frame, without approximation, even when the model layers slant (as in ztilde_gradient) or bend near the bathymetry (partial cells).
The cancellation holds only for the total acceleration. Omega’s discrete pieces \(\nabla M\) and \(p\,\nabla\alpha\) are individually coordinate-dependent and large; the design never compares those pieces to the reference, only the total. The residual between Omega’s along-layer total and the continuous total is therefore the genuine discretization error of the along-layer scheme — including the classic slant-induced pressure-gradient error that the future high-order operator is intended to reduce.

Algorithm Design: Single-column chain-rule / Leibniz expression for the reference

Date last modified: 2026/06/14

Contributors: Xylar Asay-Davis, Claude

Geometric height as a function of pseudo-height in a single column is obtained by integrating the hydrostatic relation from the surface:

\[ z(\tilde z; x) = \eta(x) + \rho_0 \int_{\tilde z_s(x)}^{\tilde z} \alpha\bigl(S_A(\tilde z', x),\,\Theta(\tilde z', x),\,p(\tilde z')\bigr) \, d\tilde z' , \]

where \(\eta(x)\) is the geometric sea-surface height, \(\tilde z_s(x)\) is the surface pseudo-height (\(\tilde z_s = -p_s/(\rho_0 g)\), zero only when the surface pressure \(p_s\) is zero), and \(p(\tilde z') = -\rho_0 g\,\tilde z'\) depends only on \(\tilde z'\) (not on \(x\)). The integral runs downward, so for \(\tilde z < \tilde z_s\) it is negative and \(z < \eta\), as required. The reference is anchored at the surface rather than the seafloor because that is the boundary the model honours: Init/Omega build the p-star column from the prescribed sea-surface height and surface pressure, so the HPGA vanishes at the surface (for zero surface slope) and grows downward. Anchoring elsewhere — e.g. at the configured seafloor pseudo-height, which is a coordinate buffer below the deepest valid layer rather than a honoured geometric boundary — would reconstruct a different \(z(\tilde z)\) and miss the model by a constant offset.

Differentiating with respect to \(x\) at fixed \(\tilde z\), the lower integration limit depends on \(x\), so the Leibniz rule contributes a boundary term:

\[ \left.\frac{\partial z}{\partial x}\right|_{\tilde z} = \eta'(x) - \rho_0\,\alpha(\tilde z_s)\,\tilde z_s'(x) + \rho_0 \int_{\tilde z_s}^{\tilde z} \left.\frac{\partial \alpha}{\partial x}\right|_{\tilde z'} d\tilde z' . \]

Because \(p\) is fixed along constant \(\tilde z'\), the integrand involves only the horizontal variation of the thermohaline state through the TEOS-10 specific volume:

\[ \left.\frac{\partial \alpha}{\partial x}\right|_{\tilde z'} = \alpha_{S_A}\,\left.\frac{\partial S_A}{\partial x}\right|_{\tilde z'} + \alpha_{\Theta}\,\left.\frac{\partial \Theta}{\partial x}\right|_{\tilde z'}, \]

where \(\alpha_{S_A} = \partial \alpha / \partial S_A\) and \(\alpha_{\Theta} = \partial \alpha / \partial \Theta\) are the first derivatives of specific volume from the equation of state, evaluated at the central column’s \((S_A, \Theta, p)\). The reference acceleration is then

\[ a(\tilde z) = -g\left.\frac{\partial z}{\partial x}\right|_{\tilde z} = -g\left[\,\eta' - \rho_0\,\alpha(\tilde z_s)\,\tilde z_s' + \rho_0 \int_{\tilde z_s}^{\tilde z} \left(\alpha_{S_A}\,\partial_x S_A + \alpha_{\Theta}\,\partial_x \Theta\right) d\tilde z'\right]. \]

Notable properties:

No geometric-height integration is required for the HPGA. Everything is expressed through the prescribed profiles \(S_A(\tilde z'), \Theta(\tilde z')\) and the analytic pressure \(p(\tilde z') = -\rho_0 g\,\tilde z'\) at the central column, plus the input slopes \(\eta'\) and \(\tilde z_s'\). The only quadrature is the smooth one-dimensional integral of \(\partial_x\alpha\). (A geometric-height integration may still be performed for optional diagnostics.)
Single column, valid to the seafloor. The expression is evaluated only at \(x = 0\) and is well defined for every \(\tilde z \in [\tilde z_b(0), \tilde z_s(0)]\). Integrating from the surface loses none of this validity — it reaches every pseudo-height down to and immediately above the bathymetry — and it never references neighboring columns, so it cannot break down above the bathymetry.
Exact in the horizontal. The horizontal gradient is analytic; there is no \(\Delta x\) truncation, satisfying the beyond-second-order fidelity requirement.

The boundary term uses \(\tilde z_s(0) = \) the surface pseudo-height (the shallowest z_tilde node) and \(\alpha(\tilde z_s)\) evaluated from the profiles at \(\tilde z_s(0)\). The slopes \(\eta'\) (from geom_ssh_grad) and \(\tilde z_s'\) (from the surface z_tilde_grad node) are read from the configuration gradients (converted from per-km to per-m and projected onto the edge normal). They are kept fully general: a nonzero sea-surface height or surface pressure makes \(\eta'\) and/or \(\tilde z_s'\) (and \(\tilde z_s\) itself) nonzero, and the formula already accounts for them. A regression test that exercises a nonzero surface pressure is deferred to a follow-up PR; the present tasks all use zero surface slope and pressure.

Algorithm Design: Horizontal derivatives of the prescribed profiles at fixed pseudo-height

Date last modified: 2026/06/14

Contributors: Xylar Asay-Davis, Claude

The factors \(\partial_x S_A|_{\tilde z'}\) and \(\partial_x \Theta|_{\tilde z'}\) are the horizontal gradients of the prescribed profiles at fixed pseudo-height. The profiles are defined by node values and node pseudo-heights, each of which varies linearly in \(x\) (configuration *_mid plus *_grad\(\cdot x\)), and a monotone PCHIP interpolant in pseudo-height (the same interpolant Init uses, so the reference and the test describe the same physical state).

The profile value at fixed \(\tilde z'\) depends on \(x\) through both the node values and, in the general case (z_tilde_grad \(\neq 0\)), the node pseudo-heights. Because PCHIP is nonlinear in its node data, the horizontal derivative is obtained by differencing the analytic input profiles: the interpolant is constructed at \(x = \pm\varepsilon\) from the exactly known linear node functions and evaluated at \(\tilde z'\), and a centered difference is taken,

\[ \left.\frac{\partial S_A}{\partial x}\right|_{\tilde z'} \approx \frac{S_A(\tilde z'; +\varepsilon) - S_A(\tilde z'; -\varepsilon)} {2\varepsilon}, \]

and similarly for \(\Theta\). Since the nodes are exactly linear in \(x\) and the interpolant is smooth in \(x\), this is accurate to round-off for a suitably chosen \(\varepsilon\) (optionally refined by Richardson extrapolation). This single construction handles all task variants uniformly, including tilted node pseudo-heights, satisfying the general-input requirement.

The specific-volume derivatives \(\alpha_{S_A}\) and \(\alpha_{\Theta}\) are obtained from the equation of state at the central column’s \((S_A(\tilde z'), \Theta(\tilde z'), p(\tilde z'))\) — concretely, from gsw.specvol_first_derivatives for TEOS-10.

Algorithm Design: Building the reference profile and the layer-averaged comparison

Date last modified: 2026/06/14

Contributors: Xylar Asay-Davis, Claude

The reference is needed as a layer mean over each model layer. Define the cumulative integral

\[ I(\tilde z) = \int_{\tilde z_s}^{\tilde z} \left(\alpha_{S_A}\,\partial_x S_A + \alpha_{\Theta}\,\partial_x \Theta\right) d\tilde z', \qquad a(\tilde z) = -g\left[\,\eta' - \rho_0\,\alpha(\tilde z_s)\,\tilde z_s' + \rho_0\,I(\tilde z)\right]. \]

The integrand is smooth, so \(I(\tilde z)\) — and therefore \(a(\tilde z)\) — is evaluated with a high-order quadrature (the existing Gauss–Legendre / adaptive machinery; see Implementation). \(a(\tilde z)\) is thus available as an accurate callable of pseudo-height.

For the comparison, each model layer \(k\) spans, at the edge, the pseudo-height interval \([\tilde z^{\text{bot}}_{k}, \tilde z^{\text{top}}_{k}]\), where the edge interfaces are the average of the two cells’ interface pseudo-heights (the edge lies midway between the columns; for linearly varying inputs this average equals the central-column value). The reference target for layer \(k\) is the layer mean

\[ \overline{a}_k = \frac{1}{\tilde z^{\text{top}}_k - \tilde z^{\text{bot}}_k} \int_{\tilde z^{\text{bot}}_k}^{\tilde z^{\text{top}}_k} a(\tilde z)\, d\tilde z , \]

evaluated by high-order quadrature within the layer. This compares the model’s layer-mean acceleration to the true layer-mean acceleration, so finite layer thickness and layer tilt are folded into the target rather than counted as error (satisfying the “on its own merits” requirement). As both resolutions refine, \(\overline{a}_k \to a(\tilde z_k)\) and the residual is the scheme’s genuine discretization error.

The model layer interfaces are taken from the test output (identical between Init and Forward), so the reference is sampled at exactly the test’s pseudo-heights and the partial-cell bottom geometry is honored automatically. As in the current analysis, the single deepest valid model layer (which abuts the bathymetry) is excluded from the error metric.

The sign/orientation of \(a\) must match the edge-normal convention used by the model output: the horizontal coordinate \(x\) increases from the first to the second cell on the internal edge. Analysis determines that orientation from the mesh cell positions and applies it consistently to the input slopes and the reference acceleration.

Algorithm Design: Removal of the separate reference step

Date last modified: 2026/06/14

Contributors: Xylar Asay-Davis, Claude

Because the reference is an inexpensive callable evaluated at the test’s pseudo-heights, it does not need to be precomputed by a dedicated step. The Reference step is removed. Analysis constructs the reference evaluator from the configuration (which fully determines the central column) and evaluates the layer-averaged target per resolution. The greatest-common-divisor reference grid, the five-column construction, the cross-column finite-difference stencil, and the exact-match sampling are all retired.

Implementation

Implementation: shared reference evaluator module

Date last modified: 2026/06/14

Contributors: Xylar Asay-Davis, Claude

Replace the Reference step in polaris/tasks/ocean/horiz_press_grad/reference.py with a Reference evaluator module in the same package (the file name may remain reference.py; it no longer defines a Step subclass). The evaluator encapsulates the central column and exposes the reference as a function of pseudo-height.

class ReferenceColumn:
    """
    High-fidelity continuous HPGA reference at the edge (x = 0) between the
    two columns, evaluated analytically via the chain-rule / Leibniz
    expansion of the along-pseudo-height gradient.
    """

    def __init__(self, config, x_sign=1.0):
        # x_sign carries the edge-normal orientation (+1 if x increases from
        # cell0 to cell1, -1 otherwise) so the reference matches the model
        # output convention.
        ...
        # Cache, at x = 0:
        #   - PCHIP interpolants for SA(z_tilde), CT(z_tilde)
        #   - the surface pseudo-height z_tilde_s (shallowest z_tilde node),
        #     the sea-surface-height slope eta' and the surface pseudo-height
        #     slope z_tilde_s' (per-metre, * x_sign); kept general so nonzero
        #     SSH / surface pressure need no further change
        #   - alpha(z_tilde_s) for the boundary term

    def specvol(self, z_tilde):
        """alpha = EOS(SA(z_tilde), CT(z_tilde), p = -RhoSw*g*z_tilde)."""

    def dalpha_dx(self, z_tilde):
        """
        alpha_SA * dSA/dx + alpha_CT * dCT/dx at fixed z_tilde, with alpha_*
        from gsw.specvol_first_derivatives and dSA/dx, dCT/dx from analytic-input
        differencing of the profiles (handles tilted nodes).
        """

    def hpga(self, z_tilde):
        """
        a(z_tilde) = -g [ eta' - RhoSw*alpha(z_tilde_s)*z_tilde_s'
                          + RhoSw * I(z_tilde) ],
        with I(z_tilde) the cumulative integral of dalpha_dx from the surface
        z_tilde_s. Vectorized over a 1-D array of target pseudo-heights.
        """

    def layer_mean_hpga(self, z_tilde_interfaces):
        """
        Given edge layer interface pseudo-heights (length nLayers+1), return
        the per-layer mean of a(z_tilde) over each [z_top, z_bot] interval by
        high-order quadrature.
        """

Reuse existing utilities rather than reimplementing them:

get_array_from_mid_grad(config, name, x) and get_pchip_interpolator(...) from column.py for node arrays and profile interpolants. Pass an array of the single point x = [0.0] (plus x = ±epsilon for the differencing).
compute_specvol from polaris.ocean.eos (TEOS-10 path) for \(\alpha\), and gsw.specvol_first_derivatives(SA, CT, p_dbar) for \(\alpha_{S_A}, \alpha_{\Theta}\) (note gsw takes pressure in dbar = Pa \(\times 10^{-4}\)).
The fixed-step Gauss–Legendre and adaptive-Simpson quadrature helpers already present in reference.py (_gauss_composite, _fixed_quadrature, _adaptive_simpson) for both \(I(\tilde z)\) and the layer-mean integral. The quadrature method remains configurable via reference_quadrature_method (default gauss4).
Gravity and RhoSw from polaris.ocean.vertical.ztilde.

Recommended evaluation strategy for hpga/layer_mean_hpga: build \(a(\tilde z)\) on a fine internal pseudo-height grid spanning \([\tilde z_b(0), 0]\) by cumulative high-order quadrature of the integrand, then integrate \(a\) over each model layer with the same quadrature using a few sub-samples per layer (the sub-sampling fineness can reuse reference_quadrature_method’s panel count). Choose the internal grid fine enough that its contribution to the error is far below the HPGA errors being measured; alternatively evaluate \(I\) on demand at the layer quadrature points to avoid any interpolation. Either approach must keep the reference fidelity well beyond the model’s.

epsilon for the profile differencing should be a small fraction of the column extent (e.g. relative to reference_horiz_res if retained for this purpose, or a fixed small value such as \(10^{-3}\) km); document the choice and optionally apply Richardson extrapolation.

Implementation: Analysis consumes the reference directly

Date last modified: 2026/06/14

Contributors: Xylar Asay-Davis, Claude

In analysis.py:

Remove the reference dependency, the reference_solution.nc input, and _sample_reference_without_interpolation.
In setup, drop the reference-step input wiring; keep the per-resolution init, forward, culled_mesh, and vert_coord inputs.
In run, construct a single ReferenceColumn from self.config (its orientation x_sign is derived once from the mesh, see below). Then for each resolution:
1. Identify the internal edge and its two cells with the existing _get_internal_edge.
2. Determine x_sign from the cells’ xCell (sign of xCell[cell1] - xCell[cell0]), so the reference gradient matches the model edge-normal convention. The current Python HPGA in Init (_compute_montgomery_and_hpga) uses x = horiz_res*[-0.5, 0.5] for cells [0, 1]; the new orientation logic must reproduce the same sign the existing omega_vs_python comparison relies on.
3. Read the edge layer-interface pseudo-heights as the average of the two cells’ ZTildeInterface from the forward output (equivalently Init).
4. Compute ref_layer_mean = reference.layer_mean_hpga(edge_interfaces).
5. Restrict to valid layers via maxLevelCell (shallowest of the two cells, as today) and drop the deepest valid layer (abuts bathymetry).
6. RMS-difference Omega’s NormalVelocityTend at the edge against ref_layer_mean; accumulate per resolution.
The power-law fit, threshold checks, plots, and omega_vs_python comparison are unchanged in structure. Optionally overlay the reference layer-mean profile on a diagnostic plot at the finest resolution.

Implementation: task wiring and configuration

Date last modified: 2026/06/14

Contributors: Xylar Asay-Davis, Claude

task.py: remove the Reference step creation and its entry in the Analysis dependencies dict ({'init': ..., 'forward': ...} only).
horiz_press_grad.cfg: remove reference_horiz_res (no neighboring columns). Keep reference_quadrature_method for the one-dimensional integral and the layer mean. Retain the omega_vs_reference_* threshold and convergence-rate options; their values will likely need re-tuning now that the reference is exact in the horizontal (see Testing).
Delete the now-unused machinery in the old reference.py: the five-column setup, _get_reference_vert_res and _fraction_gcd, the 4th-order stencil (_compute_4th_order_gradient, _check_gradient), the per-column geometric-height integration used for the cross-column difference, and the reference_solution.nc output. Keep and relocate the quadrature helpers used by ReferenceColumn.

Implementation: orientation and consistency safeguards

Date last modified: 2026/06/14

Contributors: Xylar Asay-Davis, Claude

The reference must use the same \(+x\) orientation as Init/Forward. Derive it from xCell and assert that the resulting sign reproduces the existing omega_vs_python agreement at the coarsest resolution as a guard against a flipped gradient.
The edge surface pseudo-height (average of the two cells’ surface pseudo-heights) equals \(\tilde z_s(0)\) for linearly varying inputs; assert this equality to a tolerance to catch configuration or averaging mistakes.
Keep the EOS path restricted to TEOS-10 (as the task already requires Omega), and guard dalpha_dx against requesting specvol_first_derivatives outside gsw’s valid range.

Testing

Testing and Validation: arbitrary-pseudo-height evaluation and near-bathymetry validity

Date last modified: 2026/06/14

Contributors: Xylar Asay-Davis, Claude

Add unit tests for ReferenceColumn (e.g. tests/ocean/horiz_press_grad/test_reference.py):

Constant-density column. With a configuration giving spatially constant \(\alpha\) (so \(\partial_x\alpha = 0\)), hpga must equal the closed form \(-g[\eta' - \rho_0\alpha\,\tilde z_s']\) at every pseudo-height, to round-off, for several arbitrary \(\tilde z\) including values that fall between any plausible grid points. (With a nonzero surface pseudo-height slope \(\tilde z_s'\), evaluate the boundary term at the surface, where the integral is exactly zero.)
Zero horizontal variation. With all input gradients zero, hpga must be identically zero (no spurious pressure-gradient signal) for arbitrary \(\tilde z\), including immediately above the seafloor.
Independence from sampling. Evaluate hpga on two unrelated pseudo-height sets (one of which deliberately does not align with any uniform grid) and on a shifted set emulating a nonzero surface pressure; values at coincident pseudo-heights must agree to round-off, demonstrating freedom from grid alignment.

Testing and Validation: fidelity beyond second order

Date last modified: 2026/06/14

Contributors: Xylar Asay-Davis, Claude

Brute-force gradient cross-check. For a \(\tilde z\) comfortably above the seafloor, compare ReferenceColumn.hpga to a high-accuracy finite-difference estimate formed by fully integrating \(z(\tilde z; x)\) at \(x = \pm\delta\) and differencing at constant \(\tilde z\), for a sequence of shrinking \(\delta\). The analytic value must agree with the Richardson-extrapolated finite difference to many digits, confirming the chain-rule/Leibniz implementation and that the reference carries no fixed-order horizontal truncation error.
Quadrature convergence. Confirm that refining the quadrature panel count changes hpga by amounts decreasing at the expected high order and that the result is converged well below the model HPGA error levels.

Testing and Validation: general horizontal input variation

Date last modified: 2026/06/14

Contributors: Xylar Asay-Davis, Claude

Tilted nodes. With z_tilde_grad \(\neq 0\) (tilted profile node pseudo-heights), verify dalpha_dx against an independent finite difference of \(\alpha\) in \(x\) at fixed \(\tilde z\), confirming the analytic-input differencing handles moving nodes.
Each gradient in isolation. Exercise nonzero salinity_grad, temperature_grad, geom_ssh_grad, and the surface z_tilde_grad node one at a time and confirm the corresponding boundary or integrand contribution appears with the correct sign and magnitude. The surface boundary terms (geom_ssh_grad, surface z_tilde_grad) are kept general so a follow-up nonzero-surface-pressure test needs no further changes.

Testing and Validation: layer-averaged comparison and regression behavior

Date last modified: 2026/06/14

Contributors: Xylar Asay-Davis, Claude

Layer-mean correctness. For a known smooth \(a(\tilde z)\), verify layer_mean_hpga returns the analytic layer averages, and that for vanishing layer thickness it approaches the pointwise value.
Regression suite. Run all three task variants (salinity_gradient, temperature_gradient, ztilde_gradient) in the omega_pr suite. Confirm the omega_vs_python comparison still passes at its tight tolerance (unchanged scheme), and that omega_vs_reference produces a convergence slope in the expected range. Because the reference is now exact in the horizontal, re-tune omega_vs_reference_convergence_rate_min/max, omega_vs_reference_high_res_rms_threshold, and omega_vs_reference_convergence_fit_max_resolution to the observed behavior; document the updated values and the resolution at which curvature in the convergence plot appears.
Tilted-layer case. Pay particular attention to ztilde_gradient, whose slanted layers exercise the slant-induced pressure-gradient error; confirm the layer-averaged comparison yields a clean convergence rather than a spurious floor from charging layer tilt as error.