Add finite-volume pressure gradient (phase 1 -- 2nd order) - #493
Draft
xylar wants to merge 30 commits into
Draft
Conversation
Add the first derivatives of the specific volume with respect to conservative temperature, absolute salinity and pressure to the Teos10Eos functor, as needed by the higher-order pressure gradient (PGradHighOrder.md 3.3). The GSW-C license permits redistribution only without modification, so gsw_specvol_first_derivatives cannot be adapted into Kokkos device code, and calling into the library from a device kernel is not possible in any case. The derivatives here are instead the analytic derivatives of the Roquet et al. 2015 75-term polynomial that calcPCoeffs and calcRefProfile already carry, differentiated with respect to each of the three normalized variables and converted back with the chain rule. The temperature and salinity coefficient tables are the same values already present in Teos10BruntVaisalaFreqSq::calcAlpha and calcBeta, regrouped by power of pressure so they follow the calcPCoeffs/calcDelta split. That they are the same values was verified independently before writing this code, by differentiating the 75-term table mechanically and comparing term for term: all 50 A and all 50 B coefficients agree exactly in decimal arithmetic, since A(i,j,k) = (j+1) V(i,j+1,k) and B(i,j,k) = (i+1) V(i+1,j,k). The duplication that leaves behind is removed in a later commit. The pressure derivative needs no new constants at all, being formed from the pressure coefficients calcPCoeffs already assembles, plus the derivative of the reference profile. That last term is easy to overlook and is not optional: it supplies about 99 percent of d(SpecVol)/dP. Also hoist the polynomial normalization constants to static members, so the specific volume and its derivatives are guaranteed to be evaluated at the same normalized state rather than at two spellings of it. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Give LinearEos and ConstantEos the same calcSpecVolDerivs entry point the TEOS-10 functor now has, so the higher-order pressure gradient can be run with any of the three equations of state rather than only with TEOS-10. Both are closed form. For the linear EOS the derivatives with respect to conservative temperature and absolute salinity are -DRhodT and -DRhodS times the square of the specific volume, and there is no pressure dependence at all; for the constant EOS all three vanish. The idealized Polaris cases that use these options therefore get an exactly known reference to test the pressure gradient against. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Add the array-level entry point that the higher-order pressure gradient will call, dispatching on the configured EOS to the calcSpecVolDerivs kernel of the corresponding functor and filling the specific volume and its three first derivatives in one pass. Eos owns the SpecVolDCt, SpecVolDSa and SpecVolDP arrays just as it owns SpecVol, allocating them in the constructor and registering them as fields in the Eos group so they can be written to a stream. Their valid range spans the full range of Real rather than starting at zero, since the salinity derivative is negative everywhere and the temperature derivative is negative in cold, nearly fresh water. Since computeSpecVolAndDerivs fills SpecVol as well, it replaces a call to computeSpecVol rather than accompanying one. The two are kept separate because the derivatives roughly double the TEOS-10 arithmetic per cell and layer, and only the higher-order pressure gradient needs them. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Check the new derivatives against the unmodified GSW-C library over a cross product of 216 states spanning the oceanographic range and its corners, including fresh water, sub-zero temperatures and the full pressure range. Add an array-level test over the mesh with a state that varies with depth, so the device kernel and the vertical chunking are exercised over a range of values rather than a single one, and closed-form tests for the linear and constant options. Also compare the thermal expansion and haline contraction coefficients used by the Brunt-Vaisala frequency against gsw_specvol_alpha_beta. Those two functions were previously exercised only through a single hardcoded BruntVaisalaFreqSq value at a tolerance too loose to catch a mistake in either; the check is added before the following commit rewrites them. Measured agreement with GSW-C: 2.4e-14 for the specific volume, 3.1e-14 for the temperature derivative, 8.5e-15 for the salinity derivative, and 2.4e-14 for alpha and beta. The pressure derivative agrees only to 2.3e-12 and is gated separately at 1e-10. The discrepancy is GSW-C's: its v_P comes from a table of coefficients pre-multiplied by their pressure exponents and rounded, so it departs from the exact derivative of the 75-term polynomial by about 2e-12 at 10000 dbar, growing with pressure. Evaluating the exact derivative in 60-digit arithmetic puts the Omega value within 1e-16 of it and GSW-C's at 1.8e-12, and rounding the coefficient table to 11 digits reproduces GSW-C's error pattern. The looser gate therefore bounds GSW-C's rounding rather than ours. Add a finite-difference check as well. It is redundant with the GSW-C comparison while that library is present and correct, which is the point: it pins the unit convention of the Omega interface, per degC, per (g/kg) and per Pa, without reference to GSW, and would catch a pressure derivative that silently became per dbar. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Rewrite Teos10BruntVaisalaFreqSq::calcAlpha and calcBeta over the TEOS-10 derivative helpers and delete the A and B coefficient tables, which were a second copy of the same 100 constants. Alpha is the temperature derivative of the specific volume divided by the specific volume and beta is minus the salinity derivative divided by it, so both follow directly. The coefficient assemblers are made static, since they carry no state, so the frequency functor can reach them without holding an equation of state instance. This changes answers at roundoff level: the same polynomial is evaluated in a different Horner arrangement, and the normalized salinity is now formed as sqrt((Sa + DeltaS) / SaNorm) as calcPCoeffs does it rather than from a separately rounded reciprocal. Measured against the GSW-C library over the 216 test states, the maximum relative difference moves from 2.4e-14 to 3.1e-14 for alpha and from 2.4e-14 to 3.2e-14 for beta. BruntVaisalaFreqSq is gated at 1e-10 against a hardcoded value and continues to pass unchanged. The guard test added in the previous commit is what makes this safe: alpha and beta are now pinned to gsw_specvol_alpha_beta over the full state range, so a mistake in the rewrite fails there rather than hiding inside the single frequency value. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Describe computeSpecVolAndDerivs in the developer guide: the signature, the unit convention, that the results land in Eos members registered as fields in the Eos group, that it replaces rather than accompanies computeSpecVol since it fills SpecVol too, and that the Brunt-Vaisala expansion and contraction coefficients are now derived from them. Record the licensing position, since it is the reason the implementation looks the way it does, and the measured agreement with GSW-C, including why the pressure derivative agrees less well than the other two. Add a shorter user guide section noting that the derivatives exist for all three EOS options, that they carry no configuration of their own, that they cost nothing unless a scheme that needs them is enabled, and that they are available for output like SpecVol. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Rename the single-state entry point to calcSpecVolAndDerivsAtPoint and the vertical-chunk entry point to calcSpecVolAndDerivsInChunk, so the two are told apart at the call site rather than by the presence of "And" in the name. Document what each is for.
Say explicitly that the higher-order pressure gradient will need the derivatives every time step and that computeSpecVolAndDerivs then replaces the computeSpecVol call, while configurations that do not use it keep the cheaper call. Note why there is no displaced counterpart and what adding one would take.
The Brunt-Vaisala frequency evaluates alpha and beta at the interface, from temperature, salinity and pressure averaged across the two adjacent layers, while SpecVolDCt and SpecVolDSa hold the derivatives at the layer centers. Record that in the calcAlpha and calcBeta comments so the duplication does not look accidental.
Both compare against GSW-C, but testEosTeos10Derivs covers the array-level machinery -- dispatch, chunking, layer masking, member arrays and field registration -- on one realistic profile, while checkValueGswcSpecVolDerivs covers the polynomial itself point by point at the corners of the oceanographic range. Say so in each test's comment.
AuxiliaryState::computeMomVertAux is the only caller of computeSpecVol; computeBruntVaisalaFreqSq consumes the SpecVol array rather than recomputing it, and computeSpecVolDisp is a separate evaluation. Describe the actual call graph, and give the two reasons the plain computeSpecVol is still wanted: the centered pressure gradient is the config default, and the VertMix refresh of SpecVol feeds nothing that reads the derivatives.
Add a double-precision R8Vert variable to the file IOTest writes, then read it back into a single-precision buffer and read the existing R4Vert back into a double-precision buffer. IO::readNDVar goes through PIOc_get_var/PIOc_get_vara, which fill the caller's buffer using the type of the variable in the file rather than the type of the buffer. Reading a double variable into a float buffer therefore writes eight bytes per element into four-byte slots, which overruns the buffer and corrupts the heap. This is why VertCoord::init() corrupts the heap in an OMEGA_SINGLE_PRECISION build: VertCoordMovementWeights is the only non-distributed real field read from the mesh file, and it is stored there in double precision. Both destination buffers are allocated at twice the length needed, with the unused half set to a guard value, so the defect is reported as wrong values and overwritten guards rather than as an out-of-bounds write. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
IO::readNDVar called PIOc_get_var and PIOc_get_vara, which pass PIO_NAT to SCORPIO and therefore fill the caller's buffer using the type the variable has in the file rather than the type of the buffer. Reading a double variable into a single-precision array wrote eight bytes per element into four-byte slots, overrunning the buffer. This corrupted the heap in every OMEGA_SINGLE_PRECISION build that brought up VertCoord. VertCoordMovementWeights is dimensioned on NVertLayers alone, making it the only non-distributed real field Omega reads from the mesh file, and it is stored there in double precision; the read overran its 60-element buffer by 240 bytes. The damage was silent and surfaced later as a segmentation fault in an unrelated read or an invalid free during teardown. Double-precision builds were unaffected because the two types have the same width there. Pass the type of the destination buffer to readNDVar and dispatch to the type-specific SCORPIO calls so the values are converted on read. In IOStream::readFieldData the type is set in the same switch that allocates the buffer, so the two cannot drift apart. Distributed reads already carry this information in the decomposition and are unchanged. Also fixes the frame-time read in writeStream, which read the file's time variable into an R8 and would have produced garbage for a stream written with reduced precision. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Register VERTCOORD_SINGLE_PRECISION_TEST, a single-precision build of the existing VertCoord test driver, following the pattern of TEND_PLANE_SINGLE_PRECISION_TEST. This is the coverage the readNDVar fix makes possible: VertCoord::init() reads the mesh file, which was corrupting the heap in any single-precision build, so no single-precision test could bring up a vertical coordinate. The checks were written as a fixed 1e-10 absolute tolerance, which single precision cannot meet. The expected values reach a few thousand once cell and layer indices accumulate down a column, and the sea surface height and target thicknesses are arrived at by cancelling contributions from every layer, so their error is set by the magnitude of the column rather than by the size of the result. Replace the hand-rolled comparisons with isApprox() from OceanTestCommon.h and a relative tolerance plus absolute floor, chosen per precision. Double precision is unchanged: with RTol = 0 and ATol = 1e-10, isApprox() reduces to exactly the comparison that was there before. The one difference is that isApprox() rejects NaN, which the old form silently accepted since a NaN difference is not greater than the tolerance. Also correct "Err += Err + 1" to "Err += 1" in the mid-pressure check. That doubling only affected the reported count, not pass or fail, but it overflows a signed int after about thirty failures, which single precision makes reachable. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
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. Phase 1 of 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 that both phases share. 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 Phase 1 default (HorzOrder 2, VerticalReconstruction 'linear', QuadraturePoints 2) is used. Both phases share these keys, so Phase 2 will add values rather than keys. Reject the Phase 2 values HorzOrder 4 and VerticalReconstruction 'ppm' with a clear "not implemented in Phase 1" error rather than falling back silently, which would make a Phase 2 run look like a Phase 1 pass. 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 the Phase 2 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 Phase 1 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 Phase 1 makes 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 Phase 2 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>
Author
|
This implements the design described in #448 |
Author
|
This work is built on #485 |
Author
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.
|
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
1 participant
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
In phase 1, we add the 2nd order version of the finite-volume pressure gradient. The finite-volume implementation improves on the centered implementation in that it is "consistent" to machine precision for a certain class of solutions. The 2nd-order implementation is designed to produce machine-zero horizontal pressure gradient for CT and SA fields that are linear functions of pressure even for tilted layers.
Phase 2 (not in this PR) will add a 4th-order implementation that guarantees consistency with more variation of CT and SA with pressure, and will produce less spurious flow in tilted coordinates with thin layers (like sigma and p-star squashed by ice shelves) that is desirable for longer-term work.
Checklist
Documentation:
Linting
Building
Testing
aurora, oneapi-ifx, mpich
chrysalis, oneapi-ifx, 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
Testingwith the following:have been run on and indicate that are all passing.
has passed, using the Polaris
e3sm_submodules/Omegabaseline-pfor both the baseline (Polarise3sm_submodules/Omega) and the PR buildPerformance related PRs: Please include a relevant PACE experiment link documenting performance before and after.
New tests:
Stealth Features