Skip to content

Add higher-order horizontal pressure gradient design doc - #448

Open
xylar wants to merge 22 commits into
E3SM-Project:developfrom
xylar:omega/design/add-higher-order-hpg-design
Open

Add higher-order horizontal pressure gradient design doc#448
xylar wants to merge 22 commits into
E3SM-Project:developfrom
xylar:omega/design/add-higher-order-hpg-design

Conversation

@xylar

@xylar xylar commented Jun 28, 2026

Copy link
Copy Markdown

Checklist

@xylar
xylar requested a review from sbrus89 June 28, 2026 10:29
@xylar xylar self-assigned this Jun 28, 2026
@xylar xylar added the Design label Jun 28, 2026
@xylar
xylar force-pushed the omega/design/add-higher-order-hpg-design branch from 366d4ba to 7e05a15 Compare June 28, 2026 11:12
Comment thread components/omega/doc/design/PGradHighOrder.md Outdated
@xylar
xylar force-pushed the omega/design/add-higher-order-hpg-design branch from 7e05a15 to 92e8e3b Compare July 21, 2026 12:33
xylar and others added 7 commits July 28, 2026 07:09
Convergence rate was framed as a secondary diagnostic subordinate to
absolute error. Absolute error is the metric that matters for
simulations, but the measured convergence rate verifies that the
implementation achieves its designed order. Add Requirement 2.6 for
verified order of accuracy, make the rate band in the Polaris test a
gating criterion rather than a loosened diagnostic, and note that the
resolution sweep may need retuning to resolve a fourth-order slope.

Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
The convergence test parametrized only over the legacy centered functor
and the fourth-order scheme, so the order-2 verification gate would have
measured code this design does not change and left the new scheme's
lowest-order path unverified. Add a third test configuration running the
new implementation at ReconstructionOrder 2.

Record the decision to keep PressureGradCentered as a separate functor
rather than reimplementing it as the lowest-order configuration of the
new scheme. Keeping it preserves an independent cross-check of shared
upstream state and avoids an answer-changing change to the default, at
the cost of ~40 redundant lines. Its eventual removal is follow-up work
once FiniteVolume becomes the default.

Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
The previous draft required the scheme to return zero PGF whenever layer
means are horizontally uniform, for arbitrary tilt and arbitrary vertical
profile. That cannot be delivered: under tilt each column reconstructs on
its own interfaces, so the two columns describe different water and a
nonzero PGF is the correct answer.

Replace the claim with the three conditions the implementation must
actually meet, and state for each phase which water columns cancel to
machine precision. Also record two constraints the earlier draft missed:
the EOS expansion point must be shared across an edge rather than taken
per cell, and the sloping-interface term needs the pressure-averaged
alpha*p along the interface, not the average of the two cells' values.

Flag the VertCoord geopotential as an unmet prerequisite, since
computeGeomZHeight integrates the layer-mean specific volume and so works
from a different alpha than the pressure term would.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Phase 1 delivers the consistency fix at second order: the control-volume
form, the sloping-interface integral the centered scheme omits, the
equation-of-state expansion shared across each edge, and a linear
mean-preserving reconstruction of temperature and salinity in pressure.
Phase 2 adds the parabolic reconstruction and the wide horizontal stencil
that reach fourth order.

Mark the phase each part of the formulation belongs to, and record what
each phase delivers, depends on, and leaves open.

Resolve two things the phase split exposed. The horizontal reconstruction
of section 3.6 as previously written would have interpolated temperature
and salinity to the edge and formed alpha*p there, which does not belong
to any control volume and would give up the cancellation; require instead
that the wide stencil be a weighted sum of two-column pair contributions.
And require the vertical reconstruction to recover its own degree exactly
on non-uniform layer thicknesses, without which the exact set does not
hold.

Keep constant in-layer specific volume as a verification-only setting so
the reduction to the centered scheme survives as a permanent regression
test.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Give every test an explicit statement of the requirements, table rows and
assumptions it exercises, and add a coverage summary.

Three things had no test. The adequacy of the equation-of-state expansion
across an edge now gets a horizontal-contrast sweep in the two-column
task, confirmed under dynamics by the overflow task. The vorticity
question gets a diagnostic in the seamount task. The bound on TEOS-10
evaluations gets an instrumented count that must not change with
quadrature points or stencil width.

Also correct the seamount setup description: a stratification that is a
function of pressure alone gives layer means that differ between columns
where layers tilt, which is the situation the test needs, not the
horizontally uniform means the earlier text asked for.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
The claim that sharing the equation-of-state expansion point across an
edge costs a potential-form property was wrong. In a general ALE
coordinate the layers are neither isopycnal nor isobaric, so the pressure
gradient force is not the gradient of a scalar to begin with and the
centered scheme already carries the non-gradient term explicitly. Nothing
is given up, so the discussion, the assumption it introduced, and the
test for it are all removed rather than corrected.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
@xylar
xylar force-pushed the omega/design/add-higher-order-hpg-design branch from 92e8e3b to bfe454d Compare July 28, 2026 14:12
@xylar
xylar force-pushed the omega/design/add-higher-order-hpg-design branch from e217c33 to 5b2f9ed Compare July 29, 2026 16:11
@xylar xylar changed the title Add higher-order horizontal pressure gradient design doc Add finite-volume horizontal pressure gradient design doc Jul 29, 2026
@xylar
xylar force-pushed the omega/design/add-higher-order-hpg-design branch from 5b2f9ed to 63ba646 Compare July 29, 2026 17:02
Record the final shape of the Eos derivative API in PGradHighOrder.md 4.1.2:
Eos owns the SpecVolDCt/SpecVolDSa/SpecVolDP arrays exactly as it owns
SpecVol, allocating them in its constructor and registering them as fields in
the Eos group, so any number of consumers can read them without each
allocating a copy. The names match the Ct/Sa abbreviations the TEOS-10
helpers already use, and their valid range must admit negative values. Note
also that computeSpecVolAndDerivs fills SpecVol too, so it replaces rather
than accompanies a call to computeSpecVol.

Correct the cost claim in 3.3. The first derivatives do not reuse the same
polynomial coefficients as the specific volume: the pressure derivative does,
but the temperature and salinity derivatives need coefficient sets of their
own. What is actually bounded, and what Requirement 2.2 states, is the
TEOS-10 call count.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
@xylar
xylar force-pushed the omega/design/add-higher-order-hpg-design branch from 63ba646 to 08ef74a Compare July 29, 2026 17:03
xylar and others added 2 commits July 29, 2026 14:22
These particularly relate to a correction of the geopotential at
edges needed to achieve machine-precision consisency.
Omega defines pseudo-height as z-tilde = -p/(rho0 g) with no offset, so
surfaces of constant z-tilde are surfaces of constant pressure and the
pressure-gradient term vanishes along a coordinate surface. The four
pressure terms of the inherited control-volume form therefore sum
identically to zero, and the entire horizontal PGF is the geopotential
compared at fixed pressure.

Correct that inherited form first: it is missing a p grad(alpha) term,
without which the scheme does not converge at all, and the signs on its
two metric terms are reversed. Both follow from Leibniz' rule and are
owed back to OmegaV1GoverningEqns.md.

State the reduced target, and evaluate it as PressureGradCentered plus a
remainder. The centered scheme is exactly the first-order conversion of a
height difference at fixed layer index into one at fixed pressure, so the
remainder is the centered scheme's error and setting it to zero recovers
that scheme at any tilt. Give the remainder's four sources, three local to
the layer and one a column prefix sum, and require each to be formed
without subtracting large quantities.

The reduction is a statement about coordinate surfaces, not layer
interfaces: interfaces remain tilted, time-varying and at different
pressures in neighbouring columns under any ALE layering of pseudo-height,
and every term the scheme computes is a tilt term. Say so explicitly and
bound the scope.

Consequently drop the side-wall integral, both metric integrals, the
quadrature option and the 'constant' reconstruction; drop the
pair-decomposition constraint on the phase 2 stencil, which exactness no
longer requires; and rename the functor to PressureGradFiniteVolume,
since phase 1 is second order and not higher order than the scheme it
replaces.

Add a guard test asserting the scheme responds to tilt at all, and record
the supporting measurements and their limits.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
@xylar xylar changed the title Add finite-volume horizontal pressure gradient design doc Add higher-order horizontal pressure gradient design doc Jul 30, 2026
@xylar

xylar commented Jul 30, 2026

Copy link
Copy Markdown
Author

I have made major revisions to the design that emerged in my testing. You can find the document (for now at least) on my test build:
https://web.lcrc.anl.gov/public/e3sm/diagnostic_output/xasaydavis/test_omega_docs/html/design/PGradHighOrder.html

The previous formulation built each column's hydrostatic integral separately
and reconciled the two as the centered scheme plus a remainder. A reference
implementation showed it misses exactness by ~2e-3 of the centered scheme's
answer on a profile inside the exact set, flat across three decades of tilt,
where the requirement is round-off. A brute-force evaluation of the target,
sharing no code with the closed forms, misses by the same amount, which places
the fault in the target rather than the algebra.

The cause is that the equation-of-state expansion was shared across the edge at
fixed layer *index*. Under tilt that is not fixed pressure, so the two columns
evaluated specific volume as piecewise-linear functions whose breakpoints sat at
different pressures, and integrated genuinely different functions even for a
profile both reconstructed exactly. Every other comparison in the design had
been moved to fixed pressure; this one had not.

Difference the integrand instead. With one shared coefficient set per edge layer
and each column supplying its own theta(p), S(p) from whichever of its own
layers contains p, the alpha_0 and alpha_p terms are common to both columns and
cancel, leaving the horizontal contrast in reconstructed theta and S at matched
pressure. That is identically zero pointwise on the exact set, so exactness no
longer depends on the coefficients, the quadrature, or the interfaces lining up.

Restate condition 1 accordingly: the columns must agree at every pressure from
the anchor down, not merely within a layer.

Consequences. Compressibility drops out of the horizontal gradient entirely,
which is correct physics and means the scheme needs alpha_theta and alpha_S but
not alpha_p. No large quantities are formed anywhere, so the round-off exposure
belongs to the centered scheme alone and single precision is now expected to
pass. Assumption A5 demotes from an exactness risk to an accuracy question,
since an extrapolated reconstruction still reproduces a profile it resolves.
The quadrature setting returns as a pure accuracy knob. No configuration reduces
the scheme to the centered one, so the regression test asserts the centered
identity standalone instead; that identity is confirmed to 0.5 eps over fifty
states.

Fold in five corrections from the reference implementation, including two to
this document: the claim that each component of a correction is individually
small, and the reconciliation with the accumulated centered error, which holds
only when the two columns have equal total pseudo-thickness.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
@xylar
xylar force-pushed the omega/design/add-higher-order-hpg-design branch 2 times, most recently from 249e9c4 to 43b71f3 Compare July 31, 2026 14:38
The scheme of section 3.5 has now been implemented and measured in the Polaris
two-column harness. Exactness holds: the differenced integrand is zero at every
quadrature point to 5e-21, the accumulated height difference to 1e-15 m, and the
assembled tendency to 2e-18 m/s, across the full tilt sweep at three vertical
resolutions. Several claims in this document did not survive contact, and this
commit corrects them.

Section 3.5's first consequence overreached. Exactness does not depend on the
coefficients' values, but it does depend on one set multiplying both columns;
give each its own and exactness is lost at 1e-5. The edge-shared expansion point
remains load-bearing, for sharing rather than for where the point sits.

The recorded diagnosis of the earlier formulation's failure was wrong. It said
the cause was sharing the equation-of-state expansion at fixed layer index
rather than at fixed pressure. Measured, that choice has no effect at all,
because whatever the shared set is, it multiplies a contrast that is zero. The
mechanism was never established and, the formulation being abandoned, is left
open rather than pursued. State the condition the design actually needs instead:
the per-layer version is insufficient because the column integral accumulates
from the surface.

Requirement 2.3.4 asked for the residual to shrink one order faster in layer
thickness than the centered scheme's. It does not, because the premise - that
the centered scheme is first order in layer thickness - does not hold on the
profiles measured; both converge at roughly second order on a smooth curved
resting profile, and the advantage narrows from 6.5x at 256 m layers to 2.4x at
64 m. Express the requirement as absolute error at affordable resolution, which
is what section 1 says the metric of success actually is. The first-order
behaviour that is robust is in the coordinate tilt, not in layer thickness.

The anchor was described in terms of how one test happens to initialize its sea
surface, and that had been promoted to a standing condition on the design and a
dependency on VertCoord. It is neither. The anchor is computed from the model's
sea-surface heights and surface pressures, not assumed, and it is the k=1
instance of the same fixed-pressure comparison every other interface gets. What
is a real implementation condition is that the surface difference be shifted to
a common pressure rather than taken as the raw height difference, so condition 3
says that. The test-design consequence - that a resting-state gate must be run
on a state that is exactly at rest, or must expect the residual its
initialization implies - belongs in section 5.2 and is stated there generically.

Correct the guard list to match what fires. A cell-local expansion point breaks
exactness and is the condition-1 guard; the choice of shared coefficient set
does not and is deliberately not a guard. Record prominently that the guard
which would catch a column's state being looked up by layer index instead of by
pressure cannot be made to fire on any available configuration, so an
implementation that gets this wrong will pass every check specified here and
must be pinned by direct property tests.

Keep the design to requirements and formulation. Statements about what earlier
revisions said, and the reasons a requirement changed, belong in section 6 or in
the commit history rather than in sections 1 to 5.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
@xylar
xylar force-pushed the omega/design/add-higher-order-hpg-design branch from 43b71f3 to fd9c663 Compare July 31, 2026 15:30
The scheme has been implemented and measured in Omega.  Four things the
design left open or predicted are now fixed, and the text is updated to
describe the design as it stands rather than as it was proposed.

The column scan anchors at the sea floor.  Conditioning favoured it
already, but exactness decides it: VertCoord accumulates a midpoint rule
over each column's own layers, so on an unresolved profile two columns
with different layer partitions give sea-surface heights differing at
second order in layer thickness.  Anchored at the surface that
disagreement enters the anchor directly and the machine-precision
property is lost; anchored at the sea floor it never enters.  The earlier
claim that the two ends agree to round-off held only inside the exact
set, where the two columns' sums agree term by term.

Assumption A5's rule is stated: clamp to the outermost valid layer of the
column and extrapolate that layer's reconstruction.

Single precision is measured rather than predicted.  The exactness gate
returns the same result as double precision scaled by epsilon, while
PressureGradCentered is essentially unmoved, since its error is
truncation rather than round-off.

The functor signature is corrected.  It takes the column scan's two
output arrays, not the temperature, salinity and derivative arrays: both
integrals are formed in the scan over one set of quadrature points, as
section 3.5.1 requires, so evaluating the integrand again in the functor
would double the cost of the most expensive part of the scheme.  The
scan's arrays are described with their real extents, and the Eos
evaluation counter the cost check needs is recorded.

The cost check is sharpened to what the design delivers: zero
equation-of-state evaluations in the pressure gradient, the one per cell
per layer the requirement allows being paid by AuxiliaryState.

Section 5.2 records what each guard did, including that guard (c) is
confirmed unable to fire, and states that the guards run against a
switchable assembly sharing the kernel's helper functions -- so they
establish that the three rules are load-bearing without independently
verifying the arithmetic.

Section 6.1.1 collects the Omega-side measurements.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
xylar and others added 8 commits August 2, 2026 08:49
Equation (z-increment-exact) states that VertCoord's midpoint rule is the
exact layer integral.  That is true of the linearized specific volume the
pressure gradient uses, and the derivation given -- the mean-preserving
deviations and the pressure term both integrating to zero over the layer
-- establishes exactly that.  It is not true of what VertCoord actually
accumulates, which is the exact TEOS-10 specific volume evaluated at the
layer-mean state.  Specific volume is nonlinear in pressure, so its layer
average is not its value at the midpoint state, and the increment is a
second-order-accurate approximation rather than an exact integral.
Measured on the Polaris side, the gap accumulates to 1.1e-3 m over a
column.

The conclusion the section draws is unchanged: VertCoord needs no change
and no answer-changing baseline step.  But the reason is not that its rule
is exact.  It is that the scheme never accumulates VertCoord's height,
using it once as the difference the anchor takes at the sea floor, where
over a flat floor with matching maxLevelCell it is exact prescribed input.
Where the two columns' floors differ the anchor spans a partial
accumulation and the gap enters at the level of the difference between the
columns' partial sums, which is what bathymetry_step measures.

This also corrects the mechanism given for preferring the sea-floor
anchor.  That was attributed to profiles the reconstruction does not
resolve; the disagreement in fact comes from the equation of state's
curvature in pressure, so it is present even on the exact set, where two
columns with different layer partitions still accumulate different amounts
of midpoint-rule error.  The claim that the two anchor ends agree to
round-off therefore fails on the exact set as well, not only off it.

Reported from the Polaris side.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Section 3.7.5 claimed no large quantity is formed anywhere in the scheme.
That is true of the recurrence, every increment of which is a horizontal
contrast, and false of the anchor.  Where the bathymetry steps, the height
difference at the deepest shared interface is of order the step and the
anchor's two short integrals span that step and cancel it.  Measured on a
stepped floor with a profile linear in pressure: 4.0 m cancelled to
1.9e-12 at a floor gradient of 1 m/km, 100 m to 2.6e-8 at 25 m/km, and
205.6 m to 1.7e-7 at 200 m/km -- nine to twelve digits consumed.  On the
same states the surface anchor forms nothing at all, its height difference
being zero and its correction a zero-width integral.

The two ends therefore fail in complementary ways, and the choice is a
trade-off rather than a settled preference.  The sea floor is exact where
the coordinate tilts over an unstepped floor, which is the case
Requirement 2.3.1 was written for, and the surface is not; the surface
forms nothing where the floor steps, and the sea floor does not.  Neither
gives machine precision at every tilt, thickness and bathymetry, so that
requirement is met on an unstepped floor rather than unconditionally.
Section 3.7.4 now states this as the open question Phase 1 leaves behind.

The stepped-floor residual is not quadrature error: refining from 2 to 16
Gauss points does not move it, which is conclusive rather than suggestive
since the reconstruction is linear in pressure and every rule integrates
it exactly.  What remains is the edge-shared first-order equation-of-state
expansion integrated across the whole step, some 2e6 Pa at 200 m/km.  That
attribution is a scaling argument rather than a measurement; if it holds,
the second-order expansion would reduce it, as would re-expanding about
the midpoint of each short integral, and neither has been tried.

Two things this does not disturb, both measured: the scan stays flat to
1.4e-15 m at every floor gradient, including where a column's
reconstruction is evaluated three layers below its own floor, confirming
assumption A5 on the geometry it was most doubted on; and the scheme still
beats the centered one by a factor of 2e5 on these states.

Reported from the Polaris side.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
…dition

Section 3.5.1 still described the scan as accumulating from the sea
surface and wrote the anchor at k = 1, contradicting section 3.7.4, which
settles the sea floor and gives the argument.  The implementation follows
3.7.4.  Rewrite 3.5.1 to match: the recurrence runs upward, the anchor
sits at the deepest interface the two columns share, and its short
integrals each span half the cross-edge pressure difference there.  The
surface form is kept as the alternative it is.  Also drop the description
of the anchor end as a conditioning choice, which it is not, and qualify
the conditioning argument as holding where the floor is unstepped.

Section 5.2's machine-precision gate needs a condition on how its state
was constructed, which it did not previously state.  The anchor is the one
place a state can inject a residual the scheme cannot remove, so the gate
must be run on a state whose anchor inputs vanish identically: at a
sea-floor anchor, two columns sharing both a floor depth and a bottom
pressure.  A state that pins the opposite end instead carries VertCoord's
second-order truncation difference into the anchor's input, and then no
scheme however exact returns zero -- measured on Polaris-initialized
states, whose columns share a floor depth but reach different bottom
pressures under tilt, the sea-floor anchor is -7.1e-8 m where the surface
anchor is zero to round-off.

That residual is the anchor and nothing else, since every increment of the
recurrence is still zero and the height difference is still flat down the
column, so it leaves the pointwise cancellation untouched.

Omega's gate satisfies the condition by construction: it prescribes the
pressure grid, so both columns share a surface and a bottom pressure, and
prescribes an equal bottom depth from which VertCoord derives height
upward.  Both anchor inputs are identically zero, which is why its
3.1e-18 m/s2 is a property of the scan alone; the same scheme run on a
Polaris-initialized state should return that state's anchor residual
instead, and that is the correct answer there.

Reported from the Polaris side.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Section 5.1 called for running the temperature and salinity gradient
variants at several amplitudes of the cross-edge contrast and requiring
the error to grow no faster than linearly.  That test cannot work, for
three independent reasons, each established by measurement on the Polaris
side.

The contrast axis is not independent of the resolution axis.  A two-column
initialization interpolates at plus and minus half the cell spacing, so
the state depends on the contrast and the spacing only through their
product; contrast 2 at 2 km reproduces contrast 1 at 4 km to every printed
digit.  A contrast sweep is a coarsening sweep in disguise, which is why
the convergence slope collapses from 1.46 to -0.04 across it.

A2 and the horizontal stencil's truncation error are confounded by
construction.  Both scale with the cross-edge contrast, a comparison
against the quasi-analytic reference cannot separate them, and the stencil
dominates by enough that the two schemes give identical RMS to four
significant figures at every contrast.

The initialization's shape-preserving limiter is nonlinear in the node
values, so the state is not smooth in the contrast: the second difference
of layer-mean salinity is 5.75e-2 g/kg on the shipped profile against
7.1e-15 where the limiter has no freedom, and the nonlinearity saturates
with contrast rather than scaling with it.

What replaces it is stronger.  A2 is the claim that the second-order
remainder of the equation-of-state expansion across an edge is small, and
that remainder is computable straight from the state with no reference
comparison, so it is neither confounded with the stencil nor sensitive to
how the state was built.  Measured that way it is 8.0e-6 relative on
temperature_gradient at 4 km, the hardest configuration in the family,
with a companion check confirming it vanishes when the two columns
coincide.

Section 3.7.6's A2 entry and the section 5.7 coverage table are updated to
match.

Reported from the Polaris side.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
The eight Polaris horiz_press_grad variants have been run in Omega with
both schemes over shared initial conditions.  Every prediction in
section 6.5 was reproduced, several to four or five significant figures.

Section 5.1 was still written as a specification to be implemented.  It
now records what was built and measured: the two configurations Phase 1
actually offers rather than the three an earlier draft assumed, the
per-variant and per-scheme convergence bands, and the accuracy ratios.

The headline is ztilde_gradient, which holds temperature and salinity
horizontally uniform and tilts only the coordinate, so the first-order
term of [centered-error] is exposed with nothing else mixed in: centered
converges at 1.068 there and finite volume at 1.971.  The resting-state
behaviour of section 3.7.3 and its cure, visible in a convergence rate
rather than inferred.

Two corrections to guidance this section gave.  The accuracy gate's
default must assert "not worse" rather than "better": on
temperature_gradient and salinity_gradient the schemes agree to five
figures, correctly, because level coordinates make matched-pressure and
matched-index coincide and the residual is a two-column error both
share.  And the advantage is not always a constant factor -- where the
orders genuinely differ it widens under refinement, 1.80x at 4 km to
14.8x at 0.5 km.

Add section 3.7.3's warning that the tilt exponent must not be used to
compare the schemes.  Both give q = 1.000, for unrelated reasons: the
centered scheme because alpha is effectively piecewise constant in
pressure within a layer, the finite-volume scheme because its one
remaining error term is the anchor, which inherits the tilt's own
linearity.  Reading q alone would report "no better than centered" for a
scheme up to 500 times more accurate.

Add section 6.6 with the resting-state advantage per variant, spanning
four orders of magnitude, and with bathymetry_step_linear -- a new
variant that is both inside the exact set and stepped, so it isolates
what a bathymetry step costs.  It confirms A5 on the geometry where that
assumption was least obvious: D_k stays flat to 1.4e-15 m even where the
two columns reach different maxLevelCell and the deepest edge layers
extend three layers below one column's floor.

Record that omega_vs_polaris measures 4e-19 to 5e-16 for the new scheme,
better than the legacy scheme's, and note what maintaining that
independence required.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Section 6.6 previously asserted that interior coordinate tilt
contributes essentially nothing to the global bottom-layer error and
attributed that to section 5.3, which contained no such measurement.
The claim was ahead of the evidence as the document stood; the assertion
was removed, and this restores it properly by putting the measurement
itself in section 5.3 with its numbers and provenance.

The measurement is an offline evaluation of the centered scheme's PGF on
EC30to60E2r2 with horizontally uniform temperature and salinity, which
is an exact resting state on the real mesh and bathymetry, so the whole
computed PGF is error -- section 5.2's premise at global scale.

It establishes three things.  The error is concentrated in the bottom
layer and is baroclinic there, matching the bottom-trapped character of
the KE response; it is 1.28x to 1.55x the pressure gradient the model
actually carries in that layer, which sets the bar for Requirement 2.1
at a factor of 5-10 rather than the ~2x a modest gain would give; and
the global IC's interior coordinate is 100-1000x less tilted than the
smallest tilt in the Polaris sweep, putting that mechanism some 6000x
below the observed interior error.  It also calibrates bathymetry_step
as a quantitative proxy, within 20-35% over the slope range that
dominates the mesh and over-predicting steep edges by 3-10x.

It does not settle A4, and the subsection says so at length.  Every
number is the error rather than its dynamical consequence, which is
necessary for A4 but not sufficient: it does not show the error produces
the observed bottom-layer KE rather than merely coexisting with it.  The
one forward-run observation, bottom KE 1.19x the fully forced run, is
suggestive and no more.

Section 6.6's parenthetical now cites the narrow claim the measurement
supports and keeps A4 open, rather than conflating the two.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Section 5.3.1 sets the bar for Requirement 2.1 at a factor of 5-10, from
the global resting state where the centered scheme's bottom-layer error is
1.3 to 1.6 times the pressure gradient the model carries there.  Section
6.6 reports the advantage per variant.  Neither section referred to the
other, so the requirement stayed qualitative and the ratios had nothing to
be read against.

Requirement 2.1 now states the figure and where it comes from.  Section
6.6 reads its four ratios against it: three clear the bar, and
hydrostatic_consistency does not, at 2.58x.

That is not a shortfall, and saying so matters because the number invites
the opposite reading.  Section 5.3.1 measures the interior coordinate of a
realistic global initial condition as 100 to 1000 times less tilted than
the smallest tilt in that sweep, putting interior coordinate tilt some
6000 times below the observed global interior error.  So the variant is
calibrating the scheme on a mechanism that is not the one the bar was set
from, and its ratio is not a prediction of global improvement.  The
stepped-floor variants are the quantitative proxies for that error, and
they clear the bar by two to four orders of magnitude.

Reported from the Polaris side.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
The design specified the A4 diagnostic in three places as a pre-Phase-1
check that could be run immediately with the existing centered scheme.
That cannot work, and section 3.7.3's own table is what contradicts it: a
profile linear in pressure lies in Phase 1's exact set and not in the
centered scheme's, where it stays first order in layer thickness with a
coefficient proportional to the coordinate tilt.

The diagnostic's inference requires the scheme being run to have zero
pressure gradient error on that profile.  Run with the centered scheme the
linear-in-pressure seamount is not a null configuration, only a slightly
cleaner non-null one, and a sigma coordinate over a seamount is the
largest tilt in the suite, so the contamination is worst exactly there.
The test would measure the thing it was designed to exclude and its result
would support neither conclusion.

A4 is therefore a post-Phase-1 test and became runnable for the first time
when FiniteVolume landed.  Section 4.5.3's step 2 is rewritten as what a
centered run there can legitimately establish -- that the case exercises
the pressure gradient at all, and how large the spurious velocity is --
which is worth doing early but answers a different question.

Reported from the Polaris side, after the first Omega seamount runs.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Section 5.3 asked for the seamount test and had no measurement.  It now
has one, run with TEOS-10 over six days with tracer diffusion and
hyperviscosity off so the pressure gradient is what remains.

On z-star the scheme passes the section's criterion by a wide margin: 44x
lower maximum velocity and 3e4 lower mean kinetic energy than the centered
scheme on the realistic profile, and 134x and 1780x at one hour before
anything else develops.  The centered scheme's error sits in the partial
bottom cell, matching section 5.1's bathymetry_step.  That is Requirement
2.3 under dynamics.

On sigma it does not.  The scheme develops density inversions at day 1.5
on the realistic profile and day 3 on the linear-in-pressure one, reaching
about 5e-1 m/s by day 6 where the centered scheme on the identical
configuration is clean and decaying at about 1e-2.

It is not an accuracy failure -- before the growth takes hold the new
scheme is five times more accurate than the centered one.  Three things
point at a grid-scale numerical instability: the disturbance sits at a
fixed location on the flank, the same under both equations of state, so it
follows the geometry rather than the profile; its grid-scale content rises
before any inversion appears; and the temperature anomaly reaches a third
of the profile's range.

Sigma layers follow the bathymetry, which is the steepest tilt against
isobars in the suite and is not what Omega runs, so this does not qualify
the result on the configurations Omega uses.  It does mean the scheme has
a limit on steeply tilted layers that the centered scheme does not share,
and that the limit is not understood.  The next diagnostic is a re-run
with monotonic horizontal advection.

An earlier reading attributed the sigma behaviour to a Polaris mixing
misconfiguration.  Removing that mixing fixed the centered runs and made
the finite-volume ones worse, so the attribution was wrong: the diffusion
had been damping the mode.

The former 5.3.1 becomes 5.3.2.

Reported from the Polaris side.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

Projects

None yet

Development

Successfully merging this pull request may close these issues.

2 participants