adds test for eos clamping

Author: katsmith133

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
 

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.

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;
    }

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");
    }