diff --git a/components/omega/src/ocn/HorzMesh.cpp b/components/omega/src/ocn/HorzMesh.cpp index ee47ab55614c..3e731177c58f 100644 --- a/components/omega/src/ocn/HorzMesh.cpp +++ b/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; } }); diff --git a/components/omega/src/ocn/Tendencies.cpp b/components/omega/src/ocn/Tendencies.cpp index bf4949f43949..1e58eee9b1c9 100644 --- a/components/omega/src/ocn/Tendencies.cpp +++ b/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 diff --git a/components/omega/src/ocn/TendencyTerms.cpp b/components/omega/src/ocn/TendencyTerms.cpp index 58a446560fa1..36dff5b5e258 100644 --- a/components/omega/src/ocn/TendencyTerms.cpp +++ b/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 diff --git a/components/omega/src/ocn/TendencyTerms.h b/components/omega/src/ocn/TendencyTerms.h index 3b645f497d02..0c153c8c773a 100644 --- a/components/omega/src/ocn/TendencyTerms.h +++ b/components/omega/src/ocn/TendencyTerms.h @@ -102,7 +102,7 @@ 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]; } } @@ -110,6 +110,7 @@ class PotentialVortHAdvOnEdge { 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,8 +165,9 @@ 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; } } @@ -171,6 +175,7 @@ class SSHGradOnEdge { 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 diff --git a/components/omega/src/ocn/auxiliaryVars/TracerAuxVars.cpp b/components/omega/src/ocn/auxiliaryVars/TracerAuxVars.cpp index d0224f0a54ce..b970c8136ec3 100644 --- a/components/omega/src/ocn/auxiliaryVars/TracerAuxVars.cpp +++ b/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 { diff --git a/components/omega/src/ocn/auxiliaryVars/TracerAuxVars.h b/components/omega/src/ocn/auxiliaryVars/TracerAuxVars.h index 1436dc773af2..ab8c43c4071f 100644 --- a/components/omega/src/ocn/auxiliaryVars/TracerAuxVars.h +++ b/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 diff --git a/components/omega/src/ocn/auxiliaryVars/VelocityDel2AuxVars.cpp b/components/omega/src/ocn/auxiliaryVars/VelocityDel2AuxVars.cpp index fa927d55dc73..060a380ea2a6 100644 --- a/components/omega/src/ocn/auxiliaryVars/VelocityDel2AuxVars.cpp +++ b/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) {} diff --git a/components/omega/src/ocn/auxiliaryVars/VelocityDel2AuxVars.h b/components/omega/src/ocn/auxiliaryVars/VelocityDel2AuxVars.h index e359942aa8f1..93bf91492947 100644 --- a/components/omega/src/ocn/auxiliaryVars/VelocityDel2AuxVars.h +++ b/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; }; diff --git a/components/omega/test/ocn/TendencyTermsTest.cpp b/components/omega/test/ocn/TendencyTermsTest.cpp index 29699e5c73af..f09571745408 100644 --- a/components/omega/test/ocn/TendencyTermsTest.cpp +++ b/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(