From 6b72c7d42f4b3b3afd6d36fa6c44677d99f1bb0f Mon Sep 17 00:00:00 2001 From: Hyun Kang Date: Thu, 23 Jul 2026 23:28:19 -0400 Subject: [PATCH 1/5] Add mesh scaling computation - Add mesh-density-based scaling for del2 and del4 viscosity/diffusion - Follow the MPAS-Ocean approach for computing mesh scaling --- components/omega/configs/Default.yml | 5 + components/omega/doc/userGuide/HorzMesh.md | 23 +++++ components/omega/src/ocn/HorzMesh.cpp | 98 +++++++++++++++++-- components/omega/test/ocn/HorzMeshTest.cpp | 104 +++++++++++++++++++++ 4 files changed, 223 insertions(+), 7 deletions(-) diff --git a/components/omega/configs/Default.yml b/components/omega/configs/Default.yml index 9f589f9466cf..8519e9973bdd 100644 --- a/components/omega/configs/Default.yml +++ b/components/omega/configs/Default.yml @@ -71,6 +71,11 @@ Omega: PressureGradTendencyEnable: true VelVertMixTendencyEnable: true TracerVertMixTendencyEnable: true + Hmix: + HmixScaleWithMesh: false + MaxMeshDensity: -1.0 + HmixUseRefWidth: false + HmixRefWidth: 30.0e3 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..5365cc555e9e 100644 --- a/components/omega/doc/userGuide/HorzMesh.md +++ b/components/omega/doc/userGuide/HorzMesh.md @@ -62,6 +62,29 @@ 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 `Hmix` configuration controls whether Laplacian and biharmonic mixing +coefficients vary with horizontal mesh size: + +```yaml + Hmix: + HmixScaleWithMesh: false + MaxMeshDensity: -1.0 + HmixUseRefWidth: false + HmixRefWidth: 30.0e3 +``` + +When `HmixScaleWithMesh` is false, both scaling coefficients are one. When it +is true and `HmixUseRefWidth` 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 `HmixRefWidth`, and the biharmonic scaling is the +cube of that ratio. + +When `HmixUseRefWidth` 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..6f0aa8859d56 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 { @@ -393,16 +395,98 @@ void HorzMesh::computeEdgeSign() { // equations so viscosity and diffusion scale with mesh. void HorzMesh::computeMeshScaling() { + Config *OmegaConfig = Config::getOmegaConfig(); + Config HmixConfig("Hmix"); + Error Err = OmegaConfig->get(HmixConfig); + CHECK_ERROR_ABORT(Err, + "HorzMesh: Hmix group not found in input configuration"); + + bool ScaleWithMesh; + bool UseRefWidth; + Real MaxMeshDensity; + Real RefWidth; + + Err += HmixConfig.get("HmixScaleWithMesh", ScaleWithMesh); + Err += HmixConfig.get("MaxMeshDensity", MaxMeshDensity); + Err += HmixConfig.get("HmixUseRefWidth", UseRefWidth); + Err += HmixConfig.get("HmixRefWidth", RefWidth); + 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 && UseRefWidth) { + OMEGA_REQUIRE(RefWidth > 0.0_Real, + "HorzMesh: HmixRefWidth must be positive, got {}", + RefWidth); + + OMEGA_SCOPE(o_AreaCell, AreaCell); + OMEGA_SCOPE(o_CellsOnEdge, CellsOnEdge); + OMEGA_SCOPE(o_RefWidth, RefWidth); + + // Compute an effective cell width at each edge by treating the two + // adjacent cells as circles, following MPAS-Ocean. + 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 if (ScaleWithMesh) { + // A negative configured value requests the global maximum of the mesh + // density, matching MPAS-Ocean's legacy scaling behavior. + 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()); + HmixConfig.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); + }); + + } else { + parallelFor( + {NEdgesAll}, KOKKOS_LAMBDA(int Edge) { + o_MeshScalingDel2(Edge) = 1.0_Real; + o_MeshScalingDel4(Edge) = 1.0_Real; + }); + } MeshScalingDel2H = createHostMirrorCopy(MeshScalingDel2); MeshScalingDel4H = createHostMirrorCopy(MeshScalingDel4); diff --git a/components/omega/test/ocn/HorzMeshTest.cpp b/components/omega/test/ocn/HorzMeshTest.cpp index a1849fb0cebd..574594b6a2f0 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,107 @@ 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 < LocEdges; ++Edge) { + if (Mesh->MeshScalingDel2H(Edge) != 1.0_Real || + Mesh->MeshScalingDel4H(Edge) != 1.0_Real) { + ++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 HmixConfig("Hmix"); + Error ConfigErr = OmegaConfig->get(HmixConfig); + CHECK_ERROR_ABORT(ConfigErr, + "HorzMeshTest: Hmix group not found in Config"); + + // Test scaling based on a configured reference cell width + const Real RefWidth = 30000.0_Real; + HmixConfig.set("HmixScaleWithMesh", true); + HmixConfig.set("HmixUseRefWidth", true); + HmixConfig.set("HmixRefWidth", RefWidth); + Mesh->computeMeshScaling(); + + Count = 0; + for (int Edge = 0; Edge < LocEdges; ++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. + HmixConfig.set("HmixUseRefWidth", false); + HmixConfig.set("MaxMeshDensity", -1.0); + Mesh->computeMeshScaling(); + + Real MaxMeshDensity; + ConfigErr = HmixConfig.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 < LocEdges; ++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 From e4f9ff45fd07a9330fc45c6283237012d372ddd0 Mon Sep 17 00:00:00 2001 From: Hyun Kang Date: Mon, 10 Aug 2026 09:23:51 -0400 Subject: [PATCH 2/5] Change config options from Hmix to MeshScaling --- components/omega/configs/Default.yml | 11 ++++--- components/omega/doc/userGuide/HorzMesh.md | 25 ++++++++------- components/omega/src/ocn/HorzMesh.cpp | 24 ++++++++------ components/omega/test/ocn/HorzMeshTest.cpp | 37 ++++++++++++---------- 4 files changed, 54 insertions(+), 43 deletions(-) diff --git a/components/omega/configs/Default.yml b/components/omega/configs/Default.yml index 8519e9973bdd..f04f54eb59d5 100644 --- a/components/omega/configs/Default.yml +++ b/components/omega/configs/Default.yml @@ -71,11 +71,12 @@ Omega: PressureGradTendencyEnable: true VelVertMixTendencyEnable: true TracerVertMixTendencyEnable: true - Hmix: - HmixScaleWithMesh: false - MaxMeshDensity: -1.0 - HmixUseRefWidth: false - HmixRefWidth: 30.0e3 + HorzMesh: + MeshScaling: + ScaleWithMesh: false + UseRefWidth: false + 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 5365cc555e9e..edbd1e208514 100644 --- a/components/omega/doc/userGuide/HorzMesh.md +++ b/components/omega/doc/userGuide/HorzMesh.md @@ -64,24 +64,25 @@ boolean flags in the code. ### Horizontal mixing mesh scaling -The `Hmix` configuration controls whether Laplacian and biharmonic mixing -coefficients vary with horizontal mesh size: +The `HorzMesh: MeshScaling` configuration controls whether Laplacian and +biharmonic mixing coefficients vary with horizontal mesh size: ```yaml - Hmix: - HmixScaleWithMesh: false - MaxMeshDensity: -1.0 - HmixUseRefWidth: false - HmixRefWidth: 30.0e3 + HorzMesh: + MeshScaling: + ScaleWithMesh: false + UseRefWidth: false + RefWidth: 30.0e3 + MaxMeshDensity: -1.0 ``` -When `HmixScaleWithMesh` is false, both scaling coefficients are one. When it -is true and `HmixUseRefWidth` 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 `HmixRefWidth`, and the biharmonic scaling is the +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 `HmixUseRefWidth` is false, Omega uses the legacy MPAS-Ocean scaling based +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. diff --git a/components/omega/src/ocn/HorzMesh.cpp b/components/omega/src/ocn/HorzMesh.cpp index 6f0aa8859d56..8e40f2faab14 100644 --- a/components/omega/src/ocn/HorzMesh.cpp +++ b/components/omega/src/ocn/HorzMesh.cpp @@ -396,20 +396,25 @@ void HorzMesh::computeEdgeSign() { void HorzMesh::computeMeshScaling() { Config *OmegaConfig = Config::getOmegaConfig(); - Config HmixConfig("Hmix"); - Error Err = OmegaConfig->get(HmixConfig); + Config HorzMeshConfig("HorzMesh"); + Error Err = OmegaConfig->get(HorzMeshConfig); CHECK_ERROR_ABORT(Err, - "HorzMesh: Hmix group not found in input configuration"); + "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 += HmixConfig.get("HmixScaleWithMesh", ScaleWithMesh); - Err += HmixConfig.get("MaxMeshDensity", MaxMeshDensity); - Err += HmixConfig.get("HmixUseRefWidth", UseRefWidth); - Err += HmixConfig.get("HmixRefWidth", 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"); @@ -418,8 +423,7 @@ void HorzMesh::computeMeshScaling() { if (ScaleWithMesh && UseRefWidth) { OMEGA_REQUIRE(RefWidth > 0.0_Real, - "HorzMesh: HmixRefWidth must be positive, got {}", - RefWidth); + "HorzMesh: RefWidth must be positive, got {}", RefWidth); OMEGA_SCOPE(o_AreaCell, AreaCell); OMEGA_SCOPE(o_CellsOnEdge, CellsOnEdge); @@ -454,7 +458,7 @@ void HorzMesh::computeMeshScaling() { Halo *HorzMeshHalo = Halo::get(MeshName); MaxMeshDensity = globalMaxVal(MaxMeshDensityLocal, HorzMeshHalo->getComm()); - HmixConfig.set("MaxMeshDensity", MaxMeshDensity); + MeshScalingConfig.set("MaxMeshDensity", MaxMeshDensity); } OMEGA_REQUIRE(MaxMeshDensity > 0.0_Real, diff --git a/components/omega/test/ocn/HorzMeshTest.cpp b/components/omega/test/ocn/HorzMeshTest.cpp index 574594b6a2f0..14eb9e407197 100644 --- a/components/omega/test/ocn/HorzMeshTest.cpp +++ b/components/omega/test/ocn/HorzMeshTest.cpp @@ -579,11 +579,11 @@ int main(int argc, char *argv[]) { // Test mesh scaling with scaling disabled Count = 0; - for (int Edge = 0; Edge < LocEdges; ++Edge) { - if (Mesh->MeshScalingDel2H(Edge) != 1.0_Real || - Mesh->MeshScalingDel4H(Edge) != 1.0_Real) { - ++Count; - } + 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) @@ -592,20 +592,24 @@ int main(int argc, char *argv[]) { // Retrieve mesh scaling configuration so the two MPAS-O scaling modes // can be tested independently. Config *OmegaConfig = Config::getOmegaConfig(); - Config HmixConfig("Hmix"); - Error ConfigErr = OmegaConfig->get(HmixConfig); + Config HorzMeshConfig("HorzMesh"); + Error ConfigErr = OmegaConfig->get(HorzMeshConfig); CHECK_ERROR_ABORT(ConfigErr, - "HorzMeshTest: Hmix group not found in Config"); + "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; - HmixConfig.set("HmixScaleWithMesh", true); - HmixConfig.set("HmixUseRefWidth", true); - HmixConfig.set("HmixRefWidth", RefWidth); + MeshScalingConfig.set("ScaleWithMesh", true); + MeshScalingConfig.set("UseRefWidth", true); + MeshScalingConfig.set("RefWidth", RefWidth); Mesh->computeMeshScaling(); Count = 0; - for (int Edge = 0; Edge < LocEdges; ++Edge) { + 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 = @@ -628,12 +632,13 @@ int main(int argc, char *argv[]) { // Test legacy scaling based on mesh density. A negative configured // maximum requests that the global maximum be computed from the mesh. - HmixConfig.set("HmixUseRefWidth", false); - HmixConfig.set("MaxMeshDensity", -1.0); + MeshScalingConfig.set("UseRefWidth", false); + MeshScalingConfig.set("MaxMeshDensity", -1.0); Mesh->computeMeshScaling(); Real MaxMeshDensity; - ConfigErr = HmixConfig.get("MaxMeshDensity", MaxMeshDensity); + ConfigErr = + MeshScalingConfig.get("MaxMeshDensity", MaxMeshDensity); CHECK_ERROR_ABORT( ConfigErr, "HorzMeshTest: unable to retrieve computed MaxMeshDensity"); @@ -655,7 +660,7 @@ int main(int argc, char *argv[]) { } Count = 0; - for (int Edge = 0; Edge < LocEdges; ++Edge) { + 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 = From 268f173cc526595a3406369a4cd56d7dcf12817b Mon Sep 17 00:00:00 2001 From: Hyun Kang Date: Mon, 10 Aug 2026 09:45:53 -0400 Subject: [PATCH 3/5] Update comments --- components/omega/src/ocn/HorzMesh.cpp | 90 ++++++++++++++------------- 1 file changed, 48 insertions(+), 42 deletions(-) diff --git a/components/omega/src/ocn/HorzMesh.cpp b/components/omega/src/ocn/HorzMesh.cpp index 8e40f2faab14..4410ae07d5fb 100644 --- a/components/omega/src/ocn/HorzMesh.cpp +++ b/components/omega/src/ocn/HorzMesh.cpp @@ -392,7 +392,7 @@ 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(); @@ -421,16 +421,18 @@ void HorzMesh::computeMeshScaling() { OMEGA_SCOPE(o_MeshScalingDel2, MeshScalingDel2); OMEGA_SCOPE(o_MeshScalingDel4, MeshScalingDel4); - if (ScaleWithMesh && UseRefWidth) { + 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); - OMEGA_SCOPE(o_RefWidth, RefWidth); - // Compute an effective cell width at each edge by treating the two - // adjacent cells as circles, following MPAS-Ocean. + // 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); @@ -446,43 +448,47 @@ void HorzMesh::computeMeshScaling() { Del2Scale * Del2Scale * Del2Scale; }); - } else if (ScaleWithMesh) { - // A negative configured value requests the global maximum of the mesh - // density, matching MPAS-Ocean's legacy scaling behavior. - 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; + } 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); + } - o_MeshScalingDel2(Edge) = - 1.0_Real / Kokkos::pow(DensityRatio, 0.25_Real); - o_MeshScalingDel4(Edge) = - 1.0_Real / Kokkos::pow(DensityRatio, 0.75_Real); - }); + 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( @@ -490,7 +496,7 @@ void HorzMesh::computeMeshScaling() { o_MeshScalingDel2(Edge) = 1.0_Real; o_MeshScalingDel4(Edge) = 1.0_Real; }); - } + } // if ScaleWithMesh MeshScalingDel2H = createHostMirrorCopy(MeshScalingDel2); MeshScalingDel4H = createHostMirrorCopy(MeshScalingDel4); From 54b0b66a26844163280237d4ab8a42a966a29640 Mon Sep 17 00:00:00 2001 From: Hyun Kang Date: Mon, 10 Aug 2026 09:48:24 -0400 Subject: [PATCH 4/5] Linting --- components/omega/src/ocn/HorzMesh.cpp | 53 ++++++++++------------ components/omega/test/ocn/HorzMeshTest.cpp | 16 +++---- 2 files changed, 31 insertions(+), 38 deletions(-) diff --git a/components/omega/src/ocn/HorzMesh.cpp b/components/omega/src/ocn/HorzMesh.cpp index 4410ae07d5fb..65e023d4dfdf 100644 --- a/components/omega/src/ocn/HorzMesh.cpp +++ b/components/omega/src/ocn/HorzMesh.cpp @@ -403,8 +403,8 @@ void HorzMesh::computeMeshScaling() { Config MeshScalingConfig("MeshScaling"); Err = HorzMeshConfig.get(MeshScalingConfig); - CHECK_ERROR_ABORT( - Err, "HorzMesh: MeshScaling group not found in configuration"); + CHECK_ERROR_ABORT(Err, + "HorzMesh: MeshScaling group not found in configuration"); bool ScaleWithMesh; bool UseRefWidth; @@ -415,38 +415,36 @@ void HorzMesh::computeMeshScaling() { 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"); + CHECK_ERROR_ABORT(Err, "HorzMesh: error reading mesh scaling configuration"); OMEGA_SCOPE(o_MeshScalingDel2, MeshScalingDel2); OMEGA_SCOPE(o_MeshScalingDel4, MeshScalingDel4); if (ScaleWithMesh) { if (UseRefWidth) { - OMEGA_REQUIRE(RefWidth > 0.0_Real, - "HorzMesh: RefWidth must be positive, got {}", RefWidth); + 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); + 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; - }); + // 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 @@ -479,8 +477,7 @@ void HorzMesh::computeMeshScaling() { 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)); + 0.5_Real * (o_MeshDensity(Cell0) + o_MeshDensity(Cell1)); const Real DensityRatio = AvgDensity / o_MaxMeshDensity; o_MeshScalingDel2(Edge) = diff --git a/components/omega/test/ocn/HorzMeshTest.cpp b/components/omega/test/ocn/HorzMeshTest.cpp index 14eb9e407197..a9cc37f69765 100644 --- a/components/omega/test/ocn/HorzMeshTest.cpp +++ b/components/omega/test/ocn/HorzMeshTest.cpp @@ -598,8 +598,8 @@ int main(int argc, char *argv[]) { "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"); + 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; @@ -637,8 +637,7 @@ int main(int argc, char *argv[]) { Mesh->computeMeshScaling(); Real MaxMeshDensity; - ConfigErr = - MeshScalingConfig.get("MaxMeshDensity", MaxMeshDensity); + ConfigErr = MeshScalingConfig.get("MaxMeshDensity", MaxMeshDensity); CHECK_ERROR_ABORT( ConfigErr, "HorzMeshTest: unable to retrieve computed MaxMeshDensity"); @@ -664,13 +663,10 @@ int main(int argc, char *argv[]) { 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)); + 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); + 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)) || From 52f4931615da6604c7a600b88bf66a2865cebef5 Mon Sep 17 00:00:00 2001 From: Hyun Kang Date: Mon, 10 Aug 2026 10:15:15 -0400 Subject: [PATCH 5/5] Set UseRefWidth to true as the default value --- components/omega/configs/Default.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/components/omega/configs/Default.yml b/components/omega/configs/Default.yml index f04f54eb59d5..34bba95c555e 100644 --- a/components/omega/configs/Default.yml +++ b/components/omega/configs/Default.yml @@ -74,7 +74,7 @@ Omega: HorzMesh: MeshScaling: ScaleWithMesh: false - UseRefWidth: false + UseRefWidth: true RefWidth: 30.0e3 MaxMeshDensity: -1.0 ManufacturedSolution: