Merge pull request #456 from xylar/omega/fix-del4-curl-operator

Compute the Laplacian of relative vorticity (Del2RelVortVertex) over the full
vertex valid range [MinLayerVertexTop, MaxLayerVertexBot] and clamp each edge
contribution to that edge's valid range [MinLayerEdgeTop, MaxLayerEdgeBot],
matching VorticityAuxVars::computeVarsOnVertex.

Previously Del2RelVortVertex was computed only over the narrower
[MinLayerVertexBot, MaxLayerVertexTop] range (layers where every surrounding
cell is active) and summed Del2Edge with no per-edge clamp. The biharmonic
velocity tendency (VelocityHyperDiffOnEdge) reads Del2RelVortVertex over the 
edge range up to MaxLayerEdgeTop, which for a deep edge sharing a vertex with a
shallower cell exceeds MaxLayerVertexTop. Those boundary-vertex layers were
never computed, so on partial-bottom meshes the biharmonic term was 
under-computed there. This is a latent correctness bug on its own; on the 
fill-values branch (#428) the uncomputed layers hold FillValueReal, so the term
instead blew up to ~1e45, corrupting NormalVelocity and, through the flux
divergence and vertical advection, PseudoThickness and the tracers (surfaced by
the new state validation in DRIVER_TEST).

The wider range gives boundary vertices a valid, generally non-zero value from
their active edges; inactive edges contribute zero via the per-edge clamp and 
Del2Edge's zeroed boundary band, so no fill value is read.

Author: sbrus89

components/omega/src/ocn/AuxiliaryState.cpp

@@ -117,10 +117,8 @@ void AuxiliaryState::computeMomAux(const OceanState *State,
 
    OMEGA_SCOPE(MinLayerCell, VCoord->MinLayerCell);
    OMEGA_SCOPE(MaxLayerCell, VCoord->MaxLayerCell);
-   OMEGA_SCOPE(MinLayerVertexBot, VCoord->MinLayerVertexBot);
    OMEGA_SCOPE(MinLayerVertexTop, VCoord->MinLayerVertexTop);
    OMEGA_SCOPE(MaxLayerVertexBot, VCoord->MaxLayerVertexBot);
-   OMEGA_SCOPE(MaxLayerVertexTop, VCoord->MaxLayerVertexTop);
    OMEGA_SCOPE(MinLayerEdgeTop, VCoord->MinLayerEdgeTop);
    OMEGA_SCOPE(MinLayerEdgeBot, VCoord->MinLayerEdgeBot);
    OMEGA_SCOPE(MaxLayerEdgeBot, VCoord->MaxLayerEdgeBot);
@@ -204,8 +202,12 @@ void AuxiliaryState::computeMomAux(const OceanState *State,
    parallelForOuter(
        "vertexAuxState2", {Mesh->NVerticesAll},
        KOKKOS_LAMBDA(int IVertex, const TeamMember &Team) {
-          const int KMin   = MinLayerVertexBot(IVertex);
-          const int KMax   = MaxLayerVertexTop(IVertex);
+          // Del2RelVortVertex is computed over the full vertex valid range
+          // [MinLayerVertexTop, MaxLayerVertexBot] so that boundary-vertex
+          // layers read by the biharmonic velocity tendency are valid rather
+          // than fill values (see VelocityDel2AuxVars::computeVarsOnVertex).
+          const int KMin   = MinLayerVertexTop(IVertex);
+          const int KMax   = MaxLayerVertexBot(IVertex);
           const int KRange = vertRangeChunked(KMin, KMax);
 
           parallelForInner(

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

@@ -24,8 +24,10 @@ VelocityDel2AuxVars::VelocityDel2AuxVars(const std::string &AuxStateSuffix,
       AreaTriangle(Mesh->AreaTriangle), VertexDegree(Mesh->VertexDegree),
       MinLayerEdgeBot(VCoord->MinLayerEdgeBot),
       MaxLayerEdgeTop(VCoord->MaxLayerEdgeTop),
-      MinLayerVertexBot(VCoord->MinLayerVertexBot),
-      MaxLayerVertexTop(VCoord->MaxLayerVertexTop),
+      MinLayerEdgeTop(VCoord->MinLayerEdgeTop),
+      MaxLayerEdgeBot(VCoord->MaxLayerEdgeBot),
+      MinLayerVertexTop(VCoord->MinLayerVertexTop),
+      MaxLayerVertexBot(VCoord->MaxLayerVertexBot),
       MinLayerCell(VCoord->MinLayerCell), MaxLayerCell(VCoord->MaxLayerCell) {}
 
 void VelocityDel2AuxVars::registerFields(const std::string &AuxGroupName,

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

@@ -75,25 +75,37 @@ class VelocityDel2AuxVars {
    }
 
    KOKKOS_FUNCTION void computeVarsOnVertex(int IVertex, int KChunk) const {
-      const int KStart = chunkStart(KChunk, MinLayerVertexBot(IVertex));
-      const int KLen = chunkLength(KChunk, KStart, MaxLayerVertexTop(IVertex));
+      // Compute over the full vertex valid range [MinLayerVertexTop,
+      // MaxLayerVertexBot] so that boundary-vertex layers (where only some
+      // surrounding cells are active) receive a valid value. Each edge's
+      // contribution is clamped to that edge's valid range [MinLayerEdgeTop,
+      // MaxLayerEdgeBot], where Del2Edge has been computed or zeroed; this
+      // matches VorticityAuxVars::computeVarsOnVertex and avoids reading
+      // uninitialized (fill-value) layers of Del2Edge for deeper edges.
+      const int KStartVertex = chunkStart(KChunk, MinLayerVertexTop(IVertex));
+      const int KLenVertex =
+          chunkLength(KChunk, KStartVertex, MaxLayerVertexBot(IVertex));
+      const int KEndVertex = KStartVertex + KLenVertex - 1;
 
       const Real InvAreaTriangle = 1._Real / AreaTriangle(IVertex);
 
       Real Del2RelVortVertexTmp[VecLength] = {0};
 
       for (int J = 0; J < VertexDegree; ++J) {
          const int JEdge = EdgesOnVertex(IVertex, J);
-         for (int KVec = 0; KVec < KLen; ++KVec) {
-            const int K = KStart + KVec;
+         const int KStartEdge =
+             Kokkos::max(KStartVertex, MinLayerEdgeTop(JEdge));
+         const int KEndEdge = Kokkos::min(KEndVertex, MaxLayerEdgeBot(JEdge));
+         for (int K = KStartEdge; K <= KEndEdge; ++K) {
+            const int KVec = K - KStartVertex;
             Del2RelVortVertexTmp[KVec] += InvAreaTriangle * DcEdge(JEdge) *
                                           EdgeSignOnVertex(IVertex, J) *
                                           Del2Edge(JEdge, K);
          }
       }
 
-      for (int KVec = 0; KVec < KLen; ++KVec) {
-         const int K                   = KStart + KVec;
+      for (int KVec = 0; KVec < KLenVertex; ++KVec) {
+         const int K                   = KStartVertex + KVec;
          Del2RelVortVertex(IVertex, K) = Del2RelVortVertexTmp[KVec];
       }
    }
@@ -118,8 +130,10 @@ class VelocityDel2AuxVars {
    I4 VertexDegree;
    Array1DI4 MinLayerEdgeBot;
    Array1DI4 MaxLayerEdgeTop;
-   Array1DI4 MinLayerVertexBot;
-   Array1DI4 MaxLayerVertexTop;
+   Array1DI4 MinLayerEdgeTop;
+   Array1DI4 MaxLayerEdgeBot;
+   Array1DI4 MinLayerVertexTop;
+   Array1DI4 MaxLayerVertexBot;
    Array1DI4 MinLayerCell;
    Array1DI4 MaxLayerCell;
 };