From 13b1d613b0c29f0d53649f9c850f3a64f3b24541 Mon Sep 17 00:00:00 2001 From: Alice Barthel Date: Thu, 18 Sep 2025 15:01:27 -0700 Subject: [PATCH 1/7] added p,T,S limits (with CtFreezing calc) to eos Poly75t with a check on CtFreezing calculation and updated docs --- components/omega/doc/devGuide/EOS.md | 6 +- components/omega/doc/userGuide/EOS.md | 2 + components/omega/src/ocn/Eos.h | 128 +++++++++++++++++++++++--- components/omega/test/ocn/EosTest.cpp | 29 ++++++ 4 files changed, 149 insertions(+), 16 deletions(-) diff --git a/components/omega/doc/devGuide/EOS.md b/components/omega/doc/devGuide/EOS.md index 2cef706e9382..7d3d74286193 100644 --- a/components/omega/doc/devGuide/EOS.md +++ b/components/omega/doc/devGuide/EOS.md @@ -52,7 +52,11 @@ Eos.computeSpecVolDisp(ConsrvTemp, AbsSalinity, Pressure, KDisp); where `KDisp` is the number of `k` levels you want to displace each specific volume level to. For example, to displace each level to one below, set `KDisp = 1`. -## Removal of Eos +## Bounds check (and truncation) for the state variables (under Teos-10) + +The implemented 75-term polynomial for the calculation of the specific volume under Teos-10 is considered valid for ocean states contained in the ''oceanographic funnel'' defined in [McDougall et al., 2003](https://journals.ametsoc.org/view/journals/atot/20/5/1520-0426_2003_20_730_aaceaf_2_0_co_2.xml). When using teos-10, the Eos uses member methods `calcSLimits(P)` and `calcTLimits(Sa, P)` to calculate the valid ranges of Sa and T given the pressure. The conservative temperature lower bound depends is set by the freezing temperature, using the member method `calcCtFreezing(Sa, P, SaturationFract)`. This method implements the polynomial approximation of the conservative freezing temperature (called `gsw_ct_freezing_poly` in the GSW package), which is known to produce erros in the (-5e-4 K, 6e-4 K) range. Once we calculate the upper and lower bounds of validity, the state variables are clipped to the valid range (if outside the bounds) before we run the specific volume calculation. The state fields themselves are not changed. + + ## Removal of Eos To clear the Eos instance do: diff --git a/components/omega/doc/userGuide/EOS.md b/components/omega/doc/userGuide/EOS.md index 3a7ca6b40097..59eec0206d8a 100644 --- a/components/omega/doc/userGuide/EOS.md +++ b/components/omega/doc/userGuide/EOS.md @@ -20,3 +20,5 @@ 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` is also calculated by the EOS. This calculates the density of a parcel of fluid that is adiabatically displaced by a relative `k` levels, 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 EOS, `SpecVolDisplaced` will be the same as `SpecVol` since the specific volume calculation is independent of pressure/depth. + +When using `Teos-10`, the state variables are checked against the range over which the polynomial is considered to be valid (see [Roquet et al. 2015](https://www.sciencedirect.com/science/article/pii/S1463500315000566)). If the values are outside of the accepted values, we use the valid bounds for the specific volume calculation. Note that the state variable values themselves are not modified, only that they are not used as is in the calculation. diff --git a/components/omega/src/ocn/Eos.h b/components/omega/src/ocn/Eos.h index b9eb0f949c31..e18102a62c27 100644 --- a/components/omega/src/ocn/Eos.h +++ b/components/omega/src/ocn/Eos.h @@ -26,6 +26,11 @@ enum class EosType { Teos10Eos /// Roquet et al. 2015 75 term expansion }; +struct Range { + Real lo; + Real hi; +}; + /// TEOS10 75-term Polynomial Equation of State class Teos10Eos { public: @@ -52,7 +57,7 @@ class Teos10Eos { /// Calculate the local specific volume polynomial pressure /// coefficients calcPCoeffs(LocSpecVolPCoeffs, KVec, ConservTemp(ICell, K), - AbsSalinity(ICell, K)); + AbsSalinity(ICell, K), Pressure(ICell, K)); /// Calculate the specific volume at the given pressure /// If KDisp is 0, we use the current pressure, otherwise we use the @@ -78,12 +83,17 @@ class Teos10Eos { /// TEOS-10 helpers /// Calculate pressure polynomial coefficients for TEOS-10 KOKKOS_FUNCTION void calcPCoeffs(Array2DReal SpecVolPCoeffs, const I4 K, - const Real Ct, const Real Sa) const { + const Real Ct, const Real Sa, + const Real P) const { constexpr Real SAu = 40.0 * 35.16504 / 35.0; constexpr Real CTu = 40.0; constexpr Real DeltaS = 24.0; - Real Ss = Kokkos::sqrt((Sa + DeltaS) / SAu); - Real Tt = Ct / CTu; + Range Srange = calcSLimits(P); + Range Trange = calcTLimits(Sa, P); + Real SaInFunnel = Kokkos::clamp( + Sa, Srange.lo, Srange.hi); // Salt limited to Poly75t valid range + Real Ss = Kokkos::sqrt((SaInFunnel + DeltaS) / SAu); + Real Tt = Kokkos::clamp(Ct, Trange.lo, Trange.hi) / CTu; /// Coefficients for the polynomial expansion constexpr Real V000 = 1.0769995862e-03; @@ -205,8 +215,9 @@ class Teos10Eos { KOKKOS_FUNCTION Real calcDelta(const Array2DReal &SpecVolPCoeffs, const I4 K, const Real P) const { - constexpr Real Pu = 1e4; - Real Pp = P / Pu; + constexpr Real Pu = 1e4; + constexpr Real Pmax = 8000.0; + Real Pp = Kokkos::min(P, Pmax) / Pu; // P limited to Poly75t valid range Real Delta = ((((SpecVolPCoeffs(5, K) * Pp + SpecVolPCoeffs(4, K)) * Pp + SpecVolPCoeffs(3, K)) * @@ -220,15 +231,16 @@ class Teos10Eos { } /// Calculate reference profile for TEOS-10 - KOKKOS_FUNCTION Real calcRefProfile(Real P) const { - constexpr Real Pu = 1e4; - 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 / Pu; + KOKKOS_FUNCTION Real calcRefProfile(const Real P) const { + constexpr Real Pu = 1e4; + 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; + constexpr Real Pmax = 8000.0; + Real Pp = Kokkos::min(P, Pmax) / Pu; // P limited to Poly75t valid range Real V0 = (((((V05 * Pp + V04) * Pp + V03) * Pp + V02) * Pp + V01) * Pp + V00) * @@ -236,6 +248,92 @@ class Teos10Eos { return V0; } + /// Calculate S limits of validity given pressure p + KOKKOS_FUNCTION Range calcSLimits(const Real P) const { + Real lo = 0.0; + Real hi = 42.0; + Real lo2 = P * 5e-3 - 2.5; + Real lo3 = 30.0; + + if (P >= 500.0 && P < 6500.0) { + lo = Kokkos::max(lo, lo2); + } else if (P >= 6500.0) { + lo = Kokkos::max(lo, lo3); + } + return {lo, hi}; + } + + /// Calculate T limits of validity given p,Sa + KOKKOS_FUNCTION Range calcTLimits(const Real Sa, const Real P) const { + Real lo = -15.0; + Real hi = 95.0; + + if (P < 500.0) { + lo = calcCtFreezing(Sa, P, Real{0.0}); + } else if (P < 6500.0) { + lo = calcCtFreezing(Sa, Real{500.0}, Real{0.0}); + hi = 31.66666666666667 - P * 3.333333333333334e-3; + } else { + lo = calcCtFreezing(Sa, Real{500.0}, Real{0.0}); + hi = 10.0; + } + return {lo, hi}; + } + + /// Calculates freezing CTemperature (polynomial error in [-5e-4,6e-4] K) + KOKKOS_FUNCTION Real calcCtFreezing(const Real Sa, const Real P, + const Real SaturationFract) const { + constexpr Real Sso = 35.16504; + constexpr Real c0 = 0.017947064327968736; + constexpr Real c1 = -6.076099099929818; + constexpr Real c2 = 4.883198653547851; + constexpr Real c3 = -11.88081601230542; + constexpr Real c4 = 13.34658511480257; + constexpr Real c5 = -8.722761043208607; + constexpr Real c6 = 2.082038908808201; + constexpr Real c7 = -7.389420998107497; + constexpr Real c8 = -2.110913185058476; + constexpr Real c9 = 0.2295491578006229; + constexpr Real c10 = -0.9891538123307282; + constexpr Real c11 = -0.08987150128406496; + constexpr Real c12 = 0.3831132432071728; + constexpr Real c13 = 1.054318231187074; + constexpr Real c14 = 1.065556599652796; + constexpr Real c15 = -0.7997496801694032; + constexpr Real c16 = 0.3850133554097069; + constexpr Real c17 = -2.078616693017569; + constexpr Real c18 = 0.8756340772729538; + constexpr Real c19 = -2.079022768390933; + constexpr Real c20 = 1.596435439942262; + constexpr Real c21 = 0.1338002171109174; + constexpr Real c22 = 1.242891021876471; + + // Note: a = 0.502500117621 / Sso + constexpr Real a = 0.014289763856964; + constexpr Real b = 0.057000649899720; + + Real CtFreez; + + Real Sar = Sa * 1.0e-2; + Real X = Kokkos::sqrt(Sar); + Real Pr = P * 1.0e-4; + + CtFreez = + c0 + Sar * (c1 + X * (c2 + X * (c3 + X * (c4 + X * (c5 + c6 * X))))) + + Pr * (c7 + Pr * (c8 + c9 * Pr)) + + Sar * Pr * + (c10 + Pr * (c12 + Pr * (c15 + c21 * Sar)) + + Sar * (c13 + c17 * Pr + c19 * Sar) + + X * (c11 + Pr * (c14 + c18 * Pr) + + Sar * (c16 + c20 * Pr + c22 * Sar))); + + /* Adjust for the effects of dissolved air */ + CtFreez = CtFreez - SaturationFract * (1e-3) * (2.4 - a * Sa) * + (1.0 + b * (1.0 - Sa / Sso)); + + return CtFreez; + } + private: const int NVertLayers; }; diff --git a/components/omega/test/ocn/EosTest.cpp b/components/omega/test/ocn/EosTest.cpp index 028d136e9c16..10d0ffdf4507 100644 --- a/components/omega/test/ocn/EosTest.cpp +++ b/components/omega/test/ocn/EosTest.cpp @@ -315,6 +315,34 @@ int checkValueGswcSpecVol() { return Err; } +/// Test that the external GSW-C library returns the expected specific volume +int checkValueCtFreezing() { + int Err = 0; + const Real RTol = 1e-10; + Teos10Eos TestEos(1); + // TestEos->EosChoice = EosType::Teos10Eos; + constexpr Real SaturationFrac = 0.0; + constexpr Real P = 500.0; + constexpr Real Sa = 32.0; + + /// + /// Get specific volume from GSW-C library + double CtFreezGswc = gsw_ct_freezing_poly(Sa, P, SaturationFrac); + double CtFreez = TestEos.calcCtFreezing(Sa, P, SaturationFrac); + /// Check the value against the expected TEOS-10 value + bool Check = isApprox(CtFreezGswc, CtFreez, RTol); + if (!Check) { + Err++; + LOG_ERROR( + "checkValueCtFreezing: CtFreez isApprox FAIL, expected {}, got {}", + CtFreezGswc, CtFreez); + } + if (Err == 0) { + LOG_INFO("checkValueCtFreezing: PASS"); + } + return Err; +} + // the main test (all in one to have the same log) // Single value test: // --> test calls the external GSW-C library @@ -333,6 +361,7 @@ int eosTest(const std::string &MeshFile = "OmegaMesh.nc") { LOG_INFO("Single value checks:"); Err += checkValueGswcSpecVol(); + Err += checkValueCtFreezing(); LOG_INFO("Full array checks:"); Err += testEosLinear(); From 645f10888156621d1e899b022f687f37d4ab10ff Mon Sep 17 00:00:00 2001 From: Kat Smith Date: Thu, 2 Oct 2025 14:14:21 -0700 Subject: [PATCH 2/7] adds test for eos clamping --- components/omega/doc/devGuide/EOS.md | 6 +- components/omega/doc/userGuide/EOS.md | 2 +- components/omega/src/ocn/Eos.h | 118 +++++++++++----------- components/omega/test/ocn/EosTest.cpp | 136 ++++++++++++++++++++++++-- 4 files changed, 192 insertions(+), 70 deletions(-) diff --git a/components/omega/doc/devGuide/EOS.md b/components/omega/doc/devGuide/EOS.md index 7d3d74286193..5354062eaa29 100644 --- a/components/omega/doc/devGuide/EOS.md +++ b/components/omega/doc/devGuide/EOS.md @@ -15,7 +15,7 @@ An enumeration listing all implemented schemes is provided. It needs to be exten EOS is added. It is used to identify which EOS method is to be used at run time. ```c++ -enum class EosType { Linear, Teos10Poly75t }; +enum class EosType { LinearEos, Teos10Eos }; ``` ## Initialization @@ -52,9 +52,9 @@ Eos.computeSpecVolDisp(ConsrvTemp, AbsSalinity, Pressure, KDisp); where `KDisp` is the number of `k` levels you want to displace each specific volume level to. For example, to displace each level to one below, set `KDisp = 1`. -## Bounds check (and truncation) for the state variables (under Teos-10) +## Bounds check (and truncation) for the state variables (under TEOS-10) -The implemented 75-term polynomial for the calculation of the specific volume under Teos-10 is considered valid for ocean states contained in the ''oceanographic funnel'' defined in [McDougall et al., 2003](https://journals.ametsoc.org/view/journals/atot/20/5/1520-0426_2003_20_730_aaceaf_2_0_co_2.xml). When using teos-10, the Eos uses member methods `calcSLimits(P)` and `calcTLimits(Sa, P)` to calculate the valid ranges of Sa and T given the pressure. The conservative temperature lower bound depends is set by the freezing temperature, using the member method `calcCtFreezing(Sa, P, SaturationFract)`. This method implements the polynomial approximation of the conservative freezing temperature (called `gsw_ct_freezing_poly` in the GSW package), which is known to produce erros in the (-5e-4 K, 6e-4 K) range. Once we calculate the upper and lower bounds of validity, the state variables are clipped to the valid range (if outside the bounds) before we run the specific volume calculation. The state fields themselves are not changed. +The implemented 75-term polynomial for the calculation of the specific volume under TEOS-10 is considered valid for ocean states contained in the ''oceanographic funnel'' defined in [McDougall et al., 2003](https://journals.ametsoc.org/view/journals/atot/20/5/1520-0426_2003_20_730_aaceaf_2_0_co_2.xml). When using TEOS-10, the Eos uses member methods `calcSLimits(P)` and `calcTLimits(Sa, P)` to calculate the valid ranges of Sa and T given the pressure. The conservative temperature lower bound is set by the freezing temperature, using the member method `calcCtFreezing(Sa, P, SaturationFract)`. This method implements the polynomial approximation of the conservative freezing temperature (called `gsw_ct_freezing_poly` in the GSW package), which is known to produce erros in the (-5e-4 K, 6e-4 K) range. Once we calculate the upper and lower bounds of validity, the state variables are clipped to the valid range (if outside the bounds) before we run the specific volume calculation. The state fields themselves are not changed. ## Removal of Eos diff --git a/components/omega/doc/userGuide/EOS.md b/components/omega/doc/userGuide/EOS.md index 59eec0206d8a..ad1710f1cf82 100644 --- a/components/omega/doc/userGuide/EOS.md +++ b/components/omega/doc/userGuide/EOS.md @@ -21,4 +21,4 @@ where `DRhoDT` is the thermal expansion coefficient ($\textrm{kg}/(\textrm{m}^3 In addition to `SpecVol`, the displaced specific volume `SpecVolDisplaced` is also calculated by the EOS. This calculates the density of a parcel of fluid that is adiabatically displaced by a relative `k` levels, 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 EOS, `SpecVolDisplaced` will be the same as `SpecVol` since the specific volume calculation is independent of pressure/depth. -When using `Teos-10`, the state variables are checked against the range over which the polynomial is considered to be valid (see [Roquet et al. 2015](https://www.sciencedirect.com/science/article/pii/S1463500315000566)). If the values are outside of the accepted values, we use the valid bounds for the specific volume calculation. Note that the state variable values themselves are not modified, only that they are not used as is in the calculation. +When using TEOS-10, the state variables are checked against the range over which the polynomial is considered to be valid (see [Roquet et al. 2015](https://www.sciencedirect.com/science/article/pii/S1463500315000566)). If the values are outside of the accepted values, we use the valid bounds for the specific volume calculation. Note that the state variable values themselves are not modified, only that they are not used as is in the calculation. diff --git a/components/omega/src/ocn/Eos.h b/components/omega/src/ocn/Eos.h index e18102a62c27..704b18a0e070 100644 --- a/components/omega/src/ocn/Eos.h +++ b/components/omega/src/ocn/Eos.h @@ -26,9 +26,9 @@ enum class EosType { Teos10Eos /// Roquet et al. 2015 75 term expansion }; -struct Range { - Real lo; - Real hi; +struct EosRange { + Real Lo; + Real Hi; }; /// TEOS10 75-term Polynomial Equation of State @@ -88,12 +88,12 @@ class Teos10Eos { constexpr Real SAu = 40.0 * 35.16504 / 35.0; constexpr Real CTu = 40.0; constexpr Real DeltaS = 24.0; - Range Srange = calcSLimits(P); - Range Trange = calcTLimits(Sa, P); + EosRange SRange = calcSLimits(P); + EosRange TRange = calcTLimits(Sa, P); Real SaInFunnel = Kokkos::clamp( - Sa, Srange.lo, Srange.hi); // Salt limited to Poly75t valid range + Sa, SRange.Lo, SRange.Hi); // Salt limited to Poly75t valid range Real Ss = Kokkos::sqrt((SaInFunnel + DeltaS) / SAu); - Real Tt = Kokkos::clamp(Ct, Trange.lo, Trange.hi) / CTu; + Real Tt = Kokkos::clamp(Ct, TRange.Lo, TRange.Hi) / CTu; /// Coefficients for the polynomial expansion constexpr Real V000 = 1.0769995862e-03; @@ -249,87 +249,85 @@ class Teos10Eos { } /// Calculate S limits of validity given pressure p - KOKKOS_FUNCTION Range calcSLimits(const Real P) const { - Real lo = 0.0; - Real hi = 42.0; - Real lo2 = P * 5e-3 - 2.5; - Real lo3 = 30.0; + KOKKOS_FUNCTION EosRange calcSLimits(const Real P) const { + Real Lo = 0.0; + Real Hi = 42.0; + Real Lo2 = P * 5e-3 - 2.5; + Real Lo3 = 30.0; if (P >= 500.0 && P < 6500.0) { - lo = Kokkos::max(lo, lo2); + Lo = Kokkos::max(Lo, Lo2); } else if (P >= 6500.0) { - lo = Kokkos::max(lo, lo3); + Lo = Kokkos::max(Lo, Lo3); } - return {lo, hi}; + return {Lo, Hi}; } /// Calculate T limits of validity given p,Sa - KOKKOS_FUNCTION Range calcTLimits(const Real Sa, const Real P) const { - Real lo = -15.0; - Real hi = 95.0; + KOKKOS_FUNCTION EosRange calcTLimits(const Real Sa, const Real P) const { + Real Lo = -15.0; + Real Hi = 95.0; if (P < 500.0) { - lo = calcCtFreezing(Sa, P, Real{0.0}); + Lo = calcCtFreezing(Sa, P, 0.0_Real); } else if (P < 6500.0) { - lo = calcCtFreezing(Sa, Real{500.0}, Real{0.0}); - hi = 31.66666666666667 - P * 3.333333333333334e-3; + Lo = calcCtFreezing(Sa, 500.0_Real, 0.0_Real); + Hi = 31.66666666666667 - P * 3.333333333333334e-3; } else { - lo = calcCtFreezing(Sa, Real{500.0}, Real{0.0}); - hi = 10.0; + Lo = calcCtFreezing(Sa, 500.0_Real, 0.0_Real); + Hi = 10.0; } - return {lo, hi}; + return {Lo, Hi}; } /// Calculates freezing CTemperature (polynomial error in [-5e-4,6e-4] K) KOKKOS_FUNCTION Real calcCtFreezing(const Real Sa, const Real P, const Real SaturationFract) const { constexpr Real Sso = 35.16504; - constexpr Real c0 = 0.017947064327968736; - constexpr Real c1 = -6.076099099929818; - constexpr Real c2 = 4.883198653547851; - constexpr Real c3 = -11.88081601230542; - constexpr Real c4 = 13.34658511480257; - constexpr Real c5 = -8.722761043208607; - constexpr Real c6 = 2.082038908808201; - constexpr Real c7 = -7.389420998107497; - constexpr Real c8 = -2.110913185058476; - constexpr Real c9 = 0.2295491578006229; - constexpr Real c10 = -0.9891538123307282; - constexpr Real c11 = -0.08987150128406496; - constexpr Real c12 = 0.3831132432071728; - constexpr Real c13 = 1.054318231187074; - constexpr Real c14 = 1.065556599652796; - constexpr Real c15 = -0.7997496801694032; - constexpr Real c16 = 0.3850133554097069; - constexpr Real c17 = -2.078616693017569; - constexpr Real c18 = 0.8756340772729538; - constexpr Real c19 = -2.079022768390933; - constexpr Real c20 = 1.596435439942262; - constexpr Real c21 = 0.1338002171109174; - constexpr Real c22 = 1.242891021876471; + constexpr Real C0 = 0.017947064327968736; + constexpr Real C1 = -6.076099099929818; + constexpr Real C2 = 4.883198653547851; + constexpr Real C3 = -11.88081601230542; + constexpr Real C4 = 13.34658511480257; + constexpr Real C5 = -8.722761043208607; + constexpr Real C6 = 2.082038908808201; + constexpr Real C7 = -7.389420998107497; + constexpr Real C8 = -2.110913185058476; + constexpr Real C9 = 0.2295491578006229; + constexpr Real C10 = -0.9891538123307282; + constexpr Real C11 = -0.08987150128406496; + constexpr Real C12 = 0.3831132432071728; + constexpr Real C13 = 1.054318231187074; + constexpr Real C14 = 1.065556599652796; + constexpr Real C15 = -0.7997496801694032; + constexpr Real C16 = 0.3850133554097069; + constexpr Real C17 = -2.078616693017569; + constexpr Real C18 = 0.8756340772729538; + constexpr Real C19 = -2.079022768390933; + constexpr Real C20 = 1.596435439942262; + constexpr Real C21 = 0.1338002171109174; + constexpr Real C22 = 1.242891021876471; // Note: a = 0.502500117621 / Sso - constexpr Real a = 0.014289763856964; - constexpr Real b = 0.057000649899720; - - Real CtFreez; + constexpr Real A = 0.014289763856964; + constexpr Real B = 0.057000649899720; Real Sar = Sa * 1.0e-2; Real X = Kokkos::sqrt(Sar); Real Pr = P * 1.0e-4; - CtFreez = - c0 + Sar * (c1 + X * (c2 + X * (c3 + X * (c4 + X * (c5 + c6 * X))))) + - Pr * (c7 + Pr * (c8 + c9 * Pr)) + + Real CtFreez = + C0 + Sar * (C1 + X * (C2 + X * (C3 + X * (C4 + X * (C5 + C6 * X))))) + + Pr * (C7 + Pr * (C8 + C9 * Pr)) + Sar * Pr * - (c10 + Pr * (c12 + Pr * (c15 + c21 * Sar)) + - Sar * (c13 + c17 * Pr + c19 * Sar) + - X * (c11 + Pr * (c14 + c18 * Pr) + - Sar * (c16 + c20 * Pr + c22 * Sar))); + (C10 + Pr * (C12 + Pr * (C15 + C21 * Sar)) + + Sar * (C13 + C17 * Pr + C19 * Sar) + + X * (C11 + Pr * (C14 + C18 * Pr) + + Sar * (C16 + C20 * Pr + C22 * Sar))); /* Adjust for the effects of dissolved air */ - CtFreez = CtFreez - SaturationFract * (1e-3) * (2.4 - a * Sa) * - (1.0 + b * (1.0 - Sa / Sso)); + CtFreez = CtFreez - SaturationFract * (1e-3) * (2.4 - A * Sa) * + (1.0 + B * (1.0 - Sa / Sso)); return CtFreez; } diff --git a/components/omega/test/ocn/EosTest.cpp b/components/omega/test/ocn/EosTest.cpp index 10d0ffdf4507..659faecfbe4e 100644 --- a/components/omega/test/ocn/EosTest.cpp +++ b/components/omega/test/ocn/EosTest.cpp @@ -38,6 +38,12 @@ constexpr int NVertLayers = 60; /// Published values (TEOS-10 and linear) to test against const Real TeosExpValue = 0.0009732819628; // Expected value for TEOS-10 eos +const Real TeosClampValue1 = + 0.0009714522912320203; // Expected value for TEOS-10 eos clamping test 1 +const Real TeosClampValue2 = + 0.0009662964459162306; // Expected value for TEOS-10 eos clamping test 2 +const Real TeosClampValue3 = + 0.0010086299185825206; // Expected value for TEOS-10 eos clamping test 3 const Real LinearExpValue = 0.0009784735812133072; // Expected value for Linear eos @@ -284,6 +290,117 @@ int testEosTeos10Displaced() { return Err; } +int testEosTeos10Clamping() { + int Err = 0; + const auto *Mesh = HorzMesh::getDefault(); + /// Get Eos instance to test + Eos *TestEos = Eos::getInstance(); + TestEos->EosChoice = EosType::Teos10Eos; + + /// Create and fill ocean state arrays + Array2DReal SArray = Array2DReal("SArray", Mesh->NCellsAll, NVertLevels); + Array2DReal TArray = Array2DReal("TArray", Mesh->NCellsAll, NVertLevels); + Array2DReal PArray = Array2DReal("PArray", Mesh->NCellsAll, NVertLevels); + /// Use Kokkos::deep_copy to fill the entire view with the ref value + /// Test with valid poly75t values first. + deepCopy(SArray, 35.0); + deepCopy(TArray, 5.0); + deepCopy(PArray, 400.0); + deepCopy(TestEos->SpecVol, 0.0); + + /// Compute specific volume + TestEos->computeSpecVol(TArray, SArray, PArray); + + /// Check all array values against expected value + int numMismatches = 0; + Array2DReal SpecVol = TestEos->SpecVol; + parallelReduce( + "CheckSpecVolMatrix-Teos", {Mesh->NCellsAll, NVertLevels}, + KOKKOS_LAMBDA(int i, int j, int &localCount) { + if (!isApprox(SpecVol(i, j), TeosClampValue1, RTol)) { + localCount++; + } + }, + numMismatches); + + auto SpecVolH = createHostMirrorCopy(SpecVol); + if (numMismatches != 0) { + Err++; + LOG_ERROR("EosTest: TEOS SpecVolClampingNone isApprox FAIL, " + "expected {}, got {} with {} mismatches", + TeosClampValue1, SpecVolH(1, 1), numMismatches); + } + if (Err == 0) { + LOG_INFO("EosTest SpecVolClampingNone TEOS-10: PASS"); + } + + /// Test with an ivalid poly75t salinity value second + deepCopy(SArray, 45.0); + deepCopy(TArray, 5.0); + deepCopy(PArray, 400.0); + deepCopy(TestEos->SpecVol, 0.0); + + /// Compute specific volume + TestEos->computeSpecVol(TArray, SArray, PArray); + + /// Check all array values against expected value + numMismatches = 0; + SpecVol = TestEos->SpecVol; + parallelReduce( + "CheckSpecVolMatrix-Teos", {Mesh->NCellsAll, NVertLevels}, + KOKKOS_LAMBDA(int i, int j, int &localCount) { + if (!isApprox(SpecVol(i, j), TeosClampValue2, RTol)) { + localCount++; + } + }, + numMismatches); + + SpecVolH = createHostMirrorCopy(SpecVol); + if (numMismatches != 0) { + Err++; + LOG_ERROR("EosTest: TEOS SpecVolClampingSalt isApprox FAIL, " + "expected {}, got {} with {} mismatches", + TeosClampValue2, SpecVolH(1, 1), numMismatches); + } + if (Err == 0) { + LOG_INFO("EosTest SpecVolClampingSalt TEOS-10: PASS"); + } + + /// Test with an ivalid poly75t temperature value third + deepCopy(SArray, 35.0); + deepCopy(TArray, 100.0); + deepCopy(PArray, 400.0); + deepCopy(TestEos->SpecVol, 0.0); + + /// Compute specific volume + TestEos->computeSpecVol(TArray, SArray, PArray); + + /// Check all array values against expected value + numMismatches = 0; + SpecVol = TestEos->SpecVol; + parallelReduce( + "CheckSpecVolMatrix-Teos", {Mesh->NCellsAll, NVertLevels}, + KOKKOS_LAMBDA(int i, int j, int &localCount) { + if (!isApprox(SpecVol(i, j), TeosClampValue3, RTol)) { + localCount++; + } + }, + numMismatches); + + SpecVolH = createHostMirrorCopy(SpecVol); + if (numMismatches != 0) { + Err++; + LOG_ERROR("EosTest: TEOS SpecVolClampingTemp isApprox FAIL, " + "expected {}, got {} with {} mismatches", + TeosClampValue3, SpecVolH(1, 1), numMismatches); + } + if (Err == 0) { + LOG_INFO("EosTest SpecVolClampingTemp TEOS-10: PASS"); + } + + return Err; +} + /// Finalize and clean up all test infrastructure void finalizeEosTest() { HorzMesh::clear(); @@ -296,8 +413,7 @@ void finalizeEosTest() { /// Test that the external GSW-C library returns the expected specific volume int checkValueGswcSpecVol() { - int Err = 0; - const Real RTol = 1e-10; + int Err = 0; /// Get specific volume from GSW-C library double SpecVol = gsw_specvol(Sa, Ct, P); @@ -317,8 +433,8 @@ int checkValueGswcSpecVol() { /// Test that the external GSW-C library returns the expected specific volume int checkValueCtFreezing() { - int Err = 0; - const Real RTol = 1e-10; + int Err = 0; + Teos10Eos TestEos(1); // TestEos->EosChoice = EosType::Teos10Eos; constexpr Real SaturationFrac = 0.0; @@ -344,14 +460,19 @@ int checkValueCtFreezing() { } // the main test (all in one to have the same log) -// Single value test: -// --> test calls the external GSW-C library +// Single value tests: +// --> one test calls the external GSW-C library for specific volume +// --> next test calls the external GSW-C library for freezing value // and compares the specific volume to the published value // Full array tests: // --> one tests the value on a Eos with linear option // --> next checks the value on a Eos with linear displaced option // --> next checks the value on a Eos with TEOS-10 option // --> next checks the value on a Eos with TEOS-10 displaced option +// Clamping test: +// --> test check the clamping works for values that (1) don't need +// clamping, (2) need claiming on the salinity, and (3) need clamping +// on the temperature int eosTest(const std::string &MeshFile = "OmegaMesh.nc") { int Err = initEosTest(MeshFile); if (Err != 0) { @@ -369,6 +490,9 @@ int eosTest(const std::string &MeshFile = "OmegaMesh.nc") { Err += testEosTeos10(); Err += testEosTeos10Displaced(); + LOG_INFO("Clamping check:"); + Err += testEosTeos10Clamping(); + if (Err == 0) { LOG_INFO("EosTest: Successful completion"); } From 30c56667db3d2b43c69587452675047b97f0628b Mon Sep 17 00:00:00 2001 From: Kat Smith Date: Mon, 6 Oct 2025 09:14:47 -0700 Subject: [PATCH 3/7] adds lower bounds clamping test for temperature --- components/omega/test/ocn/EosTest.cpp | 54 ++++++++++++++++++++++++--- 1 file changed, 48 insertions(+), 6 deletions(-) diff --git a/components/omega/test/ocn/EosTest.cpp b/components/omega/test/ocn/EosTest.cpp index 659faecfbe4e..53123183198d 100644 --- a/components/omega/test/ocn/EosTest.cpp +++ b/components/omega/test/ocn/EosTest.cpp @@ -44,6 +44,8 @@ const Real TeosClampValue2 = 0.0009662964459162306; // Expected value for TEOS-10 eos clamping test 2 const Real TeosClampValue3 = 0.0010086299185825206; // Expected value for TEOS-10 eos clamping test 3 +const Real TeosClampValue4 = + 0.000970887661659709; // Expected value for TEOS-10 eos clamping test 4 const Real LinearExpValue = 0.0009784735812133072; // Expected value for Linear eos @@ -302,7 +304,9 @@ int testEosTeos10Clamping() { Array2DReal TArray = Array2DReal("TArray", Mesh->NCellsAll, NVertLevels); Array2DReal PArray = Array2DReal("PArray", Mesh->NCellsAll, NVertLevels); /// Use Kokkos::deep_copy to fill the entire view with the ref value - /// Test with valid poly75t values first. + /// ------------------------------------- + /// Test with valid poly75t values first (no clamping) + /// ------------------------------------- deepCopy(SArray, 35.0); deepCopy(TArray, 5.0); deepCopy(PArray, 400.0); @@ -334,8 +338,10 @@ int testEosTeos10Clamping() { LOG_INFO("EosTest SpecVolClampingNone TEOS-10: PASS"); } + /// ------------------------------------- /// Test with an ivalid poly75t salinity value second - deepCopy(SArray, 45.0); + /// ------------------------------------- + deepCopy(SArray, 45.0); // Clamping happens at 42.0 deepCopy(TArray, 5.0); deepCopy(PArray, 400.0); deepCopy(TestEos->SpecVol, 0.0); @@ -366,9 +372,11 @@ int testEosTeos10Clamping() { LOG_INFO("EosTest SpecVolClampingSalt TEOS-10: PASS"); } - /// Test with an ivalid poly75t temperature value third + /// ------------------------------------- + /// Test with an ivalid high poly75t temperature value third + /// ------------------------------------- deepCopy(SArray, 35.0); - deepCopy(TArray, 100.0); + deepCopy(TArray, 100.0); // Clamping happens at 95.0 deepCopy(PArray, 400.0); deepCopy(TestEos->SpecVol, 0.0); @@ -390,12 +398,46 @@ int testEosTeos10Clamping() { SpecVolH = createHostMirrorCopy(SpecVol); if (numMismatches != 0) { Err++; - LOG_ERROR("EosTest: TEOS SpecVolClampingTemp isApprox FAIL, " + LOG_ERROR("EosTest: TEOS SpecVolClampingTempHi isApprox FAIL, " "expected {}, got {} with {} mismatches", TeosClampValue3, SpecVolH(1, 1), numMismatches); } if (Err == 0) { - LOG_INFO("EosTest SpecVolClampingTemp TEOS-10: PASS"); + LOG_INFO("EosTest SpecVolClampingTempHi TEOS-10: PASS"); + } + + /// ------------------------------------- + /// Test with an ivalid low poly75t temperature value last + /// ------------------------------------- + deepCopy(SArray, 35.0); + deepCopy(TArray, -5.0); // Clamping happens at -2.2160968103069774 + deepCopy(PArray, 400.0); + deepCopy(TestEos->SpecVol, 0.0); + + /// Compute specific volume + TestEos->computeSpecVol(TArray, SArray, PArray); + + /// Check all array values against expected value + numMismatches = 0; + SpecVol = TestEos->SpecVol; + parallelReduce( + "CheckSpecVolMatrix-Teos", {Mesh->NCellsAll, NVertLevels}, + KOKKOS_LAMBDA(int i, int j, int &localCount) { + if (!isApprox(SpecVol(i, j), TeosClampValue4, RTol)) { + localCount++; + } + }, + numMismatches); + + SpecVolH = createHostMirrorCopy(SpecVol); + if (numMismatches != 0) { + Err++; + LOG_ERROR("EosTest: TEOS SpecVolClampingTempLo isApprox FAIL, " + "expected {}, got {} with {} mismatches", + TeosClampValue4, SpecVolH(1, 1), numMismatches); + } + if (Err == 0) { + LOG_INFO("EosTest SpecVolClampingTempLo TEOS-10: PASS"); } return Err; From 3f00564d497d84d78a131718242f9f76a807574e Mon Sep 17 00:00:00 2001 From: Alice Barthel Date: Tue, 11 Nov 2025 10:09:49 -0800 Subject: [PATCH 4/7] switched clamping to an option, and eos limits has 2 options --- components/omega/configs/Default.yml | 2 + components/omega/src/ocn/Eos.cpp | 19 +++++- components/omega/src/ocn/Eos.h | 97 +++++++++++++++++---------- components/omega/test/ocn/EosTest.cpp | 14 ++-- 4 files changed, 90 insertions(+), 42 deletions(-) diff --git a/components/omega/configs/Default.yml b/components/omega/configs/Default.yml index b2b1a7c37b16..0c639edea9fa 100644 --- a/components/omega/configs/Default.yml +++ b/components/omega/configs/Default.yml @@ -59,6 +59,8 @@ Omega: DRhoDT: -0.2 DRhoDS: 0.8 RhoT0S0: 1000.0 + ClampingEnable: false + EosLimits: Funnel IOStreams: InitialVertCoord: UsePointerFile: false diff --git a/components/omega/src/ocn/Eos.cpp b/components/omega/src/ocn/Eos.cpp index b8d1751c24c3..d59f8c59a6ec 100644 --- a/components/omega/src/ocn/Eos.cpp +++ b/components/omega/src/ocn/Eos.cpp @@ -81,6 +81,12 @@ void Eos::init() { Err += EosConfig.get("EosType", EosTypeStr); CHECK_ERROR_ABORT(Err, "Eos::init: EosType subgroup not found in EosConfig"); + std::string EosLimStr; + Err += EosConfig.get("EosLimits", EosLimStr); + CHECK_ERROR_ABORT(Err, + "Eos::init: EosLimits subgroup not found in EosConfig"); + /// EosLimits EosLimChoice; + /// Set EosChoice and parameters based on EosTypeStr if (EosTypeStr == "Linear" or EosTypeStr == "linear") { Config EosLinConfig("Linear"); @@ -104,12 +110,23 @@ void Eos::init() { } else if ((EosTypeStr == "teos10") or (EosTypeStr == "teos-10") or (EosTypeStr == "TEOS-10")) { eos->EosChoice = EosType::Teos10Eos; + if (EosLimStr == "Funnel") { + eos->ComputeSpecVolTeos10.EosLimChoice = EosLimits::Funnel; + } else if (EosLimStr == "Cube") { + eos->ComputeSpecVolTeos10.EosLimChoice = EosLimits::Cube; + } else { + LOG_ERROR("Eos::init: Unknown EosLimits requested"); + Err += EosConfig.get("ClampingEnable", + eos->ComputeSpecVolTeos10.ClampingEnable); + CHECK_ERROR_ABORT( + Err, "Eos::init: Parameter ClampingEnable not found in EosConfig"); + } } else { LOG_ERROR("Eos::init: Unknown EosType requested"); } } // end init -/// Compute specific volume for all cells/layers (no displacement) +/// Compute specific volume for all cells/levels (no displacement) void Eos::computeSpecVol(const Array2DReal &ConservTemp, const Array2DReal &AbsSalinity, const Array2DReal &Pressure) { diff --git a/components/omega/src/ocn/Eos.h b/components/omega/src/ocn/Eos.h index 704b18a0e070..78c4894c1a15 100644 --- a/components/omega/src/ocn/Eos.h +++ b/components/omega/src/ocn/Eos.h @@ -26,6 +26,11 @@ enum class EosType { Teos10Eos /// Roquet et al. 2015 75 term expansion }; +enum class EosLimits { + Funnel, /// + Cube /// +}; + struct EosRange { Real Lo; Real Hi; @@ -35,6 +40,8 @@ struct EosRange { class Teos10Eos { public: Array2DReal SpecVolPCoeffs; + EosLimits EosLimChoice; ///< EOS clamping option in use + bool ClampingEnable; ///< EOS clamping option in use /// constructor declaration Teos10Eos(int NVertLayers); @@ -90,10 +97,14 @@ class Teos10Eos { constexpr Real DeltaS = 24.0; EosRange SRange = calcSLimits(P); EosRange TRange = calcTLimits(Sa, P); - Real SaInFunnel = Kokkos::clamp( - Sa, SRange.Lo, SRange.Hi); // Salt limited to Poly75t valid range - Real Ss = Kokkos::sqrt((SaInFunnel + DeltaS) / SAu); - Real Tt = Kokkos::clamp(Ct, TRange.Lo, TRange.Hi) / CTu; + Real Ss = Kokkos::sqrt((Sa + DeltaS) / SAu); + Real Tt = Ct / CTu; + if (ClampingEnable) { + Real SaInRange = Kokkos::clamp( + Sa, SRange.Lo, SRange.Hi); // Salt limited to Poly75t valid range + Ss = Kokkos::sqrt((SaInRange + DeltaS) / SAu); + Tt = Kokkos::clamp(Ct, TRange.Lo, TRange.Hi) / CTu; + } /// Coefficients for the polynomial expansion constexpr Real V000 = 1.0769995862e-03; @@ -215,10 +226,15 @@ class Teos10Eos { KOKKOS_FUNCTION Real calcDelta(const Array2DReal &SpecVolPCoeffs, const I4 K, const Real P) const { - constexpr Real Pu = 1e4; - constexpr Real Pmax = 8000.0; - Real Pp = Kokkos::min(P, Pmax) / Pu; // P limited to Poly75t valid range - + constexpr Real Pu = 1e4; + const Real Pmax = (EosLimChoice == EosLimits::Funnel) ? 8000.0 : 10000.0; + Real Pp = P / Pu; + if (ClampingEnable) { + Pp = Kokkos::min(P, Pmax) / Pu; // P limited to Poly75t valid range + LOG_INFO("P={} exceeds Pmax={}; Clamping the pressure", P, Pmax); + } else if (P > Pmax) { + LOG_WARN("P={} exceeds Pmax={}", P, Pmax); + } Real Delta = ((((SpecVolPCoeffs(5, K) * Pp + SpecVolPCoeffs(4, K)) * Pp + SpecVolPCoeffs(3, K)) * Pp + @@ -232,15 +248,20 @@ class Teos10Eos { /// Calculate reference profile for TEOS-10 KOKKOS_FUNCTION Real calcRefProfile(const Real P) const { - constexpr Real Pu = 1e4; - 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; - constexpr Real Pmax = 8000.0; - Real Pp = Kokkos::min(P, Pmax) / Pu; // P limited to Poly75t valid range + constexpr Real Pu = 1e4; + 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; + const Real Pmax = (EosLimChoice == EosLimits::Funnel) ? 8000.0 : 10000.0; + Real Pp = P / Pu; + if (ClampingEnable) { + Pp = Kokkos::min(P, Pmax) / Pu; // P limited to Poly75t valid range + } else if (P > Pmax) { + LOG_WARN("P={} exceeds Pmax={}", P, Pmax); + } Real V0 = (((((V05 * Pp + V04) * Pp + V03) * Pp + V02) * Pp + V01) * Pp + V00) * @@ -250,15 +271,17 @@ class Teos10Eos { /// Calculate S limits of validity given pressure p KOKKOS_FUNCTION EosRange calcSLimits(const Real P) const { - Real Lo = 0.0; - Real Hi = 42.0; - Real Lo2 = P * 5e-3 - 2.5; - Real Lo3 = 30.0; - - if (P >= 500.0 && P < 6500.0) { - Lo = Kokkos::max(Lo, Lo2); - } else if (P >= 6500.0) { - Lo = Kokkos::max(Lo, Lo3); + Real Lo = 0.0; + Real Hi = 42.0; + if (EosLimChoice == EosLimits::Funnel) { + Real Lo2 = P * 5e-3 - 2.5; + Real Lo3 = 30.0; + + if (P >= 500.0 && P < 6500.0) { + Lo = Kokkos::max(Lo, Lo2); + } else if (P >= 6500.0) { + Lo = Kokkos::max(Lo, Lo3); + } } return {Lo, Hi}; } @@ -268,14 +291,20 @@ class Teos10Eos { Real Lo = -15.0; Real Hi = 95.0; - if (P < 500.0) { - Lo = calcCtFreezing(Sa, P, 0.0_Real); - } else if (P < 6500.0) { - Lo = calcCtFreezing(Sa, 500.0_Real, 0.0_Real); - Hi = 31.66666666666667 - P * 3.333333333333334e-3; - } else { - Lo = calcCtFreezing(Sa, 500.0_Real, 0.0_Real); - Hi = 10.0; + if (EosLimChoice == EosLimits::Funnel) { + if (P < 500.0) { + Lo = calcCtFreezing(Sa, P, 0.0_Real); + } else if (P < 6500.0) { + Lo = calcCtFreezing(Sa, 500.0_Real, 0.0_Real); + Hi = 31.66666666666667 - P * 3.333333333333334e-3; + } else { + Lo = calcCtFreezing(Sa, 500.0_Real, 0.0_Real); + Hi = 10.0; + } + } + if (EosLimChoice == EosLimits::Cube) { + Lo = -2.0; + Hi = 40.0; } return {Lo, Hi}; } diff --git a/components/omega/test/ocn/EosTest.cpp b/components/omega/test/ocn/EosTest.cpp index 53123183198d..2968ef6cedb2 100644 --- a/components/omega/test/ocn/EosTest.cpp +++ b/components/omega/test/ocn/EosTest.cpp @@ -300,9 +300,9 @@ int testEosTeos10Clamping() { TestEos->EosChoice = EosType::Teos10Eos; /// Create and fill ocean state arrays - Array2DReal SArray = Array2DReal("SArray", Mesh->NCellsAll, NVertLevels); - Array2DReal TArray = Array2DReal("TArray", Mesh->NCellsAll, NVertLevels); - Array2DReal PArray = Array2DReal("PArray", Mesh->NCellsAll, NVertLevels); + Array2DReal SArray = Array2DReal("SArray", Mesh->NCellsAll, NVertLayers); + Array2DReal TArray = Array2DReal("TArray", Mesh->NCellsAll, NVertLayers); + Array2DReal PArray = Array2DReal("PArray", Mesh->NCellsAll, NVertLayers); /// Use Kokkos::deep_copy to fill the entire view with the ref value /// ------------------------------------- /// Test with valid poly75t values first (no clamping) @@ -319,7 +319,7 @@ int testEosTeos10Clamping() { int numMismatches = 0; Array2DReal SpecVol = TestEos->SpecVol; parallelReduce( - "CheckSpecVolMatrix-Teos", {Mesh->NCellsAll, NVertLevels}, + "CheckSpecVolMatrix-Teos", {Mesh->NCellsAll, NVertLayers}, KOKKOS_LAMBDA(int i, int j, int &localCount) { if (!isApprox(SpecVol(i, j), TeosClampValue1, RTol)) { localCount++; @@ -353,7 +353,7 @@ int testEosTeos10Clamping() { numMismatches = 0; SpecVol = TestEos->SpecVol; parallelReduce( - "CheckSpecVolMatrix-Teos", {Mesh->NCellsAll, NVertLevels}, + "CheckSpecVolMatrix-Teos", {Mesh->NCellsAll, NVertLayers}, KOKKOS_LAMBDA(int i, int j, int &localCount) { if (!isApprox(SpecVol(i, j), TeosClampValue2, RTol)) { localCount++; @@ -387,7 +387,7 @@ int testEosTeos10Clamping() { numMismatches = 0; SpecVol = TestEos->SpecVol; parallelReduce( - "CheckSpecVolMatrix-Teos", {Mesh->NCellsAll, NVertLevels}, + "CheckSpecVolMatrix-Teos", {Mesh->NCellsAll, NVertLayers}, KOKKOS_LAMBDA(int i, int j, int &localCount) { if (!isApprox(SpecVol(i, j), TeosClampValue3, RTol)) { localCount++; @@ -421,7 +421,7 @@ int testEosTeos10Clamping() { numMismatches = 0; SpecVol = TestEos->SpecVol; parallelReduce( - "CheckSpecVolMatrix-Teos", {Mesh->NCellsAll, NVertLevels}, + "CheckSpecVolMatrix-Teos", {Mesh->NCellsAll, NVertLayers}, KOKKOS_LAMBDA(int i, int j, int &localCount) { if (!isApprox(SpecVol(i, j), TeosClampValue4, RTol)) { localCount++; From 4ee8e9a5cf3edd257598459f86eee3ef00a99327 Mon Sep 17 00:00:00 2001 From: Alice Barthel Date: Thu, 13 Nov 2025 10:30:29 -0800 Subject: [PATCH 5/7] added a setClamping method to toggle the Clamping option for the relevant tests --- components/omega/src/ocn/Eos.h | 4 ++++ components/omega/test/ocn/EosTest.cpp | 1 + 2 files changed, 5 insertions(+) diff --git a/components/omega/src/ocn/Eos.h b/components/omega/src/ocn/Eos.h index 78c4894c1a15..333d9f10b46e 100644 --- a/components/omega/src/ocn/Eos.h +++ b/components/omega/src/ocn/Eos.h @@ -46,6 +46,9 @@ class Teos10Eos { /// constructor declaration Teos10Eos(int NVertLayers); + KOKKOS_INLINE_FUNCTION + void setClamping(bool Flag) { ClampingEnable = Flag; } + // 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 @@ -402,6 +405,7 @@ class Eos { /// Destroy instance (frees Kokkos views) static void destroyInstance(); + void setTeosClamping(bool Flag) { ComputeSpecVolTeos10.setClamping(Flag); } EosType EosChoice; ///< Current EOS type in use Array2DReal SpecVol; ///< Specific volume field Array2DReal SpecVolDisplaced; ///< Displaced specific volume field diff --git a/components/omega/test/ocn/EosTest.cpp b/components/omega/test/ocn/EosTest.cpp index 2968ef6cedb2..13e83e439424 100644 --- a/components/omega/test/ocn/EosTest.cpp +++ b/components/omega/test/ocn/EosTest.cpp @@ -298,6 +298,7 @@ int testEosTeos10Clamping() { /// Get Eos instance to test Eos *TestEos = Eos::getInstance(); TestEos->EosChoice = EosType::Teos10Eos; + TestEos->setTeosClamping(true); /// Create and fill ocean state arrays Array2DReal SArray = Array2DReal("SArray", Mesh->NCellsAll, NVertLayers); From 0ab6d7c0378feaad60deb1140e238ea905ba24f2 Mon Sep 17 00:00:00 2001 From: Alice Barthel Date: Mon, 17 Nov 2025 09:00:45 -0800 Subject: [PATCH 6/7] updated the documentation and added T,S warnings --- components/omega/doc/devGuide/EOS.md | 2 +- components/omega/doc/userGuide/EOS.md | 6 ++++-- components/omega/src/ocn/Eos.h | 20 +++++++++++++++++--- 3 files changed, 22 insertions(+), 6 deletions(-) diff --git a/components/omega/doc/devGuide/EOS.md b/components/omega/doc/devGuide/EOS.md index 5354062eaa29..6754d27fcf9d 100644 --- a/components/omega/doc/devGuide/EOS.md +++ b/components/omega/doc/devGuide/EOS.md @@ -54,7 +54,7 @@ For example, to displace each level to one below, set `KDisp = 1`. ## Bounds check (and truncation) for the state variables (under TEOS-10) -The implemented 75-term polynomial for the calculation of the specific volume under TEOS-10 is considered valid for ocean states contained in the ''oceanographic funnel'' defined in [McDougall et al., 2003](https://journals.ametsoc.org/view/journals/atot/20/5/1520-0426_2003_20_730_aaceaf_2_0_co_2.xml). When using TEOS-10, the Eos uses member methods `calcSLimits(P)` and `calcTLimits(Sa, P)` to calculate the valid ranges of Sa and T given the pressure. The conservative temperature lower bound is set by the freezing temperature, using the member method `calcCtFreezing(Sa, P, SaturationFract)`. This method implements the polynomial approximation of the conservative freezing temperature (called `gsw_ct_freezing_poly` in the GSW package), which is known to produce erros in the (-5e-4 K, 6e-4 K) range. Once we calculate the upper and lower bounds of validity, the state variables are clipped to the valid range (if outside the bounds) before we run the specific volume calculation. The state fields themselves are not changed. +The implemented 75-term polynomial for the calculation of the specific volume under TEOS-10 has been evaluated for ocean states in the ''cube" (-2-40 C ;0-42 g/kg; 0-10,000 dbar) and the ''oceanographic funnel'' defined in [McDougall et al., 2003](https://journals.ametsoc.org/view/journals/atot/20/5/1520-0426_2003_20_730_aaceaf_2_0_co_2.xml). When using TEOS-10, the Eos uses member methods `calcSLimits(P)` and `calcTLimits(Sa, P)` to calculate the valid ranges of Sa and T. When using the `Funnel` opion of `EosLimits`, the salinity limits are calculated as a function of pressure and the temperature limits as a function of pressure and salinity. The conservative temperature lower bound is set by the freezing temperature, using the member method `calcCtFreezing(Sa, P, SaturationFract)`. This method implements the polynomial approximation of the conservative freezing temperature (called `gsw_ct_freezing_poly` in the GSW package), which is known to produce erros in the (-5e-4 K, 6e-4 K) range. When using the `Cube` option of `EosLimits`, the bounds are constant. Once we calculate the upper and lower bounds of validity, warnings are issued when the state variables are outside the validity bounds. If `ClampingEnable` is true, the state variables are clipped to the valid range (if outside the bounds) before we run the specific volume calculation. The state fields themselves are not changed. If `ClampingEnable` is false, the warnings are issued but the specific volume is calculated based on the unchanged state variables. ## Removal of Eos diff --git a/components/omega/doc/userGuide/EOS.md b/components/omega/doc/userGuide/EOS.md index ad1710f1cf82..662dbc7f88e4 100644 --- a/components/omega/doc/userGuide/EOS.md +++ b/components/omega/doc/userGuide/EOS.md @@ -15,10 +15,12 @@ Eos: DRhoDT: -0.2 DRhoDS: 0.8 RhoT0S0: 1000.0 + ClampingEnable: false + EosLimits: Funnel ``` -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$). +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$). The `ClampingEnable` flag restricts the state variables to the ranges specified by the `EosLimits` options if turned on, or it issues a warning if off. Note that `ClampingEnable` and `EosLimits` are currently used only in the Teos-10 option. In addition to `SpecVol`, the displaced specific volume `SpecVolDisplaced` is also calculated by the EOS. This calculates the density of a parcel of fluid that is adiabatically displaced by a relative `k` levels, 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 EOS, `SpecVolDisplaced` will be the same as `SpecVol` since the specific volume calculation is independent of pressure/depth. -When using TEOS-10, the state variables are checked against the range over which the polynomial is considered to be valid (see [Roquet et al. 2015](https://www.sciencedirect.com/science/article/pii/S1463500315000566)). If the values are outside of the accepted values, we use the valid bounds for the specific volume calculation. Note that the state variable values themselves are not modified, only that they are not used as is in the calculation. +When using TEOS-10, the state variables are checked against the range over which the polynomial is considered to be valid (see [Roquet et al. 2015](https://www.sciencedirect.com/science/article/pii/S1463500315000566)). The possible ranges are refered to as `Funnel` (by default) or `Cube` (wider range). If the values are outside of the accepted values, the code issues a warning. When `ClampingEnable=true`, the code uses the valid bounds for the specific volume calculation. Note that in that case, the state variable values themselves are not modified, only that they are not used as is in the calculation. diff --git a/components/omega/src/ocn/Eos.h b/components/omega/src/ocn/Eos.h index 333d9f10b46e..1d0f0a3482cd 100644 --- a/components/omega/src/ocn/Eos.h +++ b/components/omega/src/ocn/Eos.h @@ -109,6 +109,19 @@ class Teos10Eos { Tt = Kokkos::clamp(Ct, TRange.Lo, TRange.Hi) / CTu; } + if (Sa > SRange.Hi) { + LOG_WARN("In calcPCoeffs: S={} exceeds Smax={}", Sa, SRange.Hi); + } + if (Sa < SRange.Lo) { + LOG_WARN("In calcPCoeffs: S={} lower than Smin={}", Sa, SRange.Lo); + } + if (Ct > TRange.Hi) { + LOG_WARN("In calcPCoeffs: T={} exceeds Tmax={}", Ct, TRange.Hi); + } + if (Ct < TRange.Lo) { + LOG_WARN("In calcPCoeffs: T={} lower than Tmin={}", Ct, TRange.Lo); + } + /// Coefficients for the polynomial expansion constexpr Real V000 = 1.0769995862e-03; constexpr Real V100 = -3.1038981976e-04; @@ -234,9 +247,10 @@ class Teos10Eos { Real Pp = P / Pu; if (ClampingEnable) { Pp = Kokkos::min(P, Pmax) / Pu; // P limited to Poly75t valid range - LOG_INFO("P={} exceeds Pmax={}; Clamping the pressure", P, Pmax); + LOG_INFO("In calcDelta: P={} exceeds Pmax={}; Clamping the pressure", + P, Pmax); } else if (P > Pmax) { - LOG_WARN("P={} exceeds Pmax={}", P, Pmax); + LOG_WARN("In calcDelta: P={} exceeds Pmax={}", P, Pmax); } Real Delta = ((((SpecVolPCoeffs(5, K) * Pp + SpecVolPCoeffs(4, K)) * Pp + SpecVolPCoeffs(3, K)) * @@ -263,7 +277,7 @@ class Teos10Eos { if (ClampingEnable) { Pp = Kokkos::min(P, Pmax) / Pu; // P limited to Poly75t valid range } else if (P > Pmax) { - LOG_WARN("P={} exceeds Pmax={}", P, Pmax); + LOG_WARN("In calcRefProfile: P={} exceeds Pmax={}", P, Pmax); } Real V0 = From 3a84e53a5525b2130a4d04b3cd457e0e2f1b6900 Mon Sep 17 00:00:00 2001 From: Alice Barthel Date: Mon, 17 Nov 2025 10:48:28 -0800 Subject: [PATCH 7/7] added documentation on internal helper functions for developers --- components/omega/doc/devGuide/EOS.md | 28 ++++++++++++++++++++++++---- 1 file changed, 24 insertions(+), 4 deletions(-) diff --git a/components/omega/doc/devGuide/EOS.md b/components/omega/doc/devGuide/EOS.md index 6754d27fcf9d..2af9e6a2733d 100644 --- a/components/omega/doc/devGuide/EOS.md +++ b/components/omega/doc/devGuide/EOS.md @@ -35,10 +35,10 @@ OMEGA::Eos* DefEos = OMEGA::Eos::getInstance(); ## Computation of Eos -To compute `SpecVol` for a particular set of temperature, salinity, and pressure arrays, do +To compute `SpecVol` for a particular set of temperature `Ct`, salinity `Sa`, and pressure `P` arrays, do ```c++ -Eos.computeSpecVol(ConsrvTemp, AbsSalinity, Pressure); +Eos.computeSpecVol(Ct, Sa, P); ``` `SpecVolDisplaced` is calculated using local temperature and salinity values, but a pressure @@ -46,7 +46,7 @@ value at `K + KDisp`. To compute `SpecVolDisplaced` for a particular set of temp and pressure arrays and displaced vertical index level, do ```c++ -Eos.computeSpecVolDisp(ConsrvTemp, AbsSalinity, Pressure, KDisp); +Eos.computeSpecVolDisp(Ct, Sa, P, KDisp); ``` where `KDisp` is the number of `k` levels you want to displace each specific volume level to. @@ -54,7 +54,27 @@ For example, to displace each level to one below, set `KDisp = 1`. ## Bounds check (and truncation) for the state variables (under TEOS-10) -The implemented 75-term polynomial for the calculation of the specific volume under TEOS-10 has been evaluated for ocean states in the ''cube" (-2-40 C ;0-42 g/kg; 0-10,000 dbar) and the ''oceanographic funnel'' defined in [McDougall et al., 2003](https://journals.ametsoc.org/view/journals/atot/20/5/1520-0426_2003_20_730_aaceaf_2_0_co_2.xml). When using TEOS-10, the Eos uses member methods `calcSLimits(P)` and `calcTLimits(Sa, P)` to calculate the valid ranges of Sa and T. When using the `Funnel` opion of `EosLimits`, the salinity limits are calculated as a function of pressure and the temperature limits as a function of pressure and salinity. The conservative temperature lower bound is set by the freezing temperature, using the member method `calcCtFreezing(Sa, P, SaturationFract)`. This method implements the polynomial approximation of the conservative freezing temperature (called `gsw_ct_freezing_poly` in the GSW package), which is known to produce erros in the (-5e-4 K, 6e-4 K) range. When using the `Cube` option of `EosLimits`, the bounds are constant. Once we calculate the upper and lower bounds of validity, warnings are issued when the state variables are outside the validity bounds. If `ClampingEnable` is true, the state variables are clipped to the valid range (if outside the bounds) before we run the specific volume calculation. The state fields themselves are not changed. If `ClampingEnable` is false, the warnings are issued but the specific volume is calculated based on the unchanged state variables. +The implemented 75-term polynomial for the calculation of the specific volume under TEOS-10 has been evaluated for ocean states in the ''cube" (-2-40 C; 0-42 g/kg; 0-10,000 dbar) and the ''oceanographic funnel'' defined in [McDougall et al., 2003](https://journals.ametsoc.org/view/journals/atot/20/5/1520-0426_2003_20_730_aaceaf_2_0_co_2.xml). When using TEOS-10, the Eos uses member methods `calcSLimits(P)` and `calcTLimits(Sa, P)` to calculate the valid ranges of Sa and Ct. When using the `Funnel` opion of `EosLimits`, the salinity limits are calculated as a function of pressure and the temperature limits as a function of pressure and salinity. The conservative temperature lower bound is set by the freezing temperature, using the member method `calcCtFreezing(Sa, P, SaturationFract)`. This method implements the polynomial approximation of the conservative freezing temperature (called `gsw_ct_freezing_poly` in the GSW package), which is known to produce errors in the (-5e-4 K, 6e-4 K) range. When using the `Cube` option of `EosLimits`, the bounds are constant. Once we calculate the upper and lower bounds of validity, warnings are issued when the state variables are outside the validity bounds. If `ClampingEnable` is true, the state variables are clipped to the valid range (if outside the bounds) before we run the specific volume calculation. The state fields themselves are not changed. If `ClampingEnable` is false, the warnings are issued but the specific volume is calculated based on the unchanged state variables. + +## Underlying helper functions when using TEOS-10 + +The computation of the TEOS-10 specific volume relies on 3 underlying member functions: + +```c++ +calcPCoeffs(SpecVolPCoeffs, K, Ct, Sa, P); +``` +which calculates the coefficents that will be applied to the pressure. These coefficients are dependent on the temperature and salinity variables but not the pressure. The pressure argument present in the function call is only used to set the bounds of Ct,Sa validity. This function is where the bulk of the polynomial calculation takes place and thus it is advised to reuse the `SpecVolPCoeffs` which are a class data member if possible to reduce computational expense. + +```c++ +calcRefProfile(Pressure); +``` +which calculates the ocean reference profile as a function of pressure in the layers. This is a simple 6th order polynomial in P with constant coefficients. It could be reused within a timestep provided that the layer pressures do not change but is cheap to recalculate. + +```c++ +calcDelta(SpecVolPCoeffs, K, Pressure); +``` +which applies the `SpecVolPCoeffs` calculated above to the pressure state variable. The resulting "delta" in specific volume is added to the reference profile calculated above to form the specific volume. This step is a simple 5th order polynomial in pressure, which combines the effects of T,S (pre-calculated in `calcPCoeffs`) and the pressure. + ## Removal of Eos