Use EdgeMask to handle boundary edges

components/omega/src/ocn/HorzMesh.cpp

@@ -635,12 +635,18 @@ void HorzMesh::setMasks(int NVertLevels) {
 
    EdgeMask = Array2DReal("EdgeMask", NEdgesSize, NVertLevels);
 
-   OMEGA_SCOPE(O_EdgeMask, EdgeMask);
+   OMEGA_SCOPE(LocEdgeMask, EdgeMask);
+   OMEGA_SCOPE(LocCellsOnEdge, CellsOnEdge);
+   OMEGA_SCOPE(LocNCellsAll, NCellsAll);
 
+   deepCopy(EdgeMask, 1.0);
    parallelFor(
-       {NEdgesAll}, KOKKOS_LAMBDA(int Edge) {
-          for (int K = 0; K < NVertLevels; ++K) {
-             O_EdgeMask(Edge, K) = 1.0;
+       {NEdgesAll, NVertLevels}, 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;
           }
        });
 

components/omega/src/ocn/Tendencies.cpp

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

components/omega/src/ocn/TendencyTerms.cpp

@@ -24,13 +24,15 @@ ThicknessFluxDivOnCell::ThicknessFluxDivOnCell(const HorzMesh *Mesh)
 
 PotentialVortHAdvOnEdge::PotentialVortHAdvOnEdge(const HorzMesh *Mesh)
     : NEdgesOnEdge(Mesh->NEdgesOnEdge), EdgesOnEdge(Mesh->EdgesOnEdge),
-      WeightsOnEdge(Mesh->WeightsOnEdge) {}
+      WeightsOnEdge(Mesh->WeightsOnEdge), EdgeMask(Mesh->EdgeMask) {}
 
 KEGradOnEdge::KEGradOnEdge(const HorzMesh *Mesh)
-    : CellsOnEdge(Mesh->CellsOnEdge), DcEdge(Mesh->DcEdge) {}
+    : CellsOnEdge(Mesh->CellsOnEdge), DcEdge(Mesh->DcEdge),
+      EdgeMask(Mesh->EdgeMask) {}
 
 SSHGradOnEdge::SSHGradOnEdge(const HorzMesh *Mesh)
-    : CellsOnEdge(Mesh->CellsOnEdge), DcEdge(Mesh->DcEdge) {}
+    : CellsOnEdge(Mesh->CellsOnEdge), DcEdge(Mesh->DcEdge),
+      EdgeMask(Mesh->EdgeMask) {}
 
 VelocityDiffusionOnEdge::VelocityDiffusionOnEdge(const HorzMesh *Mesh)
     : CellsOnEdge(Mesh->CellsOnEdge), VerticesOnEdge(Mesh->VerticesOnEdge),
@@ -42,26 +44,30 @@ VelocityHyperDiffOnEdge::VelocityHyperDiffOnEdge(const HorzMesh *Mesh)
       DcEdge(Mesh->DcEdge), DvEdge(Mesh->DvEdge),
       MeshScalingDel4(Mesh->MeshScalingDel4), EdgeMask(Mesh->EdgeMask) {}
 
+WindForcingOnEdge::WindForcingOnEdge(const HorzMesh *Mesh)
+    : Enabled(false), EdgeMask(Mesh->EdgeMask) {}
+
 BottomDragOnEdge::BottomDragOnEdge(const HorzMesh *Mesh)
     : Enabled(false), Coeff(0), CellsOnEdge(Mesh->CellsOnEdge),
-      NVertLevels(Mesh->NVertLevels) {}
+      NVertLevels(Mesh->NVertLevels), EdgeMask(Mesh->EdgeMask) {}
 
 TracerHorzAdvOnCell::TracerHorzAdvOnCell(const HorzMesh *Mesh)
     : NEdgesOnCell(Mesh->NEdgesOnCell), EdgesOnCell(Mesh->EdgesOnCell),
       CellsOnEdge(Mesh->CellsOnEdge), EdgeSignOnCell(Mesh->EdgeSignOnCell),
-      DvEdge(Mesh->DvEdge), AreaCell(Mesh->AreaCell) {}
+      DvEdge(Mesh->DvEdge), AreaCell(Mesh->AreaCell), EdgeMask(Mesh->EdgeMask) {
+}
 
 TracerDiffOnCell::TracerDiffOnCell(const HorzMesh *Mesh)
     : NEdgesOnCell(Mesh->NEdgesOnCell), EdgesOnCell(Mesh->EdgesOnCell),
       CellsOnEdge(Mesh->CellsOnEdge), EdgeSignOnCell(Mesh->EdgeSignOnCell),
       DvEdge(Mesh->DvEdge), DcEdge(Mesh->DcEdge), AreaCell(Mesh->AreaCell),
-      MeshScalingDel2(Mesh->MeshScalingDel2) {}
+      MeshScalingDel2(Mesh->MeshScalingDel2), EdgeMask(Mesh->EdgeMask) {}
 
 TracerHyperDiffOnCell::TracerHyperDiffOnCell(const HorzMesh *Mesh)
     : NEdgesOnCell(Mesh->NEdgesOnCell), EdgesOnCell(Mesh->EdgesOnCell),
       CellsOnEdge(Mesh->CellsOnEdge), EdgeSignOnCell(Mesh->EdgeSignOnCell),
       DvEdge(Mesh->DvEdge), DcEdge(Mesh->DcEdge), AreaCell(Mesh->AreaCell),
-      MeshScalingDel4(Mesh->MeshScalingDel4) {}
+      MeshScalingDel4(Mesh->MeshScalingDel4), EdgeMask(Mesh->EdgeMask) {}
 
 } // end namespace OMEGA
 

components/omega/src/ocn/TendencyTerms.h

@@ -102,14 +102,15 @@ class PotentialVortHAdvOnEdge {
 
       for (int KVec = 0; KVec < VecLength; ++KVec) {
          const I4 K = KStart + KVec;
-         Tend(IEdge, K) += VortTmp[KVec];
+         Tend(IEdge, K) += EdgeMask(IEdge, K) * VortTmp[KVec];
       }
    }
 
  private:
    Array1DI4 NEdgesOnEdge;
    Array2DI4 EdgesOnEdge;
    Array2DReal WeightsOnEdge;
+   Array2DReal EdgeMask;
 };
 
 /// Gradient of kinetic energy defined on edges, for momentum equation
@@ -132,13 +133,15 @@ class KEGradOnEdge {
 
       for (int KVec = 0; KVec < VecLength; ++KVec) {
          const I4 K = KStart + KVec;
-         Tend(IEdge, K) -= (KECell(JCell1, K) - KECell(JCell0, K)) * InvDcEdge;
+         Tend(IEdge, K) -= EdgeMask(IEdge, K) *
+                           (KECell(JCell1, K) - KECell(JCell0, K)) * InvDcEdge;
       }
    }
 
  private:
    Array2DI4 CellsOnEdge;
    Array1DReal DcEdge;
+   Array2DReal EdgeMask;
 };
 
 /// Gradient of sea surface height defined on edges multipled by gravitational
@@ -162,15 +165,17 @@ class SSHGradOnEdge {
 
       for (int KVec = 0; KVec < VecLength; ++KVec) {
          const I4 K = KStart + KVec;
-         Tend(IEdge, K) -=
-             Grav * (SshCell(ICell1, K) - SshCell(ICell0, K)) * InvDcEdge;
+         Tend(IEdge, K) -= EdgeMask(IEdge, K) * Grav *
+                           (SshCell(ICell1, K) - SshCell(ICell0, K)) *
+                           InvDcEdge;
       }
    }
 
  private:
    Real Grav = 9.80665_Real;
    Array2DI4 CellsOnEdge;
    Array1DReal DcEdge;
+   Array2DReal EdgeMask;
 };
 
 /// Laplacian horizontal mixing, for momentum equation
@@ -277,6 +282,9 @@ class WindForcingOnEdge {
    bool Enabled;
    Real SaltWaterDensity;
 
+   /// constructor declaration
+   WindForcingOnEdge(const HorzMesh *Mesh);
+
    /// The functor takes the edge index, vertical chunk index, and arrays for
    /// normal wind stress and edge layer thickness, outputs tendency array
    KOKKOS_FUNCTION void operator()(const Array2DReal &Tend, I4 IEdge, I4 KChunk,
@@ -286,10 +294,13 @@ class WindForcingOnEdge {
          const I4 K = 0;
 
          const Real InvThickEdge = 1._Real / LayerThickEdge(IEdge, K);
-         Tend(IEdge, K) +=
-             InvThickEdge * NormalStressEdge(IEdge) / SaltWaterDensity;
+         Tend(IEdge, K) += EdgeMask(IEdge, K) * InvThickEdge *
+                           NormalStressEdge(IEdge) / SaltWaterDensity;
       }
    }
+
+ private:
+   Array2DReal EdgeMask;
 };
 
 /// Bottom drag
@@ -317,13 +328,14 @@ class BottomDragOnEdge {
           Kokkos::sqrt(KECell(JCell0, KBot) + KECell(JCell1, KBot));
 
       const Real InvThickEdge = 1._Real / LayerThickEdge(IEdge, KBot);
-      Tend(IEdge, KBot) -=
-          Coeff * VelNormEdge * InvThickEdge * NormalVelEdge(IEdge, KBot);
+      Tend(IEdge, KBot) -= EdgeMask(IEdge, KBot) * Coeff * VelNormEdge *
+                           InvThickEdge * NormalVelEdge(IEdge, KBot);
    }
 
  private:
    I4 NVertLevels;
    Array2DI4 CellsOnEdge;
+   Array2DReal EdgeMask;
 };
 
 // Tracer horizontal advection term
@@ -347,7 +359,8 @@ class TracerHorzAdvOnCell {
 
          for (int KVec = 0; KVec < VecLength; ++KVec) {
             const I4 K = KStart + KVec;
-            HAdvTmp[KVec] -= DvEdge(JEdge) * EdgeSignOnCell(ICell, J) *
+            HAdvTmp[KVec] -= EdgeMask(JEdge, K) * DvEdge(JEdge) *
+                             EdgeSignOnCell(ICell, J) *
                              HTracersOnEdge(L, JEdge, K) *
                              NormVelEdge(JEdge, K) * InvAreaCell;
          }
@@ -365,6 +378,7 @@ class TracerHorzAdvOnCell {
    Array2DReal EdgeSignOnCell;
    Array1DReal DvEdge;
    Array1DReal AreaCell;
+   Array2DReal EdgeMask;
 };
 
 // Tracer horizontal diffusion term
@@ -400,8 +414,8 @@ class TracerDiffOnCell {
             const Real TracerGrad =
                 (TracerCell(L, JCell1, K) - TracerCell(L, JCell0, K));
 
-            DiffTmp[KVec] -= EdgeSignOnCell(ICell, J) * RTemp *
-                             MeanLayerThickEdge(JEdge, K) * TracerGrad;
+            DiffTmp[KVec] -= EdgeMask(JEdge, K) * EdgeSignOnCell(ICell, J) *
+                             RTemp * MeanLayerThickEdge(JEdge, K) * TracerGrad;
          }
       }
       for (int KVec = 0; KVec < VecLength; ++KVec) {
@@ -419,6 +433,7 @@ class TracerDiffOnCell {
    Array1DReal DcEdge;
    Array1DReal AreaCell;
    Array1DReal MeshScalingDel2;
+   Array2DReal EdgeMask;
 };
 
 // Tracer biharmonic horizontal mixing term
@@ -453,7 +468,8 @@ class TracerHyperDiffOnCell {
             const Real Del2TrGrad =
                 (TrDel2Cell(L, JCell1, K) - TrDel2Cell(L, JCell0, K));
 
-            HypTmp[KVec] -= EdgeSignOnCell(ICell, J) * RTemp * Del2TrGrad;
+            HypTmp[KVec] -= EdgeMask(JEdge, K) * EdgeSignOnCell(ICell, J) *
+                            RTemp * Del2TrGrad;
          }
       }
       for (int KVec = 0; KVec < VecLength; ++KVec) {
@@ -471,6 +487,7 @@ class TracerHyperDiffOnCell {
    Array1DReal DcEdge;
    Array1DReal AreaCell;
    Array1DReal MeshScalingDel4;
+   Array2DReal EdgeMask;
 };
 
 } // namespace OMEGA

components/omega/src/ocn/auxiliaryVars/TracerAuxVars.cpp

@@ -14,7 +14,8 @@ TracerAuxVars::TracerAuxVars(const std::string &AuxStateSuffix,
                       Mesh->NCellsSize, NVertLevels),
       NEdgesOnCell(Mesh->NEdgesOnCell), EdgesOnCell(Mesh->EdgesOnCell),
       CellsOnEdge(Mesh->CellsOnEdge), EdgeSignOnCell(Mesh->EdgeSignOnCell),
-      DcEdge(Mesh->DcEdge), DvEdge(Mesh->DvEdge), AreaCell(Mesh->AreaCell) {}
+      DcEdge(Mesh->DcEdge), DvEdge(Mesh->DvEdge), AreaCell(Mesh->AreaCell),
+      EdgeMask(Mesh->EdgeMask) {}
 
 void TracerAuxVars::registerFields(const std::string &AuxGroupName,
                                    const std::string &MeshName) const {

components/omega/src/ocn/auxiliaryVars/TracerAuxVars.h

@@ -79,7 +79,8 @@ class TracerAuxVars {
          for (int KVec = 0; KVec < VecLength; ++KVec) {
             const int K           = KStart + KVec;
             const Real TracerGrad = TrCell(L, JCell1, K) - TrCell(L, JCell0, K);
-            Del2TrCellTmp[KVec] -= EdgeSignOnCell(ICell, J) * DvDcEdge *
+            Del2TrCellTmp[KVec] -= EdgeMask(JEdge, K) *
+                                   EdgeSignOnCell(ICell, J) * DvDcEdge *
                                    LayerThickEdgeMean(JEdge, K) * TracerGrad;
          }
       }
@@ -101,6 +102,7 @@ class TracerAuxVars {
    Array1DReal DcEdge;
    Array1DReal DvEdge;
    Array1DReal AreaCell;
+   Array2DReal EdgeMask;
 };
 
 } // namespace OMEGA

components/omega/src/ocn/auxiliaryVars/VelocityDel2AuxVars.cpp

@@ -17,7 +17,7 @@ VelocityDel2AuxVars::VelocityDel2AuxVars(const std::string &AuxStateSuffix,
       EdgeSignOnCell(Mesh->EdgeSignOnCell), DcEdge(Mesh->DcEdge),
       DvEdge(Mesh->DvEdge), AreaCell(Mesh->AreaCell),
       EdgesOnVertex(Mesh->EdgesOnVertex), CellsOnEdge(Mesh->CellsOnEdge),
-      VerticesOnEdge(Mesh->VerticesOnEdge),
+      VerticesOnEdge(Mesh->VerticesOnEdge), EdgeMask(Mesh->EdgeMask),
       EdgeSignOnVertex(Mesh->EdgeSignOnVertex),
       AreaTriangle(Mesh->AreaTriangle), VertexDegree(Mesh->VertexDegree) {}
 

components/omega/src/ocn/auxiliaryVars/VelocityDel2AuxVars.h

@@ -40,7 +40,7 @@ class VelocityDel2AuxVars {
          const Real CurlVort =
              -(RelVortVertex(JVertex1, K) - RelVortVertex(JVertex0, K)) *
              InvDvEdge;
-         Del2Edge(IEdge, K) = GradDiv + CurlVort;
+         Del2Edge(IEdge, K) = EdgeMask(IEdge, K) * GradDiv + CurlVort;
       }
    }
 
@@ -104,6 +104,7 @@ class VelocityDel2AuxVars {
    Array2DI4 VerticesOnEdge;
    Array2DReal EdgeSignOnVertex;
    Array1DReal AreaTriangle;
+   Array2DReal EdgeMask;
    I4 VertexDegree;
 };
 

components/omega/test/ocn/TendencyTermsTest.cpp

@@ -736,7 +736,7 @@ int testWindForcing(int NVertLevels) {
    // Compute numerical result
    Array2DReal NumWindForcing("NumWindForcing", Mesh->NEdgesOwned, NVertLevels);
 
-   WindForcingOnEdge WindForcingOnE;
+   WindForcingOnEdge WindForcingOnE(Mesh);
    WindForcingOnE.SaltWaterDensity = SaltWaterDensity;
 
    parallelFor(