added a new Eos constant option and ctest

components/omega/doc/devGuide/EOS.md

@@ -3,10 +3,13 @@
 # Equation of State (EOS)
 
 Omega includes an `Eos` class that provides functions that compute `SpecVol`, `SpecVolDisplaced`,
-and `BruntVaisalaFreqSq`. Current EOS options are a linear EOS or an EOS computed using the TEOS-10
-75 term expansion from [Roquet et al. 2015](https://www.sciencedirect.com/science/article/pii/S1463500315000566).
-If `SpecVolDisplaced` is calculated with the linear EOS option, it will be equal to `SpecVol` as there
-is no pressure/depth dependence for the linear EOS. `SpecVolDisplaced` computes specific volume
+and `BruntVaisalaFreqSq`. Current EOS options are a linear EOS, a constant EOS,
+or an EOS computed using the TEOS-10 75 term expansion from
+[Roquet et al. 2015](https://www.sciencedirect.com/science/article/pii/S1463500315000566).
+If `SpecVolDisplaced` is calculated with the linear or constant EOS option,
+it will be equal to `SpecVol` as there is no pressure/depth dependence for
+those EOS options. For the constant EOS option, `SpecVol` is set to `1/RhoSw`
+for all active cells/layers. `SpecVolDisplaced` computes specific volume
 adiabatically displaced to `K + KDisp` (where `K` counted positive downward, ie `K+1` is one layer below `K`). Note: `SpecVol` must be calculated before `BruntVaisalaFreqSq`, as
 `SpecVol` is an input for the `BruntVaisalaFreqSq` calculation. If the linear EOS option is used, then the `BruntVaisalaFreqSq`
 is calculated using linear coefficients. If the TEOS-10 option is used, the `BruntVaisalaFreqSq` is calculated with non-linear
@@ -15,14 +18,15 @@ for the `BruntVaisalaFreqSq` TEOS-10 option that differ from how it is calculate
 (1) gravity is assumed to be constant and not a function of depth and latitude, and (2) the interface value of the specific volume is
 calculated as the average between two layer values, rather than being recalculated using the interface values of temperature,
 salinity, and pressure. Both of these assumptions incur less than a 1% error.
+For the constant EOS option, `BruntVaisalaFreqSq` is identically zero.
 
 ## Eos type
 
 An enumeration listing all implemented schemes is provided. It needs to be extended every time an
 EOS is added. It is used to identify which EOS method is to be used at run time.
 
 ```c++
-enum class EosType { LinearEos, Teos10Eos };
+enum class EosType { LinearEos, Teos10Eos, ConstantEos };
 ```
 
 ## Initialization
@@ -72,3 +76,4 @@ To clear the Eos instance do:
 
 ```c++
 OMEGA::Eos::destroyInstance();
+```

components/omega/doc/userGuide/EOS.md

@@ -4,13 +4,13 @@
 
 The equation of state (EOS) for the ocean describes the relationship between specific volume of seawater (in $\textrm{m}^3/\textrm{kg}$; the reciprocal of density) and temperature (in $^{\circ}\textrm{C}$), salinity (in $\textrm{g/kg}$), and pressure (in $\textrm{dbar}$). Through the hydrostatic balance (which relates density/specific volume gradients to pressure gradients), the equation of state provides a connection between active tracers (temperature and salinity) and the fluid dynamics.
 
-Two choices of EOS are provided by Omega: a linear EOS and a TEOS-10 EOS. The linear EOS simplifies the relationship by excluding the influence of pressure and using constant expansion/contraction coefficients, making the specific volume a simple linear function of temperature and salinity. However, this option should only be used for simpler idealized test cases as it is not appropriate for realistic ocean simulations. The TEOS-10 EOS is a 75-term polynomial expression from [Roquet et al. 2015](https://www.sciencedirect.com/science/article/pii/S1463500315000566) that approximates the [Thermodynamic Equation of Seawater 2010](https://www.teos-10.org/pubs/TEOS-10_Manual.pdf), but in a less complex and more computationally efficient manner, and is the preferred EOS for real ocean simulations in Omega.
+Three choices of EOS are provided by Omega: a linear EOS, a constant EOS, and a TEOS-10 EOS. The linear EOS simplifies the relationship by excluding the influence of pressure and using constant expansion/contraction coefficients, making the specific volume a simple linear function of temperature and salinity. The constant EOS is intended for idealized/debug-style configurations and sets the specific volume to `1/RhoSw` for all active cells/layers (with `RhoSw` set to `1026.0` in `GlobalConstants.h`). However, both simplified options should only be used for simpler idealized test cases and are not appropriate for realistic ocean simulations. The TEOS-10 EOS is a 75-term polynomial expression from [Roquet et al. 2015](https://www.sciencedirect.com/science/article/pii/S1463500315000566) that approximates the [Thermodynamic Equation of Seawater 2010](https://www.teos-10.org/pubs/TEOS-10_Manual.pdf), but in a less complex and more computationally efficient manner, and is the preferred EOS for real ocean simulations in Omega.
 
-The user-configurable options are: `EosType` (choose either `Linear` or `Teos-10`), as well as the parameters needed for the linear EOS.
+The user-configurable option `EosType` can be set to `linear`, `teos10`, or `constant`. Parameters in the `Linear` subsection are used only when `EosType` is set to `linear`.
 
 ```yaml
 Eos:
-   EosType : teos10
+   EosType: teos10
    Linear:
       DRhoDT: -0.2
       DRhoDS: 0.8
@@ -23,8 +23,8 @@ In addition to `SpecVol`, the displaced specific volume `SpecVolDisplaced` and `
 
 ## 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` EOS option, `SpecVolDisplaced` will be the same as `SpecVol` since the specific volume calculation is independent of pressure/depth.
+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.
 
 ## Squared Brunt Vaisala Frequency
 
-The `Eos` class also calculates the squared Brunt Vaisala Frequency, which is used by the `VertMix` class to determine water column stability for both the convective adjustment and shear-instability driven mixing. When using the `Linear` EOS option, the `BruntVaisalaFreqSq` is calculated using linear coefficients. When using the `Teos-10` EOS option, the `BruntVaisalaFreqSq` is calculated with non-linear coefficients according to the TOES-10. For both options, `SpecVol` must be calculated prior to calculating `BruntVaisalaFreqSq`, as it is an input to the `BruntVaisalaFreqSq` calculation.
+The `Eos` class also calculates the squared Brunt Vaisala Frequency, which is used by the `VertMix` class to determine water column stability for both the convective adjustment and shear-instability driven mixing. When using the `Linear` EOS option, the `BruntVaisalaFreqSq` is calculated using linear coefficients. When using the `Teos-10` EOS option, the `BruntVaisalaFreqSq` is calculated with non-linear coefficients according to the TOES-10. When using the `constant` EOS option, `BruntVaisalaFreqSq` is zero in this implementation. For all options, `SpecVol` must be calculated prior to calculating `BruntVaisalaFreqSq`, as it is an input to the `BruntVaisalaFreqSq` calculation.

components/omega/src/ocn/Eos.cpp

@@ -21,6 +21,10 @@ Teos10Eos::Teos10Eos(const VertCoord *VCoord)
 LinearEos::LinearEos(const VertCoord *VCoord)
     : MinLayerCell(VCoord->MinLayerCell), MaxLayerCell(VCoord->MaxLayerCell) {}
 
+/// Constructor for ConstantEos
+ConstantEos::ConstantEos(const VertCoord *VCoord)
+    : MinLayerCell(VCoord->MinLayerCell), MaxLayerCell(VCoord->MaxLayerCell) {}
+
 /// Constructor for Teos10 squared Brunt-Vaisala frequency
 Teos10BruntVaisalaFreqSq::Teos10BruntVaisalaFreqSq(const VertCoord *VCoord)
     : MinLayerCell(VCoord->MinLayerCell), MaxLayerCell(VCoord->MaxLayerCell) {}
@@ -36,7 +40,7 @@ Eos::Eos(const std::string &Name, ///< [in] Name for eos object
          const VertCoord *VCoord  ///< [in] Vertical coordinate
          )
     : ComputeSpecVolLinear(VCoord), ComputeSpecVolTeos10(VCoord),
-      ComputeBruntVaisalaFreqSqLinear(VCoord),
+      ComputeSpecVolConstant(VCoord), ComputeBruntVaisalaFreqSqLinear(VCoord),
       ComputeBruntVaisalaFreqSqTeos10(VCoord), Name(Name), Mesh(Mesh),
       VCoord(VCoord) {
    SpecVol = Array2DReal("SpecVol", Mesh->NCellsSize, VCoord->NVertLayers);
@@ -45,6 +49,9 @@ Eos::Eos(const std::string &Name, ///< [in] Name for eos object
    BruntVaisalaFreqSq =
        Array2DReal("BruntVaisalaFreqSq", Mesh->NCellsSize, VCoord->NVertLayers);
 
+   deepCopy(SpecVol, 1.0_Real / RhoSw);
+   deepCopy(SpecVolDisplaced, 1.0_Real / RhoSw);
+
    defineFields();
 }
 
@@ -114,6 +121,8 @@ void Eos::init() {
    } else if ((EosTypeStr == "teos10") or (EosTypeStr == "teos-10") or
               (EosTypeStr == "TEOS-10")) {
       eos->EosChoice = EosType::Teos10Eos;
+   } else if ((EosTypeStr == "constant") or (EosTypeStr == "Constant")) {
+      eos->EosChoice = EosType::ConstantEos;
    } else {
       ABORT_ERROR("Eos::init: Unknown EosType requested");
    }
@@ -128,6 +137,8 @@ void Eos::computeSpecVol(const Array2DReal &ConservTemp,
                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
    OMEGA_SCOPE(MinLayerCell, VCoord->MinLayerCell);
    OMEGA_SCOPE(MaxLayerCell, VCoord->MaxLayerCell);
 
@@ -165,6 +176,20 @@ void Eos::computeSpecVol(const Array2DReal &ConservTemp,
                                             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);
+                 });
+          });
    }
 }
 
@@ -178,6 +203,8 @@ void Eos::computeSpecVolDisp(const Array2DReal &ConservTemp,
                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
    OMEGA_SCOPE(MinLayerCell, VCoord->MinLayerCell);
    OMEGA_SCOPE(MaxLayerCell, VCoord->MaxLayerCell);
 
@@ -214,6 +241,19 @@ void Eos::computeSpecVolDisp(const Array2DReal &ConservTemp,
                                             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);
+                 });
+          });
    }
 }
 

components/omega/src/ocn/Eos.h

@@ -23,8 +23,9 @@
 namespace OMEGA {
 
 enum class EosType {
-   LinearEos, /// Linear equation of state
-   Teos10Eos  /// Roquet et al. 2015 75 term expansion
+   LinearEos,  /// Linear equation of state
+   Teos10Eos,  /// Roquet et al. 2015 75 term expansion
+   ConstantEos /// Constant specific volume equation of state
 };
 
 /// TEOS10 75-term Polynomial Equation of State
@@ -278,6 +279,35 @@ class LinearEos {
    Array1DI4 MaxLayerCell;
 };
 
+/// Constant Equation of State
+class ConstantEos {
+ public:
+   /// 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
+   //   value for all active layers.
+   KOKKOS_FUNCTION void operator()(Array2DReal SpecVol, I4 ICell, I4 KChunk,
+                                   const Array2DReal &ConservTemp,
+                                   const Array2DReal &AbsSalinity) const {
+
+      const I4 KStart = chunkStart(KChunk, MinLayerCell(ICell));
+      const I4 KLen   = chunkLength(KChunk, KStart, 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;
+      }
+   }
+
+ private:
+   Array1DI4 MinLayerCell;
+   Array1DI4 MaxLayerCell;
+};
+
 /// Functor for calculating the squared Brunt-Vaisala frequency using TEOS-10
 class Teos10BruntVaisalaFreqSq {
  public:
@@ -589,8 +619,9 @@ class Eos {
    const HorzMesh *Mesh;    ///< Horizontal mesh
    const VertCoord *VCoord; ///< Vertical coordinate
 
-   Teos10Eos ComputeSpecVolTeos10; ///< TEOS-10 specific volume calculator
-   LinearEos ComputeSpecVolLinear; ///< Linear specific volume calculator
+   Teos10Eos ComputeSpecVolTeos10;     ///< TEOS-10 specific volume calculator
+   LinearEos ComputeSpecVolLinear;     ///< Linear specific volume calculator
+   ConstantEos ComputeSpecVolConstant; ///< Constant specific volume calculator
    Teos10BruntVaisalaFreqSq
        ComputeBruntVaisalaFreqSqTeos10; ///< TEOS-10 squared Brunt-Vaisala
                                         ///< calculator

components/omega/src/ocn/OceanInit.cpp

@@ -136,8 +136,8 @@ int initOmegaModules(MPI_Comm Comm) {
    Tracers::init();
    VertAdv::init();
    AuxiliaryState::init();
-   PressureGrad::init();
    Eos::init();
+   PressureGrad::init();
    Tendencies::init();
    TimeStepper::init2();
 

components/omega/test/ocn/EosTest.cpp

@@ -41,6 +41,8 @@ const Real TeosSVExpValue =
     0.0009732819628; // Expected value for TEOS-10 specific volume
 const Real LinearExpValue =
     0.0009784735812133072; // Expected value for Linear specific volume
+const Real ConstantExpValue =
+    1.0_Real / RhoSw; // Expected value for constant specific volume
 const Real TeosBVFExpValue =
     0.020913834194283325; // Expected value for TEOS-10 squared Brunt-Vaisala
                           // frequency
@@ -242,6 +244,75 @@ void testEosLinearDisplaced() {
    return;
 }
 
+/// Test Constant EOS calculation for all cells/layers
+void testEosConstant() {
+   /// Get mesh and coordinate info
+   const auto Mesh     = HorzMesh::getDefault();
+   const auto VCoord   = VertCoord::getDefault();
+   VCoord->NVertLayers = NVertLayers;
+   I4 NCellsAll        = Mesh->NCellsAll;
+   /// Get Eos instance to test
+   Eos *TestEos       = Eos::getInstance();
+   TestEos->EosChoice = EosType::ConstantEos;
+
+   /// Create and fill ocean state arrays
+   Array2DReal SArray = Array2DReal("SArray", NCellsAll, NVertLayers);
+   Array2DReal TArray = Array2DReal("TArray", NCellsAll, NVertLayers);
+   Array2DReal PArray = Array2DReal("PArray", NCellsAll, NVertLayers);
+   deepCopy(SArray, Sa);
+   deepCopy(TArray, Ct);
+   deepCopy(PArray, P);
+   deepCopy(TestEos->SpecVol, 0.0);
+
+   /// Compute specific volume
+   TestEos->computeSpecVol(TArray, SArray, PArray);
+
+   const auto &MinLayerCell = VCoord->MinLayerCell;
+   const auto &MaxLayerCell = VCoord->MaxLayerCell;
+
+   /// Check all active layers against expected constant value
+   int NumMismatches   = 0;
+   Array2DReal SpecVol = TestEos->SpecVol;
+   parallelReduceOuter(
+       "CheckSpecVolMatrix-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)) {
+                    InnerCount++;
+                 }
+              },
+              NumMismatchesCol);
+
+          Kokkos::single(PerTeam(Team),
+                         [&]() { OuterCount += NumMismatchesCol; });
+       },
+       NumMismatches);
+
+   // If test fails, print bad values and abort
+   if (NumMismatches != 0) {
+      auto SpecVolH = createHostMirrorCopy(SpecVol);
+      for (int I = 0; I < NCellsAll; ++I) {
+         for (int K = 0; K < NVertLayers; ++K) {
+            if (!isApprox(SpecVolH(I, K), ConstantExpValue, RTol))
+               LOG_ERROR("EosTest: SpecVol Constant Bad Value: "
+                         "SpecVol({},{}) = {}; Expected {}",
+                         I, K, SpecVolH(I, K), ConstantExpValue);
+         }
+      }
+      ABORT_ERROR("EosTest: SpecVol Constant FAIL with {} bad values",
+                  NumMismatches);
+   }
+
+   return;
+}
+
 /// Test linear squared Brunt-Vaisala frequency calculation for all cells/layers
 void testBruntVaisalaFreqSqLinear() {
    /// Get mesh and coordinate info
@@ -712,6 +783,7 @@ void eosTest(const std::string &MeshFile = "OmegaMesh.nc") {
 
    testEosLinear();
    testEosLinearDisplaced();
+   testEosConstant();
    testBruntVaisalaFreqSqLinear();
    testEosTeos10();
    testEosTeos10Displaced();