diff --git a/components/omega/configs/Default.yml b/components/omega/configs/Default.yml index 9f589f9466cf..34bba95c555e 100644 --- a/components/omega/configs/Default.yml +++ b/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 diff --git a/components/omega/doc/userGuide/HorzMesh.md b/components/omega/doc/userGuide/HorzMesh.md index 475ff844b07a..edbd1e208514 100644 --- a/components/omega/doc/userGuide/HorzMesh.md +++ b/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 diff --git a/components/omega/src/ocn/HorzMesh.cpp b/components/omega/src/ocn/HorzMesh.cpp index bb8b909aad42..65e023d4dfdf 100644 --- a/components/omega/src/ocn/HorzMesh.cpp +++ b/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); diff --git a/components/omega/test/ocn/HorzMeshTest.cpp b/components/omega/test/ocn/HorzMeshTest.cpp index a1849fb0cebd..a9cc37f69765 100644 --- a/components/omega/test/ocn/HorzMeshTest.cpp +++ b/components/omega/test/ocn/HorzMeshTest.cpp @@ -25,6 +25,9 @@ #include "TimeMgr.h" #include "mpi.h" +#include +#include + 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