diff --git a/components/omega/doc/devGuide/EOS.md b/components/omega/doc/devGuide/EOS.md index 4ba3c27e294b..778a71afec45 100644 --- a/components/omega/doc/devGuide/EOS.md +++ b/components/omega/doc/devGuide/EOS.md @@ -70,6 +70,97 @@ volume arrays, do Eos.computeBruntVaisalaFreqSq(ConservTemp, AbsSalinity, Pressure, SpecVol); ``` +## First derivatives of specific volume + +The `Eos` class can also compute the first derivatives of the specific volume +with respect to conservative temperature, absolute salinity, and pressure, +together with the specific volume itself: + +```c++ +Eos.computeSpecVolAndDerivs(ConservTemp, AbsSalinity, Pressure); +``` + +`Pressure` is the relative pressure (gauge pressure in Pa) as elsewhere in +`Eos`, and the derivatives are returned per `degC`, per `(g/kg)`, and per `Pa` +respectively. Note the pressure derivative is per Pascal, not per decibar. + +The results are stored in the `SpecVolDCt`, `SpecVolDSa` and `SpecVolDP` +members alongside `SpecVol`, and all three are registered as fields in the +`Eos` group so they can be written to a stream. Because `SpecVol` is computed +here as well, `computeSpecVolAndDerivs` replaces a call to `computeSpecVol` +rather than accompanying one; calling both would evaluate the equation of state +twice. The valid range of the derivative fields 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. + +The two methods are kept separate rather than always computing the derivatives +because the derivatives roughly double the TEOS-10 arithmetic per cell and +layer, and not every call needs them. `AuxiliaryState::computeMomVertAux` is +the only place that calls `computeSpecVol`; everything else consumes the +`Eos::SpecVol` array rather than recomputing it, including +`computeBruntVaisalaFreqSq`, which takes the specific volume as an argument. +That one call site is reached once per time stepper stage through +`computeMomAux` and `computeAll` in the tendency calculation, and once more +from `VertMix::VertMixImplicit`, which refreshes the pressure and specific +volume before the vertical mixing coefficients are formed. + +A run using the higher-order pressure gradient therefore needs the derivatives +at every time step, and at those call sites `computeSpecVolAndDerivs` takes the +place of the `computeSpecVol` call that would otherwise be made, leaving one +evaluation of the equation of state where there was one before. Two things +still call for the plain `computeSpecVol`. First, `PressureGradType` is a +runtime option that defaults to `Centered`, so a run may never need the +derivatives at all. Second, even with the higher-order pressure gradient +selected, the `VertMix::VertMixImplicit` update feeds only +`computeGeomZHeight` and `computeBruntVaisalaFreqSq`, neither of which reads +the derivatives, so computing them there would be work that nothing consumes. +Which method to call is thus a decision for each call site, not one the `Eos` +class should make for it. + +There is no displaced counterpart to `computeSpecVolAndDerivs`. The pressure +gradient needs the derivatives at the in-situ pressure of the layer, whereas +`computeSpecVolDisp` exists to evaluate the specific volume at the pressure of +a displaced layer. Nothing about the derivatives prevents an adiabatically +displaced version: it would take the same `KDisp` argument as +`computeSpecVolDisp` and evaluate the same coefficients at the displaced +pressure, with no new polynomial. It is left out here only because no caller +needs it yet, and it would mean carrying three more model-sized arrays and +fields. + +All four values come from a single pass over the equation of state. For +`EosType::Teos10Eos` the derivatives are the analytic derivatives of the same +75-term polynomial used for the specific volume, evaluated at the same +normalized state, so no second call to the equation of state is made. The +pressure derivative reuses the pressure coefficients already assembled for the +specific volume; the temperature and salinity derivatives need coefficient sets +of their own but share the normalization and the square root. For +`EosType::LinearEos` the derivatives are `-DRhoDT` and `-DRhoDS` times the +square of the specific volume, with no pressure dependence, and for +`EosType::ConstantEos` all three vanish. + +The thermal expansion and haline contraction coefficients used by the +`BruntVaisalaFreqSq` calculation are formed from these same derivatives, +`alpha = SpecVolDCt / SpecVol` and `beta = -SpecVolDSa / SpecVol`, so the +polynomial coefficients exist in only one place. + +### A note on GSW-C + +The GSW toolbox may be redistributed only without modification, so the +derivative routines in GSW-C are not ported or adapted here; they also could +not be called from a Kokkos device kernel. The implementation instead +differentiates the published Roquet et al. 2015 polynomial that `Teos10Eos` +already carries. GSW-C is used unmodified, through its public API, as an +independent check in the unit test. + +That test compares against `gsw_specvol_first_derivatives` over a range of +states and finds agreement of order `1e-14` for the temperature and salinity +derivatives. The pressure derivative agrees only to about `2e-12`, and the +difference is on the GSW-C side: its `v_P` is evaluated from coefficients that +have been pre-multiplied by their pressure exponents and rounded, whereas the +Omega implementation differentiates the full-precision coefficients and matches +the exact derivative to roughly `1e-16`. + ## Helper functions for conversion The TEOS-10 implementation includes helper functions for temperature diff --git a/components/omega/doc/userGuide/EOS.md b/components/omega/doc/userGuide/EOS.md index df7247cac47e..d7aa1e38df63 100644 --- a/components/omega/doc/userGuide/EOS.md +++ b/components/omega/doc/userGuide/EOS.md @@ -19,7 +19,7 @@ Eos: where `DRhoDT` is the thermal expansion coefficient ($\textrm{kg}/(\textrm{m}^3 \cdot ^{\circ}\textrm{C})$), `DRhoDS` is the saline contraction coefficient ($\textrm{kg}/\textrm{m}^3$), and `RhoT0S0` is the reference density at (T,S)=(0,0) (in $\textrm{kg}/\textrm{m}^3$). -In addition to `SpecVol`, the displaced specific volume `SpecVolDisplaced` and `BruntVaisalaFreqSq` are also calculated by the EOS. +In addition to `SpecVol`, the displaced specific volume `SpecVolDisplaced`, the squared Brunt-Vaisala frequency `BruntVaisalaFreqSq` and the first derivatives of specific volume `SpecVolDCt`, `SpecVolDSa` and `SpecVolDP` are also calculated by the EOS. ## TEOS-10 Helper Conversions @@ -37,6 +37,16 @@ These helper methods are available through the EOS implementation but do not replace the standard `computeSpecVol`, `computeSpecVolDisp`, or `computeBruntVaisalaFreqSq` calculations. +## First Derivatives of Specific Volume + +The `Eos` class can also compute the first derivatives of the specific volume with respect to conservative temperature (in $\textrm{m}^3\textrm{kg}^{-1}\,^{\circ}\textrm{C}^{-1}$), absolute salinity (in $\textrm{m}^3\textrm{g}^{-1}$), and pressure (in $\textrm{m}^3\textrm{kg}^{-1}\textrm{Pa}^{-1}$). These are needed by the higher-order horizontal pressure gradient, which expands the specific volume about a reference state within each layer instead of evaluating the full equation of state at every quadrature point. + +There is no user-configurable option associated with the derivatives. They are computed on request by the parts of the model that need them, in the same pass that computes the specific volume, so selecting `teos10` does not make the model slower unless a scheme that uses them is enabled. They are available for all three `EosType` choices: for `teos10` they are the analytic derivatives of the same 75-term polynomial, for `linear` they follow from the configured `DRhoDT` and `DRhoDS` and have no pressure dependence, and for `constant` they are zero. + +The derivatives are stored in the `SpecVolDCt`, `SpecVolDSa` and `SpecVolDP` fields of the `Eos` field group and can be requested in a stream's contents just like `SpecVol`. + +The thermal expansion and haline contraction coefficients that enter the squared Brunt-Vaisala frequency are computed from these same derivatives. + ## Displaced Specific Volume The `Eos` class calculates the density of a parcel of fluid that is adiabatically displaced by a relative `k` levels (`k` counted positive downward), capturing the effects of pressure/depth changes. This is primarily used to calculate quantities for determining the water column stability (i.e. the stratification) and the vertical mixing coefficients (viscosity and diffusivity). Note: when using the `Linear` or `constant` EOS option, `SpecVolDisplaced` will be the same as `SpecVol` since the specific volume calculation is independent of pressure/depth. diff --git a/components/omega/src/ocn/Eos.cpp b/components/omega/src/ocn/Eos.cpp index 16ea4da9f5a2..8c2d8721181a 100644 --- a/components/omega/src/ocn/Eos.cpp +++ b/components/omega/src/ocn/Eos.cpp @@ -48,6 +48,11 @@ Eos::Eos(const std::string &Name, ///< [in] Name for eos object Array2DReal("SpecVolDisplaced", Mesh->NCellsSize, VCoord->NVertLayers); BruntVaisalaFreqSq = Array2DReal("BruntVaisalaFreqSq", Mesh->NCellsSize, VCoord->NVertLayersP1); + SpecVolDCt = + Array2DReal("SpecVolDCt", Mesh->NCellsSize, VCoord->NVertLayers); + SpecVolDSa = + Array2DReal("SpecVolDSa", Mesh->NCellsSize, VCoord->NVertLayers); + SpecVolDP = Array2DReal("SpecVolDP", Mesh->NCellsSize, VCoord->NVertLayers); defineFields(); } @@ -136,8 +141,6 @@ void Eos::computeSpecVol(const Array2DReal &ConservTemp, ComputeSpecVolTeos10); /// Local view for TEOS-10 computation OMEGA_SCOPE(LocComputeSpecVolConstant, ComputeSpecVolConstant); /// Local view for constant computation - OMEGA_SCOPE(MinLayerCell, VCoord->MinLayerCell); - OMEGA_SCOPE(MaxLayerCell, VCoord->MaxLayerCell); I4 KDisp = 0; /// No displacement in this case @@ -146,44 +149,22 @@ void Eos::computeSpecVol(const Array2DReal &ConservTemp, parallelForOuter( "eos-linear", {Mesh->NCellsAll}, KOKKOS_LAMBDA(I4 ICell, const TeamMember &Team) { - const int KMin = MinLayerCell(ICell); - const int KMax = MaxLayerCell(ICell); - const int KRange = vertRangeChunked(KMin, KMax); - - parallelForInner( - Team, KRange, INNER_LAMBDA(int KChunk) { - LocComputeSpecVolLinear(LocSpecVol, ICell, KChunk, - ConservTemp, AbsSalinity); - }); + LocComputeSpecVolLinear(LocSpecVol, Team, ICell, ConservTemp, + AbsSalinity); }); } else if (EosChoice == EosType::Teos10Eos) { parallelForOuter( "eos-teos10", {Mesh->NCellsAll}, KOKKOS_LAMBDA(I4 ICell, const TeamMember &Team) { - const int KMin = MinLayerCell(ICell); - const int KMax = MaxLayerCell(ICell); - const int KRange = vertRangeChunked(KMin, KMax); - - parallelForInner( - Team, KRange, INNER_LAMBDA(int KChunk) { - LocComputeSpecVolTeos10(LocSpecVol, ICell, KChunk, - ConservTemp, AbsSalinity, Pressure, - KDisp); - }); + LocComputeSpecVolTeos10(LocSpecVol, Team, ICell, ConservTemp, + AbsSalinity, Pressure, KDisp); }); } else if (EosChoice == EosType::ConstantEos) { parallelForOuter( "eos-constant", {Mesh->NCellsAll}, KOKKOS_LAMBDA(I4 ICell, const TeamMember &Team) { - const int KMin = MinLayerCell(ICell); - const int KMax = MaxLayerCell(ICell); - const int KRange = vertRangeChunked(KMin, KMax); - - parallelForInner( - Team, KRange, INNER_LAMBDA(int KChunk) { - LocComputeSpecVolConstant(LocSpecVol, ICell, KChunk, - ConservTemp, AbsSalinity); - }); + LocComputeSpecVolConstant(LocSpecVol, Team, ICell, ConservTemp, + AbsSalinity); }); } } @@ -200,8 +181,6 @@ void Eos::computeSpecVolDisp(const Array2DReal &ConservTemp, ComputeSpecVolTeos10); /// Local view for TEOS-10 computation OMEGA_SCOPE(LocComputeSpecVolConstant, ComputeSpecVolConstant); /// Local view for constant computation - OMEGA_SCOPE(MinLayerCell, VCoord->MinLayerCell); - OMEGA_SCOPE(MaxLayerCell, VCoord->MaxLayerCell); /// Dispatch to the correct EOS calculation /// If EosChoice is Linear, the displaced specific @@ -210,41 +189,65 @@ void Eos::computeSpecVolDisp(const Array2DReal &ConservTemp, parallelForOuter( "eos-linear", {Mesh->NCellsAll}, KOKKOS_LAMBDA(I4 ICell, const TeamMember &Team) { - const int KMin = MinLayerCell(ICell); - const int KMax = MaxLayerCell(ICell); - const int KRange = vertRangeChunked(KMin, KMax); - parallelForInner( - Team, KRange, INNER_LAMBDA(int KChunk) { - LocComputeSpecVolLinear(LocSpecVolDisplaced, ICell, KChunk, - ConservTemp, AbsSalinity); - }); + LocComputeSpecVolLinear(LocSpecVolDisplaced, Team, ICell, + ConservTemp, AbsSalinity); }); } else if (EosChoice == EosType::Teos10Eos) { parallelForOuter( "eos-teos10", {Mesh->NCellsAll}, KOKKOS_LAMBDA(I4 ICell, const TeamMember &Team) { - const int KMin = MinLayerCell(ICell); - const int KMax = MaxLayerCell(ICell); - const int KRange = vertRangeChunked(KMin, KMax); - parallelForInner( - Team, KRange, INNER_LAMBDA(int KChunk) { - LocComputeSpecVolTeos10(LocSpecVolDisplaced, ICell, KChunk, - ConservTemp, AbsSalinity, Pressure, - KDisp); - }); + LocComputeSpecVolTeos10(LocSpecVolDisplaced, Team, ICell, + ConservTemp, AbsSalinity, Pressure, KDisp); }); } else if (EosChoice == EosType::ConstantEos) { parallelForOuter( "eos-constant", {Mesh->NCellsAll}, KOKKOS_LAMBDA(I4 ICell, const TeamMember &Team) { - const int KMin = MinLayerCell(ICell); - const int KMax = MaxLayerCell(ICell); - const int KRange = vertRangeChunked(KMin, KMax); - parallelForInner( - Team, KRange, INNER_LAMBDA(int KChunk) { - LocComputeSpecVolConstant(LocSpecVolDisplaced, ICell, - KChunk, ConservTemp, AbsSalinity); - }); + LocComputeSpecVolConstant(LocSpecVolDisplaced, Team, ICell, + ConservTemp, AbsSalinity); + }); + } +} + +/// Compute specific volume and its first derivatives for all cells/layers +void Eos::computeSpecVolAndDerivs(const Array2DReal &ConservTemp, + const Array2DReal &AbsSalinity, + const Array2DReal &Pressure) { + OMEGA_SCOPE(LocSpecVol, SpecVol); /// Local views for computation + OMEGA_SCOPE(LocSpecVolDCt, SpecVolDCt); /// Temperature derivative + OMEGA_SCOPE(LocSpecVolDSa, SpecVolDSa); /// Salinity derivative + OMEGA_SCOPE(LocSpecVolDP, SpecVolDP); /// Pressure derivative + OMEGA_SCOPE(LocComputeSpecVolLinear, + ComputeSpecVolLinear); /// Local view for linear EOS computation + OMEGA_SCOPE(LocComputeSpecVolTeos10, + ComputeSpecVolTeos10); /// Local view for TEOS-10 computation + OMEGA_SCOPE(LocComputeSpecVolConstant, + ComputeSpecVolConstant); /// Local view for constant computation + + /// Dispatch to the correct EOS calculation + if (EosChoice == EosType::LinearEos) { + parallelForOuter( + "eos-derivs-linear", {Mesh->NCellsAll}, + KOKKOS_LAMBDA(I4 ICell, const TeamMember &Team) { + LocComputeSpecVolLinear.calcSpecVolAndDerivsOnCells( + LocSpecVol, LocSpecVolDCt, LocSpecVolDSa, LocSpecVolDP, Team, + ICell, ConservTemp, AbsSalinity); + }); + } else if (EosChoice == EosType::Teos10Eos) { + parallelForOuter( + "eos-derivs-teos10", {Mesh->NCellsAll}, + KOKKOS_LAMBDA(I4 ICell, const TeamMember &Team) { + LocComputeSpecVolTeos10.calcSpecVolAndDerivsOnCells( + LocSpecVol, LocSpecVolDCt, LocSpecVolDSa, LocSpecVolDP, Team, + ICell, ConservTemp, AbsSalinity, Pressure); + }); + } else if (EosChoice == EosType::ConstantEos) { + parallelForOuter( + "eos-derivs-constant", {Mesh->NCellsAll}, + KOKKOS_LAMBDA(I4 ICell, const TeamMember &Team) { + LocComputeSpecVolConstant.calcSpecVolAndDerivsOnCells( + LocSpecVol, LocSpecVolDCt, LocSpecVolDSa, LocSpecVolDP, Team, + ICell, ConservTemp, AbsSalinity); }); } } @@ -272,14 +275,11 @@ void Eos::computeBruntVaisalaFreqSq(const Array2DReal &ConservTemp, "bvf-linear", {Mesh->NCellsAll}, KOKKOS_LAMBDA(I4 ICell, const TeamMember &Team) { // Compute Brunt-Vaisala frequency at interior vertical interfaces - const int KMin = MinLayerCell(ICell) + 1; - const int KMax = MaxLayerCell(ICell); - const int KRange = vertRangeChunked(KMin, KMax); - parallelForInner( - Team, KRange, INNER_LAMBDA(int KChunk) { - LocComputeBruntVaisalaFreqSqLinear(LocBruntVaisalaFreqSq, - ICell, KChunk, SpecVol); - }); + const int KMin = MinLayerCell(ICell) + 1; + const int KMax = MaxLayerCell(ICell); + + LocComputeBruntVaisalaFreqSqLinear(LocBruntVaisalaFreqSq, Team, + ICell, SpecVol); teamBarrier(Team); @@ -301,15 +301,12 @@ void Eos::computeBruntVaisalaFreqSq(const Array2DReal &ConservTemp, "bvf-teos10", {Mesh->NCellsAll}, KOKKOS_LAMBDA(I4 ICell, const TeamMember &Team) { // Compute Brunt-Vaisala frequency at interior vertical interfaces - const int KMin = MinLayerCell(ICell) + 1; - const int KMax = MaxLayerCell(ICell); - const int KRange = vertRangeChunked(KMin, KMax); - parallelForInner( - Team, KRange, INNER_LAMBDA(int KChunk) { - LocComputeBruntVaisalaFreqSqTeos10( - LocBruntVaisalaFreqSq, ICell, KChunk, ConservTemp, - AbsSalinity, Pressure, SpecVol); - }); + const int KMin = MinLayerCell(ICell) + 1; + const int KMax = MaxLayerCell(ICell); + + LocComputeBruntVaisalaFreqSqTeos10(LocBruntVaisalaFreqSq, Team, + ICell, ConservTemp, AbsSalinity, + Pressure, SpecVol); teamBarrier(Team); @@ -350,10 +347,16 @@ void Eos::defineFields() { SpecVolFldName = "SpecVol"; SpecVolDisplacedFldName = "SpecVolDisplaced"; BruntVaisalaFreqSqFldName = "BruntVaisalaFreqSq"; + SpecVolDCtFldName = "SpecVolDCt"; + SpecVolDSaFldName = "SpecVolDSa"; + SpecVolDPFldName = "SpecVolDP"; if (Name != "Default") { SpecVolFldName.append(Name); SpecVolDisplacedFldName.append(Name); BruntVaisalaFreqSqFldName.append(Name); + SpecVolDCtFldName.append(Name); + SpecVolDSaFldName.append(Name); + SpecVolDPFldName.append(Name); } /// Create fields for state variables @@ -386,6 +389,45 @@ void Eos::defineFields() { DimNames // Dimension names ); + /// The specific volume derivatives are legitimately negative, so their + /// valid range spans the full range of Real rather than starting at zero + auto SpecVolDCtField = Field::create( + SpecVolDCtFldName, // Field name + "Derivative of specific volume with respect to conservative " + "temperature", // Long Name + "m3 kg-1 degC-1", // Units + // CF-ish Name + "sea_water_specific_volume_derivative_wrt_conservative_temperature", + std::numeric_limits::lowest(), // Min valid value + std::numeric_limits::max(), // Max valid value + NDims, // Number of dimensions + DimNames // Dimension names + ); + + auto SpecVolDSaField = Field::create( + SpecVolDSaFldName, // Field name + "Derivative of specific volume with respect to absolute " + "salinity", // Long Name + "m3 g-1", // Units + // CF-ish Name + "sea_water_specific_volume_derivative_wrt_absolute_salinity", + std::numeric_limits::lowest(), // Min valid value + std::numeric_limits::max(), // Max valid value + NDims, // Number of dimensions + DimNames // Dimension names + ); + + auto SpecVolDPField = Field::create( + SpecVolDPFldName, // Field name + "Derivative of specific volume with respect to pressure", // Long Name + "m3 kg-1 Pa-1", // Units + "sea_water_specific_volume_derivative_wrt_pressure", // CF-ish Name + std::numeric_limits::lowest(), // Min valid value + std::numeric_limits::max(), // Max valid value + NDims, // Num dimensions + DimNames // Dimension names + ); + // Brunt-Vaisala frequency is located at interfaces DimNames[1] = "NVertLayersP1"; @@ -412,11 +454,17 @@ void Eos::defineFields() { EosGroup->addField(SpecVolDisplacedFldName); EosGroup->addField(SpecVolFldName); EosGroup->addField(BruntVaisalaFreqSqFldName); + EosGroup->addField(SpecVolDCtFldName); + EosGroup->addField(SpecVolDSaFldName); + EosGroup->addField(SpecVolDPFldName); // Attach Kokkos views to the fields SpecVolDisplacedField->attachData(SpecVolDisplaced); SpecVolField->attachData(SpecVol); BruntVaisalaFreqSqField->attachData(BruntVaisalaFreqSq); + SpecVolDCtField->attachData(SpecVolDCt); + SpecVolDSaField->attachData(SpecVolDSa); + SpecVolDPField->attachData(SpecVolDP); } // end defineIOFields diff --git a/components/omega/src/ocn/Eos.h b/components/omega/src/ocn/Eos.h index 5e4fd89cae21..462103c2c3f2 100644 --- a/components/omega/src/ocn/Eos.h +++ b/components/omega/src/ocn/Eos.h @@ -34,60 +34,68 @@ class Teos10Eos { /// constructor declaration Teos10Eos(const VertCoord *VCoord); + /// Normalization used by the Roquet et al. 2015 75-term polynomial. The + /// polynomial is written in the normalized variables + /// Ss = sqrt((Sa + DeltaS) / SaNorm), Tt = Ct / CtNorm, Pp = P * PNorm + /// with P the relative pressure in dbar. These are shared by the specific + /// volume and by its derivatives, so that both are evaluated at exactly + /// the same normalized state. + static constexpr Real SaNorm = 40.0 * SS0 / 35.0; ///< salinity scale (g/kg) + static constexpr Real CtNorm = 40.0; ///< temperature scale (degC) + static constexpr Real DeltaS = 24.0; ///< salinity offset (g/kg) + static constexpr Real PNorm = 1.0e-4; ///< pressure scale (1/dbar) + // The functor takes the full arrays of specific volume (inout), - // the indices ICell and KChunk, and the ocean tracers (conservative) - // temperature, (absolute) salinity, and relative pressure (gauge pressure, - // i.e., absolute pressure minus the standard atmosphere) as inputs, and - // outputs the specific volume according to the Roquet et al. 2015 75 term - // expansion. - KOKKOS_FUNCTION void operator()(Array2DReal SpecVol, I4 ICell, I4 KChunk, - const Array2DReal &ConservTemp, + // the team member and the cell index ICell, and the ocean tracers + // (conservative) temperature, (absolute) salinity, and relative pressure + // (gauge pressure, i.e., absolute pressure minus the standard atmosphere) + // as inputs, and outputs the specific volume according to the Roquet et + // al. 2015 75 term expansion. + KOKKOS_FUNCTION void operator()(Array2DReal SpecVol, const TeamMember &Team, + I4 ICell, const Array2DReal &ConservTemp, const Array2DReal &AbsSalinity, const Array2DReal &Pressure, I4 KDisp) const { - Real SpecVolPCoeffs[6 * VecLength]; - const I4 KStart = chunkStart(KChunk, MinLayerCell(ICell)); - const I4 KLen = chunkLength(KChunk, KStart, MaxLayerCell(ICell)); - - for (int KVec = 0; KVec < KLen; ++KVec) { - const I4 K = KStart + KVec; - /// Calculate the local specific volume polynomial pressure - /// coefficients with cell center values - calcPCoeffs(SpecVolPCoeffs, KVec, ConservTemp(ICell, K), - AbsSalinity(ICell, K)); - - /// Calculate the specific volume at the given pressure - /// If KDisp is 0, we use the current pressure, otherwise - /// we use the displaced pressure (K + KDisp) - /// Note: KDisp is only used for TEOS-10, for Linear EOS it - /// is always 0. - if (KDisp == 0) { - // No displacement - SpecVol(ICell, K) = - calcRefProfile(Pressure(ICell, K) * Pa2Db) + - calcDelta(SpecVolPCoeffs, KVec, Pressure(ICell, K) * Pa2Db); - } else { - // Displacement, use the displaced pressure - I4 KTmp = Kokkos::min(K + KDisp, MaxLayerCell(ICell)); - KTmp = Kokkos::max(MinLayerCell(ICell), KTmp); - SpecVol(ICell, K) = - calcRefProfile(Pressure(ICell, KTmp) * Pa2Db) + - calcDelta(SpecVolPCoeffs, KVec, Pressure(ICell, KTmp) * Pa2Db); - } - } + const I4 KMin = MinLayerCell(ICell); + const I4 KMax = MaxLayerCell(ICell); + + parallelForInner( + Team, Range{KMin, KMax}, INNER_LAMBDA(int K) { + Real SpecVolPCoeffs[6]; + + /// Calculate the local specific volume polynomial pressure + /// coefficients with cell center values + calcPCoeffs(SpecVolPCoeffs, ConservTemp(ICell, K), + AbsSalinity(ICell, K)); + + /// Calculate the specific volume at the given pressure + /// If KDisp is 0, we use the current pressure, otherwise + /// we use the displaced pressure (K + KDisp) + /// Note: KDisp is only used for TEOS-10, for Linear EOS it + /// is always 0. + if (KDisp == 0) { + // No displacement + SpecVol(ICell, K) = + calcRefProfile(Pressure(ICell, K) * Pa2Db) + + calcDelta(SpecVolPCoeffs, Pressure(ICell, K) * Pa2Db); + } else { + // Displacement, use the displaced pressure + I4 KTmp = Kokkos::min(K + KDisp, KMax); + KTmp = Kokkos::max(KMin, KTmp); + SpecVol(ICell, K) = + calcRefProfile(Pressure(ICell, KTmp) * Pa2Db) + + calcDelta(SpecVolPCoeffs, Pressure(ICell, KTmp) * Pa2Db); + } + }); } /// TEOS-10 helpers /// Calculate pressure polynomial coefficients for TEOS-10 - KOKKOS_FUNCTION void calcPCoeffs(Real (&SpecVolPCoeffs)[6 * VecLength], - const I4 KVec, const Real Ct, + KOKKOS_FUNCTION void calcPCoeffs(Real (&SpecVolPCoeffs)[6], const Real Ct, const Real Sa) const { - constexpr Real SaNorm = 40.0 * 35.16504 / 35.0; - constexpr Real CtNorm = 40.0; - constexpr Real DeltaS = 24.0; - Real Ss = Kokkos::sqrt((Sa + DeltaS) / SaNorm); - Real Tt = Ct / CtNorm; + Real Ss = Kokkos::sqrt((Sa + DeltaS) / SaNorm); + Real Tt = Ct / CtNorm; /// Coefficients for the polynomial expansion constexpr Real V000 = 1.0769995862e-03; @@ -165,18 +173,18 @@ class Teos10Eos { constexpr Real V014 = 3.1454099902e-07; constexpr Real V005 = 4.2369007180e-09; - SpecVolPCoeffs[5 + 6 * KVec] = V005; - SpecVolPCoeffs[4 + 6 * KVec] = V014 * Tt + V104 * Ss + V004; - SpecVolPCoeffs[3 + 6 * KVec] = + SpecVolPCoeffs[5] = V005; + SpecVolPCoeffs[4] = V014 * Tt + V104 * Ss + V004; + SpecVolPCoeffs[3] = (V023 * Tt + V113 * Ss + V013) * Tt + (V203 * Ss + V103) * Ss + V003; - SpecVolPCoeffs[2 + 6 * KVec] = - (((V042 * Tt + V132 * Ss + V032) * Tt + (V222 * Ss + V122) * Ss + - V022) * - Tt + - ((V312 * Ss + V212) * Ss + V112) * Ss + V012) * - Tt + - (((V402 * Ss + V302) * Ss + V202) * Ss + V102) * Ss + V002; - SpecVolPCoeffs[1 + 6 * KVec] = + SpecVolPCoeffs[2] = (((V042 * Tt + V132 * Ss + V032) * Tt + + (V222 * Ss + V122) * Ss + V022) * + Tt + + ((V312 * Ss + V212) * Ss + V112) * Ss + V012) * + Tt + + (((V402 * Ss + V302) * Ss + V202) * Ss + V102) * Ss + + V002; + SpecVolPCoeffs[1] = ((((V051 * Tt + V141 * Ss + V041) * Tt + (V231 * Ss + V131) * Ss + V031) * Tt + @@ -186,7 +194,7 @@ class Teos10Eos { Tt + ((((V501 * Ss + V401) * Ss + V301) * Ss + V201) * Ss + V101) * Ss + V001; - SpecVolPCoeffs[0 + 6 * KVec] = + SpecVolPCoeffs[0] = (((((V060 * Tt + V150 * Ss + V050) * Tt + (V240 * Ss + V140) * Ss + V040) * Tt + @@ -204,36 +212,32 @@ class Teos10Eos { } /// Evaluate pressure polynomial delta for TEOS-10 - KOKKOS_FUNCTION Real calcDelta(const Real (&SpecVolPCoeffs)[6 * VecLength], - const I4 KVec, const Real P) const { + KOKKOS_FUNCTION Real calcDelta(const Real (&SpecVolPCoeffs)[6], + const Real P) const { - constexpr Real PNorm = 1e-4; - Real Pp = P * PNorm; + Real Pp = P * PNorm; - Real Delta = ((((SpecVolPCoeffs[5 + 6 * KVec] * Pp + - SpecVolPCoeffs[4 + 6 * KVec]) * - Pp + - SpecVolPCoeffs[3 + 6 * KVec]) * + Real Delta = ((((SpecVolPCoeffs[5] * Pp + SpecVolPCoeffs[4]) * Pp + + SpecVolPCoeffs[3]) * Pp + - SpecVolPCoeffs[2 + 6 * KVec]) * + SpecVolPCoeffs[2]) * Pp + - SpecVolPCoeffs[1 + 6 * KVec]) * + SpecVolPCoeffs[1]) * Pp + - SpecVolPCoeffs[0 + 6 * KVec]; + SpecVolPCoeffs[0]; return Delta; } /// Calculate reference profile for TEOS-10 KOKKOS_FUNCTION Real calcRefProfile(Real P) const { - constexpr Real PNorm = 1e-4; - constexpr Real V00 = -4.4015007269e-05; - constexpr Real V01 = 6.9232335784e-06; - constexpr Real V02 = -7.5004675975e-07; - constexpr Real V03 = 1.7009109288e-08; - constexpr Real V04 = -1.6884162004e-08; - constexpr Real V05 = 1.9613503930e-09; - Real Pp = P * PNorm; + constexpr Real V00 = -4.4015007269e-05; + constexpr Real V01 = 6.9232335784e-06; + constexpr Real V02 = -7.5004675975e-07; + constexpr Real V03 = 1.7009109288e-08; + constexpr Real V04 = -1.6884162004e-08; + constexpr Real V05 = 1.9613503930e-09; + Real Pp = P * PNorm; Real V0 = (((((V05 * Pp + V04) * Pp + V03) * Pp + V02) * Pp + V01) * Pp + V00) * @@ -241,6 +245,307 @@ class Teos10Eos { return V0; } + /// Calculate pressure polynomial coefficients for the derivative of the + /// TEOS-10 specific volume with respect to the normalized temperature Tt. + /// + /// The coefficients are the analytic Tt-derivative of the 75-term + /// polynomial in calcPCoeffs above: A(i,j,k) = (j+1) * V(i,j+1,k), grouped + /// here by power of pressure to match the calcPCoeffs/calcDelta split. The + /// Tt-derivative is one degree lower in pressure than the specific volume + /// itself, so there are five coefficients rather than six. + /// + /// This is the only copy of these coefficients: the thermal expansion + /// coefficient used by the Brunt-Vaisala frequency is derived from here. + KOKKOS_FUNCTION static void calcPCoeffsDTt(Real (&DTtPCoeffs)[5], + const Real Ss, const Real Tt) { + + constexpr Real A000 = -1.56497346750e-5; + constexpr Real A001 = 1.85057654290e-5; + constexpr Real A002 = -1.17363867310e-6; + constexpr Real A003 = -3.65270065530e-7; + constexpr Real A004 = 3.14540999020e-7; + constexpr Real A010 = 5.55242129680e-5; + constexpr Real A011 = -2.34332137060e-5; + constexpr Real A012 = 4.26100574800e-6; + constexpr Real A013 = 5.73918103180e-7; + constexpr Real A020 = -4.95634777770e-5; + constexpr Real A021 = 2.37838968519e-5; + constexpr Real A022 = -1.38397620111e-6; + constexpr Real A030 = 2.76445290808e-5; + constexpr Real A031 = -1.36408749928e-5; + constexpr Real A032 = -2.53411666056e-7; + constexpr Real A040 = -4.02698077700e-6; + constexpr Real A041 = 2.53683834070e-6; + constexpr Real A050 = 1.23258565608e-6; + constexpr Real A100 = 3.50095997640e-5; + constexpr Real A101 = -9.56770881560e-6; + constexpr Real A102 = -5.56991545570e-6; + constexpr Real A103 = -2.72956962370e-7; + constexpr Real A110 = -7.48716846880e-5; + constexpr Real A111 = -4.73566167220e-7; + constexpr Real A112 = 7.82747741600e-7; + constexpr Real A120 = 7.24244384490e-5; + constexpr Real A121 = -1.03676320965e-5; + constexpr Real A122 = 2.32856664276e-8; + constexpr Real A130 = -3.50383492616e-5; + constexpr Real A131 = 5.18268711320e-6; + constexpr Real A140 = -1.65263794500e-6; + constexpr Real A200 = -4.35926785610e-5; + constexpr Real A201 = 1.11008347650e-5; + constexpr Real A202 = 5.46207488340e-6; + constexpr Real A210 = 7.18156455200e-5; + constexpr Real A211 = 5.85666925900e-6; + constexpr Real A212 = -1.31462208134e-6; + constexpr Real A220 = -4.30608991440e-5; + constexpr Real A221 = 9.49659182340e-7; + constexpr Real A230 = 1.74814722392e-5; + constexpr Real A300 = 3.45324618280e-5; + constexpr Real A301 = -9.84471178440e-6; + constexpr Real A302 = -1.35441856270e-6; + constexpr Real A310 = -3.73971683740e-5; + constexpr Real A311 = -9.76522784000e-7; + constexpr Real A320 = 6.85899736680e-6; + constexpr Real A400 = -1.19594097880e-5; + constexpr Real A401 = 2.59092252600e-6; + constexpr Real A410 = 7.71906784880e-6; + constexpr Real A500 = 1.38645945810e-6; + + DTtPCoeffs[4] = A004; + DTtPCoeffs[3] = A013 * Tt + A103 * Ss + A003; + DTtPCoeffs[2] = ((A032 * Tt + A122 * Ss + A022) * Tt + + (A212 * Ss + A112) * Ss + A012) * + Tt + + ((A302 * Ss + A202) * Ss + A102) * Ss + A002; + DTtPCoeffs[1] = (((A041 * Tt + A131 * Ss + A031) * Tt + + (A221 * Ss + A121) * Ss + A021) * + Tt + + ((A311 * Ss + A211) * Ss + A111) * Ss + A011) * + Tt + + (((A401 * Ss + A301) * Ss + A201) * Ss + A101) * Ss + + A001; + DTtPCoeffs[0] = + ((((A050 * Tt + A140 * Ss + A040) * Tt + (A230 * Ss + A130) * Ss + + A030) * + Tt + + ((A320 * Ss + A220) * Ss + A120) * Ss + A020) * + Tt + + (((A410 * Ss + A310) * Ss + A210) * Ss + A110) * Ss + A010) * + Tt + + ((((A500 * Ss + A400) * Ss + A300) * Ss + A200) * Ss + A100) * Ss + + A000; + } + + /// Calculate pressure polynomial coefficients for the derivative of the + /// TEOS-10 specific volume with respect to the normalized salinity Ss. + /// + /// As above, these are the analytic Ss-derivative of the 75-term + /// polynomial: B(i,j,k) = (i+1) * V(i+1,j,k), and this is the only copy of + /// them; the haline contraction coefficient is derived from here. + KOKKOS_FUNCTION static void calcPCoeffsDSs(Real (&DSsPCoeffs)[5], + const Real Ss, const Real Tt) { + + constexpr Real B000 = -3.10389819760e-4; + constexpr Real B001 = 2.42624687470e-5; + constexpr Real B002 = -5.84844329840e-7; + constexpr Real B003 = 3.63101885150e-7; + constexpr Real B004 = -1.11471254230e-7; + constexpr Real B010 = 3.50095997640e-5; + constexpr Real B011 = -9.56770881560e-6; + constexpr Real B012 = -5.56991545570e-6; + constexpr Real B013 = -2.72956962370e-7; + constexpr Real B020 = -3.74358423440e-5; + constexpr Real B021 = -2.36783083610e-7; + constexpr Real B022 = 3.91373870800e-7; + constexpr Real B030 = 2.41414794830e-5; + constexpr Real B031 = -3.45587736550e-6; + constexpr Real B032 = 7.76188880920e-9; + constexpr Real B040 = -8.75958731540e-6; + constexpr Real B041 = 1.29567177830e-6; + constexpr Real B050 = -3.30527589000e-7; + constexpr Real B100 = 1.33856134076e-3; + constexpr Real B101 = -6.95849219480e-5; + constexpr Real B102 = -9.62445031940e-6; + constexpr Real B103 = 3.34926075600e-8; + constexpr Real B110 = -8.71853571220e-5; + constexpr Real B111 = 2.22016695300e-5; + constexpr Real B112 = 1.09241497668e-5; + constexpr Real B120 = 7.18156455200e-5; + constexpr Real B121 = 5.85666925900e-6; + constexpr Real B122 = -1.31462208134e-6; + constexpr Real B130 = -2.87072660960e-5; + constexpr Real B131 = 6.33106121560e-7; + constexpr Real B140 = 8.74073611960e-6; + constexpr Real B200 = -2.55143801811e-3; + constexpr Real B201 = 1.12412331915e-4; + constexpr Real B202 = 1.47789320994e-5; + constexpr Real B210 = 1.03597385484e-4; + constexpr Real B211 = -2.95341353532e-5; + constexpr Real B212 = -4.06325568810e-6; + constexpr Real B220 = -5.60957525610e-5; + constexpr Real B221 = -1.46478417600e-6; + constexpr Real B230 = 6.85899736680e-6; + constexpr Real B300 = 2.32344279772e-3; + constexpr Real B301 = -6.92888744480e-5; + constexpr Real B302 = -7.12478989080e-6; + constexpr Real B310 = -4.78376391520e-5; + constexpr Real B311 = 1.03636901040e-5; + constexpr Real B320 = 1.54381356976e-5; + constexpr Real B400 = -1.05461852535e-3; + constexpr Real B401 = 1.54637136265e-5; + constexpr Real B410 = 6.93229729050e-6; + constexpr Real B500 = 1.91594743830e-4; + + DSsPCoeffs[4] = B004; + DSsPCoeffs[3] = B013 * Tt + B103 * Ss + B003; + DSsPCoeffs[2] = ((B032 * Tt + B122 * Ss + B022) * Tt + + (B212 * Ss + B112) * Ss + B012) * + Tt + + ((B302 * Ss + B202) * Ss + B102) * Ss + B002; + DSsPCoeffs[1] = (((B041 * Tt + B131 * Ss + B031) * Tt + + (B221 * Ss + B121) * Ss + B021) * + Tt + + ((B311 * Ss + B211) * Ss + B111) * Ss + B011) * + Tt + + (((B401 * Ss + B301) * Ss + B201) * Ss + B101) * Ss + + B001; + DSsPCoeffs[0] = + ((((B050 * Tt + B140 * Ss + B040) * Tt + (B230 * Ss + B130) * Ss + + B030) * + Tt + + ((B320 * Ss + B220) * Ss + B120) * Ss + B020) * + Tt + + (((B410 * Ss + B310) * Ss + B210) * Ss + B110) * Ss + B010) * + Tt + + ((((B500 * Ss + B400) * Ss + B300) * Ss + B200) * Ss + B100) * Ss + + B000; + } + + /// Evaluate one of the degree-4 derivative pressure polynomials assembled + /// by calcPCoeffsDTt or calcPCoeffsDSs. P is relative pressure in dbar. + KOKKOS_FUNCTION static Real calcDeltaDeriv(const Real (&DPCoeffs)[5], + const Real P) { + + Real Pp = P * PNorm; + + Real DDelta = + (((DPCoeffs[4] * Pp + DPCoeffs[3]) * Pp + DPCoeffs[2]) * Pp + + DPCoeffs[1]) * + Pp + + DPCoeffs[0]; + + return DDelta; + } + + /// Evaluate the derivative of the TEOS-10 pressure polynomial with respect + /// to pressure, from the coefficients calcPCoeffs has already assembled for + /// the specific volume itself. P is relative pressure in dbar and the + /// result is per dbar. + KOKKOS_FUNCTION Real calcDeltaDP(const Real (&SpecVolPCoeffs)[6], + const Real P) const { + + Real Pp = P * PNorm; + + Real DDelta = + (((5.0_Real * SpecVolPCoeffs[5] * Pp + 4.0_Real * SpecVolPCoeffs[4]) * + Pp + + 3.0_Real * SpecVolPCoeffs[3]) * + Pp + + 2.0_Real * SpecVolPCoeffs[2]) * + Pp + + SpecVolPCoeffs[1]; + + return DDelta * PNorm; + } + + /// Calculate the derivative of the TEOS-10 reference profile with respect + /// to pressure. P is relative pressure in dbar and the result is per dbar. + /// The reference profile does not depend on temperature or salinity, so it + /// contributes only to the pressure derivative -- but it must not be left + /// out of that one. + KOKKOS_FUNCTION Real calcRefProfileDP(Real P) const { + constexpr Real V00 = -4.4015007269e-05; + constexpr Real V01 = 6.9232335784e-06; + constexpr Real V02 = -7.5004675975e-07; + constexpr Real V03 = 1.7009109288e-08; + constexpr Real V04 = -1.6884162004e-08; + constexpr Real V05 = 1.9613503930e-09; + Real Pp = P * PNorm; + + Real DV0 = + (((((6.0_Real * V05 * Pp + 5.0_Real * V04) * Pp + 4.0_Real * V03) * + Pp + + 3.0_Real * V02) * + Pp + + 2.0_Real * V01) * + Pp + + V00); + + return DV0 * PNorm; + } + + /// Calculate the TEOS-10 specific volume and its three first derivatives at + /// a single state. P is the relative pressure (gauge pressure in Pa, i.e. + /// absolute pressure minus the standard atmosphere). The derivatives are + /// returned per degC, per (g/kg), and per Pa respectively. + /// + /// This is the point-wise entry point: it takes scalars and returns scalars, + /// with no cell or layer indexing. It is what the unit tests call to check + /// the polynomial against GSW-C and against finite differences at chosen + /// states, and it is the form to use anywhere a single state needs to be + /// evaluated. It is also the single implementation of the derivatives: the + /// array-level path over the mesh, calcSpecVolAndDerivsOnCells below, calls + /// this at each cell and layer. + KOKKOS_FUNCTION void + calcSpecVolAndDerivsAtPoint(const Real Ct, const Real Sa, const Real P, + Real &SpecVol, Real &SpecVolDCt, + Real &SpecVolDSa, Real &SpecVolDP) const { + + Real SpecVolPCoeffs[6]; + Real DTtPCoeffs[5]; + Real DSsPCoeffs[5]; + + const Real Ss = Kokkos::sqrt((Sa + DeltaS) / SaNorm); + const Real Tt = Ct / CtNorm; + const Real Pdb = P * Pa2Db; + + calcPCoeffs(SpecVolPCoeffs, Ct, Sa); + calcPCoeffsDTt(DTtPCoeffs, Ss, Tt); + calcPCoeffsDSs(DSsPCoeffs, Ss, Tt); + + SpecVol = calcRefProfile(Pdb) + calcDelta(SpecVolPCoeffs, Pdb); + + /// Chain rule from the normalized variables of the polynomial to the + /// physical ones: dTt/dCt = 1 / CtNorm, dSs/dSa = 1 / (2 SaNorm Ss), + /// and dPdb/dP = Pa2Db. + SpecVolDCt = calcDeltaDeriv(DTtPCoeffs, Pdb) / CtNorm; + SpecVolDSa = calcDeltaDeriv(DSsPCoeffs, Pdb) * 0.5_Real / (SaNorm * Ss); + SpecVolDP = + (calcRefProfileDP(Pdb) + calcDeltaDP(SpecVolPCoeffs, Pdb)) * Pa2Db; + } + + /// Calculate the TEOS-10 specific volume and its three first derivatives + /// over the active layers of a cell. Pressure is the relative pressure in + /// Pa; the derivatives are per degC, per (g/kg), and per Pa. This is the + /// array-level counterpart of calcSpecVolAndDerivsAtPoint above and is what + /// Eos::computeSpecVolAndDerivs calls for each cell. + KOKKOS_FUNCTION void calcSpecVolAndDerivsOnCells( + Array2DReal SpecVol, Array2DReal SpecVolDCt, Array2DReal SpecVolDSa, + Array2DReal SpecVolDP, const TeamMember &Team, I4 ICell, + const Array2DReal &ConservTemp, const Array2DReal &AbsSalinity, + const Array2DReal &Pressure) const { + + const I4 KMin = MinLayerCell(ICell); + const I4 KMax = MaxLayerCell(ICell); + + parallelForInner( + Team, Range{KMin, KMax}, INNER_LAMBDA(int K) { + calcSpecVolAndDerivsAtPoint( + ConservTemp(ICell, K), AbsSalinity(ICell, K), + Pressure(ICell, K), SpecVol(ICell, K), SpecVolDCt(ICell, K), + SpecVolDSa(ICell, K), SpecVolDP(ICell, K)); + }); + } + /// Calculate 2nd derivative of Gibbs wrt pot temp at ref P for TEOS-10 KOKKOS_FUNCTION Real calcGibbsDerivPt0Pt0(Real Sa, Real P) const { Real x2 = Sfac * Sa; @@ -426,22 +731,48 @@ class LinearEos { LinearEos(const VertCoord *VCoord); // The functor takes the full arrays of specific volume (inout), - // the indices ICell and KChunk, and the ocean tracers (conservative) - // temperature, and (absolute) salinity as inputs, and outputs the - // linear specific volume. - KOKKOS_FUNCTION void operator()(Array2DReal SpecVol, I4 ICell, I4 KChunk, - const Array2DReal &ConservTemp, + // the team member and the cell index ICell, and the ocean tracers + // (conservative) temperature, and (absolute) salinity as inputs, and + // outputs the linear specific volume. + KOKKOS_FUNCTION void operator()(Array2DReal SpecVol, const TeamMember &Team, + I4 ICell, const Array2DReal &ConservTemp, const Array2DReal &AbsSalinity) const { - const I4 KStart = chunkStart(KChunk, MinLayerCell(ICell)); - const I4 KLen = chunkLength(KChunk, KStart, MaxLayerCell(ICell)); + const I4 KMin = MinLayerCell(ICell); + const I4 KMax = MaxLayerCell(ICell); + + parallelForInner( + Team, Range{KMin, KMax}, INNER_LAMBDA(int K) { + SpecVol(ICell, K) = + 1.0_Real / (RhoT0S0 + (DRhodT * ConservTemp(ICell, K) + + DRhodS * AbsSalinity(ICell, K))); + }); + } - for (int KVec = 0; KVec < KLen; ++KVec) { - const I4 K = KStart + KVec; - SpecVol(ICell, K) = - 1.0_Real / (RhoT0S0 + (DRhodT * ConservTemp(ICell, K) + - DRhodS * AbsSalinity(ICell, K))); - } + /// Calculate the linear specific volume and its three first derivatives + /// over the active layers of a cell. With + /// SpecVol = 1 / (RhoT0S0 + DRhodT * Ct + DRhodS * Sa) + /// the derivatives are -DRhodT * SpecVol^2 and -DRhodS * SpecVol^2, and + /// the linear EOS has no pressure dependence at all. + KOKKOS_FUNCTION void calcSpecVolAndDerivsOnCells( + Array2DReal SpecVol, Array2DReal SpecVolDCt, Array2DReal SpecVolDSa, + Array2DReal SpecVolDP, const TeamMember &Team, I4 ICell, + const Array2DReal &ConservTemp, const Array2DReal &AbsSalinity) const { + + const I4 KMin = MinLayerCell(ICell); + const I4 KMax = MaxLayerCell(ICell); + + parallelForInner( + Team, Range{KMin, KMax}, INNER_LAMBDA(int K) { + const Real Sv = + 1.0_Real / (RhoT0S0 + (DRhodT * ConservTemp(ICell, K) + + DRhodS * AbsSalinity(ICell, K))); + + SpecVol(ICell, K) = Sv; + SpecVolDCt(ICell, K) = -DRhodT * Sv * Sv; + SpecVolDSa(ICell, K) = -DRhodS * Sv * Sv; + SpecVolDP(ICell, K) = 0.0_Real; + }); } private: @@ -455,22 +786,43 @@ class ConstantEos { /// constructor declaration ConstantEos(const VertCoord *VCoord); - // The functor takes the full arrays of specific volume (inout), - // the indices ICell and KChunk, and returns a constant specific volume + // The functor takes the full arrays of specific volume (inout), the team + // member and the cell index ICell, and returns a constant specific volume // value for all active layers. - KOKKOS_FUNCTION void operator()(Array2DReal SpecVol, I4 ICell, I4 KChunk, - const Array2DReal &ConservTemp, + KOKKOS_FUNCTION void operator()(Array2DReal SpecVol, const TeamMember &Team, + I4 ICell, const Array2DReal &ConservTemp, const Array2DReal &AbsSalinity) const { - const I4 KStart = chunkStart(KChunk, MinLayerCell(ICell)); - const I4 KLen = chunkLength(KChunk, KStart, MaxLayerCell(ICell)); + const I4 KMin = MinLayerCell(ICell); + const I4 KMax = MaxLayerCell(ICell); + (void)ConservTemp; + (void)AbsSalinity; + + parallelForInner( + Team, Range{KMin, KMax}, + INNER_LAMBDA(int K) { SpecVol(ICell, K) = 1.0_Real / RhoSw; }); + } + + /// Calculate the constant specific volume and its three first derivatives + /// over the active layers of a cell. The specific volume does not depend on + /// temperature, salinity or pressure, so all three derivatives vanish. + KOKKOS_FUNCTION void calcSpecVolAndDerivsOnCells( + Array2DReal SpecVol, Array2DReal SpecVolDCt, Array2DReal SpecVolDSa, + Array2DReal SpecVolDP, const TeamMember &Team, I4 ICell, + const Array2DReal &ConservTemp, const Array2DReal &AbsSalinity) const { + + const I4 KMin = MinLayerCell(ICell); + const I4 KMax = MaxLayerCell(ICell); (void)ConservTemp; (void)AbsSalinity; - for (int KVec = 0; KVec < KLen; ++KVec) { - const I4 K = KStart + KVec; - SpecVol(ICell, K) = 1.0_Real / RhoSw; - } + parallelForInner( + Team, Range{KMin, KMax}, INNER_LAMBDA(int K) { + SpecVol(ICell, K) = 1.0_Real / RhoSw; + SpecVolDCt(ICell, K) = 0.0_Real; + SpecVolDSa(ICell, K) = 0.0_Real; + SpecVolDP(ICell, K) = 0.0_Real; + }); } private: @@ -485,196 +837,91 @@ class Teos10BruntVaisalaFreqSq { Teos10BruntVaisalaFreqSq(const VertCoord *VCoord); // The functor takes the full arrays of squared Brunt-Vaisala frequency - // (inout) the index ICell, and the ocean tracers (conservative) - // temperature, (absolute) salinity, relative pressure (gauge pressure, - // i.e., absolute pressure minus the standard atmosphere), and specific - // volume as inputs, and outputs the squared Brunt-Vaisala frequency. - KOKKOS_FUNCTION void operator()(Array2DReal BruntVaisalaFreqSq, I4 ICell, - I4 KChunk, const Array2DReal &ConservTemp, - const Array2DReal &AbsSalinity, - const Array2DReal &Pressure, - const Array2DReal &SpecVol) const { - - const I4 KStart = chunkStart(KChunk, MinLayerCell(ICell) + 1); - const I4 KLen = chunkLength(KChunk, KStart, MaxLayerCell(ICell)); - - for (int KVec = 0; KVec < KLen; ++KVec) { - const I4 K = KStart + KVec; - // Calculate squared Brunt-Vaisala frequency - Real CtInt = - 0.5_Real * (ConservTemp(ICell, K) + ConservTemp(ICell, K - 1)); - Real SaInt = - 0.5_Real * (AbsSalinity(ICell, K) + AbsSalinity(ICell, K - 1)); - Real PInt = 0.5_Real * (Pressure(ICell, K) + Pressure(ICell, K - 1)); - Real SpInt = 0.5_Real * (SpecVol(ICell, K) + SpecVol(ICell, K - 1)); - Real AlphaInt = calcAlpha(SaInt, CtInt, PInt * Pa2Db, SpInt); - Real BetaInt = calcBeta(SaInt, CtInt, PInt * Pa2Db, SpInt); - Real DSa = AbsSalinity(ICell, K) - AbsSalinity(ICell, K - 1); - Real DCt = ConservTemp(ICell, K) - ConservTemp(ICell, K - 1); - Real DP = Pressure(ICell, K) - Pressure(ICell, K - 1); - - BruntVaisalaFreqSq(ICell, K) = Gravity * Gravity * - (BetaInt * DSa - AlphaInt * DCt) / - (SpInt * DP); - } + // (inout) the team member and the cell index ICell, and the ocean tracers + // (conservative) temperature, (absolute) salinity, relative pressure + // (gauge pressure, i.e., absolute pressure minus the standard + // atmosphere), and specific volume as inputs, and outputs the squared + // Brunt-Vaisala frequency. + KOKKOS_FUNCTION void + operator()(Array2DReal BruntVaisalaFreqSq, const TeamMember &Team, I4 ICell, + const Array2DReal &ConservTemp, const Array2DReal &AbsSalinity, + const Array2DReal &Pressure, const Array2DReal &SpecVol) const { + + // Compute at interior vertical interfaces only + const I4 KMin = MinLayerCell(ICell) + 1; + const I4 KMax = MaxLayerCell(ICell); + + parallelForInner( + Team, Range{KMin, KMax}, INNER_LAMBDA(int K) { + // Calculate squared Brunt-Vaisala frequency + Real CtInt = + 0.5_Real * (ConservTemp(ICell, K) + ConservTemp(ICell, K - 1)); + Real SaInt = + 0.5_Real * (AbsSalinity(ICell, K) + AbsSalinity(ICell, K - 1)); + Real PInt = + 0.5_Real * (Pressure(ICell, K) + Pressure(ICell, K - 1)); + Real SpInt = + 0.5_Real * (SpecVol(ICell, K) + SpecVol(ICell, K - 1)); + Real AlphaInt = calcAlpha(SaInt, CtInt, PInt * Pa2Db, SpInt); + Real BetaInt = calcBeta(SaInt, CtInt, PInt * Pa2Db, SpInt); + Real DSa = AbsSalinity(ICell, K) - AbsSalinity(ICell, K - 1); + Real DCt = ConservTemp(ICell, K) - ConservTemp(ICell, K - 1); + Real DP = Pressure(ICell, K) - Pressure(ICell, K - 1); + + BruntVaisalaFreqSq(ICell, K) = Gravity * Gravity * + (BetaInt * DSa - AlphaInt * DCt) / + (SpInt * DP); + }); } - /// Calculate alpha values for the squared Brunt-Vaisala frequency + /// Calculate alpha (the thermal expansion coefficient) for the squared + /// Brunt-Vaisala frequency. Alpha is the temperature derivative of the + /// specific volume divided by the specific volume, so it is formed from the + /// TEOS-10 derivative helpers rather than from a second copy of the same + /// polynomial coefficients. P is relative pressure in dbar and Sp is the + /// specific volume. + /// + /// The coefficients are assembled here rather than taken from the + /// Eos::SpecVolDCt array because the two are not evaluated at the same + /// state: this is called at the interface, with the temperature, salinity + /// and pressure averaged from the two adjacent layers, while SpecVolDCt + /// holds the derivative at the layer centers. Averaging the layer-center + /// derivatives instead would be a different approximation and would change + /// answers. The stored derivatives are also filled only when + /// computeSpecVolAndDerivs is called, which the Brunt-Vaisala calculation + /// cannot assume. KOKKOS_FUNCTION Real calcAlpha(Real Sa, Real Ct, Real P, Real Sp) const { - constexpr Real Factor = 0.0248826675584615; - constexpr Real Offset = 5.971840214030754e-1; - constexpr Real PNorm = 1.0e-4; - Real Ss = Kokkos::sqrt(Factor * Sa + Offset); - Real Tt = 0.025_Real * Ct; - Real Pp = P * PNorm; - constexpr Real A000 = -1.56497346750e-5; - constexpr Real A001 = 1.85057654290e-5; - constexpr Real A002 = -1.17363867310e-6; - constexpr Real A003 = -3.65270065530e-7; - constexpr Real A004 = 3.14540999020e-7; - constexpr Real A010 = 5.55242129680e-5; - constexpr Real A011 = -2.34332137060e-5; - constexpr Real A012 = 4.26100574800e-6; - constexpr Real A013 = 5.73918103180e-7; - constexpr Real A020 = -4.95634777770e-5; - constexpr Real A021 = 2.37838968519e-5; - constexpr Real A022 = -1.38397620111e-6; - constexpr Real A030 = 2.76445290808e-5; - constexpr Real A031 = -1.36408749928e-5; - constexpr Real A032 = -2.53411666056e-7; - constexpr Real A040 = -4.02698077700e-6; - constexpr Real A041 = 2.53683834070e-6; - constexpr Real A050 = 1.23258565608e-6; - constexpr Real A100 = 3.50095997640e-5; - constexpr Real A101 = -9.56770881560e-6; - constexpr Real A102 = -5.56991545570e-6; - constexpr Real A103 = -2.72956962370e-7; - constexpr Real A110 = -7.48716846880e-5; - constexpr Real A111 = -4.73566167220e-7; - constexpr Real A112 = 7.82747741600e-7; - constexpr Real A120 = 7.24244384490e-5; - constexpr Real A121 = -1.03676320965e-5; - constexpr Real A122 = 2.32856664276e-8; - constexpr Real A130 = -3.50383492616e-5; - constexpr Real A131 = 5.18268711320e-6; - constexpr Real A140 = -1.65263794500e-6; - constexpr Real A200 = -4.35926785610e-5; - constexpr Real A201 = 1.11008347650e-5; - constexpr Real A202 = 5.46207488340e-6; - constexpr Real A210 = 7.18156455200e-5; - constexpr Real A211 = 5.85666925900e-6; - constexpr Real A212 = -1.31462208134e-6; - constexpr Real A220 = -4.30608991440e-5; - constexpr Real A221 = 9.49659182340e-7; - constexpr Real A230 = 1.74814722392e-5; - constexpr Real A300 = 3.45324618280e-5; - constexpr Real A301 = -9.84471178440e-6; - constexpr Real A302 = -1.35441856270e-6; - constexpr Real A310 = -3.73971683740e-5; - constexpr Real A311 = -9.76522784000e-7; - constexpr Real A320 = 6.85899736680e-6; - constexpr Real A400 = -1.19594097880e-5; - constexpr Real A401 = 2.59092252600e-6; - constexpr Real A410 = 7.71906784880e-6; - constexpr Real A500 = 1.38645945810e-6; + Real DTtPCoeffs[5]; + + const Real Ss = + Kokkos::sqrt((Sa + Teos10Eos::DeltaS) / Teos10Eos::SaNorm); + const Real Tt = Ct / Teos10Eos::CtNorm; - Real Rval = - A000 + - Ss * (A100 + Ss * (A200 + Ss * (A300 + Ss * (A400 + A500 * Ss)))) + - Tt * (A010 + Ss * (A110 + Ss * (A210 + Ss * (A310 + A410 * Ss))) + - Tt * (A020 + Ss * (A120 + Ss * (A220 + A320 * Ss)) + - Tt * (A030 + Ss * (A130 + A230 * Ss) + - Tt * (A040 + A140 * Ss + A050 * Tt)))) + - Pp * (A001 + Ss * (A101 + Ss * (A201 + Ss * (A301 + A401 * Ss))) + - Tt * (A011 + Ss * (A111 + Ss * (A211 + A311 * Ss)) + - Tt * (A021 + Ss * (A121 + A221 * Ss) + - Tt * (A031 + A131 * Ss + A041 * Tt))) + - Pp * (A002 + Ss * (A102 + Ss * (A202 + A302 * Ss)) + - Tt * (A012 + Ss * (A112 + A212 * Ss) + - Tt * (A022 + A122 * Ss + A032 * Tt)) + - Pp * (A003 + A103 * Ss + A013 * Tt + A004 * Pp))); - - return 0.025_Real * Rval / Sp; + Teos10Eos::calcPCoeffsDTt(DTtPCoeffs, Ss, Tt); + + return Teos10Eos::calcDeltaDeriv(DTtPCoeffs, P) / + (Teos10Eos::CtNorm * Sp); } - /// Calculate beta values for the squared Brunt-Vaisala frequency + /// Calculate beta (the haline contraction coefficient) for the squared + /// Brunt-Vaisala frequency. Beta is minus the salinity derivative of the + /// specific volume divided by the specific volume. P is relative pressure + /// in dbar and Sp is the specific volume. As for calcAlpha above, this is + /// evaluated at the interface state and so cannot reuse the layer-center + /// Eos::SpecVolDSa array. KOKKOS_FUNCTION Real calcBeta(Real Sa, Real Ct, Real P, Real Sp) const { - constexpr Real Factor = 0.0248826675584615; - constexpr Real Offset = 5.971840214030754e-1; - constexpr Real PNorm = 1.0e-4; - Real Ss = Kokkos::sqrt(Factor * Sa + Offset); - Real Tt = 0.025_Real * Ct; - Real Pp = P * PNorm; - constexpr Real B000 = -3.10389819760e-4; - constexpr Real B003 = 3.63101885150e-7; - constexpr Real B004 = -1.11471254230e-7; - constexpr Real B010 = 3.50095997640e-5; - constexpr Real B013 = -2.72956962370e-7; - constexpr Real B020 = -3.74358423440e-5; - constexpr Real B030 = 2.41414794830e-5; - constexpr Real B040 = -8.75958731540e-6; - constexpr Real B050 = -3.30527589000e-7; - constexpr Real B100 = 1.33856134076e-3; - constexpr Real B103 = 3.34926075600e-8; - constexpr Real B110 = -8.71853571220e-5; - constexpr Real B120 = 7.18156455200e-5; - constexpr Real B130 = -2.87072660960e-5; - constexpr Real B140 = 8.74073611960e-6; - constexpr Real B200 = -2.55143801811e-3; - constexpr Real B210 = 1.03597385484e-4; - constexpr Real B220 = -5.60957525610e-5; - constexpr Real B230 = 6.85899736680e-6; - constexpr Real B300 = 2.32344279772e-3; - constexpr Real B310 = -4.78376391520e-5; - constexpr Real B320 = 1.54381356976e-5; - constexpr Real B400 = -1.05461852535e-3; - constexpr Real B410 = 6.93229729050e-6; - constexpr Real B500 = 1.91594743830e-4; - constexpr Real B001 = 2.42624687470e-5; - constexpr Real B011 = -9.56770881560e-6; - constexpr Real B021 = -2.36783083610e-7; - constexpr Real B031 = -3.45587736550e-6; - constexpr Real B041 = 1.29567177830e-6; - constexpr Real B101 = -6.95849219480e-5; - constexpr Real B111 = 2.22016695300e-5; - constexpr Real B121 = 5.85666925900e-6; - constexpr Real B131 = 6.33106121560e-7; - constexpr Real B201 = 1.12412331915e-4; - constexpr Real B211 = -2.95341353532e-5; - constexpr Real B221 = -1.46478417600e-6; - constexpr Real B301 = -6.92888744480e-5; - constexpr Real B311 = 1.03636901040e-5; - constexpr Real B401 = 1.54637136265e-5; - constexpr Real B002 = -5.84844329840e-7; - constexpr Real B012 = -5.56991545570e-6; - constexpr Real B022 = 3.91373870800e-7; - constexpr Real B032 = 7.76188880920e-9; - constexpr Real B102 = -9.62445031940e-6; - constexpr Real B112 = 1.09241497668e-5; - constexpr Real B122 = -1.31462208134e-6; - constexpr Real B202 = 1.47789320994e-5; - constexpr Real B212 = -4.06325568810e-6; - constexpr Real B302 = -7.12478989080e-6; + Real DSsPCoeffs[5]; - Real Rval = - B000 + - Ss * (B100 + Ss * (B200 + Ss * (B300 + Ss * (B400 + B500 * Ss)))) + - Tt * (B010 + Ss * (B110 + Ss * (B210 + Ss * (B310 + B410 * Ss))) + - Tt * (B020 + Ss * (B120 + Ss * (B220 + B320 * Ss)) + - Tt * (B030 + Ss * (B130 + B230 * Ss) + - Tt * (B040 + B140 * Ss + B050 * Tt)))) + - Pp * (B001 + Ss * (B101 + Ss * (B201 + Ss * (B301 + B401 * Ss))) + - Tt * (B011 + Ss * (B111 + Ss * (B211 + B311 * Ss)) + - Tt * (B021 + Ss * (B121 + B221 * Ss) + - Tt * (B031 + B131 * Ss + B041 * Tt))) + - Pp * (B002 + Ss * (B102 + Ss * (B202 + B302 * Ss)) + - Tt * (B012 + Ss * (B112 + B212 * Ss) + - Tt * (B022 + B122 * Ss + B032 * Tt)) + - Pp * (B003 + B103 * Ss + B013 * Tt + B004 * Pp))); - - return -0.5_Real * Rval * Factor / (Sp * Ss); + const Real Ss = + Kokkos::sqrt((Sa + Teos10Eos::DeltaS) / Teos10Eos::SaNorm); + const Real Tt = Ct / Teos10Eos::CtNorm; + + Teos10Eos::calcPCoeffsDSs(DSsPCoeffs, Ss, Tt); + + return -0.5_Real * Teos10Eos::calcDeltaDeriv(DSsPCoeffs, P) / + (Teos10Eos::SaNorm * Ss * Sp); } private: @@ -689,26 +936,28 @@ class LinearBruntVaisalaFreqSq { LinearBruntVaisalaFreqSq(const VertCoord *VCoord); // The functor takes the full arrays of squared Brunt-Vaisala frequency - // (inout), the index ICell, and the specific volume and pseudo-thickness - // as inputs, and outputs the squared Brunt-Vaisala frequency. - KOKKOS_FUNCTION void operator()(Array2DReal BruntVaisalaFreqSq, I4 ICell, - I4 KChunk, + // (inout), the team member and the cell index ICell, and the specific + // volume and pseudo-thickness as inputs, and outputs the squared + // Brunt-Vaisala frequency. + KOKKOS_FUNCTION void operator()(Array2DReal BruntVaisalaFreqSq, + const TeamMember &Team, I4 ICell, const Array2DReal &SpecVol) const { - const I4 KStart = chunkStart(KChunk, MinLayerCell(ICell) + 1); - const I4 KLen = chunkLength(KChunk, KStart, MaxLayerCell(ICell)); - - for (int KVec = 0; KVec < KLen; ++KVec) { - const I4 K = KStart + KVec; - /// Calculate squared Brunt-Vaisala frequency at mid-point between - /// K-1 and K Do not need to use displaced specific volume here - /// since only the linear EOS is used with this BVF formulation. - BruntVaisalaFreqSq(ICell, K) = - -(Gravity / RhoSw) * - ((1.0_Real / SpecVol(ICell, K - 1)) - - (1.0_Real / SpecVol(ICell, K))) / - (GeomZMid(ICell, K - 1) - GeomZMid(ICell, K)); - } + // Compute at interior vertical interfaces only + const I4 KMin = MinLayerCell(ICell) + 1; + const I4 KMax = MaxLayerCell(ICell); + + parallelForInner( + Team, Range{KMin, KMax}, INNER_LAMBDA(int K) { + /// Calculate squared Brunt-Vaisala frequency at mid-point between + /// K-1 and K Do not need to use displaced specific volume here + /// since only the linear EOS is used with this BVF formulation. + BruntVaisalaFreqSq(ICell, K) = + -(Gravity / RhoSw) * + ((1.0_Real / SpecVol(ICell, K - 1)) - + (1.0_Real / SpecVol(ICell, K))) / + (GeomZMid(ICell, K - 1) - GeomZMid(ICell, K)); + }); } private: @@ -730,14 +979,20 @@ class Eos { Array2DReal SpecVol; ///< Specific volume field at level centers Array2DReal SpecVolDisplaced; ///< Displaced specific volume field Array2DReal BruntVaisalaFreqSq; ///< Squared Brunt-Vaisala frequency field + Array2DReal SpecVolDCt; ///< d(SpecVol)/d(ConservTemp), per degC + Array2DReal SpecVolDSa; ///< d(SpecVol)/d(AbsSalinity), per (g/kg) + Array2DReal SpecVolDP; ///< d(SpecVol)/d(Pressure), per Pa std::string SpecVolFldName; ///< Field name for specific volume std::string SpecVolDisplacedFldName; ///< Field name for displaced specific volume std::string BruntVaisalaFreqSqFldName; ///< Field name for squared ///< Brunt-Vaisala frequency - std::string EosGroupName; ///< EOS group name (for config) - std::string Name; ///< Name of this EOS instance + std::string SpecVolDCtFldName; ///< Field name for temperature derivative + std::string SpecVolDSaFldName; ///< Field name for salinity derivative + std::string SpecVolDPFldName; ///< Field name for pressure derivative + std::string EosGroupName; ///< EOS group name (for config) + std::string Name; ///< Name of this EOS instance /// Compute specific volume for all cells/layers void computeSpecVol(const Array2DReal &ConservTemp, @@ -749,6 +1004,17 @@ class Eos { const Array2DReal &AbsSalinity, const Array2DReal &Pressure, I4 KDisp); + /// Compute specific volume together with its first derivatives with respect + /// to conservative temperature, absolute salinity and pressure, in a single + /// pass over the equation of state. Pressure is the relative pressure in Pa + /// and the derivatives are returned per degC, per (g/kg) and per Pa. The + /// results are stored in the SpecVol, SpecVolDCt, SpecVolDSa and SpecVolDP + /// members. Since SpecVol is computed here too, this replaces rather than + /// accompanies a call to computeSpecVol. + void computeSpecVolAndDerivs(const Array2DReal &ConservTemp, + const Array2DReal &AbsSalinity, + const Array2DReal &Pressure); + /// Compute squared Brunt-Vaisala frequency for all cells/layers void computeBruntVaisalaFreqSq(const Array2DReal &ConservTemp, const Array2DReal &AbsSalinity, diff --git a/components/omega/test/ocn/EosTest.cpp b/components/omega/test/ocn/EosTest.cpp index 0e7eb2c5be76..ed12201a6750 100644 --- a/components/omega/test/ocn/EosTest.cpp +++ b/components/omega/test/ocn/EosTest.cpp @@ -54,6 +54,14 @@ const Real LinearBVFExpValue = const Real GswBVFExpValue = 0.02081197958166906; // Expected value from GSW-C library +/// Linear EOS coefficients, matching the Linear subsection of the Eos group in +/// the test configuration. The expected derivatives follow from +/// SpecVol = 1 / (RhoT0S0 + DRhodT * Ct + DRhodS * Sa). +const Real LinearDRhodT = -0.2; +const Real LinearDRhodS = 0.8; +const Real LinearDCtExpValue = -LinearDRhodT * LinearExpValue * LinearExpValue; +const Real LinearDSaExpValue = -LinearDRhodS * LinearExpValue * LinearExpValue; + /// Test input values const Real Sa = 30.0; // Absolute Salinity in g/kg const Real Ct = 10.0; // Conservative Temperature in degC @@ -62,6 +70,48 @@ const Real P = 1000.0 * Db2Pa; // Pressure in Pa const I4 KDisp = 1; // Displate parcel to K=1 for TEOS-10 eos const Real RTol = 1e-10; // Relative tolerance for isApprox checks +/// States spanning the oceanographic range and its corners, used by the +/// specific volume derivative checks. The fresh end is included because that +/// is where the normalized salinity is smallest and the salinity derivative +/// worst conditioned, and the full pressure range because the reference +/// profile supplies almost all of the pressure derivative. +constexpr int NSaTest = 6; +constexpr int NCtTest = 6; +constexpr int NPTest = 6; +const Real SaTest[NSaTest] = {0.0, 5.0, 20.0, 30.0, 35.0, 38.5}; // g/kg +const Real CtTest[NCtTest] = {-2.0, 0.0, 4.0, 10.0, 25.0, 35.0}; // degC +const Real PTest[NPTest] = {0.0, 100.0, 1000.0, + 4000.0, 8000.0, 10000.0}; // dbar + +/// Relative tolerance for the specific volume derivative checks against the +/// GSW-C library. Omega and GSW-C evaluate the same polynomial in different +/// arrangements, so a few ulp of disagreement is expected; anything larger +/// than this means a term has been dropped or mis-scaled. +const Real DerivRTol = 1e-12; + +/// The pressure derivative needs a looser tolerance, and the reason is on the +/// GSW-C side rather than ours. GSW-C evaluates its v_P from a table of +/// coefficients that have been pre-multiplied by their pressure exponents and +/// rounded, so its result departs from the exact derivative of the 75-term +/// polynomial by about 2e-12 relative at 10000 dbar, growing with pressure. +/// The Omega implementation differentiates the full-precision coefficients and +/// agrees with the exact derivative, evaluated in 60-digit arithmetic, to +/// around 1e-16. This tolerance therefore bounds GSW-C's rounding, not ours; +/// it is still far tighter than any real mistake would produce, and the +/// finite-difference check below pins the value independently. +const Real DerivDPRTol = 1e-10; + +/// The temperature derivative passes through zero near the density maximum of +/// nearly fresh water, where a relative tolerance means nothing. This absolute +/// floor sits well above the roundoff of the polynomial sum that forms it +/// (terms of order 1e-5, so roundoff of order 1e-21) and far below any value +/// of physical interest (order 1e-7). +const Real DerivDCtATol = 1e-19; + +/// Likewise for the thermal expansion coefficient, which is the temperature +/// derivative divided by the specific volume (order 1e-3). +const Real AlphaATol = 1e-16; + /// The initialization routine for Eos testing. It calls various /// init routines, including the creation of the default decomposition. void initEosTest(const std::string &mesh) { @@ -706,6 +756,339 @@ void testBruntVaisalaFreqSqTeos10() { return; } +/// Test the array-level TEOS-10 specific volume derivatives over the mesh. +/// +/// The state varies with depth rather than being uniform, so that the check +/// covers a range of values and exercises the vertical chunking, and the +/// expected values are obtained layer by layer from GSW-C on the host. +/// +/// This is a test of the plumbing, not of the polynomial. It runs +/// Eos::computeSpecVolAndDerivs on the device over the whole mesh and so +/// covers the dispatch on EosChoice, the chunked vertical loop and its chunk +/// boundaries, the MinLayerCell/MaxLayerCell masking, the writing of the four +/// results into the Eos member arrays, and the registration of the derivative +/// fields in the Eos group. The math itself is covered point by point, over a +/// much wider range of states, by checkValueGswcSpecVolDerivs below; GSW-C +/// appears here only as a convenient source of expected values. +void testEosTeos10Derivs() { + /// Get mesh and coordinate info + const auto Mesh = HorzMesh::getDefault(); + const auto VCoord = VertCoord::getDefault(); + VCoord->NVertLayers = NVertLayers; + I4 NCellsSize = Mesh->NCellsSize; + /// Get Eos instance to test + Eos *TestEos = Eos::getInstance(); + TestEos->EosChoice = EosType::Teos10Eos; + + /// Create the ocean state arrays + Array2DReal SArray = Array2DReal("SArray", NCellsSize, NVertLayers); + Array2DReal TArray = Array2DReal("TArray", NCellsSize, NVertLayers); + Array2DReal PArray = Array2DReal("PArray", NCellsSize, NVertLayers); + deepCopy(TestEos->SpecVol, 0.0); + deepCopy(TestEos->SpecVolDCt, 0.0); + deepCopy(TestEos->SpecVolDSa, 0.0); + deepCopy(TestEos->SpecVolDP, 0.0); + + /// A state that gets saltier, colder and deeper with depth, spanning a + /// realistic part of the oceanographic range over the column + parallelFor( + "populateDerivArrays", {Mesh->NCellsAll, NVertLayers}, + KOKKOS_LAMBDA(I4 ICell, I4 K) { + SArray(ICell, K) = Sa + 0.1_Real * K; + TArray(ICell, K) = Ct - 0.15_Real * K; + PArray(ICell, K) = (100.0_Real + 150.0_Real * K) * Db2Pa; + }); + + /// Compute specific volume and its derivatives + TestEos->computeSpecVolAndDerivs(TArray, SArray, PArray); + + /// Take local handles on the Eos members for the reduction kernels + Array2DReal SpecVol = TestEos->SpecVol; + Array2DReal SpecVolDCt = TestEos->SpecVolDCt; + Array2DReal SpecVolDSa = TestEos->SpecVolDSa; + Array2DReal SpecVolDP = TestEos->SpecVolDP; + + /// Expected values per layer from GSW-C, computed on the host and copied to + /// the device for the comparison + HostArray1DReal ExpSpecVolH("ExpSpecVolH", NVertLayers); + HostArray1DReal ExpDCtH("ExpDCtH", NVertLayers); + HostArray1DReal ExpDSaH("ExpDSaH", NVertLayers); + HostArray1DReal ExpDPH("ExpDPH", NVertLayers); + + for (int K = 0; K < NVertLayers; ++K) { + const double SaVal = Sa + 0.1 * K; + const double CtVal = Ct - 0.15 * K; + const double PDb = 100.0 + 150.0 * K; + + double GswDSa, GswDCt, GswDP; + gsw_specvol_first_derivatives(SaVal, CtVal, PDb, &GswDSa, &GswDCt, + &GswDP); + + ExpSpecVolH(K) = gsw_specvol(SaVal, CtVal, PDb); + ExpDCtH(K) = GswDCt; + ExpDSaH(K) = GswDSa; + ExpDPH(K) = GswDP; + } + + auto ExpSpecVol = createDeviceMirrorCopy(ExpSpecVolH); + auto ExpDCt = createDeviceMirrorCopy(ExpDCtH); + auto ExpDSa = createDeviceMirrorCopy(ExpDSaH); + auto ExpDP = createDeviceMirrorCopy(ExpDPH); + + const auto &MinLayerCell = VCoord->MinLayerCell; + const auto &MaxLayerCell = VCoord->MaxLayerCell; + + /// Check all active cells and layers against the expected values + int NumMismatches = 0; + parallelReduceOuter( + "CheckSpecVolDerivs-Teos", {Mesh->NCellsAll}, + KOKKOS_LAMBDA(int ICell, const TeamMember &Team, int &OuterCount) { + int NumMismatchesCol; + const int KMin = MinLayerCell(ICell); + const int KMax = MaxLayerCell(ICell); + const int KRange = vertRange(KMin, KMax); + parallelReduceInner( + Team, KRange, + INNER_LAMBDA(int KOff, int &InnerCount) { + const int K = KMin + KOff; + if (!isApprox(SpecVol(ICell, K), ExpSpecVol(K), DerivRTol) or + !isApprox(SpecVolDCt(ICell, K), ExpDCt(K), DerivRTol, + DerivDCtATol) or + !isApprox(SpecVolDSa(ICell, K), ExpDSa(K), DerivRTol) or + !isApprox(SpecVolDP(ICell, K), ExpDP(K), DerivDPRTol)) { + InnerCount++; + } + }, + NumMismatchesCol); + + Kokkos::single(PerTeam(Team), + [&]() { OuterCount += NumMismatchesCol; }); + }, + NumMismatches); + + // If test fails, print bad values and abort + if (NumMismatches != 0) { + auto SpecVolH = createHostMirrorCopy(SpecVol); + auto SpecVolDCtH = createHostMirrorCopy(SpecVolDCt); + auto SpecVolDSaH = createHostMirrorCopy(SpecVolDSa); + auto SpecVolDPH = createHostMirrorCopy(SpecVolDP); + for (int I = 0; I < Mesh->NCellsAll; ++I) { + for (int K = 0; K < NVertLayers; ++K) { + if (!isApprox(SpecVolH(I, K), ExpSpecVolH(K), DerivRTol)) + LOG_ERROR("EosTest: SpecVol Deriv Bad Value: " + "SpecVol({},{}) = {}; Expected {}", + I, K, SpecVolH(I, K), ExpSpecVolH(K)); + if (!isApprox(SpecVolDCtH(I, K), ExpDCtH(K), DerivRTol, + DerivDCtATol)) + LOG_ERROR("EosTest: SpecVolDCt Bad Value: " + "SpecVolDCt({},{}) = {}; Expected {}", + I, K, SpecVolDCtH(I, K), ExpDCtH(K)); + if (!isApprox(SpecVolDSaH(I, K), ExpDSaH(K), DerivRTol)) + LOG_ERROR("EosTest: SpecVolDSa Bad Value: " + "SpecVolDSa({},{}) = {}; Expected {}", + I, K, SpecVolDSaH(I, K), ExpDSaH(K)); + if (!isApprox(SpecVolDPH(I, K), ExpDPH(K), DerivDPRTol)) + LOG_ERROR("EosTest: SpecVolDP Bad Value: " + "SpecVolDP({},{}) = {}; Expected {}", + I, K, SpecVolDPH(I, K), ExpDPH(K)); + } + } + ABORT_ERROR("EosTest: SpecVol Derivs TEOS FAIL with {} bad values", + NumMismatches); + } + + /// Check that each derivative is registered as a field in the Eos group + /// with the member array attached, so that it can be written to a stream + const std::string DerivFldNames[3] = {TestEos->SpecVolDCtFldName, + TestEos->SpecVolDSaFldName, + TestEos->SpecVolDPFldName}; + const Array2DReal DerivArrays[3] = {SpecVolDCt, SpecVolDSa, SpecVolDP}; + + for (int IFld = 0; IFld < 3; ++IFld) { + const std::string &FldName = DerivFldNames[IFld]; + + if (!Field::exists(FldName)) { + ABORT_ERROR("EosTest: SpecVol Derivs field {} does not exist", + FldName); + } + + if (!FieldGroup::isFieldInGroup(FldName, TestEos->EosGroupName)) { + ABORT_ERROR("EosTest: SpecVol Derivs field {} not in group {}", + FldName, TestEos->EosGroupName); + } + + auto DerivField = Field::get(FldName); + auto FieldData = DerivField->getDataArray(); + if (FieldData.data() != DerivArrays[IFld].data()) { + ABORT_ERROR("EosTest: SpecVol Derivs field {} does not alias the " + "Eos member array", + FldName); + } + } + + return; +} + +/// Test the array-level linear EOS specific volume derivatives, which are +/// known in closed form from the configured linear coefficients +void testEosLinearDerivs() { + /// Get mesh and coordinate info + const auto Mesh = HorzMesh::getDefault(); + const auto VCoord = VertCoord::getDefault(); + VCoord->NVertLayers = NVertLayers; + I4 NCellsSize = Mesh->NCellsSize; + /// Get Eos instance to test + Eos *TestEos = Eos::getInstance(); + TestEos->EosChoice = EosType::LinearEos; + + /// Create and fill ocean state arrays + Array2DReal SArray = Array2DReal("SArray", NCellsSize, NVertLayers); + Array2DReal TArray = Array2DReal("TArray", NCellsSize, NVertLayers); + Array2DReal PArray = Array2DReal("PArray", NCellsSize, NVertLayers); + deepCopy(SArray, Sa); + deepCopy(TArray, Ct); + deepCopy(PArray, P); + deepCopy(TestEos->SpecVol, 0.0); + deepCopy(TestEos->SpecVolDCt, 0.0); + deepCopy(TestEos->SpecVolDSa, 0.0); + deepCopy(TestEos->SpecVolDP, 0.0); + + TestEos->computeSpecVolAndDerivs(TArray, SArray, PArray); + + /// Take local handles on the Eos members for the reduction kernels + Array2DReal SpecVol = TestEos->SpecVol; + Array2DReal SpecVolDCt = TestEos->SpecVolDCt; + Array2DReal SpecVolDSa = TestEos->SpecVolDSa; + Array2DReal SpecVolDP = TestEos->SpecVolDP; + + const auto &MinLayerCell = VCoord->MinLayerCell; + const auto &MaxLayerCell = VCoord->MaxLayerCell; + + /// Check all array values against the expected values + int NumMismatches = 0; + parallelReduceOuter( + "CheckSpecVolDerivs-linear", {Mesh->NCellsAll}, + KOKKOS_LAMBDA(int ICell, const TeamMember &Team, int &OuterCount) { + int NumMismatchesCol; + const int KMin = MinLayerCell(ICell); + const int KMax = MaxLayerCell(ICell); + const int KRange = vertRange(KMin, KMax); + parallelReduceInner( + Team, KRange, + INNER_LAMBDA(int KOff, int &InnerCount) { + const int K = KMin + KOff; + if (!isApprox(SpecVol(ICell, K), LinearExpValue, RTol) or + !isApprox(SpecVolDCt(ICell, K), LinearDCtExpValue, RTol) or + !isApprox(SpecVolDSa(ICell, K), LinearDSaExpValue, RTol) or + SpecVolDP(ICell, K) != 0.0_Real) { + InnerCount++; + } + }, + NumMismatchesCol); + + Kokkos::single(PerTeam(Team), + [&]() { OuterCount += NumMismatchesCol; }); + }, + NumMismatches); + + // If test fails, print bad values and abort + if (NumMismatches != 0) { + auto SpecVolDCtH = createHostMirrorCopy(SpecVolDCt); + auto SpecVolDSaH = createHostMirrorCopy(SpecVolDSa); + auto SpecVolDPH = createHostMirrorCopy(SpecVolDP); + for (int I = 0; I < Mesh->NCellsAll; ++I) { + for (int K = 0; K < NVertLayers; ++K) { + if (!isApprox(SpecVolDCtH(I, K), LinearDCtExpValue, RTol)) + LOG_ERROR("EosTest: SpecVolDCt Linear Bad Value: " + "SpecVolDCt({},{}) = {}; Expected {}", + I, K, SpecVolDCtH(I, K), LinearDCtExpValue); + if (!isApprox(SpecVolDSaH(I, K), LinearDSaExpValue, RTol)) + LOG_ERROR("EosTest: SpecVolDSa Linear Bad Value: " + "SpecVolDSa({},{}) = {}; Expected {}", + I, K, SpecVolDSaH(I, K), LinearDSaExpValue); + if (SpecVolDPH(I, K) != 0.0_Real) + LOG_ERROR("EosTest: SpecVolDP Linear Bad Value: " + "SpecVolDP({},{}) = {}; Expected 0", + I, K, SpecVolDPH(I, K)); + } + } + ABORT_ERROR("EosTest: SpecVol Derivs Linear FAIL with {} bad values", + NumMismatches); + } + + return; +} + +/// Test the array-level constant EOS specific volume derivatives, all of which +/// must be identically zero +void testEosConstantDerivs() { + /// Get mesh and coordinate info + const auto Mesh = HorzMesh::getDefault(); + const auto VCoord = VertCoord::getDefault(); + VCoord->NVertLayers = NVertLayers; + I4 NCellsSize = Mesh->NCellsSize; + /// Get Eos instance to test + Eos *TestEos = Eos::getInstance(); + TestEos->EosChoice = EosType::ConstantEos; + + /// Create and fill ocean state arrays + Array2DReal SArray = Array2DReal("SArray", NCellsSize, NVertLayers); + Array2DReal TArray = Array2DReal("TArray", NCellsSize, NVertLayers); + Array2DReal PArray = Array2DReal("PArray", NCellsSize, NVertLayers); + deepCopy(SArray, Sa); + deepCopy(TArray, Ct); + deepCopy(PArray, P); + deepCopy(TestEos->SpecVol, 0.0); + deepCopy(TestEos->SpecVolDCt, 0.0); + deepCopy(TestEos->SpecVolDSa, 0.0); + deepCopy(TestEos->SpecVolDP, 0.0); + + TestEos->computeSpecVolAndDerivs(TArray, SArray, PArray); + + /// Take local handles on the Eos members for the reduction kernels + Array2DReal SpecVol = TestEos->SpecVol; + Array2DReal SpecVolDCt = TestEos->SpecVolDCt; + Array2DReal SpecVolDSa = TestEos->SpecVolDSa; + Array2DReal SpecVolDP = TestEos->SpecVolDP; + + const auto &MinLayerCell = VCoord->MinLayerCell; + const auto &MaxLayerCell = VCoord->MaxLayerCell; + + /// Check all array values against the expected values + int NumMismatches = 0; + parallelReduceOuter( + "CheckSpecVolDerivs-Constant", {Mesh->NCellsAll}, + KOKKOS_LAMBDA(int ICell, const TeamMember &Team, int &OuterCount) { + int NumMismatchesCol; + const int KMin = MinLayerCell(ICell); + const int KMax = MaxLayerCell(ICell); + const int KRange = vertRange(KMin, KMax); + parallelReduceInner( + Team, KRange, + INNER_LAMBDA(int KOff, int &InnerCount) { + const int K = KMin + KOff; + if (!isApprox(SpecVol(ICell, K), ConstantExpValue, RTol) or + SpecVolDCt(ICell, K) != 0.0_Real or + SpecVolDSa(ICell, K) != 0.0_Real or + SpecVolDP(ICell, K) != 0.0_Real) { + InnerCount++; + } + }, + NumMismatchesCol); + + Kokkos::single(PerTeam(Team), + [&]() { OuterCount += NumMismatchesCol; }); + }, + NumMismatches); + + if (NumMismatches != 0) { + ABORT_ERROR("EosTest: SpecVol Derivs Constant FAIL with {} bad values", + NumMismatches); + } + + return; +} + /// Finalize and clean up all test infrastructure void finalizeEosTest() { Eos::destroyInstance(); @@ -788,6 +1171,261 @@ void checkValueCtFreezing() { return; } +/// Relative difference between two values, zero when both vanish +Real relDiff(Real X, Real Y) { + const Real Scale = std::max(std::abs(X), std::abs(Y)); + return Scale > 0.0 ? std::abs(X - Y) / Scale : 0.0; +} + +/// Test the TEOS-10 specific volume derivatives against the GSW-C library over +/// a range of states. +/// +/// GSW-C is used here unmodified and through its public API, as an independent +/// oracle. The Omega implementation does not derive from it: the derivatives +/// are the analytic derivatives of the Roquet et al. 2015 polynomial that the +/// Teos10Eos functor already carries. +/// +/// This is the test of the polynomial itself, as opposed to +/// testEosTeos10Derivs above, which tests the array-level machinery. It calls +/// the point-wise calcSpecVolAndDerivsAtPoint on the host at every combination +/// of the salinity, temperature and pressure values in SaTest, CtTest and +/// PTest, which reach the corners of the oceanographic range -- fresh and +/// salty, freezing and warm, surface and 10000 dbar -- rather than the single +/// realistic profile the mesh test uses. A dropped or mis-scaled term shows up +/// here, and the tolerances are tight enough to say so. +void checkValueGswcSpecVolDerivs() { + + Teos10Eos TestEos(VertCoord::getDefault()); + + int NumBad = 0; + Real WorstSv = 0.0; + Real WorstDCt = 0.0; + Real WorstDSa = 0.0; + Real WorstDP = 0.0; + int NumChecked = 0; + + for (int ISa = 0; ISa < NSaTest; ++ISa) { + for (int ICt = 0; ICt < NCtTest; ++ICt) { + for (int IP = 0; IP < NPTest; ++IP) { + + const Real SaVal = SaTest[ISa]; + const Real CtVal = CtTest[ICt]; + const Real PDb = PTest[IP]; + + Real SpecVol, SpecVolDCt, SpecVolDSa, SpecVolDP; + TestEos.calcSpecVolAndDerivsAtPoint(CtVal, SaVal, PDb * Db2Pa, + SpecVol, SpecVolDCt, SpecVolDSa, + SpecVolDP); + + /// GSW-C takes pressure in dbar and returns the derivatives per + /// (g/kg), per degC, and per Pa + double GswDSa, GswDCt, GswDP; + gsw_specvol_first_derivatives(SaVal, CtVal, PDb, &GswDSa, &GswDCt, + &GswDP); + const double GswSpecVol = gsw_specvol(SaVal, CtVal, PDb); + + WorstSv = std::max(WorstSv, relDiff(SpecVol, GswSpecVol)); + WorstDCt = std::max(WorstDCt, relDiff(SpecVolDCt, GswDCt)); + WorstDSa = std::max(WorstDSa, relDiff(SpecVolDSa, GswDSa)); + WorstDP = std::max(WorstDP, relDiff(SpecVolDP, GswDP)); + ++NumChecked; + + if (!isApprox(SpecVol, GswSpecVol, DerivRTol)) { + LOG_ERROR("EosTest: SpecVol Bad Value at Sa={}, Ct={}, " + "P={} dbar: expected {}, got {}", + SaVal, CtVal, PDb, GswSpecVol, SpecVol); + ++NumBad; + } + if (!isApprox(SpecVolDCt, GswDCt, DerivRTol, DerivDCtATol)) { + LOG_ERROR("EosTest: SpecVolDCt Bad Value at Sa={}, Ct={}, " + "P={} dbar: expected {}, got {}", + SaVal, CtVal, PDb, GswDCt, SpecVolDCt); + ++NumBad; + } + if (!isApprox(SpecVolDSa, GswDSa, DerivRTol)) { + LOG_ERROR("EosTest: SpecVolDSa Bad Value at Sa={}, Ct={}, " + "P={} dbar: expected {}, got {}", + SaVal, CtVal, PDb, GswDSa, SpecVolDSa); + ++NumBad; + } + if (!isApprox(SpecVolDP, GswDP, DerivDPRTol)) { + LOG_ERROR("EosTest: SpecVolDP Bad Value at Sa={}, Ct={}, " + "P={} dbar: expected {}, got {}", + SaVal, CtVal, PDb, GswDP, SpecVolDP); + ++NumBad; + } + } + } + } + + LOG_INFO("EosTest: TEOS-10 derivatives vs GSW-C over {} states, max " + "relative difference: SpecVol {}, d/dCt {}, d/dSa {}, d/dP {}", + NumChecked, WorstSv, WorstDCt, WorstDSa, WorstDP); + + if (NumBad != 0) { + ABORT_ERROR("EosTest: SpecVol derivatives vs GSW-C FAIL with {} bad " + "values", + NumBad); + } + + return; +} + +/// Check the specific volume derivatives against centered finite differences of +/// the Omega specific volume itself. +/// +/// This overlaps the GSW-C comparison above whenever that library is present +/// and correct, and that is deliberate: it pins the unit convention of the +/// Omega interface -- per degC, per (g/kg), and per Pa -- without reference to +/// GSW. A change that silently made the pressure derivative per dbar instead, +/// a factor of 1e4, would be caught here as well as there. +void checkFiniteDiffSpecVolDerivs() { + + Teos10Eos TestEos(VertCoord::getDefault()); + + /// Step sizes and tolerance are set by the finite difference itself: large + /// enough that the difference of two specific volumes is not lost to + /// roundoff, small enough that the truncation error stays below the + /// tolerance. + const Real DCtStep = 1.0e-2; // degC + const Real DSaStep = 1.0e-2; // g/kg + const Real DPStep = 1.0e5; // Pa (10 dbar) + const Real FDRTol = 1.0e-5; // limited by finite difference truncation + const Real FDCtATol = 1.0e-14; // finite difference noise floor + const Real FDSaATol = 1.0e-14; + const Real FDPATol = 1.0e-20; + + /// Evaluate only the specific volume at a perturbed state + auto SpecVolAt = [&TestEos](Real CtVal, Real SaVal, Real PPa) { + Real SpecVol, DCt, DSa, DP; + TestEos.calcSpecVolAndDerivsAtPoint(CtVal, SaVal, PPa, SpecVol, DCt, DSa, + DP); + return SpecVol; + }; + + int NumBad = 0; + + for (int ISa = 0; ISa < NSaTest; ++ISa) { + for (int ICt = 0; ICt < NCtTest; ++ICt) { + for (int IP = 0; IP < NPTest; ++IP) { + + const Real SaVal = SaTest[ISa]; + const Real CtVal = CtTest[ICt]; + const Real PPa = PTest[IP] * Db2Pa; + + Real SpecVol, SpecVolDCt, SpecVolDSa, SpecVolDP; + TestEos.calcSpecVolAndDerivsAtPoint( + CtVal, SaVal, PPa, SpecVol, SpecVolDCt, SpecVolDSa, SpecVolDP); + + const Real FDCt = (SpecVolAt(CtVal + DCtStep, SaVal, PPa) - + SpecVolAt(CtVal - DCtStep, SaVal, PPa)) / + (2.0_Real * DCtStep); + const Real FDSa = (SpecVolAt(CtVal, SaVal + DSaStep, PPa) - + SpecVolAt(CtVal, SaVal - DSaStep, PPa)) / + (2.0_Real * DSaStep); + const Real FDP = (SpecVolAt(CtVal, SaVal, PPa + DPStep) - + SpecVolAt(CtVal, SaVal, PPa - DPStep)) / + (2.0_Real * DPStep); + + if (!isApprox(SpecVolDCt, FDCt, FDRTol, FDCtATol)) { + LOG_ERROR("EosTest: SpecVolDCt disagrees with finite difference " + "at Sa={}, Ct={}, P={} dbar: {} vs {}", + SaVal, CtVal, PTest[IP], SpecVolDCt, FDCt); + ++NumBad; + } + if (!isApprox(SpecVolDSa, FDSa, FDRTol, FDSaATol)) { + LOG_ERROR("EosTest: SpecVolDSa disagrees with finite difference " + "at Sa={}, Ct={}, P={} dbar: {} vs {}", + SaVal, CtVal, PTest[IP], SpecVolDSa, FDSa); + ++NumBad; + } + if (!isApprox(SpecVolDP, FDP, FDRTol, FDPATol)) { + LOG_ERROR("EosTest: SpecVolDP disagrees with finite difference " + "at Sa={}, Ct={}, P={} dbar: {} vs {}", + SaVal, CtVal, PTest[IP], SpecVolDP, FDP); + ++NumBad; + } + } + } + } + + if (NumBad != 0) { + ABORT_ERROR("EosTest: SpecVol derivatives vs finite differences FAIL " + "with {} bad values", + NumBad); + } + + return; +} + +/// Test the thermal expansion and haline contraction coefficients used by the +/// TEOS-10 Brunt-Vaisala frequency against the GSW-C library. +/// +/// These are the specific volume derivatives divided by the specific volume, +/// so this covers the same polynomial from the other side. Before this check +/// existed, calcAlpha and calcBeta were exercised only through the single +/// hardcoded BruntVaisalaFreqSq value below, which is too loose to catch a +/// mistake in either of them. +void checkValueGswcAlphaBeta() { + + Teos10Eos TestEos(VertCoord::getDefault()); + Teos10BruntVaisalaFreqSq TestBvf(VertCoord::getDefault()); + + int NumBad = 0; + Real WorstAlpha = 0.0; + Real WorstBeta = 0.0; + + for (int ISa = 0; ISa < NSaTest; ++ISa) { + for (int ICt = 0; ICt < NCtTest; ++ICt) { + for (int IP = 0; IP < NPTest; ++IP) { + + const Real SaVal = SaTest[ISa]; + const Real CtVal = CtTest[ICt]; + const Real PDb = PTest[IP]; + + Real SpecVol, SpecVolDCt, SpecVolDSa, SpecVolDP; + TestEos.calcSpecVolAndDerivsAtPoint(CtVal, SaVal, PDb * Db2Pa, + SpecVol, SpecVolDCt, SpecVolDSa, + SpecVolDP); + + const Real Alpha = TestBvf.calcAlpha(SaVal, CtVal, PDb, SpecVol); + const Real Beta = TestBvf.calcBeta(SaVal, CtVal, PDb, SpecVol); + + double GswSpecVol, GswAlpha, GswBeta; + gsw_specvol_alpha_beta(SaVal, CtVal, PDb, &GswSpecVol, &GswAlpha, + &GswBeta); + + WorstAlpha = std::max(WorstAlpha, relDiff(Alpha, GswAlpha)); + WorstBeta = std::max(WorstBeta, relDiff(Beta, GswBeta)); + + if (!isApprox(Alpha, GswAlpha, DerivRTol, AlphaATol)) { + LOG_ERROR("EosTest: Alpha Bad Value at Sa={}, Ct={}, " + "P={} dbar: expected {}, got {}", + SaVal, CtVal, PDb, GswAlpha, Alpha); + ++NumBad; + } + if (!isApprox(Beta, GswBeta, DerivRTol)) { + LOG_ERROR("EosTest: Beta Bad Value at Sa={}, Ct={}, " + "P={} dbar: expected {}, got {}", + SaVal, CtVal, PDb, GswBeta, Beta); + ++NumBad; + } + } + } + } + + LOG_INFO("EosTest: alpha and beta vs GSW-C, max relative difference: " + "alpha {}, beta {}", + WorstAlpha, WorstBeta); + + if (NumBad != 0) { + ABORT_ERROR("EosTest: alpha and beta vs GSW-C FAIL with {} bad values", + NumBad); + } + + return; +} + /// Test that the Eos CT-from-PT helper matches GSW-C void checkValueGswcCtFromPt() { Eos *TestEos = Eos::getInstance(); @@ -822,6 +1460,12 @@ void checkValueGswcPtFromCt() { // Single value test: // --> test calls the external GSW-C library // and compares the specific volume to the published value +// --> next compares the TEOS-10 specific volume and its three first +// derivatives against GSW-C over a range of states +// --> next checks those derivatives against centered finite differences of the +// Omega specific volume, which pins the unit conventions independently of GSW +// --> next compares the thermal expansion and haline contraction coefficients +// against GSW-C over the same range of states // Full array tests: // --> one tests the value on a Eos with linear option // --> next checks the value on a Eos with linear displaced option @@ -829,8 +1473,10 @@ void checkValueGswcPtFromCt() { // calculation // --> next checks the value on a Eos with TEOS-10 option // --> next checks the value on a Eos with TEOS-10 displaced option -// --> last checks the value of the TOES-10 squared Brunt Vaisala Freq. +// --> next checks the value of the TOES-10 squared Brunt Vaisala Freq. // calculation +// --> last checks the specific volume derivatives for each of the three EOS +// options over the whole mesh void eosTest(const std::string &MeshFile = "OmegaMesh.nc") { initEosTest(MeshFile); const auto &Mesh = HorzMesh::getDefault(); @@ -840,6 +1486,9 @@ void eosTest(const std::string &MeshFile = "OmegaMesh.nc") { checkValueCtFreezing(); checkValueGswcCtFromPt(); checkValueGswcPtFromCt(); + checkValueGswcSpecVolDerivs(); + checkFiniteDiffSpecVolDerivs(); + checkValueGswcAlphaBeta(); testEosLinear(); testEosLinearDisplaced(); @@ -848,6 +1497,9 @@ void eosTest(const std::string &MeshFile = "OmegaMesh.nc") { testEosTeos10(); testEosTeos10Displaced(); testBruntVaisalaFreqSqTeos10(); + testEosTeos10Derivs(); + testEosLinearDerivs(); + testEosConstantDerivs(); finalizeEosTest();