Vertical coordinate updates and fixes

components/omega/doc/devGuide/VertCoord.md

@@ -4,21 +4,36 @@
 
 ### Initialization
 
-The default `VertCoord` instance is  created by the `init` method, which assumes that `Decomp` has already been initialized.
+The default `VertCoord` instance is created by the `init` method, which assumes that `Decomp` has already been initialized.
 ```
 Decomp::init();
 VertCoord::init();
 ```
+The `init` method accepts two optional arguments with default values:
+```
+VertCoord::init(const bool ReadStream = true, const int NVertLayers = 0)
+```
+These arguments provide flexibility, particularly for unit testing. When `init(false)` is used, the method skips reading
+the `InitialVertCoord` stream. This is useful for unit tests that rely on meshes lacking the fields required by that stream.
+In this case, the min/max layer arrays are initialized with default values based on the number of vertical layers, and
+`bottomDepth` remains uninitialized. Some tests may utilize a specific number of vertical layers that differs from what is
+defined in the mesh file. To handle this, the `VertCoord` can be initialized with `init(false, LocNVertLayers)` to explicitly
+set the number of vertical layers. An argument for `ReadStream` must be provided in order to explicitly set `NVertLayers`.
+For example, `init(42)` is invalid, `init(false, 42)` must be called instead.
+
 The default instance can be retrieved by:
 ```
 auto *DefVertCoord = VertCoord::getDefault();
 ```
 
 Additional instances can be created by calling the `create` method.
 ```
-VertCoord *VertCoord::create(const std::string &Name,  ///< [in] Name for vertical coordinate
-                             const Decomp *MeshDecomp, ///< [in] Decomp for mesh
-                             Config *Options)          ///< [in] Confiuration options
+VertCoord *VertCoord::create(const std::string &Name,      ///< [in] Name for vertical coordinate
+                             const Decomp *MeshDecomp,     ///< [in] Decomp for mesh
+                             Config *Options,              ///< [in] Confiuration options
+                             const bool ReadStream = true, ///< [in] optional logical to read stream
+                             const int NVertLayers = 0     ///< [in] optional int to set vertical dim
+)
 ```
 This calls the constructor and places the instance in a map that can be used to retrieve the instance by name:
 ```

components/omega/src/infra/IOStream.cpp

@@ -1578,11 +1578,19 @@ Error IOStream::readFieldData(
    // lower case
    std::string OldFieldName = FieldName;
    OldFieldName[0]          = std::tolower(OldFieldName[0]);
-   bool OnHost              = FieldPtr->isOnHost();
-   bool IsDistributed       = FieldPtr->isDistributed();
-   bool IsTimeDependent     = FieldPtr->isTimeDependent();
-   ArrayDataType MyType     = FieldPtr->getType();
-   int NDims                = FieldPtr->getNumDims();
+   // Check for "Layer" in field name, backwards compatibility requires
+   // replacing "Layer" with "Level"
+   std::string OmegaSubStr = "Layer";
+   std::string MPASSubStr  = "Level";
+   size_t pos              = OldFieldName.find(OmegaSubStr);
+   if (pos != std::string::npos) {
+      OldFieldName.replace(pos, OmegaSubStr.length(), MPASSubStr);
+   }
+   bool OnHost          = FieldPtr->isOnHost();
+   bool IsDistributed   = FieldPtr->isDistributed();
+   bool IsTimeDependent = FieldPtr->isTimeDependent();
+   ArrayDataType MyType = FieldPtr->getType();
+   int NDims            = FieldPtr->getNumDims();
    if (NDims < 0)
       ABORT_ERROR("IOStream readFieldData: "
                   "Invalid number of dimensions for Field {}",

components/omega/src/ocn/AuxiliaryState.cpp

@@ -18,14 +18,15 @@ static std::string stripDefault(const std::string &Name) {
 // Constructor. Constructs the member auxiliary variables and registers their
 // fields with IOStreams
 AuxiliaryState::AuxiliaryState(const std::string &Name, const HorzMesh *Mesh,
-                               Halo *MeshHalo, int NVertLayers, int NTracers)
-    : Mesh(Mesh), MeshHalo(MeshHalo), Name(stripDefault(Name)),
+                               const VertCoord *VCoord, Halo *MeshHalo,
+                               int NVertLayers, int NTracers)
+    : Mesh(Mesh), VCoord(VCoord), MeshHalo(MeshHalo), Name(stripDefault(Name)),
       KineticAux(stripDefault(Name), Mesh, NVertLayers),
       LayerThicknessAux(stripDefault(Name), Mesh, NVertLayers),
       VorticityAux(stripDefault(Name), Mesh, NVertLayers),
-      VelocityDel2Aux(stripDefault(Name), Mesh, NVertLayers),
+      VelocityDel2Aux(stripDefault(Name), Mesh, VCoord, NVertLayers),
       WindForcingAux(stripDefault(Name), Mesh),
-      TracerAux(stripDefault(Name), Mesh, NVertLayers, NTracers) {
+      TracerAux(stripDefault(Name), Mesh, VCoord, NVertLayers, NTracers) {
 
    GroupName = "AuxiliaryState";
    if (Name != "Default") {
@@ -192,7 +193,8 @@ void AuxiliaryState::computeAll(const OceanState *State,
 
 // Create a non-default auxiliary state
 AuxiliaryState *AuxiliaryState::create(const std::string &Name,
-                                       const HorzMesh *Mesh, Halo *MeshHalo,
+                                       const HorzMesh *Mesh,
+                                       const VertCoord *VCoord, Halo *MeshHalo,
                                        int NVertLayers, const int NTracers) {
    if (AllAuxStates.find(Name) != AllAuxStates.end()) {
       LOG_ERROR("Attempted to create a new AuxiliaryState with name {} but it "
@@ -202,7 +204,7 @@ AuxiliaryState *AuxiliaryState::create(const std::string &Name,
    }
 
    auto *NewAuxState =
-       new AuxiliaryState(Name, Mesh, MeshHalo, NVertLayers, NTracers);
+       new AuxiliaryState(Name, Mesh, VCoord, MeshHalo, NVertLayers, NTracers);
    AllAuxStates.emplace(Name, NewAuxState);
 
    return NewAuxState;
@@ -211,14 +213,15 @@ AuxiliaryState *AuxiliaryState::create(const std::string &Name,
 // Create the default auxiliary state. Assumes that HorzMesh, VertCoord and
 // Halo have been initialized.
 void AuxiliaryState::init() {
-   const HorzMesh *DefMesh = HorzMesh::getDefault();
-   Halo *DefHalo           = Halo::getDefault();
+   const HorzMesh *DefMesh    = HorzMesh::getDefault();
+   const VertCoord *DefVCoord = VertCoord::getDefault();
+   Halo *DefHalo              = Halo::getDefault();
 
    int NVertLayers = VertCoord::getDefault()->NVertLayers;
    int NTracers    = Tracers::getNumTracers();
 
    AuxiliaryState::DefaultAuxState = AuxiliaryState::create(
-       "Default", DefMesh, DefHalo, NVertLayers, NTracers);
+       "Default", DefMesh, DefVCoord, DefHalo, NVertLayers, NTracers);
 
    Config *OmegaConfig = Config::getOmegaConfig();
    DefaultAuxState->readConfigOptions(OmegaConfig);

components/omega/src/ocn/AuxiliaryState.h

@@ -50,7 +50,8 @@ class AuxiliaryState {
 
    // Create a non-default auxiliary state
    static AuxiliaryState *create(const std::string &Name, const HorzMesh *Mesh,
-                                 Halo *MeshHalo, int NVertLayers, int NTracers);
+                                 const VertCoord *VCoord, Halo *MeshHalo,
+                                 int NVertLayers, int NTracers);
 
    /// Get the default auxiliary state
    static AuxiliaryState *getDefault();
@@ -82,13 +83,15 @@ class AuxiliaryState {
                    int TimeLevel) const;
 
  private:
-   AuxiliaryState(const std::string &Name, const HorzMesh *Mesh, Halo *MeshHalo,
-                  int NVertLayers, int NTracers);
+   AuxiliaryState(const std::string &Name, const HorzMesh *Mesh,
+                  const VertCoord *VCoord, Halo *MeshHalo, int NVertLayers,
+                  int NTracers);
 
    AuxiliaryState(const AuxiliaryState &) = delete;
    AuxiliaryState(AuxiliaryState &&)      = delete;
 
    const HorzMesh *Mesh;
+   const VertCoord *VCoord;
    Halo *MeshHalo;
    static AuxiliaryState *DefaultAuxState;
    static std::map<std::string, std::unique_ptr<AuxiliaryState>> AllAuxStates;

components/omega/src/ocn/HorzMesh.cpp

@@ -37,28 +37,22 @@ void HorzMesh::init() {
    // Retrieve the default decomposition
    Decomp *DefDecomp = Decomp::getDefault();
 
-   I4 NVertLayers = VertCoord::getDefault()->NVertLayers;
-
    // Create the default mesh and set pointer to it
-   HorzMesh::DefaultHorzMesh = create("Default", DefDecomp, NVertLayers);
+   HorzMesh::DefaultHorzMesh = create("Default", DefDecomp);
 }
 
 //------------------------------------------------------------------------------
 // Construct a new local mesh given a decomposition
 
 HorzMesh::HorzMesh(const std::string &Name, //< [in] Name for new mesh
-                   Decomp *MeshDecomp,      //< [in] Decomp for the new mesh
-                   I4 InNVertLayers         //< [in} num vertical layers
+                   Decomp *MeshDecomp       //< [in] Decomp for the new mesh
 ) {
 
    MeshName = Name;
 
    // Retrieve mesh files name from Decomp
    MeshFileName = MeshDecomp->MeshFileName;
 
-   // Set NVertLayers
-   NVertLayers = InNVertLayers;
-
    // Retrieve mesh cell/edge/vertex totals from Decomp
    NCellsHalo  = MeshDecomp->NCellsHalo;
    NCellsHaloH = MeshDecomp->NCellsHaloH;
@@ -141,9 +135,6 @@ HorzMesh::HorzMesh(const std::string &Name, //< [in] Name for new mesh
    // Compute EdgeSignOnCells and EdgeSignOnVertex
    computeEdgeSign();
 
-   // TODO: implement setMasks during Mesh constructor
-   setMasks(NVertLayers);
-
    // set mesh scaling coefficients
    setMeshScaling();
 
@@ -152,8 +143,7 @@ HorzMesh::HorzMesh(const std::string &Name, //< [in] Name for new mesh
 /// Creates a new mesh by calling the constructor and puts it in the
 /// AllHorzMeshes map
 HorzMesh *HorzMesh::create(const std::string &Name, //< [in] Name for new mesh
-                           Decomp *MeshDecomp, //< [in] Decomp for the new mesh
-                           I4 InNVertLayers    //< [in] num vertical layers
+                           Decomp *MeshDecomp //< [in] Decomp for the new mesh
 ) {
    // Check to see if a mesh of the same name already exists and
    // if so, exit with an error
@@ -166,7 +156,7 @@ HorzMesh *HorzMesh::create(const std::string &Name, //< [in] Name for new mesh
 
    // create a new mesh on the heap and put it in a map of
    // unique_ptrs, which will manage its lifetime
-   auto *NewHorzMesh = new HorzMesh(Name, MeshDecomp, InNVertLayers);
+   auto *NewHorzMesh = new HorzMesh(Name, MeshDecomp);
    AllHorzMeshes.emplace(Name, NewHorzMesh);
 
    return NewHorzMesh;
@@ -574,33 +564,6 @@ void HorzMesh::computeEdgeSign() {
    EdgeSignOnVertexH = createHostMirrorCopy(EdgeSignOnVertex);
 } // end computeEdgeSign
 
-//------------------------------------------------------------------------------
-// set computational masks for mesh elements
-// TODO: this is just a placeholder, implement actual masks for edges, cells,
-// and vertices
-void HorzMesh::setMasks(int NVertLayers) {
-
-   EdgeMask = Array2DReal("EdgeMask", NEdgesSize, NVertLayers);
-
-   OMEGA_SCOPE(LocEdgeMask, EdgeMask);
-   OMEGA_SCOPE(LocCellsOnEdge, CellsOnEdge);
-   OMEGA_SCOPE(LocNCellsAll, NCellsAll);
-
-   deepCopy(EdgeMask, 1.0);
-   parallelFor(
-       {NEdgesAll, NVertLayers}, KOKKOS_LAMBDA(int Edge, int K) {
-          const I4 Cell1 = LocCellsOnEdge(Edge, 0);
-          const I4 Cell2 = LocCellsOnEdge(Edge, 1);
-          if (!(Cell1 >= 0 and Cell1 < LocNCellsAll) or
-              !(Cell2 >= 0 and Cell2 < LocNCellsAll)) {
-             LocEdgeMask(Edge, K) = 0.0;
-          }
-       });
-
-   EdgeMaskH = createHostMirrorCopy(EdgeMask);
-
-} // end setMasks
-
 //------------------------------------------------------------------------------
 // Set mesh scaling coefficients for mixing terms in momentum and tracer
 // equations so viscosity and diffusion scale with mesh.

components/omega/src/ocn/HorzMesh.h

@@ -14,7 +14,6 @@
 #include "Decomp.h"
 #include "MachEnv.h"
 #include "OmegaKokkos.h"
-#include "VertCoord.h"
 
 #include <memory>
 #include <string>
@@ -68,8 +67,7 @@ class HorzMesh {
 
    /// Construct a new local mesh for a given decomposition
    HorzMesh(const std::string &Name, ///< [in] Name for mesh
-            Decomp *Decomp,          ///< [in] Decomposition for mesh
-            I4 InNVertLayers         ///< [in] num vertical layers
+            Decomp *Decomp           ///< [in] Decomposition for mesh
    );
 
    // Forbid copy and move construction
@@ -81,8 +79,6 @@ class HorzMesh {
    // private function.
    void computeEdgeSign();
 
-   void setMasks(int NVertLayers);
-
    void setMeshScaling();
 
    // Variables
@@ -97,8 +93,6 @@ class HorzMesh {
    // Note that all sizes are actual counts (1-based) so that loop extents
    // should always use the 0:NCellsXX-1 form.
 
-   I4 NVertLayers; ///< number of vertical layers
-
    Array1DI4 NCellsHalo;      ///< num cells owned+halo for halo layer
    HostArray1DI4 NCellsHaloH; ///< num cells owned+halo for halo layer
    I4 NCellsOwned;            ///< Number of cells owned by this task
@@ -272,8 +266,7 @@ class HorzMesh {
    /// Creates a new mesh by calling the constructor and puts it in the
    /// AllHorzMeshes map
    static HorzMesh *create(const std::string &Name, ///< [in] Name for mesh
-                           Decomp *Decomp,  ///< [in] Decomposition for mesh
-                           I4 InNVertLayers ///< [in] num vertival layers
+                           Decomp *Decomp ///< [in] Decomposition for mesh
    );
 
    /// Destructor - deallocates all memory and deletes a HorzMesh

components/omega/src/ocn/OceanInit.cpp

@@ -128,9 +128,8 @@ int initOmegaModules(MPI_Comm Comm) {
       ABORT_ERROR("ocnInit: Error initializing default halo");
    }
 
-   VertCoord::init1();
    HorzMesh::init();
-   VertCoord::init2();
+   VertCoord::init();
    Tracers::init();
    AuxiliaryState::init();
    Tendencies::init();

components/omega/src/ocn/Tendencies.cpp

@@ -221,11 +221,12 @@ Tendencies::Tendencies(const std::string &Name, ///< [in] Name for tendencies
                        Config *Options,         ///< [in] Configuration options
                        CustomTendencyType InCustomThicknessTend,
                        CustomTendencyType InCustomVelocityTend)
-    : ThicknessFluxDiv(Mesh), PotientialVortHAdv(Mesh), KEGrad(Mesh),
-      SSHGrad(Mesh), VelocityDiffusion(Mesh), VelocityHyperDiff(Mesh),
-      WindForcing(Mesh), BottomDrag(Mesh, VCoord), TracerHorzAdv(Mesh),
-      TracerDiffusion(Mesh), TracerHyperDiff(Mesh),
-      CustomThicknessTend(InCustomThicknessTend),
+    : ThicknessFluxDiv(Mesh), PotientialVortHAdv(Mesh, VCoord),
+      KEGrad(Mesh, VCoord), SSHGrad(Mesh, VCoord),
+      VelocityDiffusion(Mesh, VCoord), VelocityHyperDiff(Mesh, VCoord),
+      WindForcing(Mesh, VCoord), BottomDrag(Mesh, VCoord),
+      TracerHorzAdv(Mesh, VCoord), TracerDiffusion(Mesh, VCoord),
+      TracerHyperDiff(Mesh, VCoord), CustomThicknessTend(InCustomThicknessTend),
       CustomVelocityTend(InCustomVelocityTend) {
 
    // Tendency arrays