Skip to content

Fix the FiniteVolume sea-floor anchor over steep slopes - #635

Draft
xylar wants to merge 27 commits into
E3SM-Project:developfrom
xylar:omega/fix-fv-hpg-sloping-floor-anchor
Draft

xylar wants to merge 27 commits into
E3SM-Project:developfrom
xylar:omega/fix-fv-hpg-sloping-floor-anchor

Conversation

@xylar

@xylar xylar commented Oct 25, 2026 •

Copy link
Copy Markdown

Fix the FiniteVolume pressure gradient's sea-floor anchor where the two columns' floors are more than a layer apart in pressure, as with terrain-following layers over a steep slope. The anchor now compares the shallower floor with the deeper column's own interface just below it, so its shift to a common pressure never spans more than one layer. Answers are unchanged wherever the floors are within one layer, which includes flat floors and p-star with partial bottom cells. Builds on #493; only the last three commits are new.


This is the anchor error shown in this comment on #493. It does not cure the sigma seamount instability (details).

The Polaris reference needs the same change, in E3SM-Project/polaris#900.

Checklist

  • Documentation:

  • Linting

  • Building

    • CMake build does not produce any new warnings from changes in this PR
  • Testing

    aurora, intel, mpich

    • CTests Pass
    • Polaris omega_pr Pass

    aurora, intelgpu, mpich

    • CTests Pass
    • Polaris omega_pr Pass

    chrysalis, intel, openmpi

    • CTests Pass
    • Polaris omega_pr Pass

    frontier, craygnu, mpich

    • CTests Pass
    • Polaris omega_pr Pass

    frontier, craygnu-mphipcc, mpich

    • CTests Pass
    • Polaris omega_pr Pass

    pm-cpu, gnu, mpich

    • CTests Pass
    • Polaris omega_pr Pass

    pm-gpu, gnugpu, mpich

    • CTests Pass
    • Polaris omega_pr Pass
  • Provide relevant details in a comment to the PR titled Testing

  • New tests:

    • CTest unit tests for new features have been added per the approved design.

Posted by Claude Code on @xylar's behalf. The wording above is AI-authored; please check it accordingly.

🤖 Generated with Claude Code

xylar and others added 26 commits October 12, 2026 01:59
Rename PressureGradType::HighOrder1 to FiniteVolume and drop the
unimplemented HighOrder2, leaving a commented stub for a future variant.
Rename the functor class PressureGradHighOrder to
PressureGradFiniteVolume and the PressureGrad member HighOrderPGrad to
FiniteVolumePGrad.  Accept 'FiniteVolume' or 'finiteVolume' in place of
'HighOrder1' when parsing PressureGradType, and update the PGrad dev and
user guides accordingly.

The new scheme uses the same two-cell horizontal stencil as the centered
scheme and is second order, so naming the class for an order it does
not have would misdescribe what ships.  FiniteVolume names the
layer-integrated control-volume form instead.

This is answer-neutral: a pure rename plus a config-facing string change.
PressureGradCentered is untouched.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Add HorzOrder, VerticalReconstruction and QuadraturePoints to the
PressureGrad configuration group, parsed and validated in the
PressureGrad constructor and cached on the PressureGradFiniteVolume
functor.  All three are optional, so a configuration written before they
existed still parses; where a key is absent the default
(HorzOrder 2, VerticalReconstruction 'linear', QuadraturePoints 2) is
used.

Reject values that are not implemented, such as HorzOrder 4 and
VerticalReconstruction 'ppm', with a clear "not implemented" error
rather than falling back silently, which would make a run requesting
them look like a run that used them.
QuadraturePoints is an accuracy knob only: the matched-pressure integrand
is zero pointwise on the exact set, so no quadrature rule can break the
robustness property.  There is no 'constant' reconstruction and no
setting that reduces the scheme to Centered.

Add accessors for the configured type and options, and extend PGradTest
with a parse and dispatch check that creates a FiniteVolume instance.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Widen PressureGrad::computePressureGrad and
PressureGradFiniteVolume::operator() to take the layer-mean conservative
temperature and absolute salinity and the Eos specific volume
derivatives.  PressureGradCentered is untouched: its operator() keeps its
existing signature and the Centered dispatch branch passes the same
arguments it always did.  Update the Tendencies call site to supply Theta
and S from the tracer array and the Eos instance it already holds.

Switch AuxiliaryState::computeMomVertAux to call
Eos::computeSpecVolAndDerivs in place of Eos::computeSpecVol when the
FiniteVolume scheme is selected.  computeSpecVolAndDerivs fills SpecVol
as well, so it replaces rather than accompanies computeSpecVol; the
branch keeps a Centered run from paying for the derivative arithmetic.

The FiniteVolume path is still a no-op, so Centered answers are
unchanged.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Assert design section 3.9's identity directly in PGradTest: with the
tidal and self-attraction-and-loading potentials zero,
PressureGradCentered is exactly -g/d_e times the fixed-pressure shift

   S = 1/2 (dZ_k + dZ_k+1) + alphaBar / (2 g) (dq_k + dq_k+1)

with dZ and dq the cross-edge differences of GeomZInterface and
PressureInterface.  The check runs over a sweep of interface tilts and is
scaled by the size of the hydrostatic terms the centered functor forms
and cancels, since agreement is to round-off rather than bit-for-bit.

Factor the two-column state construction out of the refinement loop into
setupTwoColumnState so the sweep can rebuild it at each tilt.  The
refinement loop is unchanged in what it computes.

This is the permanent regression test of design section 5.5, not a
transitional one.  Because the two expressions read the mesh, VertCoord
and Eos state through independently written code, their agreement tests
the shared upstream state -- edge masks, interface indexing, VertCoord
conventions -- and not just the pressure gradient arithmetic.  A guard
asserts that dropping the pressure half of S misses by far more than the
tolerance, so a pass cannot come from both sides being trivially equal.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Make an unrecognized PressureGradType a fatal error instead of falling
back to the centered scheme with an informational message.  A silent
fallback turns a typo, or a configuration naming a scheme that no longer
exists such as the former 'HighOrder1', into a run that looks like a
passing centered run.  This is the same hazard unimplemented
sub-option values are already rejected for.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Add PGradRecon.h, holding the mean-preserving linear reconstruction of
conservative temperature and absolute salinity in pressure that the
FiniteVolume pressure gradient uses.  The slope is a centered difference
of the layer means with respect to mid-layer pressure, one-sided in the
shallowest and deepest valid layer of a column and zero where a column
has a single valid layer.

Differencing against the actual mid-layer pressures rather than against
layer index is what makes the estimator exact for a profile linear in
pressure on Omega's non-uniform vertical grid: the mean of a linear
function over a layer is its value at the layer's midpoint, so the layer
means lie exactly on the line as a function of mid-layer pressure.  A
formula that assumed uniform thickness would pass a uniform-thickness
test and fail on a real grid.

The deviation integrates to zero over the layer because PressureMid is
the exact arithmetic midpoint of the two interface pressures, which is
also what makes VertCoord's midpoint rule the exact layer integral of a
linear reconstruction and so spares VertCoord any change.  There is no
limiter.

The functions take scalars so they can be unit tested without a mesh.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Check the mean-preserving linear reconstruction on a single column,
without a mesh.  Layer means are sampled from a profile exactly linear in
pressure, so the recovered slope must match the profile's slope to
round-off in every layer, including the two that use a one-sided
difference.

The layer thicknesses are deliberately non-uniform, which is the case the
test exists for: a formula that assumed uniform thickness would pass the
uniform control that runs alongside and fail here.  Running both means a
failure says which of the two it is.

Also asserted are the two properties the scheme leans on: that the
reconstruction reproduces the layer mean at mid-layer pressure, and that
its deviation integrates to zero over the layer, by two-point Gauss
quadrature which is exact for a linear integrand.  A single-valid-layer
column must give a zero slope, a constant being the only mean-preserving
reconstruction available there.

This is worth testing separately because it is the most likely place for
the exactness gate to fail, and it localizes such a failure immediately.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Add PGradFiniteVolume.h with the first-order equation-of-state expansion
shared across an edge layer, built by averaging the two adjacent cells'
coefficients and reference states.  One set multiplies both columns,
which is what the robustness property rests on: give each column its own
expansion point and the alpha0 and alphaP terms no longer cancel in the
matched-pressure difference, and the scheme generates spurious flow out
of nothing but its own equation-of-state approximation.  Which set is
used is an ordinary accuracy question and cannot break exactness.

Add findLayerForPress to PGradRecon.h, the lookup that makes the scheme a
fixed-pressure comparison rather than a fixed-layer-index one.  At a
pressure in edge layer K each column supplies its state from whichever of
its own layers contains that pressure, which under tilt is generally not
layer K.  A pressure outside the column clamps to the outermost valid
layer and extrapolates its reconstruction, which is the rule where the
edge control volume extends past a column's own floor; exactness survives
extrapolation because an extrapolated reconstruction still reproduces a
resolved profile.

Pin the lookup with direct property tests, on a fabricated column pair
without a mesh.  These are not optional: no answer-level check anywhere
in the test plan can distinguish this from a layer-index lookup, since on
a profile linear in pressure every layer's reconstruction is the same
line and looking up the wrong layer costs nothing.  The tests assert that
the returned layer brackets the pressure, that it differs from the edge
layer index by up to three layers so the lookup is demonstrably not an
index lookup, that the answer is independent of the starting hint so the
column scan's incremented cursors cannot disagree with a search, and that
a pressure outside the column clamps to the outermost valid layer.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Add matchedPressSpecVolDiff, the central quantity of the scheme: the
difference in specific volume between an edge's two columns, evaluated at
one pressure with one shared expansion.  Because both columns use the
same expansion, the SpecVol0 and SpecVolDP terms are identical and cancel,
leaving a coefficient times the horizontal contrast in reconstructed
temperature and salinity at matched pressure.  Compressibility drops out
entirely, which is correct physics rather than an approximation: a
horizontally uniform compressibility exerts no horizontal pressure
gradient.

Add the Gauss-Legendre rule for 1 to 4 points, and move
MaxPGradQuadPoints alongside it.  The number of points is an accuracy
setting only, since the integrand is zero at every point on a resolved
profile.

Unit test on a fabricated column pair, without a mesh, exercising the
reconstruction, the pressure lookup and the shared expansion together.
Both columns' layer means are the exact layer averages of one prescribed
continuous profile, so under the offset between the columns those means
genuinely differ.  On a profile linear in pressure the integrand is zero
at every quadrature point, not merely in the integral, measured at
5.8e-21 m3/kg or 0.9 eps of the terms that had to cancel.  A quadratic
profile runs as a control and gives 5.3e-10 m3/kg, eleven orders larger,
so the linear result is not a statement that everything in the test is
zero.  The tolerance is scaled by the cancelling terms rather than being
an absolute constant.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Accumulate the fixed-pressure height difference down each edge's column,
into a per-edge array of interface values.  This is a prefix sum with
edge-dependent coefficients, so it is not expressible as an independent
per-vertical-chunk operation and cannot live in the functor; it is the one
structural addition beyond the per-edge, per-chunk pattern the centered
scheme uses.  It runs as a parallelForOuter over edges with a
parallelScanInner down the column, in the same shape as
VertCoord::computeGeomZHeight.

The scan runs in two passes.  The first evaluates the matched-pressure
integrand at each edge layer's quadrature points and forms two integrals
of it over the same points: the increment of the recurrence, and its
first moment about the layer's top interface, which the layer mean will
need.  The second anchors the column and turns the increments into the
height difference at every interface.  The per-cell reconstruction slopes
are computed once and reused across each cell's edges, so the per-edge
work is polynomial arithmetic on cached values.

Anchor at the sea floor, per design section 3.7.4, which leaves the end
open and prefers this one on conditioning grounds.  There turns out to be
a second and stronger reason.  VertCoord builds geometric height by
accumulating a midpoint rule over each column's own layers, and on a
curved profile two columns with different layer partitions give sums that
differ at second order in layer thickness.  Anchored at the surface that
discrepancy enters the height difference directly and the exactness gate
would not be machine-zero; anchored at the sea floor over a flat floor the
height difference is exact input and vanishes identically, and the
discrepancy never enters, because the scheme integrates its own
reconstruction rather than accumulating VertCoord's height.  This is a
prediction to be checked against the measurement.

The anchor is computed, not assumed: each column's height at the deepest
shared interface is shifted to the common pressure by integrating its own
reconstruction over half the cross-edge pressure difference, and both
short integrals vanish where the two columns' interface pressures agree.

Test that the height difference is zero at every interface, not merely in
the layer mean, on a two-column state whose layer means are the exact
layer averages of one profile linear in pressure -- so that under the
offset between the columns those means genuinely differ.  The shape of a
failure is reported as well as its size, since a residual growing with
depth points at the recurrence while one flat with depth points at the
anchor.  The tolerance is provisional pending the first measurement.

Allocate the working arrays only when the FiniteVolume scheme is
selected, so a Centered run pays no memory for them.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Turn the column scan's output into the tendency: form the layer mean of
the fixed-pressure height difference from its value at the layer's bottom
interface and the first moment of the integrand over the layer, scale by
-g/d_e, subtract the tidal and self-attraction-and-loading difference, and
accumulate into Tend with EdgeMask exactly as the centered functor does.
The edge control volume's pressure thickness is exactly the edge average
of the two columns' own.

This is the first point at which the FiniteVolume scheme produces a
number.

Narrow the functor's arguments to the two arrays the column scan fills.
Design section 4.1.3's illustrative signature lists the temperature,
salinity and specific volume derivative arrays, but it also omits the
scan's own array while describing the functor as reading it "like any
other input", so it is not self-consistent.  With both integrals formed
in the scan over one set of quadrature points, which section 3.5.1
requires, the functor reads none of the state the scan consumed; passing
it would only invite evaluating the integrand a second time.

Report the tendency on the exact-set state alongside PressureGradCentered's
on the same state, and require the centered value to be nonzero, so that
the FiniteVolume result is a cancellation against a known nonzero
reference rather than a comparison with zero.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
The gating test of design section 5.2, in three groups.

Exactness.  On a two-column state whose layer means are the exact layer
averages of one profile linear in pressure, with a horizontally uniform
surface pressure, the tendency and the fixed-pressure height difference
must both be zero to machine precision.  Measured at commit 10: tendency
1.6e-18 m/s2 against PressureGradCentered's 1.6e-05 m/s2 on the same
state, and height difference 6.7e-16 m, both 0.027 eps of the terms that
had to cancel.  The gates replace commit 9's provisional value with four
eps, which leaves over a hundredfold margin while still catching drift,
and they track Real's epsilon and the size of the hydrostatic terms
rather than a physical tolerance.  Absolute tendencies are reported
alongside the centered scheme's, so a pass is a cancellation of thirteen
orders against a known nonzero reference.  The maximum is the right norm
here, since the claim is a zero at every edge and layer.

Convergence.  Quadratic and cubic profiles, which the linear
reconstruction does not resolve, over 15, 30 and 60 layers at fixed tilt.
The residual must shrink at second order in layer thickness; a rate below
that means one of the scheme's conditions has been broken somewhere and
is a bug to find, not a tolerance to widen.  The rate is measured on the
RMS rather than the maximum: a maximum is set by whichever layer happens
to be worst, and that layer moves under refinement, which makes the rate
noisy and non-monotone.  RMS is also the norm the design's measured rates
were taken in, so the two are comparable.

Guards.  Guard (a), tilt sensitivity, runs on the kernel and is the only
one that can fire where every other check passes.  Guards (b), (c) and
(d) need the scheme's rules switched off one at a time, so they run
against a host reference assembled from the same helper functions the
kernel uses.  That reference is only trusted after it reproduces the
kernel, and the fidelity check that has content is on a curved profile
where both sides are large -- agreeing at 1e-18 on the exact set would
also be satisfied by a harness that computed nothing.  Note that the
reference shares the helper functions, so the guards establish that the
three rules are load-bearing rather than independently verifying the
arithmetic; that verification comes from the centered identity and from
the Omega-versus-Polaris comparison.

Guard (b), a cell-local expansion point, must fire: two expansion points
mean the SpecVol0 and SpecVolDP terms no longer cancel.  Guard (d), the
anchor taken as the raw height difference, must fire where the two
columns' end pressures differ and must be flat with depth, which is what
distinguishes it from (b); at a sea-floor anchor that means a surface
pressure contrast reaching the bottom interface.

Guard (c), a layer-index lookup in place of the pressure lookup, is
recorded rather than required to fire.  It cannot fire on the exact set,
because the profile is a single line in pressure and every layer's
reconstruction is that same line, so writing it as a guard that must fire
would assert something false.  It is logged either way, and if it ever
does fire an answer-level check on the lookup has become available.  Until
then the lookup is pinned only by the property tests.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Add PGRAD_SINGLE_PRECISION_TEST, which builds PGradTest against the
single-precision library.  Design section 3.7.5 predicts the exactness
gate should pass there: the scheme differences the integrand before
integrating it, so every quantity in the column scan is already a
horizontal contrast and no large quantity is formed anywhere.  Running
PressureGradCentered alongside, which does form and cancel quantities of
order a hundred metres, is what makes the result interpretable.

The gate is loosened from four epsilons to sixteen in single precision.
A build with fewer guard digits accumulates more round-off down the
column, so the same measurement in epsilons is expected to be somewhat
larger; the tolerance moves with that rather than the claim changing.

Skip the convergence rates and guards (b) to (d) in single precision.
Those measure truncation errors, and this build cannot resolve them: the
round-off floor of the tendency is about 3e-8 m/s2, which swamps the
quadratic profile's residual at 60 layers and leaves only a factor of ten
between guard (b)'s signal and the noise.  They are statements about the
discretization rather than about precision, so measuring them once in
double precision is enough.  The exactness gate and guard (a) run in
both, which is what design section 5.2 asks for.

Make the matched-pressure integrand's control a ratio of the quadratic
result to the linear one rather than a count of epsilons.  In single
precision the linear result rises with the round-off floor while the
quadratic one, being a truncation error, does not, so an absolute epsilon
count would fail for a reason that has nothing to do with the property
being tested.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Instrument Eos with a count of specific-volume evaluations, one per cell
per active layer per call, and assert against it in PGradTest.

Requirement 2.2 of the design bounds the number of equation-of-state
evaluations at about one per cell per layer per step, independent of the
reconstruction order, the stencil width and the quadrature.  Nothing else
in the test suite would notice that being violated: an evaluation inside
the quadrature loop would change run time without changing any answer, so
every accuracy gate would still pass.

The check is sharper than the requirement.  The pressure gradient
performs no equation-of-state evaluations at all -- the one per cell per
layer the requirement allows is paid once by AuxiliaryState, before the
tendency is computed, and the scheme works from the specific volume and
derivatives that call leaves behind.  So the count across a call to
computePressureGrad must be exactly zero, and it is asserted at
QuadraturePoints 1 through 4, the quadrature being where a call would
most naturally creep in.

A count of zero would also be produced by broken instrumentation, so the
test additionally confirms that one call to computeSpecVolAndDerivs
registers exactly one evaluation per cell per layer.

This is a counter comparison rather than a timing measurement, so it is
deterministic and suitable for continuous integration.  The wall-time
comparison against PressureGradCentered remains a separate, non-CI
measurement for the pull request description.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Bring the PGrad dev and user guides up to what ships.

The user guide gains the three sub-options and their defaults, and says
plainly that QuadraturePoints is an accuracy setting only, since the
integrand is zero at every point on a resolved profile and no quadrature
rule can affect the exactness.  It records that unimplemented values
are rejected rather than falling back, that an unrecognized PressureGradType
aborts, and that no setting reduces the finite-volume scheme to the
centered one.

The dev guide gains four things the code alone does not make obvious: why
the per-edge column scan cannot live in the functor, that it anchors at
the sea floor and why, that each column's state is looked up by pressure
rather than by layer index -- with the warning that no answer-level test
can catch getting that wrong, so the property tests pinning it are not
optional -- and that the equation-of-state expansion is shared across each
edge, with what is and is not load-bearing about that.

Also correct the computePressureGrad call signature in the dev guide,
which described arguments the method has never taken.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
The dispatch check in PGradTest asserted only that the tendency came out
finite, which is weak enough to pass on garbage: a specific volume
derivative array still holding the fill value its field was attached
with, 9.969e36, propagates through the arithmetic without overflowing a
double.

Require instead that the tendency be nonzero and no more than three
orders of magnitude larger than PressureGradCentered's on the same state.
The bound is one-sided on purpose.  The two schemes agree to second order
in the cross-edge pressure difference, which limits how much larger the
finite-volume answer can be but says nothing about how much smaller; on a
state whose profile is close to linear in pressure it is legitimately far
smaller, which is the whole point of the scheme.  A fill value misses the
upper bound by thirty orders, in either precision.

Also check the result of every PressureGrad::create in this test.  create
returns null if an instance of that name already exists, and
dereferencing that gives a segmentation fault with no message rather than
a diagnosable failure.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Co-Authored-By: Claude Opus 5.5 (1M context) <noreply@anthropic.com>
Make computeReconSlopes and computeColumnScan public.  nvcc rejects an
extended host-device lambda whose enclosing function is private or
protected in its class, which broke the build on GPUs.

Co-Authored-By: Claude Opus 5.5 (1M context) <noreply@anthropic.com>
Template buildEdgeEos on the array type, as findLayerForPress already
is.  PGradTest's host-side reference scan calls it with host arrays,
which are the same type as device arrays on CPUs but not on GPUs, so
PGradTest would not have compiled with nvcc.

Co-Authored-By: Claude Opus 5.5 (1M context) <noreply@anthropic.com>
The counter only saw host-side calls to the Eos compute methods, which
the pressure gradient cannot make through its const Eos pointer. An
evaluation inside the quadrature loop would be a device-side call the
counter never sees, so the cost check could not catch the regression it
was written for. Remove the counter, the check and the references to it.

Co-Authored-By: Claude Opus 5.5 (1M context) <noreply@anthropic.com>
Launch the reconstruction slopes and the matched-pressure integrals as
team loops over each column's active layers, as the other pressure
gradient kernels do, instead of flat loops over every layer that skip
the inactive ones. The integrals now share the column scan's per-edge
team, with a barrier before the scan, so the two passes run as one
kernel. Nothing reads these arrays outside the active layers, so the
inactive entries no longer need zeroing.

Co-Authored-By: Claude Opus 5.5 (1M context) <noreply@anthropic.com>
The pressure gradient is exact for linear-in-pressure water only where
the two columns' sea floors are at the same pressure. Over a sloping
floor the sea-floor anchor leaves a residual that grows with the cube
of the floor pressure difference, and on terrain-following layers it
approaches the centered scheme's error.

Co-Authored-By: Claude Opus 5.5 (1M context) <noreply@anthropic.com>
Co-Authored-By: Claude Opus 5.5 (1M context) <noreply@anthropic.com>
Where the two columns' floors are more than a layer apart in pressure, as
with terrain-following layers over a steep slope, compare the shallower
floor with the deeper column's own interface just below it, then carry
the anchor down to the edge's deepest interface with the matched-pressure
integrand. Shifting both floors to their mean pressure with one
first-order expansion left an error that grows with the cube of the
offset and is the same at every depth.

Where the floors are within one layer of each other, on a flat floor or
where partial bottom cells cut a sloping one, the anchor is unchanged.

Co-Authored-By: Claude Opus 5.5 (1M context) <noreply@anthropic.com>
A resting, linear-in-pressure state on terrain-following layers whose
floors are 6.4 layers apart. The host reference mirrors the new anchor
and must reproduce the kernel, the tendency must stay well below the
centered scheme's, and the mean-pressure anchor must fail that gate.

Co-Authored-By: Claude Opus 5.5 (1M context) <noreply@anthropic.com>
Co-Authored-By: Claude Opus 5.5 (1M context) <noreply@anthropic.com>
@xylar

xylar commented Oct 25, 2026

Copy link
Copy Markdown
Author

Testing

Chrysalis, Intel, OpenMPI, built at d81ca64 (the tip adds docs only): CTests 59/59, including the new sloping-floor test in double and single precision. With E3SM-Project/polaris#900, all 8 horiz_press_grad tasks pass. Against #493, 220 of 224 horiz_press_grad forward steps and both seamount z-star runs are bit-for-bit. The 4 steps that change are the 50 m/km hydrostatic_consistency tilts on 128 m and 64 m layers, where the error drops by up to 6x. omega_pr has not been run.

/lcrc/group/e3sm/ac.xylar/polaris_1.2/chrysalis/test_20261025/fv-anchor-fix
/lcrc/group/e3sm/ac.xylar/polaris_1.2/chrysalis/test_20261025/fv-anchor-fix-polaris

Posted by Claude Code on @xylar's behalf. The testing, analysis and wording above are AI-authored; please check them accordingly.

@xylar xylar added the bug Something isn't working label Oct 25, 2026
@xylar xylar added the Omega label Oct 25, 2026
@xylar

xylar commented Oct 25, 2026

Copy link
Copy Markdown
Author

Why this is separate from #493

This came up on 2026-10-25, while making figures for #493's reviewers, after #493 was already in review. It is kept separate so that the code under review in #493 doesn't change underneath its reviewers. Folding it in is easy if reviewers prefer: it is three commits on top of #493, plus one on top of E3SM-Project/polaris#675 (as E3SM-Project/polaris#900) for the matching Polaris reference.

It is also safe to review on its own. It changes no answers where the two columns' sea floors are within one layer of each other, which includes flat floors and p-star with partial bottom cells: 220 of 224 Polaris horiz_press_grad forward steps and the seamount z-star runs are bit-for-bit with #493. It matters for terrain-following layers over steep slopes, where #493's anchor stretches one first-order equation-of-state expansion across a floor offset of several layers. It does not cure the sigma seamount instability.


Posted by Claude Code on @xylar's behalf. The wording above is AI-authored; please check it accordingly.

@xylar

xylar commented Oct 25, 2026

Copy link
Copy Markdown
Author

What #635 changes

one edge, terrain-following layers

On one edge with terrain-following layers and temperature and salinity linear in pressure, #493's error climbs steeply with the slope and nearly matches Centered's at 100 m/km. With #635 it stays at least 29 times below Centered at every slope, and at gentle slopes the two are identical.

Omega, hydrostatic_consistency_linear

In Omega, the Polaris hydrostatic_consistency_linear sweep is unchanged except at the steepest tilt on the thinner layers, where FiniteVolume's error falls 1.8 times (128 m layers) and 6.1 times (64 m layers).

sigma seamount

On the sigma seamount, the linear-in-pressure run starts about 260 times lower in kinetic energy, but the instability still grows at the same rate. The exponential profile, whose error comes from the reconstruction rather than the anchor, barely changes.

Scripts: https://web.lcrc.anl.gov/public/e3sm/diagnostic_output/xasaydavis/omega_pr635_fv_anchor_20261025/scripts/


Posted by Claude Code on @xylar's behalf. The figures, analysis and wording above are AI-authored; please check them accordingly.

The functor's inputs now match the design's signature, so drop the claim
that they differ and keep only the reason. The sea-floor anchor is chosen
for exactness, not conditioning, as design section 3.7.4 now says.

Co-Authored-By: Claude Opus 5.5 (1M context) <noreply@anthropic.com>

This branch has not been deployed

No deployments
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

bug Something isn't working Omega

Projects

None yet

Development

Successfully merging this pull request may close these issues.

1 participant