Repository navigation
Conversation
|
This implements the design described in #448 |
|
This work is built on #485 |
Where the finite-volume pressure gradient standsWhat the centered scheme misses. Over sloping bathymetry it produces spurious currents in the bottom layer, strong enough to drive numerical instability in realistic runs. The cause is that it compares two neighbouring columns layer by layer — but when layers tilt, the same layer sits at different pressures in the two columns, so it is comparing water at different depths. What this branch does instead. It compares the two columns at the same pressure. For a resting ocean with a simple stratification the answer is then exactly zero rather than merely small, no matter how steeply the layers tilt or how rough the bathymetry. Evidence it works.
Caveats.
|
d9cc591 to
17e6406
Compare
Where things stand on 2026-09-11The branch is rebased onto develop (16 commits; #485 and the VertCoord single-precision fix have merged). The rebase changed no answers: CTests and the Polaris pressure-gradient and seamount tasks pass, and the sigma seamount runs reproduce the August numbers to three digits. The sigma instability is unchanged and still unexplained. With explicit mixing now off, Centered on z-star overturns from day 0.6, so the August 44x / 3e4 comparison is superseded: at one hour, before either run has done anything else, FiniteVolume has 92x lower max |u| on the realistic profile and is at round-off on the linear-in-pressure one. Outstanding: the omega_pr suite, other machines, a wall-time comparison with Centered, and a seamount run initialized on p-star. Posted by Claude Code on @xylar's behalf. The testing, analysis and wording above are AI-authored; please check them accordingly. |
TestingCTests on chrysalis with intel and openmpi: 50/50 pass in two builds of 17e6406, rebased on develop at d086d49. The Polaris horiz_press_grad tasks (8) and seamount/nonlinear tasks (4), from Polaris branches Omega and the Polaris reference agree to 1.6e-13 m/s2 or better for both schemes. The omega_pr suite has not been run. Posted by Claude Code on @xylar's behalf. The testing, analysis and wording above are AI-authored; please check them accordingly. |
17e6406 to
c26a9b4
Compare
TestingI ran the Polaris Area-mean bottom-layer KE (m² s⁻²):
Day 1 and 10 plots, daily statistics and the plotting script are on the portal. Both runs took about 3 minutes on 6 nodes. The pressure-gradient tendency cost 8.0 s with finite volume and 1.2 s with centered. The work directory is: The task runs 1 day by default, so I set the 10-day duration and the scheme by editing each forward step's Posted by Claude Code on @xylar's behalf. The testing, analysis and wording above are AI-authored; please check them accordingly. |
c26a9b4 to
8e86c05
Compare
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>
8e86c05 to
4d59331
Compare
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>
Testing: aurora, intelgpu, mpichCTests: 59 of 59 passed. Tested the merge of #493 ( Polaris
|
Testing: pm-gpu, gnugpu, mpichCTests: 59 of 59 passed. Tested the merge of #493 ( Polaris
|
Testing: evidence for making FiniteVolume the defaultThis PR now makes FiniteVolume the default pressure gradient. In the configuration Omega will run in E3SM (p-star, TEOS-10, no ice shelf), it removes the grid-scale bottom-layer kinetic energy that Centered produces across the open ocean in a realistic global run. In idealized tests, it matches a separate Python implementation to round-off and converges at second order where Centered does not. The weaknesses of the second-order scheme show up on tilted coordinates: on a tilted p-star coordinate, the error at rest is still linear in the tilt, and on sigma the seamount overturns within a few days. Both point to higher order. These come from Polaris Realistic globalI repeated the 10-day Bottom-layer KE (m² s⁻²):
Day 1 and 10 maps, daily statistics and the plotting script are on the portal. The pressure-gradient tendency cost 9.0 s with FV and 1.0 s with Centered, out of under 3 minutes per run on 6 nodes. p-star, TEOS-10The
In Limits of second orderTilted p-star coordinate. In Sigma. Sigma, like the coordinate under an ice shelf, is far from p-star. FV overturns in all four sigma seamount configurations, while Centered stays stable. At 1 h, FV is still 1.6x to 10x better than Centered, so the problem grows over time. On z-star, by contrast, FV is 17x to 7e7x better with a linear EOS, and with TEOS-10 it stays stable where Centered overturns.
|
|
@sbrus89 and @hyungyukang, this needs to be reviewed in concert with: |
|
|
||
| // Pass one: the two integrals of the matched-pressure integrand over each | ||
| // edge layer, from the same quadrature points. | ||
| parallelFor( |
There was a problem hiding this comment.
Do we want to use hierarchical parallelism here and only compute over active layers?
There was a problem hiding this comment.
Done in 88bba35. computeReconSlopes and the integrals now use parallelForOuter with parallelForInner over the active layers, and the integrals share the column scan's per-edge team with a teamBarrier, so the two passes are one kernel. On Chrysalis (intel), CTests pass and the horiz_press_grad, seamount and overflow tasks are bit-for-bit with 480626b, and the FV pressure-gradient time dropped by about a third. pm-gpu testing will follow.
Posted by Claude Code on @xylar's behalf. The testing, analysis and wording above are AI-authored; please check them accordingly.
| const I4 KTop = LocMinLayerEdgeBot(IEdge); | ||
| const I4 KBot = LocMaxLayerEdgeTop(IEdge); | ||
| if (K < KTop || K > KBot) | ||
| return; |
There was a problem hiding this comment.
The inactive layers are being skipped here, but it may be better to use hierarchical parallelism instead.
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>
|
@sbrus89, good call! I'm retesting with the fixes on Chrysalis and pm-gpu. |
Testing: chrysalis, intel, openmpiCTests on 88bba35: 59 of 59 pass. The 8 Posted by Claude Code on @xylar's behalf. The testing, analysis and wording above are AI-authored; please check them accordingly. |
Testing: pm-gpu, gnugpu, mpichCTests on 88bba35: 59 of 59 pass, including
All 18 pass in both builds, and all 740 output NetCDF files are bit-for-bit identical between them. Omega and the Python reference agree to 5.4e-15 m/s2 (FV) and 9.5e-14 m/s2 (centered), unchanged. The FV pressure-gradient time (
The seamount mesh is too small on a GPU for the kernels to matter much. Total run-loop times are the same within noise. The PR branch predates the fix for #599, so a plain standalone pm-gpu build links Posted by Claude Code on @xylar's behalf. The testing, analysis and wording above are AI-authored; please check them accordingly. |
It seems like the agent made a mistake and that merging with |
|
Testing on Aurora (with a fix) and Frontier now. |
Testing: frontier, craygnu and craygnu-mphipcc, mpichCTests on 966aa04 pass 62 of 62 with both compilers (
I set up the new build again with the old work directory as its baseline. On both rows, all 362 steps with baseline variables pass, with zero norms.
Work directories ( Unrelated failure
Posted by Claude Code on @xylar's behalf. The testing, analysis and wording above are AI-authored; please check them accordingly. |
Testing: aurora, intel and intelgpu, mpichCTests on 966aa04 (#493 at 88bba35 merged into
Work directoriesUnrelated failure
Posted by Claude Code on @xylar's behalf. The testing, analysis and wording above are AI-authored; please check them accordingly. |
How FiniteVolume differs from Centered, on one edgeThese figures come from a Python model of one edge that mirrors Centered compares layer k with layer k, which sit at different pressures when layers tilt. FiniteVolume compares the two columns at the same pressure. At each pressure, each column supplies its own linear reconstruction of T and S. For water linear in pressure the two agree exactly and there is nothing to integrate. For curved water they differ slightly, and that difference is the error. A parabolic reconstruction, planned for fourth order, would make quadratic water exact too. With cells 10 km apart and interfaces offset by up to 150 m, FiniteVolume is at round-off for uniform or linear water at every layer thickness. For curved water it is 3 to 800 times below Centered. The exception is terrain-following layers over a sloping floor. The column scan starts from the sea floor and integrates the first-order equation of state across the whole bottom-pressure difference, so FiniteVolume's advantage is gone by about 100 m/km. p-star with partial bottom cells keeps that difference within the bottom cell, and there Polaris Posted by Claude Code on @xylar's behalf. The figures, analysis and wording above are AI-authored; please check them accordingly. |
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>
Sigma seamount instability: not the sea-floor anchorI tested whether the sea-floor anchor error in the last figure above causes the sigma seamount instability noted in the description. It does not. An anchor that avoids that error starts the linear-in-pressure sigma run about 700 times lower in kinetic energy, but the same mode then grows at the same rate, about 2.6 per day, to the same state. Also turning off extrapolation below a column's floor makes things worse. z-star stays clean throughout, and the monotonic tracer limiter was already on. The anchor error itself is fixable by anchoring at the deeper column's interface just below the shallower floor, so the anchor never spans more than one layer. In the Python model that is bit-for-bit identical to the current anchor on flat floors and on p-star with partial bottom cells, and removes the error on sigma. The anchor fix and the instability are follow-ups. If either is unsolved when this PR and the Polaris tests merge, we'll open issues for them. Posted by Claude Code on @xylar's behalf. The testing, 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>
Add a FiniteVolume horizontal pressure gradient, phase 1 of the design in #448, and make it the default. It compares the two columns sharing an edge at the same pressure rather than the same layer index, so the gradient is zero to machine precision for a resting ocean with temperature and salinity linear in pressure, whatever the coordinate tilt, provided the two columns' sea floors are at the same pressure. Over a sloping floor a residual remains, largest with terrain-following layers. Centered is unchanged and can still be selected. Builds on #485.
Behaviour changes:
PressureGradTypedefaults to FiniteVolume, so simulations with horizontal density gradients are no longer bit-for-bit with Centered.PressureGrad:HorzOrder,VerticalReconstructionandQuadraturePoints. All are optional; the phase 2 values are rejected with an error.PressureGradTypeis now an error instead of falling back to Centered.Known limitation: on the seamount test's sigma coordinate the new scheme develops a grid-scale instability that Centered does not; the cause is not yet understood. z-star, the nearer analogue of Omega's p-star, is clean.
Not in this PR: phase 2, and the Polaris tasks that exercise the scheme (E3SM-Project/polaris#675 and E3SM-Project/polaris#681). The sea-floor anchor error where the two columns' floors are more than a layer apart is fixed in draft #635. If the instability on the seamount's sigma coordinate is still unsolved once this PR and the Polaris PRs merge, open an issue for it (see the 2026-10-25 comments).
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