Repository navigation
Conversation
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>
TestingChrysalis, 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 Posted by Claude Code on @xylar's behalf. The testing, analysis and wording above are AI-authored; please check them accordingly. |
Why this is separate from #493This 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 Posted by Claude Code on @xylar's behalf. The wording above is AI-authored; please check it accordingly. |
What #635 changesOn 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. In Omega, the Polaris 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. 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>



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
Testing
aurora, intel, mpich
aurora, intelgpu, mpich
chrysalis, intel, openmpi
frontier, craygnu, mpich
frontier, craygnu-mphipcc, mpich
pm-cpu, gnu, mpich
pm-gpu, gnugpu, mpich
Provide relevant details in a comment to the PR titled
TestingNew tests:
Posted by Claude Code on @xylar's behalf. The wording above is AI-authored; please check it accordingly.
🤖 Generated with Claude Code