Add TEOS-10 first derivatives of specific volume - #485
Conversation
|
This work is needed for the finite-volume pressure-gradient approach described in #448. |
TestingCTest unit tests:
Polaris
|
|
@sbrus89, please have a look and see if the organization is what you wanted or if you have suggestions. @katsmith133, could you have a look, since this makes some modifications to your previous work? I'd appreciate feedback on the organization and testing. |
|
@sbrus89 and @katsmith133, hold off for now on looking this over. I made a bad call in the design phase. I thought |
494e8fd to
0f80640
Compare
Retesting after reorgCTest unit tests:
Polaris
|
|
@katsmith133 and @sbrus89, sorry for the moving target. The reorganization is now in. |
|
I can review this on Sunday |
ab84f77 to
222301e
Compare
katsmith133
left a comment
katsmith133
left a comment
There was a problem hiding this comment.
This is great work @xylar!
I think this PR is almost there. I think most of my comments are probably due to me not understanding something, so just some explanation/clarification would be helpful.
I am working on testing on Frontier. Ctests for both craygnu and craygnu-mphipcc have passed. I am working on getting results for the omega_pr suite. Do you want me to check them off on the PR description as I go, or do you prefer to not check them off until we are testing the final code implementation?
Also, thanks for adding in more rigorous testing and making sure we don't violate GSW-C's terms. Are there any other currently applied functions that you noticed we need to make this modification for? I think some of the frazil helper functions that I am implementing will need close consideration as they are not explicitly in Roquet et al., 2014.
Testing on FrontierCTest unit tests:
Polaris
|
|
@katsmith133, thanks so much for your review and testing, particularly over the weekend! I'll reply to your comments shortly. |
Good question. I think there are 4 existing functions where we might need to re-implement them without directly referencing our using the GSW code: They're untouched by this PR and it wouldn't make sense to address them here (you weren't suggesting it was). But an issue would probably be a good idea. I'll try to find time. For the frazil helpers, the pattern I'd suggest is the one used here: implement from the published equations rather than from the GSW-C source, cite the paper in the comment (the GSW documentation names the source reference for each routine, which makes this straightforward), and use GSW-C unmodified through its public API for testing. That keeps following their license. |
Ok, great, those were the ones I was thinking we would need to re-do. We touch some of them with the thermo-coupling/frazil codes, so perhaps I can address some of them there. Though I think just getting those PRs in is the higher priority. So an issue and a later PR to address them all might be good. Beyond that, I think you have addressed all of my comments and questions. Is there anything else you'd like for me to do to help with this PR? Testing on PM, re-testing on Frontier, etc? Those are easy enough that I should be able to do them while at the meeting this week. |
c508b3a to
48ce21e
Compare
RetestingCTest unit tests:
Polaris
|
|
@katsmith133 and @sbrus89, anything else I can do to help this along? |
|
Sounds good @katsmith133, thanks! |
|
I tested a local merge of this PR on Frontier. There were machine-precision level diffs with the CTests pass with both compilers. Polaris
|
|
@sbrus89 even the small diffs are unexpected. I will have to think about whether there's an explanation for them. |
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.
Vertical chunking is on its way out of Omega (PR E3SM-Project#473 removes it from the auxiliary variables and tendency terms), so this branch should not add more of it. The EOS is not touched by that PR, so remove chunking here across the board rather than leaving the new specific volume derivative code to be converted later. The functors now follow the same pattern as the de-chunked auxiliary variables: each takes a TeamMember and a cell index and loops over the active layers with parallelForInner(Team, Range{KMin, KMax}, INNER_LAMBDA(int K) { ... }); so the callers in Eos.cpp reduce to a single call inside parallelForOuter, with no vertRangeChunked and no inner loop over chunks. Dropping the chunk loop also removes VecLength from the TEOS-10 polynomial helpers. calcPCoeffs, calcPCoeffsDTt, calcPCoeffsDSs, calcDelta, calcDeltaDeriv and calcDeltaDP took arrays sized 6 * VecLength or 5 * VecLength together with a KVec index into them; they now take plain [6] and [5] arrays and no index. That storage was never shared between layers -- the coefficients are recomputed for every layer -- so nothing is lost, and the point-wise calcSpecVolAndDerivsAtPoint and the calcAlpha and calcBeta helpers no longer allocate VecLength times more stack than they use. With the chunk loop gone, the array-level derivative routine is exactly calcSpecVolAndDerivsAtPoint evaluated at each cell and layer, so it now calls it instead of repeating the polynomial evaluation. That leaves a single implementation of the TEOS-10 derivatives. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
The InChunk suffix named an implementation detail that no longer exists, and it was flagged in review as unusual and potentially confusing. The auxiliary variables name these routines for the mesh element they act on rather than for the vertical loop structure -- computeVarsOnCell, computeVarsOnEdge, computeVarsOnVertex -- and keep those names through the removal of chunking, so follow that convention here. The point-wise entry point keeps the name calcSpecVolAndDerivsAtPoint; the pair now reads as scalars in and out versus arrays over the mesh. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
48ce21e to
fa2ffdf
Compare
|
@sbrus89, I'm also seeing the same diffs as you. I imagine some order of operations changed on craygnu on Frontier. I think we should bless the change since it's clearly at machine round-off. |
Testing on Frontier (CPU)Polaris
|
|
Ctests pass on PM-CPU and PM-GPU. Waiting for LCRC to come back up (down for maintenance today) before I can finish the |
Testing on Perlmutter (CPU)Polaris
|
|
Thank you @katsmith133! |
|
Local merge passes CTests and Polaris
|
|
@sbrus89, thank you so much for your testing and review, and coordinating these Omega merges. |
- Brings in non-B4B behavior on Frontier (craygnu) due to the TEOS10 derivative PR: E3SM-Project/Omega#485
3 participants
Add the first derivatives of specific volume with respect to conservative temperature, absolute salinity and pressure, for all three equation of state options. These are needed by the higher-order horizontal pressure gradient, which expands specific volume about a reference state in each layer rather than calling the equation of state at every quadrature point. Nothing calls them yet; this adds the capability and its tests.
The GSW toolbox may only be redistributed without modification, so its derivative routine is not ported or adapted. The derivatives are instead obtained by differentiating the Roquet et al. 2015 polynomial that Omega already carries in Eos.h. GSW-C is used unmodified, through its public API, as an independent check in the unit test.
One commit is answer-changing at roundoff level: the thermal expansion and haline contraction coefficients used by the Brunt-Vaisala frequency are now derived from the new derivatives instead of from a second copy of the same coefficients.
What this adds
Eos::computeSpecVolAndDerivs, which fills specific volume and its three first derivatives in a single pass over the equation of state.SpecVolDCt,SpecVolDSaandSpecVolDPmembers and their fields in theEosgroup, owned and allocated byEosexactly asSpecVolis.Teos10Eos,LinearEosandConstantEos, so the higher-order pressure gradient will work with any EOS option.Units: pressure in is relative pressure in Pa; derivatives out are per
degC, per(g/kg), and per Pa (not per dbar).Eosowns and registers the three derivative arrays alongsideSpecVol, so they are available for output like any other EOS field. Nothing computes them yet;PressureGradHighOrderwill be the first consumer.On the GSW-C license
gsw_specvol_first_derivativesinto a Kokkos kernel would be redistribution with modification, so that is not done here. Calling into the library from device code is not an option either.calcPCoeffsandcalcRefProfilealready contain the published Roquet et al. 2015 coefficients. Differentiating that polynomial analytically is elementary calculus on material Omega already has.EosTest, exactly as it already is forgsw_specvolandgsw_ct_freezing_poly.No new numerical constants
Worth knowing before reading the diff, because it makes the review much cheaper.
Teos10BruntVaisalaFreqSq::calcAlpha/calcBetaalready held the differentiated coefficient tables under different names.A***and all 50B***coefficients match exactly in decimal arithmetic, sinceA(i,j,k) = (j+1)·V(i,j+1,k)andB(i,j,k) = (i+1)·V(i+1,j,k).calcPCoeffsalready assembles.Verification
Maximum relative difference against GSW-C over 216 states (salinity 0–38.5 g/kg, temperature −2 to 35 °C, pressure 0–10000 dbar):
SpecVold(SpecVol)/dCtd(SpecVol)/dSad(SpecVol)/dPalpha,betaWhy the pressure derivative is gated more loosely
The initial 1e-12 gate failed for
d(SpecVol)/dPabove 8000 dbar. The error is GSW-C's, not ours, and it is worth stating explicitly so the looser tolerance does not read as a fudge:v_Pis evaluated from coefficients pre-multiplied by their pressure exponents and rounded. Omega differentiates the full-precision ones.The finite-difference check is deliberately redundant with the GSW comparison: it pins the unit convention without reference to GSW, so a pressure derivative that silently became per dbar would fail there too.
The one answer-changing commit
Compute Brunt-Vaisala coefficients from the specific volume derivativesrewritescalcAlpha/calcBetaover the new derivatives and deletes the duplicate tables. Same polynomial, different Horner arrangement, and the normalized salinity is now formed the waycalcPCoeffsforms it, so answers move at roundoff:BruntVaisalaFreqSqis gated at 1e-10 and passes unchangedVERTMIX_TEST, which consumes it, passesA baseline comparison is appropriate rather than a bit-for-bit check.
Note the commit ordering here: the guard test comparing
calcAlpha/calcBetaagainstgsw_specvol_alpha_betais added in the preceding commit, on purpose. Those two functions were previously exercised only through a single hardcodedBruntVaisalaFreqSqvalue at a tolerance too loose to catch a mistake in either, so the rewrite lands against a real oracle rather than against itself.Reading the commits
Each is self-contained and they are meant to be read in order:
Eos::computeSpecVolAndDerivs— the dispatch, plus the three arrays it owns and their field registration.Testing
omega_prsuite from Polaris.Out of scope
Everything else in Phase 1 of the higher-order pressure gradient design: the
PressureGradHighOrderfunctor, the mean-preserving vertical reconstruction, the sloping-interface integral, and theVertCoordgeopotential prerequisite. Second derivatives are Phase 2.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 buildNew tests:
Stealth Features