Compute MeshScaling for variable-resolution meshes

components/omega/configs/Default.yml

@@ -71,6 +71,12 @@ Omega:
     PressureGradTendencyEnable: true
     VelVertMixTendencyEnable: true
     TracerVertMixTendencyEnable: true
+  HorzMesh:
+    MeshScaling:
+      ScaleWithMesh: false
+      UseRefWidth: true
+      RefWidth: 30.0e3
+      MaxMeshDensity: -1.0
   ManufacturedSolution:
     WavelengthX: 5.0e6
     WavelengthY: 4.33013e6

components/omega/doc/userGuide/HorzMesh.md

@@ -62,6 +62,30 @@ variables. These include ``OnSphere`` (or ``on_a_sphere`` for backcompatibility)
 OnSphere and IsPeriodic are stored as YES/NO strings in the metadata but as
 boolean flags in the code.
 
+### Horizontal mixing mesh scaling
+
+The `HorzMesh: MeshScaling` configuration controls whether Laplacian and
+biharmonic mixing coefficients vary with horizontal mesh size:
+
+```yaml
+  HorzMesh:
+    MeshScaling:
+      ScaleWithMesh: false
+      UseRefWidth: false
+      RefWidth: 30.0e3
+      MaxMeshDensity: -1.0
+```
+
+When `ScaleWithMesh` is false, both scaling coefficients are one. When it is
+true and `UseRefWidth` is true, Omega computes the effective cell width at
+each edge from the areas of its two adjacent cells. The Laplacian scaling is
+the ratio of this width to `RefWidth`, and the biharmonic scaling is the
+cube of that ratio.
+
+When `UseRefWidth` is false, Omega uses the legacy MPAS-Ocean scaling based
+on `MeshDensity`. A negative `MaxMeshDensity` requests that Omega compute its
+global value from the mesh during initialization.
+
 In the future, the Mesh class will optionally compute the mesh variables that
 are dependent on the Cartesian mesh coordinates internally.
 This includes the various areas, lengths, angles, and weights needed for the

components/omega/src/ocn/HorzMesh.cpp

@@ -10,6 +10,7 @@
 //===----------------------------------------------------------------------===//
 
 #include "HorzMesh.h"
+#include "Config.h"
 #include "DataTypes.h"
 #include "Decomp.h"
 #include "Dimension.h"
@@ -19,6 +20,7 @@
 #include "Halo.h"
 #include "IOStream.h"
 #include "OmegaKokkos.h"
+#include "Reductions.h"
 
 namespace OMEGA {
 
@@ -390,19 +392,108 @@ void HorzMesh::computeEdgeSign() {
 
 //------------------------------------------------------------------------------
 // Set mesh scaling coefficients for mixing terms in momentum and tracer
-// equations so viscosity and diffusion scale with mesh.
+// equations so viscosity and diffusion scale with mesh, following MPAS-Ocean
 void HorzMesh::computeMeshScaling() {
 
+   Config *OmegaConfig = Config::getOmegaConfig();
+   Config HorzMeshConfig("HorzMesh");
+   Error Err = OmegaConfig->get(HorzMeshConfig);
+   CHECK_ERROR_ABORT(Err,
+                     "HorzMesh: HorzMesh group not found in configuration");
+
+   Config MeshScalingConfig("MeshScaling");
+   Err = HorzMeshConfig.get(MeshScalingConfig);
+   CHECK_ERROR_ABORT(Err,
+                     "HorzMesh: MeshScaling group not found in configuration");
+
+   bool ScaleWithMesh;
+   bool UseRefWidth;
+   Real MaxMeshDensity;
+   Real RefWidth;
+
+   Err += MeshScalingConfig.get("ScaleWithMesh", ScaleWithMesh);
+   Err += MeshScalingConfig.get("UseRefWidth", UseRefWidth);
+   Err += MeshScalingConfig.get("RefWidth", RefWidth);
+   Err += MeshScalingConfig.get("MaxMeshDensity", MaxMeshDensity);
+   CHECK_ERROR_ABORT(Err, "HorzMesh: error reading mesh scaling configuration");
+
    OMEGA_SCOPE(o_MeshScalingDel2, MeshScalingDel2);
    OMEGA_SCOPE(o_MeshScalingDel4, MeshScalingDel4);
 
-   // TODO: implement mesh scaling by cell area, only no scaling
-   // option for now
-   parallelFor(
-       {NEdgesAll}, KOKKOS_LAMBDA(int Edge) {
-          o_MeshScalingDel2(Edge) = 1.0;
-          o_MeshScalingDel4(Edge) = 1.0;
-       });
+   if (ScaleWithMesh) {
+      if (UseRefWidth) {
+         OMEGA_REQUIRE(RefWidth > 0.0_Real,
+                       "HorzMesh: RefWidth must be positive, got {}", RefWidth);
+
+         OMEGA_SCOPE(o_RefWidth, RefWidth);
+         OMEGA_SCOPE(o_AreaCell, AreaCell);
+         OMEGA_SCOPE(o_CellsOnEdge, CellsOnEdge);
+
+         // Mesh scaling is computed using CellWidth derived from AreaCell
+         // and the input reference width, RefWidth (see Eqs. (1) and (2) of
+         // Hoch et al. 2020, JAMES).
+         parallelFor(
+             {NEdgesAll}, KOKKOS_LAMBDA(int Edge) {
+                const int Cell0 = o_CellsOnEdge(Edge, 0);
+                const int Cell1 = o_CellsOnEdge(Edge, 1);
+                const Real CellWidth =
+                    2.0_Real *
+                    Kokkos::sqrt((o_AreaCell(Cell0) + o_AreaCell(Cell1)) /
+                                 (2.0_Real * Pi));
+                const Real Del2Scale = CellWidth / o_RefWidth;
+
+                o_MeshScalingDel2(Edge) = Del2Scale;
+                o_MeshScalingDel4(Edge) = Del2Scale * Del2Scale * Del2Scale;
+             });
+
+      } else {
+         // Mesh scaling is set by MeshDensity. This is both confusing
+         // and inconvenient, as the flags like ViscDel2 need
+         // to be reset for every resolution. It is kept for backwards
+         // compatibility, but should become defunct.
+         // Del2 scales as Dc**1, Del4 scales as Dc**3
+         if (MaxMeshDensity < 0.0_Real) {
+            Real MaxMeshDensityLocal = 0.0_Real;
+            for (int Cell = 0; Cell < NCellsOwned; ++Cell)
+               MaxMeshDensityLocal =
+                   std::max(MaxMeshDensityLocal, MeshDensityH(Cell));
+
+            Halo *HorzMeshHalo = Halo::get(MeshName);
+            MaxMeshDensity =
+                globalMaxVal(MaxMeshDensityLocal, HorzMeshHalo->getComm());
+            MeshScalingConfig.set("MaxMeshDensity", MaxMeshDensity);
+         }
+
+         OMEGA_REQUIRE(MaxMeshDensity > 0.0_Real,
+                       "HorzMesh: MaxMeshDensity must be positive, got {}",
+                       MaxMeshDensity);
+
+         OMEGA_SCOPE(o_CellsOnEdge, CellsOnEdge);
+         OMEGA_SCOPE(o_MeshDensity, MeshDensity);
+         OMEGA_SCOPE(o_MaxMeshDensity, MaxMeshDensity);
+
+         parallelFor(
+             {NEdgesAll}, KOKKOS_LAMBDA(int Edge) {
+                const int Cell0 = o_CellsOnEdge(Edge, 0);
+                const int Cell1 = o_CellsOnEdge(Edge, 1);
+                const Real AvgDensity =
+                    0.5_Real * (o_MeshDensity(Cell0) + o_MeshDensity(Cell1));
+                const Real DensityRatio = AvgDensity / o_MaxMeshDensity;
+
+                o_MeshScalingDel2(Edge) =
+                    1.0_Real / Kokkos::pow(DensityRatio, 0.25_Real);
+                o_MeshScalingDel4(Edge) =
+                    1.0_Real / Kokkos::pow(DensityRatio, 0.75_Real);
+             });
+      } // if UseRefWidth
+
+   } else {
+      parallelFor(
+          {NEdgesAll}, KOKKOS_LAMBDA(int Edge) {
+             o_MeshScalingDel2(Edge) = 1.0_Real;
+             o_MeshScalingDel4(Edge) = 1.0_Real;
+          });
+   } // if ScaleWithMesh
 
    MeshScalingDel2H = createHostMirrorCopy(MeshScalingDel2);
    MeshScalingDel4H = createHostMirrorCopy(MeshScalingDel4);

components/omega/test/ocn/HorzMeshTest.cpp

@@ -25,6 +25,9 @@
 #include "TimeMgr.h"
 #include "mpi.h"
 
+#include <algorithm>
+#include <cmath>
+
 using namespace OMEGA;
 
 //------------------------------------------------------------------------------
@@ -574,6 +577,108 @@ int main(int argc, char *argv[]) {
       if (Count > 0)
          ABORT_ERROR("HorzMeshTest: edgeSignOnVertex test FAIL");
 
+      // Test mesh scaling with scaling disabled
+      Count = 0;
+      for (int Edge = 0; Edge < Mesh->NEdgesAll; Edge++) {
+         if (abs(Mesh->MeshScalingDel2H(Edge) - 1.0) > Tol)
+            Count++;
+         if (abs(Mesh->MeshScalingDel4H(Edge) - 1.0) > Tol)
+            Count++;
+      }
+
+      if (Count > 0)
+         ABORT_ERROR("HorzMeshTest: disabled mesh scaling test FAIL");
+
+      // Retrieve mesh scaling configuration so the two MPAS-O scaling modes
+      // can be tested independently.
+      Config *OmegaConfig = Config::getOmegaConfig();
+      Config HorzMeshConfig("HorzMesh");
+      Error ConfigErr = OmegaConfig->get(HorzMeshConfig);
+      CHECK_ERROR_ABORT(ConfigErr,
+                        "HorzMeshTest: HorzMesh group not found in Config");
+      Config MeshScalingConfig("MeshScaling");
+      ConfigErr = HorzMeshConfig.get(MeshScalingConfig);
+      CHECK_ERROR_ABORT(ConfigErr,
+                        "HorzMeshTest: MeshScaling group not found in Config");
+
+      // Test scaling based on a configured reference cell width
+      const Real RefWidth = 30000.0_Real;
+      MeshScalingConfig.set("ScaleWithMesh", true);
+      MeshScalingConfig.set("UseRefWidth", true);
+      MeshScalingConfig.set("RefWidth", RefWidth);
+      Mesh->computeMeshScaling();
+
+      Count = 0;
+      for (int Edge = 0; Edge < Mesh->NEdgesAll; Edge++) {
+         const int Cell0 = Mesh->CellsOnEdgeH(Edge, 0);
+         const int Cell1 = Mesh->CellsOnEdgeH(Edge, 1);
+         const Real CellWidth =
+             2.0_Real *
+             std::sqrt((Mesh->AreaCellH(Cell0) + Mesh->AreaCellH(Cell1)) /
+                       (2.0_Real * Pi));
+         const Real RefDel2 = CellWidth / RefWidth;
+         const Real RefDel4 = RefDel2 * RefDel2 * RefDel2;
+
+         if (abs(Mesh->MeshScalingDel2H(Edge) - RefDel2) >
+                 Tol * std::max(1.0_Real, std::abs(RefDel2)) ||
+             abs(Mesh->MeshScalingDel4H(Edge) - RefDel4) >
+                 Tol * std::max(1.0_Real, std::abs(RefDel4))) {
+            ++Count;
+         }
+      }
+
+      if (Count > 0)
+         ABORT_ERROR("HorzMeshTest: reference-width mesh scaling test FAIL");
+
+      // Test legacy scaling based on mesh density. A negative configured
+      // maximum requests that the global maximum be computed from the mesh.
+      MeshScalingConfig.set("UseRefWidth", false);
+      MeshScalingConfig.set("MaxMeshDensity", -1.0);
+      Mesh->computeMeshScaling();
+
+      Real MaxMeshDensity;
+      ConfigErr = MeshScalingConfig.get("MaxMeshDensity", MaxMeshDensity);
+      CHECK_ERROR_ABORT(
+          ConfigErr,
+          "HorzMeshTest: unable to retrieve computed MaxMeshDensity");
+
+      Real RefMaxMeshDensityLocal = 0.0_Real;
+      for (int Cell = 0; Cell < LocCells; ++Cell) {
+         RefMaxMeshDensityLocal =
+             std::max(RefMaxMeshDensityLocal, Mesh->MeshDensityH(Cell));
+      }
+      Real RefMaxMeshDensity;
+      Err = MPI_Allreduce(&RefMaxMeshDensityLocal, &RefMaxMeshDensity, 1,
+                          MPI_RealKind, MPI_MAX, Comm);
+      if (Err != MPI_SUCCESS)
+         ABORT_ERROR("HorzMeshTest: MPI error finding max mesh density");
+
+      if (abs(MaxMeshDensity - RefMaxMeshDensity) >
+          Tol * std::max(1.0_Real, std::abs(RefMaxMeshDensity))) {
+         ABORT_ERROR("HorzMeshTest: computed MaxMeshDensity test FAIL");
+      }
+
+      Count = 0;
+      for (int Edge = 0; Edge < Mesh->NEdgesAll; Edge++) {
+         const int Cell0 = Mesh->CellsOnEdgeH(Edge, 0);
+         const int Cell1 = Mesh->CellsOnEdgeH(Edge, 1);
+         const Real AvgDensity =
+             0.5_Real * (Mesh->MeshDensityH(Cell0) + Mesh->MeshDensityH(Cell1));
+         const Real DensityRatio = AvgDensity / RefMaxMeshDensity;
+         const Real RefDel2      = 1.0_Real / std::pow(DensityRatio, 0.25_Real);
+         const Real RefDel4      = 1.0_Real / std::pow(DensityRatio, 0.75_Real);
+
+         if (abs(Mesh->MeshScalingDel2H(Edge) - RefDel2) >
+                 Tol * std::max(1.0_Real, std::abs(RefDel2)) ||
+             abs(Mesh->MeshScalingDel4H(Edge) - RefDel4) >
+                 Tol * std::max(1.0_Real, std::abs(RefDel4))) {
+            ++Count;
+         }
+      }
+
+      if (Count > 0)
+         ABORT_ERROR("HorzMeshTest: mesh-density scaling test FAIL");
+
       // Test cell halo values
       // Perform halo exhange on owned cell only array and compare
       // read values