From 00e51764fa255739e1b989160dc060e888e2f48b Mon Sep 17 00:00:00 2001 From: Xylar Asay-Davis Date: Fri, 11 Sep 2026 16:49:38 -0500 Subject: [PATCH 01/12] Use the mesh sphere radius in operators and prescribed velocities SecondDerivativeOnCell and the prescribed nondivergent and divergent velocities used the REarth constant to scale spherical geometry, while the rest of the mesh geometry comes from the input file. Take the radius from HorzMesh::SphereRadius instead, so a single value describes the sphere the mesh lies on; HorzMesh checks it against REarth on read. DetermineSphericalPatchGeometry now receives the radius as an argument, and the unit test that projects a synthetic hexagonal patch onto a sphere passes REarth explicitly, since it builds its own sphere and also runs on the planar mesh. Also correct the comment on the divergent-flow period (12 days, not 14). Co-Authored-By: Claude Opus 5 (1M context) --- components/omega/src/ocn/HorzOperators.cpp | 13 +++++++------ components/omega/src/ocn/HorzOperators.h | 14 ++++++++------ components/omega/src/timeStepping/TimeStepper.cpp | 12 +++++++----- components/omega/test/ocn/HorzOperatorsTest.cpp | 15 +++++++++------ 4 files changed, 31 insertions(+), 23 deletions(-) diff --git a/components/omega/src/ocn/HorzOperators.cpp b/components/omega/src/ocn/HorzOperators.cpp index a530e8254dd1..4a01cd320c8d 100644 --- a/components/omega/src/ocn/HorzOperators.cpp +++ b/components/omega/src/ocn/HorzOperators.cpp @@ -30,12 +30,13 @@ InterpCellToEdge::InterpCellToEdge(const HorzMesh *Mesh) VertexDegree(Mesh->VertexDegree) {} SecondDerivativeOnCell::SecondDerivativeOnCell(HorzMesh const *Mesh) - : OnSphere(Mesh->OnSphere), NCellsAll(Mesh->NCellsAll), - MaxEdges(1 + Mesh->MaxEdges), NEdgesOnCell(Mesh->NEdgesOnCell), - EdgesOnCell(Mesh->EdgesOnCell), CellsOnCell(Mesh->CellsOnCell), - CellsOnEdge(Mesh->CellsOnEdge), VerticesOnEdge(Mesh->VerticesOnEdge), - XCell(Mesh->XCell), YCell(Mesh->YCell), ZCell(Mesh->ZCell), - DvEdge(Mesh->DvEdge), DcEdge(Mesh->DcEdge), AngleEdge(Mesh->AngleEdge), + : OnSphere(Mesh->OnSphere), SphereRadius(Mesh->SphereRadius), + NCellsAll(Mesh->NCellsAll), MaxEdges(1 + Mesh->MaxEdges), + NEdgesOnCell(Mesh->NEdgesOnCell), EdgesOnCell(Mesh->EdgesOnCell), + CellsOnCell(Mesh->CellsOnCell), CellsOnEdge(Mesh->CellsOnEdge), + VerticesOnEdge(Mesh->VerticesOnEdge), XCell(Mesh->XCell), + YCell(Mesh->YCell), ZCell(Mesh->ZCell), DvEdge(Mesh->DvEdge), + DcEdge(Mesh->DcEdge), AngleEdge(Mesh->AngleEdge), AreaCell(Mesh->AreaCell), EdgeSignOnCell(Mesh->EdgeSignOnCell), XVertex(Mesh->XVertex), YVertex(Mesh->YVertex), ZVertex(Mesh->ZVertex), ThetaAbs("SphereAngle", NCellsAll), XPCell("XP", NCellsAll, MaxEdges), diff --git a/components/omega/src/ocn/HorzOperators.h b/components/omega/src/ocn/HorzOperators.h index bcb96b12dcab..2bd395cd774b 100644 --- a/components/omega/src/ocn/HorzOperators.h +++ b/components/omega/src/ocn/HorzOperators.h @@ -235,9 +235,10 @@ class SecondDerivativeOnCell { if (OnSphere) { const Array1DI4 edgesOnCell = Kokkos::subview(EdgesOnCell, ICell, Kokkos::ALL); - DetermineSphericalPatchGeometry( - NEdges, edgesOnCell, VerticesOnEdge, XCell, YCell, ZCell, XVertex, - YVertex, ZVertex, CellList, XP, YP, Angle2D, ThetaAbs[ICell]); + DetermineSphericalPatchGeometry(NEdges, edgesOnCell, VerticesOnEdge, + XCell, YCell, ZCell, XVertex, YVertex, + ZVertex, CellList, SphereRadius, XP, + YP, Angle2D, ThetaAbs[ICell]); } else { // On an x-y plane const Array1DI4 edgesOnCell = Kokkos::subview(EdgesOnCell, ICell, Kokkos::ALL); @@ -270,6 +271,7 @@ class SecondDerivativeOnCell { static constexpr R8 Pii = 3.141592653589793_Real; const bool OnSphere; + const Real SphereRadius; ///< radius (m) of the sphere the mesh lies on const I4 NCellsAll; const I4 MaxEdges; Array1DI4 NEdgesOnCell; @@ -339,10 +341,10 @@ class SecondDerivativeOnCell { const Array2DI4 VerticesOnEdge, const Array1DReal XCell, const Array1DReal YCell, const Array1DReal ZCell, const Array1DReal XVertex, const Array1DReal YVertex, - const Array1DReal ZVertex, const Array1DI4 CellList, Array1DReal XP, - Array1DReal YP, Array1DReal Angle2D, Real &ThetaAbs) { + const Array1DReal ZVertex, const Array1DI4 CellList, + const Real sphereRadius, Array1DReal XP, Array1DReal YP, + Array1DReal Angle2D, Real &ThetaAbs) { const Real length_scale = 1._Real; - const Real sphereRadius = REarth; Real XC[MaxMaxEdges] = {}; Real YC[MaxMaxEdges] = {}; Real ZC[MaxMaxEdges] = {}; diff --git a/components/omega/src/timeStepping/TimeStepper.cpp b/components/omega/src/timeStepping/TimeStepper.cpp index 759fdaae3198..58e12908208a 100644 --- a/components/omega/src/timeStepping/TimeStepper.cpp +++ b/components/omega/src/timeStepping/TimeStepper.cpp @@ -639,6 +639,7 @@ void TimeStepper::prescribeVelocity(OceanState *State1, int TimeLevel1, OMEGA_SCOPE(LatEdge, Mesh->LatEdge); OMEGA_SCOPE(LonEdge, Mesh->LonEdge); OMEGA_SCOPE(AngleEdge, Mesh->AngleEdge); + const Real SphereRadius = Mesh->SphereRadius; OMEGA_SCOPE(MinLayerEdgeBot, VCoord->MinLayerEdgeBot); OMEGA_SCOPE(MaxLayerEdgeTop, VCoord->MaxLayerEdgeTop); @@ -664,8 +665,8 @@ void TimeStepper::prescribeVelocity(OceanState *State1, int TimeLevel1, 2.0 * Pi * cos(LatEdge(IEdge))); const R8 v = (10.0 / Tau) * sin(2.0 * lon_p) * cos(LatEdge(IEdge)) * cos(Pi * TSim / Tau); - const R8 normalVel = REarth * (u * cos(AngleEdge(IEdge)) + - v * sin(AngleEdge(IEdge))); + const R8 normalVel = SphereRadius * (u * cos(AngleEdge(IEdge)) + + v * sin(AngleEdge(IEdge))); parallelForInner( Team, KRange, INNER_LAMBDA(int KChunk) { @@ -680,6 +681,7 @@ void TimeStepper::prescribeVelocity(OceanState *State1, int TimeLevel1, OMEGA_SCOPE(LatEdge, Mesh->LatEdge); OMEGA_SCOPE(LonEdge, Mesh->LonEdge); OMEGA_SCOPE(AngleEdge, Mesh->AngleEdge); + const Real SphereRadius = Mesh->SphereRadius; OMEGA_SCOPE(MinLayerEdgeBot, VCoord->MinLayerEdgeBot); OMEGA_SCOPE(MaxLayerEdgeTop, VCoord->MaxLayerEdgeTop); @@ -688,7 +690,7 @@ void TimeStepper::prescribeVelocity(OceanState *State1, int TimeLevel1, TimeInterval ElapsedTimeInterval = SimTime - ModelClock->getStartTime(); ElapsedTimeInterval.get(ElapsedTimeSec, TimeUnits::Seconds); - const R8 Tau = 12. * Day2Sec; // 14 days in seconds + const R8 Tau = 12. * Day2Sec; // 12 days in seconds const R8 TSim = ElapsedTimeSec; parallelForOuter( @@ -708,8 +710,8 @@ void TimeStepper::prescribeVelocity(OceanState *State1, int TimeLevel1, const R8 v = ((2.5 / Tau) * sin(lon_p) * Kokkos::pow(cos(LatEdge(IEdge)), 3) * cos(Pi * TSim / Tau)); - const R8 normalVel = REarth * (u * cos(AngleEdge(IEdge)) + - v * sin(AngleEdge(IEdge))); + const R8 normalVel = SphereRadius * (u * cos(AngleEdge(IEdge)) + + v * sin(AngleEdge(IEdge))); parallelForInner( Team, KRange, INNER_LAMBDA(int KChunk) { diff --git a/components/omega/test/ocn/HorzOperatorsTest.cpp b/components/omega/test/ocn/HorzOperatorsTest.cpp index f4af88fd369c..347747359c24 100644 --- a/components/omega/test/ocn/HorzOperatorsTest.cpp +++ b/components/omega/test/ocn/HorzOperatorsTest.cpp @@ -751,8 +751,9 @@ int testsecondderivativeoncellDetermineSphericalPatchGeometry(Real RTol) { {20, 21, 25, 29, 28, 24}, {23, 24, 28, 32, 31, 27}, {26, 27, 31, 35, 34, 30}, {28, 29, 33, 37, 36, 32}}; const int CtoC = 1039230; - const double R = OMEGA::REarth; - const R8 Pii = 3.141592653589793_Real; + // synthetic patch on its own sphere, independent of the test mesh + const double R = OMEGA::REarth; + const R8 Pii = 3.141592653589793_Real; // Project coordinates to sphere. @@ -899,11 +900,13 @@ int testsecondderivativeoncellDetermineSphericalPatchGeometry(Real RTol) { const Array2DI4 VerticiesOnEdge, const Array1DReal XCell, const Array1DReal YCell, const Array1DReal ZCell, const Array1DReal XVertex, const Array1DReal YVertex, - const Array1DReal ZVertex, const Array1DI4 CellList, Array1DReal XP, - Array1DReal YP, Array1DReal Angle2D, Real &ThetaAbs) { + const Array1DReal ZVertex, const Array1DI4 CellList, + const Real SphereRadius, Array1DReal XP, Array1DReal YP, + Array1DReal Angle2D, Real &ThetaAbs) { SecondDerivativeOnCell::DetermineSphericalPatchGeometry( NEdges, EdgesOnCell, VerticiesOnEdge, XCell, YCell, ZCell, XVertex, - YVertex, ZVertex, CellList, XP, YP, Angle2D, ThetaAbs); + YVertex, ZVertex, CellList, SphereRadius, XP, YP, Angle2D, + ThetaAbs); } }; { @@ -925,7 +928,7 @@ int testsecondderivativeoncellDetermineSphericalPatchGeometry(Real RTol) { Real ThetaAbs = {}; SecondDerivativeOnCellTest::Test( NEdges, edgesOnCell, VerticesOnEdge, XCell, YCell, ZCell, - XVertex, YVertex, ZVertex, CellList, XP, YP, Angle2D, + XVertex, YVertex, ZVertex, CellList, R, XP, YP, Angle2D, ThetaAbs); }); auto XPH = createHostMirrorCopy(XP); From aa5cbe42b2cc818a64783acb4b6f1906f1e47d0d Mon Sep 17 00:00:00 2001 From: Xylar Asay-Davis Date: Fri, 11 Sep 2026 16:50:21 -0500 Subject: [PATCH 02/12] Take the test sphere radius from the mesh The spherical test setups evaluated their analytic fields with REarth while the mesh they compared against could lie on a slightly different sphere; the CTest icosahedral mesh was 3.3e-5 larger, which contributed a systematic error of that size to every spherical operator test. Read the radius from the default HorzMesh instead, which resolves the long-standing TODO in three of the setups. Co-Authored-By: Claude Opus 5 (1M context) --- components/omega/test/ocn/AuxiliaryStateTest.cpp | 3 ++- components/omega/test/ocn/AuxiliaryVarsTest.cpp | 3 ++- components/omega/test/ocn/HorzOperatorsTest.cpp | 6 ++---- components/omega/test/ocn/TendenciesTest.cpp | 3 ++- components/omega/test/ocn/TendencyTermsTest.cpp | 3 +-- 5 files changed, 9 insertions(+), 9 deletions(-) diff --git a/components/omega/test/ocn/AuxiliaryStateTest.cpp b/components/omega/test/ocn/AuxiliaryStateTest.cpp index f9bd3633bfae..d8afbb701c93 100644 --- a/components/omega/test/ocn/AuxiliaryStateTest.cpp +++ b/components/omega/test/ocn/AuxiliaryStateTest.cpp @@ -26,7 +26,8 @@ using namespace OMEGA; struct TestSetup { - Real Radius = REarth; + // radius of spherical mesh + Real Radius = HorzMesh::getDefault()->SphereRadius; KOKKOS_FUNCTION Real pseudoThickness(Real Lon, Real Lat) const { return (2 + std::cos(Lon) * std::pow(std::cos(Lat), 4)); diff --git a/components/omega/test/ocn/AuxiliaryVarsTest.cpp b/components/omega/test/ocn/AuxiliaryVarsTest.cpp index 4f916d998c7d..98340ceace4f 100644 --- a/components/omega/test/ocn/AuxiliaryVarsTest.cpp +++ b/components/omega/test/ocn/AuxiliaryVarsTest.cpp @@ -146,7 +146,8 @@ struct TestSetupPlane { struct TestSetupSphere { - Real Radius = REarth; + // radius of spherical mesh + Real Radius = HorzMesh::getDefault()->SphereRadius; ErrorMeasures ExpectedKineticEnergyErrors = {0.0143579382532765844, 0.00681096618897046764}; diff --git a/components/omega/test/ocn/HorzOperatorsTest.cpp b/components/omega/test/ocn/HorzOperatorsTest.cpp index 347747359c24..9a43ac780e6c 100644 --- a/components/omega/test/ocn/HorzOperatorsTest.cpp +++ b/components/omega/test/ocn/HorzOperatorsTest.cpp @@ -80,8 +80,7 @@ struct TestSetupPlane { struct TestSetupSphere1 { // radius of spherical mesh - // TODO: get this from the mesh - Real Radius = REarth; + Real Radius = HorzMesh::getDefault()->SphereRadius; ErrorMeasures ExpectedDivErrors = {0.013652414501664885, 0.0036904315983599676}; @@ -132,8 +131,7 @@ struct TestSetupSphere1 { struct TestSetupSphere2 { // radius of spherical mesh - // TODO: get this from the mesh - Real Radius = REarth; + Real Radius = HorzMesh::getDefault()->SphereRadius; ErrorMeasures ExpectedDivErrors = {1.37734693033362766e-10, 0.000484370621558727582}; diff --git a/components/omega/test/ocn/TendenciesTest.cpp b/components/omega/test/ocn/TendenciesTest.cpp index cf53a858d2eb..1451b12301a4 100644 --- a/components/omega/test/ocn/TendenciesTest.cpp +++ b/components/omega/test/ocn/TendenciesTest.cpp @@ -31,7 +31,8 @@ using namespace OMEGA; struct TestSetup { - Real Radius = REarth; + // radius of spherical mesh + Real Radius = HorzMesh::getDefault()->SphereRadius; KOKKOS_FUNCTION Real pseudoThickness(Real Lon, Real Lat) const { return (2 + std::cos(Lon) * std::pow(std::cos(Lat), 4)); diff --git a/components/omega/test/ocn/TendencyTermsTest.cpp b/components/omega/test/ocn/TendencyTermsTest.cpp index aa4698b87b5f..5a359fb8b675 100644 --- a/components/omega/test/ocn/TendencyTermsTest.cpp +++ b/components/omega/test/ocn/TendencyTermsTest.cpp @@ -203,8 +203,7 @@ struct TestSetupPlane { struct TestSetupSphere { // radius of spherical mesh - // TODO: get this from the mesh - Real Radius = REarth; + Real Radius = HorzMesh::getDefault()->SphereRadius; ErrorMeasures ExpectedDivErrors = {0.013652414501664885, 0.0036904315983599676}; From 4c3084cddf0f0b0e7eb538b5f61e291eea52dbf4 Mon Sep 17 00:00:00 2001 From: Xylar Asay-Davis Date: Fri, 11 Sep 2026 16:51:00 -0500 Subject: [PATCH 03/12] Check the mesh sphere radius against coordinates and REarth HorzMeshTest verified that cell centers are equidistant from the origin but never compared that distance with the SphereRadius attribute read from the file, nor with the REarth constant the model assumes. Add both checks, with a tolerance that follows the build precision since the coordinates are stored as Real. Co-Authored-By: Claude Opus 5 (1M context) --- components/omega/test/ocn/HorzMeshTest.cpp | 19 +++++++++++++++++++ 1 file changed, 19 insertions(+) diff --git a/components/omega/test/ocn/HorzMeshTest.cpp b/components/omega/test/ocn/HorzMeshTest.cpp index 101048ee74c8..d88c56cbfec8 100644 --- a/components/omega/test/ocn/HorzMeshTest.cpp +++ b/components/omega/test/ocn/HorzMeshTest.cpp @@ -15,6 +15,7 @@ #include "Dimension.h" #include "Error.h" #include "Field.h" +#include "GlobalConstants.h" #include "Halo.h" #include "IO.h" #include "IOStream.h" @@ -191,6 +192,24 @@ int main(int argc, char *argv[]) { if (Count > 0) ABORT_ERROR("HorzMeshTest: Cell sphere radius test FAIL"); + // Test that the sphere radius read from the mesh metadata agrees with + // the radius implied by the cell coordinates and with the Earth radius + // from the Physical Constants Dictionary. The mesh reader aborts when + // the metadata disagrees with REarth, so this mostly confirms that the + // coordinates were generated on that same sphere. Coordinates are + // stored as Real, so the tolerance depends on the build precision. + R8 RadiusRelTol = sizeof(Real) == 4 ? 1.e-6 : 1.e-10; + if (!Mesh->OnSphere) + ABORT_ERROR("HorzMeshTest: OnSphere flag test FAIL"); + if (abs((Mesh->SphereRadius - SphereRadius) / SphereRadius) > + RadiusRelTol) + ABORT_ERROR("HorzMeshTest: SphereRadius vs coordinates test FAIL " + "{} {}", + Mesh->SphereRadius, SphereRadius); + if (abs((Mesh->SphereRadius - REarth) / REarth) > RadiusRelTol) + ABORT_ERROR("HorzMeshTest: SphereRadius vs REarth test FAIL {} {}", + Mesh->SphereRadius, REarth); + // Test lon/lat coordinates of cell centers // Convert Cartesian coordinates to lon/lat and check these agree with the // values that have been read in From 97d8e475db5fa6e43c02fea487564eafa107cea1 Mon Sep 17 00:00:00 2001 From: Xylar Asay-Davis Date: Fri, 11 Sep 2026 16:51:25 -0500 Subject: [PATCH 04/12] Tighten the mesh sphere radius check to 1e-8 The relative tolerance on the input mesh's sphere radius was 1e-4 so that the CTest icosahedral mesh, generated with an older Earth radius, would still load. That mesh has been regenerated on the sphere of the Physical Constants Dictionary, so require agreement to 1e-8. Compare the file value with pcd::mean_radius in double precision, before the cast to Real, so that a single-precision build applies the same test. Closes #346 Co-Authored-By: Claude Opus 5 (1M context) --- components/omega/src/ocn/HorzMesh.cpp | 18 ++++++++++++------ 1 file changed, 12 insertions(+), 6 deletions(-) diff --git a/components/omega/src/ocn/HorzMesh.cpp b/components/omega/src/ocn/HorzMesh.cpp index 78f041584e1d..703be96cbd87 100644 --- a/components/omega/src/ocn/HorzMesh.cpp +++ b/components/omega/src/ocn/HorzMesh.cpp @@ -186,12 +186,18 @@ HorzMesh::HorzMesh(const std::string &Name, //< [in] Name for new mesh } else { ABORT_ERROR("Mesh is on sphere but sphere radius either missing or 0"); } - // This tolerance should be tightened once we have appropriate test - // mesh files that were generated with the same Earth radius - if (std::abs((REarth - SphereRadius) / REarth) > 1.e-4) - ABORT_ERROR("Input mesh has inaccurate earth radius: " - "REarth = {} input SphereRadius = {}", - REarth, SphereRadius); + // The mesh must have been generated on the same sphere as the Earth + // radius from the Physical Constants Dictionary. Compare in double + // precision (before the cast to Real) so that the tolerance does not + // depend on the build precision; the tolerance leaves room for a mesh + // file that stores the radius with fewer significant digits. + constexpr R8 RadiusRelTol = 1.e-8; + if (std::abs((pcd::mean_radius - TmpRadius) / pcd::mean_radius) > + RadiusRelTol) + ABORT_ERROR("Input mesh sphere radius {} does not match the Earth " + "radius {} from the Physical Constants Dictionary " + "(relative tolerance {})", + TmpRadius, pcd::mean_radius, RadiusRelTol); IsPeriodic = false; XPeriod = 0.0; YPeriod = 0.0; From 150e5faf55c6f809f5ab4161577ce35c40dd66b1 Mon Sep 17 00:00:00 2001 From: Xylar Asay-Davis Date: Fri, 11 Sep 2026 16:52:12 -0500 Subject: [PATCH 05/12] Make the spherical tracer diffusion test field smooth at the poles TestSetupSphere::scalarB was 2 + cos(Lon) sin(Lat), which has no limit at the poles: its value there depends on the direction of approach. The tracer diffusion test differentiates it, so the exact tendency at a polar cell depended on the arbitrary longitude assigned to that cell's center, through a term 2 cos^3(Lon)/R. On the old CTest mesh the polar pentagons happened to sit near Lon = 58 and 257 degrees, where the term is small (LInf 0.049); on the regenerated mesh they sit near 20 and 191 degrees and the LInf error rose to 0.235, all of it in those two cells. Use 2 + cos(Lon) cos(Lat) sin(Lat) = 2 + x z / R^2 instead, which is smooth on the sphere, and update the exact divergence accordingly. The error is then 0.0156 on either mesh, with its maximum at 30 degrees latitude as ordinary truncation error. The other uses of scalarB on the sphere are pointwise at edges, away from the poles, and unaffected. Co-Authored-By: Claude Opus 5 (1M context) --- .../omega/test/ocn/TendencyTermsTest.cpp | 26 ++++++++++++++----- 1 file changed, 19 insertions(+), 7 deletions(-) diff --git a/components/omega/test/ocn/TendencyTermsTest.cpp b/components/omega/test/ocn/TendencyTermsTest.cpp index 5a359fb8b675..7b3efe716535 100644 --- a/components/omega/test/ocn/TendencyTermsTest.cpp +++ b/components/omega/test/ocn/TendencyTermsTest.cpp @@ -315,17 +315,29 @@ struct TestSetupSphere { return Radius * std::pow(std::sin(Lon), 2) * std::pow(std::cos(Lat), 2); } + // scalarB is differentiated by the tracer diffusion test, so it must be + // smooth on the sphere: cos(Lon)*cos(Lat)*sin(Lat) = x*z/R^2 vanishes at + // the poles, whereas cos(Lon)*sin(Lat) has no limit there and makes the + // exact diffusion at a polar cell depend on the arbitrary longitude of + // that cell's center. KOKKOS_FUNCTION Real scalarB(Real Lon, Real Lat) const { - return 2. + std::cos(Lon) * std::sin(Lat); + return 2. + std::cos(Lon) * std::cos(Lat) * std::sin(Lat); } + // div(scalarB grad(scalarA)) on the sphere KOKKOS_FUNCTION Real tracerDiff(Real Lon, Real Lat, Real EddyDiff2) const { - return EddyDiff2 * - (4 * std::pow(std::cos(Lon), 2) - - 2 * (1. + 3 * std::cos(2 * Lat)) * std::pow(std::sin(Lon), 2) + - 2 * std::pow(std::cos(Lon), 3) * std::sin(Lat) - - 8 * std::cos(Lon) * std::pow(std::cos(Lat), 2) * - std::pow(std::sin(Lon), 2) * std::sin(Lat)) / + const Real CosLon = std::cos(Lon); + const Real SinLon = std::sin(Lon); + const Real CosLat = std::cos(Lat); + const Real SinLat = std::sin(Lat); + return 2 * EddyDiff2 * + (2 * std::cos(2 * Lon) + + CosLat * SinLat * CosLon * + (CosLon * CosLon - 2 * SinLon * SinLon) - + SinLon * SinLon * + (2 * (CosLat * CosLat - 2 * SinLat * SinLat) + + CosLon * CosLat * SinLat * + (2 * CosLat * CosLat - 3 * SinLat * SinLat))) / Radius; } From d09c546e47d7167c3f82f482dee75c2754026a33 Mon Sep 17 00:00:00 2001 From: Xylar Asay-Davis Date: Fri, 11 Sep 2026 16:52:35 -0500 Subject: [PATCH 06/12] Update expected errors for the regenerated spherical CTest mesh The icosahedral CTest mesh, cosine_bell_icos480.omega_vars.260911.nc, was regenerated by Polaris on the sphere with the Physical Constants Dictionary radius. The mesh has the same connectivity as before but JIGSAW places the cells up to 7e-4 (relative) differently, and the analytic test fields are now evaluated with the mesh radius, so the error norms of the spherical operator, auxiliary variable, tendency and forcing tests move by up to a few tenths of a percent (2 percent for those with a radius-scaled analytic solution). Record the new values. Co-Authored-By: Claude Opus 5 (1M context) --- .../omega/test/ocn/AuxiliaryVarsTest.cpp | 60 +++++++++---------- components/omega/test/ocn/ForcingTest.cpp | 4 +- .../omega/test/ocn/HorzOperatorsTest.cpp | 28 ++++----- .../omega/test/ocn/TendencyTermsTest.cpp | 40 ++++++------- 4 files changed, 66 insertions(+), 66 deletions(-) diff --git a/components/omega/test/ocn/AuxiliaryVarsTest.cpp b/components/omega/test/ocn/AuxiliaryVarsTest.cpp index 98340ceace4f..0f34953a99cc 100644 --- a/components/omega/test/ocn/AuxiliaryVarsTest.cpp +++ b/components/omega/test/ocn/AuxiliaryVarsTest.cpp @@ -149,39 +149,39 @@ struct TestSetupSphere { // radius of spherical mesh Real Radius = HorzMesh::getDefault()->SphereRadius; - ErrorMeasures ExpectedKineticEnergyErrors = {0.0143579382532765844, - 0.00681096618897046764}; - ErrorMeasures ExpectedVelocityDivErrors = {0.013652414501663887, - 0.00369043159835992}; - - ErrorMeasures ExpectedFluxThickErrors = {0.0159821090867812224, - 0.010364511516135164}; - ErrorMeasures ExpectedMeanThickErrors = {0.000800109287518277435, - 0.000406527457820634436}; - - ErrorMeasures ExpectedRelVortVertexErrors = {0.0271404735181343393, - 0.0252023166109219786}; - ErrorMeasures ExpectedNormRelVortVertexErrors = {0.0348741350737879693, - 0.0259506101504540822}; - ErrorMeasures ExpectedNormPlanetVortVertexErrors = {0.00451268952953497778, - 0.00101771171197261793}; - - ErrorMeasures ExpectedNormRelVortEdgeErrors = {0.01255832726488633, - 0.003093726580373857}; - ErrorMeasures ExpectedNormPlanetVortEdgeErrors = {0.00495174534686814403, - 0.000855432390947949515}; - - ErrorMeasures ExpectedDel2Errors = {0.0036327578533323366, - 0.003161584491352198}; - ErrorMeasures ExpectedDel2DivErrors = {0.0177782108439020134, - 0.0075764886471247385}; - ErrorMeasures ExpectedDel2RelVortErrors = {0.09148975765556343, - 0.0246736311927726465}; + ErrorMeasures ExpectedKineticEnergyErrors = {0.014358161778181201, + 0.0068116910597281113}; + ErrorMeasures ExpectedVelocityDivErrors = {0.013659556526123424, + 0.0036698023569595927}; + + ErrorMeasures ExpectedFluxThickErrors = {0.015981055597025, + 0.010364584620791126}; + ErrorMeasures ExpectedMeanThickErrors = {0.00079985423942030413, + 0.00040654771018863028}; + + ErrorMeasures ExpectedRelVortVertexErrors = {0.027146182264485331, + 0.025202560969432641}; + ErrorMeasures ExpectedNormRelVortVertexErrors = {0.034882855311945148, + 0.025950288113940474}; + ErrorMeasures ExpectedNormPlanetVortVertexErrors = {0.0045134270290698875, + 0.0010179210697515548}; + + ErrorMeasures ExpectedNormRelVortEdgeErrors = {0.012528744550304859, + 0.0030749833896684914}; + ErrorMeasures ExpectedNormPlanetVortEdgeErrors = {0.0049522329721762544, + 0.0008554077431009539}; + + ErrorMeasures ExpectedDel2Errors = {0.0036017312925143057, + 0.003134053694265801}; + ErrorMeasures ExpectedDel2DivErrors = {0.017767406015402661, + 0.007518774954098988}; + ErrorMeasures ExpectedDel2RelVortErrors = {0.091446507985509698, + 0.024667998533570906}; ErrorMeasures ExpectedHTracerErrors = {0.01603249913425972, 0.00546762028673672059}; - ErrorMeasures ExpectedDel2TracerErrors = {0.0081206665417422382, - 0.004969575978774801}; + ErrorMeasures ExpectedDel2TracerErrors = {0.0081209686316016313, + 0.0049181207204490136}; KOKKOS_FUNCTION Real pseudoThickness(Real Lon, Real Lat) const { return (2 + std::cos(Lon) * std::pow(std::cos(Lat), 4)); diff --git a/components/omega/test/ocn/ForcingTest.cpp b/components/omega/test/ocn/ForcingTest.cpp index 2ac027553f05..15afb6eddd1a 100644 --- a/components/omega/test/ocn/ForcingTest.cpp +++ b/components/omega/test/ocn/ForcingTest.cpp @@ -49,8 +49,8 @@ struct TestSetupPlane { }; struct TestSetupSphere { - ErrorMeasures ExpectedNormalStressErrors = {0.0038588958862868362, - 0.003813760171030077}; + ErrorMeasures ExpectedNormalStressErrors = {0.003704100692552216, + 0.0037705223292748575}; KOKKOS_FUNCTION Real sfcStressX(Real Lon, Real Lat) const { return -4 * std::sin(Lon) * std::cos(Lon) * std::pow(std::cos(Lat), 3) * std::sin(Lat); diff --git a/components/omega/test/ocn/HorzOperatorsTest.cpp b/components/omega/test/ocn/HorzOperatorsTest.cpp index 9a43ac780e6c..c29dff8291b1 100644 --- a/components/omega/test/ocn/HorzOperatorsTest.cpp +++ b/components/omega/test/ocn/HorzOperatorsTest.cpp @@ -82,20 +82,20 @@ struct TestSetupSphere1 { // radius of spherical mesh Real Radius = HorzMesh::getDefault()->SphereRadius; - ErrorMeasures ExpectedDivErrors = {0.013652414501664885, - 0.0036904315983599676}; - ErrorMeasures ExpectedGradErrors = {0.0019094381714837498, - 0.0015218320661105687}; - ErrorMeasures ExpectedCurlErrors = {0.02713957370128636, - 0.025202095212756463}; - ErrorMeasures ExpectedReconErrors = {0.0206375134079833517, - 0.00692590524910695858}; - ErrorMeasures ExpectedVectorReconErrors = {0.003308078126087204, - 0.0029425911574823007}; - ErrorMeasures ExpectedAnisoInterpErrors = {0.0024015775047603197, - 0.0018490649516209202}; - ErrorMeasures ExpectedIsoInterpErrors = {0.007438367234983312, - 0.0029921955942401697}; + ErrorMeasures ExpectedDivErrors = {0.013659556526126423, + 0.0036698023569596279}; + ErrorMeasures ExpectedGradErrors = {0.0018790366180557886, + 0.0014984647360648147}; + ErrorMeasures ExpectedCurlErrors = {0.027146182264483853, + 0.025202560969432571}; + ErrorMeasures ExpectedReconErrors = {0.020633072221392729, + 0.0069262670519179576}; + ErrorMeasures ExpectedVectorReconErrors = {0.0033076912646329529, + 0.0029392950172411601}; + ErrorMeasures ExpectedAnisoInterpErrors = {0.0024007805796189899, + 0.0018491570972366891}; + ErrorMeasures ExpectedIsoInterpErrors = {0.0074341565605974018, + 0.0029922776661517052}; KOKKOS_FUNCTION Real exactScalar(Real Lon, Real Lat) const { return Radius * std::cos(Lon) * std::pow(std::cos(Lat), 4); diff --git a/components/omega/test/ocn/TendencyTermsTest.cpp b/components/omega/test/ocn/TendencyTermsTest.cpp index 7b3efe716535..525540afb25a 100644 --- a/components/omega/test/ocn/TendencyTermsTest.cpp +++ b/components/omega/test/ocn/TendencyTermsTest.cpp @@ -205,16 +205,16 @@ struct TestSetupSphere { // radius of spherical mesh Real Radius = HorzMesh::getDefault()->SphereRadius; - ErrorMeasures ExpectedDivErrors = {0.013652414501664885, - 0.0036904315983599676}; - ErrorMeasures ExpectedPVErrors = {0.0219217796608757037, - 0.0122537418367830303}; - ErrorMeasures ExpectedGradErrors = {0.0019094381714837498, - 0.0015218320661105702}; - ErrorMeasures ExpectedLaplaceErrors = {0.28193638497826856, - 0.270546491554748}; - ErrorMeasures ExpectedTrHAdvErrors = {0.013259410329645643, - 0.004094907022292395}; + ErrorMeasures ExpectedDivErrors = {0.013659556526126423, + 0.0036698023569596279}; + ErrorMeasures ExpectedPVErrors = {0.021937128290005589, + 0.012253081309866708}; + ErrorMeasures ExpectedGradErrors = {0.0018790366180557886, + 0.0014984647360648147}; + ErrorMeasures ExpectedLaplaceErrors = {0.2819223346929331, + 0.27053032782204006}; + ErrorMeasures ExpectedTrHAdvErrors = {0.013259084476801913, + 0.004075236217689366}; ErrorMeasures ExpectedFCTErrors = {0.00, 0.00}; ErrorMeasures ExpectedFCTHProv = {3.0542657508680904e-05, 1.0779233406323438e-06}; @@ -226,18 +226,18 @@ struct TestSetupSphere { ErrorMeasures ExpectedFCT_Low = {0.2958192441638326, 0.09980000200239414}; ErrorMeasures ExpectedFCTToNonFCT = {1.5120472343676013e-15, 6.215666716489436e-16}; - ErrorMeasures ExpectedTrDel2Errors = {0.04865718541236144, - 0.005105510870642706}; - ErrorMeasures ExpectedTrDel4Errors = {0.0008646345116716073, - 0.0007118574326665881}; + ErrorMeasures ExpectedTrDel2Errors = {0.015620987792697782, + 0.0032365878896264888}; + ErrorMeasures ExpectedTrDel4Errors = {0.00081985237645818541, + 0.00064699727012728107}; ErrorMeasures ExpectedSurfTrRestErrors = {0, 0}; ErrorMeasures ExpectedSfcStressForcingErrors = {0, 0}; - ErrorMeasures ExpectedBottomDragErrors = {0.0015333449035655053, - 0.0014897009917655022}; - ErrorMeasures ExpectedCoriolis2DErrors = {0.017830942909137566, - 0.00958613271059214}; - ErrorMeasures ExpectedCoriolis1DErrors = {0.01756710962800044, - 0.011526527317437694}; + ErrorMeasures ExpectedBottomDragErrors = {0.0015343843060931499, + 0.0014861805971789941}; + ErrorMeasures ExpectedCoriolis2DErrors = {0.017808748900211516, + 0.0095867629560083354}; + ErrorMeasures ExpectedCoriolis1DErrors = {0.017532349880532885, + 0.011526365807152578}; KOKKOS_FUNCTION Real velEdgeX(Real X, Real Y) const { return 2; } From 24c03492ba0161752b9f284db157c64f761e7728 Mon Sep 17 00:00:00 2001 From: Xylar Asay-Davis Date: Fri, 11 Sep 2026 16:53:15 -0500 Subject: [PATCH 07/12] Document the sphere radius requirement and update CTest mesh links The MPAS mesh page said there was no requirement on the sphere radius, which has not been true since the model started checking it. State that spherical meshes must use the Physical Constants Dictionary Earth radius, and point to the Polaris converter for older MPAS-Ocean meshes. Mention the check on the horizontal mesh page, and update the quick start links to the regenerated icosahedral mesh and the current planar mesh. Co-Authored-By: Claude Opus 5 (1M context) --- components/omega/doc/devGuide/QuickStart.md | 12 +++++++----- components/omega/doc/userGuide/HorzMesh.md | 6 +++++- components/omega/doc/userGuide/MpasMesh.md | 13 ++++++++++--- 3 files changed, 22 insertions(+), 9 deletions(-) diff --git a/components/omega/doc/devGuide/QuickStart.md b/components/omega/doc/devGuide/QuickStart.md index 38dd8e2475af..ef6166fa96b5 100644 --- a/components/omega/doc/devGuide/QuickStart.md +++ b/components/omega/doc/devGuide/QuickStart.md @@ -122,16 +122,18 @@ At a high level: Some tests require a valid Omega mesh file. Different tests require different meshes. At the moment, mesh files need to be copied or linked to specifically -named files under the `test` directory. Appropriate mesh files can be +named files under the `test` directory. Spherical meshes must lie on a sphere +with the Earth radius from the Physical Constants Dictionary (see +[Horizontal Mesh](#omega-user-horz-mesh)). Appropriate mesh files can be downloaded from: - [Ocean Mesh](https://web.lcrc.anl.gov/public/e3sm/polaris/ocean/omega_ctest/ocean.QU.240km.omega_vars.260807.nc) -- [Global Mesh](https://web.lcrc.anl.gov/public/e3sm/polaris/ocean/omega_ctest/cosine_bell_icos480.omega_vars.260807.nc) -- [Planar Mesh](https://web.lcrc.anl.gov/public/e3sm/polaris/ocean/omega_ctest/PlanarPeriodic48x48.omega_vars.260720.nc) +- [Global Mesh](https://web.lcrc.anl.gov/public/e3sm/polaris/ocean/omega_ctest/cosine_bell_icos480.omega_vars.260911.nc) +- [Planar Mesh](https://web.lcrc.anl.gov/public/e3sm/polaris/ocean/omega_ctest/PlanarPeriodic48x48.omega_vars.260825.nc) ```sh cd test wget -O OmegaMesh.nc https://web.lcrc.anl.gov/public/e3sm/polaris/ocean/omega_ctest/ocean.QU.240km.omega_vars.260807.nc -wget -O OmegaSphereMesh.nc https://web.lcrc.anl.gov/public/e3sm/polaris/ocean/omega_ctest/cosine_bell_icos480.omega_vars.260807.nc -wget -O OmegaPlanarMesh.nc https://web.lcrc.anl.gov/public/e3sm/polaris/ocean/omega_ctest/PlanarPeriodic48x48.omega_vars.260720.nc +wget -O OmegaSphereMesh.nc https://web.lcrc.anl.gov/public/e3sm/polaris/ocean/omega_ctest/cosine_bell_icos480.omega_vars.260911.nc +wget -O OmegaPlanarMesh.nc https://web.lcrc.anl.gov/public/e3sm/polaris/ocean/omega_ctest/PlanarPeriodic48x48.omega_vars.260825.nc cd .. ``` diff --git a/components/omega/doc/userGuide/HorzMesh.md b/components/omega/doc/userGuide/HorzMesh.md index 4cd0ddc26405..8ffb2c482728 100644 --- a/components/omega/doc/userGuide/HorzMesh.md +++ b/components/omega/doc/userGuide/HorzMesh.md @@ -43,7 +43,11 @@ variables. These include ``OnSphere`` (or ``on_a_sphere`` for backcompatibility) ``IsPeriodic`` (``is_periodic``), ``SphereRadius`` (``sphere_radius``), ``XPeriod`` (``x_period``), and ``YPeriod`` (``y_period``). The two flags OnSphere and IsPeriodic are stored as YES/NO strings in the metadata but as -boolean flags in the code. +boolean flags in the code. For a spherical mesh, ``SphereRadius`` must match +the Earth radius from the Physical Constants Dictionary (``REarth`` in +``GlobalConstants.h``) to a relative tolerance of $10^{-8}$; the model aborts +otherwise. See the {ref}`MPAS mesh specification ` for how to +bring older meshes into compliance. In the future, the Mesh class will optionally compute internally those fields the specification marks as *could be computed internally* — the areas, lengths, diff --git a/components/omega/doc/userGuide/MpasMesh.md b/components/omega/doc/userGuide/MpasMesh.md index 95ef34023e38..5d8654efded6 100644 --- a/components/omega/doc/userGuide/MpasMesh.md +++ b/components/omega/doc/userGuide/MpasMesh.md @@ -43,8 +43,15 @@ geometry; cells, edges, and vertices are constrained to lie on the surface of the sphere, and the sphere is centered at the origin, $(x,y,z) = (0,0,0)$. In all cases, the coordinate systems assumed by MPAS meshes are right-handed. {numref}`fig-voronoi-examples` illustrates VT meshes on the Cartesian plane and -on the sphere. There are no requirements for the radius of the sphere, however -an attribute `sphere_radius` is required that defines the radius of the sphere. +on the sphere. An attribute `sphere_radius` is required that defines the +radius of the sphere. The MPAS mesh specification places no requirement on +that radius, but Omega does: a spherical mesh must lie on a sphere with the +Earth radius from the Physical Constants Dictionary (`mean_radius`, +6371008.7714 m), and Omega aborts if `sphere_radius` differs from it by more +than a relative tolerance of $10^{-8}$. Meshes generated by Polaris satisfy +this; older MPAS-Ocean meshes (typically 6371229 m) must be rescaled, which +Polaris' `utils/omega/convert_mpaso_ic_to_omega.py` does along with the +variable renaming. ```{figure} images/meshPlanarSpherical.png :name: fig-voronoi-examples @@ -385,7 +392,7 @@ meshes should contain these global attributes. | Attribute | Type | Valid values | Description | |-----------------|-----------|-----------------------------------------------|-------------| | `on_a_sphere` | character | `"YES"` or `"NO"` | Defines whether the mesh describes points that lie on the surface of a sphere. | -| `sphere_radius` | double | any non-negative real value | Radius of the sphere the points are defined on. Set to 0 for a planar mesh. | +| `sphere_radius` | double | any non-negative real value | Radius of the sphere the points are defined on. Set to 0 for a planar mesh. Omega requires the PCD Earth radius (6371008.7714 m) for spherical meshes. | | `is_periodic` | character | `"YES"` or `"NO"` | Defines whether the mesh has any periodic boundaries. Only meaningful when `on_a_sphere` is `"NO"`. | | `mesh_spec` | character | any valid version of the MPAS mesh specification | Version of the MPAS mesh specification the mesh conforms to. | | `file_id` | character | any combination of lowercase letters and numbers | Random string used for tracking mesh provenance. | From 8e031f06308b32f82b3a44485c1cdbf4b37f62d7 Mon Sep 17 00:00:00 2001 From: Xylar Asay-Davis Date: Tue, 15 Sep 2026 03:53:27 -0500 Subject: [PATCH 08/12] Use REarth rather than the mesh radius in operators and tests With the input mesh now required to lie on the sphere of radius REarth to 1e-8, carrying HorzMesh::SphereRadius into SecondDerivativeOnCell, the prescribed velocities and the test setups adds plumbing without adding information. Use the REarth constant directly, as suggested in review. Co-Authored-By: Claude Opus 5 (1M context) --- components/omega/src/ocn/HorzOperators.cpp | 13 ++++++------- components/omega/src/ocn/HorzOperators.h | 15 +++++++-------- .../omega/src/timeStepping/TimeStepper.cpp | 10 ++++------ .../omega/test/ocn/AuxiliaryStateTest.cpp | 2 +- .../omega/test/ocn/AuxiliaryVarsTest.cpp | 2 +- .../omega/test/ocn/HorzOperatorsTest.cpp | 19 ++++++++----------- components/omega/test/ocn/TendenciesTest.cpp | 2 +- .../omega/test/ocn/TendencyTermsTest.cpp | 2 +- 8 files changed, 29 insertions(+), 36 deletions(-) diff --git a/components/omega/src/ocn/HorzOperators.cpp b/components/omega/src/ocn/HorzOperators.cpp index 4a01cd320c8d..a530e8254dd1 100644 --- a/components/omega/src/ocn/HorzOperators.cpp +++ b/components/omega/src/ocn/HorzOperators.cpp @@ -30,13 +30,12 @@ InterpCellToEdge::InterpCellToEdge(const HorzMesh *Mesh) VertexDegree(Mesh->VertexDegree) {} SecondDerivativeOnCell::SecondDerivativeOnCell(HorzMesh const *Mesh) - : OnSphere(Mesh->OnSphere), SphereRadius(Mesh->SphereRadius), - NCellsAll(Mesh->NCellsAll), MaxEdges(1 + Mesh->MaxEdges), - NEdgesOnCell(Mesh->NEdgesOnCell), EdgesOnCell(Mesh->EdgesOnCell), - CellsOnCell(Mesh->CellsOnCell), CellsOnEdge(Mesh->CellsOnEdge), - VerticesOnEdge(Mesh->VerticesOnEdge), XCell(Mesh->XCell), - YCell(Mesh->YCell), ZCell(Mesh->ZCell), DvEdge(Mesh->DvEdge), - DcEdge(Mesh->DcEdge), AngleEdge(Mesh->AngleEdge), + : OnSphere(Mesh->OnSphere), NCellsAll(Mesh->NCellsAll), + MaxEdges(1 + Mesh->MaxEdges), NEdgesOnCell(Mesh->NEdgesOnCell), + EdgesOnCell(Mesh->EdgesOnCell), CellsOnCell(Mesh->CellsOnCell), + CellsOnEdge(Mesh->CellsOnEdge), VerticesOnEdge(Mesh->VerticesOnEdge), + XCell(Mesh->XCell), YCell(Mesh->YCell), ZCell(Mesh->ZCell), + DvEdge(Mesh->DvEdge), DcEdge(Mesh->DcEdge), AngleEdge(Mesh->AngleEdge), AreaCell(Mesh->AreaCell), EdgeSignOnCell(Mesh->EdgeSignOnCell), XVertex(Mesh->XVertex), YVertex(Mesh->YVertex), ZVertex(Mesh->ZVertex), ThetaAbs("SphereAngle", NCellsAll), XPCell("XP", NCellsAll, MaxEdges), diff --git a/components/omega/src/ocn/HorzOperators.h b/components/omega/src/ocn/HorzOperators.h index 2bd395cd774b..a78eba11fe5d 100644 --- a/components/omega/src/ocn/HorzOperators.h +++ b/components/omega/src/ocn/HorzOperators.h @@ -235,10 +235,9 @@ class SecondDerivativeOnCell { if (OnSphere) { const Array1DI4 edgesOnCell = Kokkos::subview(EdgesOnCell, ICell, Kokkos::ALL); - DetermineSphericalPatchGeometry(NEdges, edgesOnCell, VerticesOnEdge, - XCell, YCell, ZCell, XVertex, YVertex, - ZVertex, CellList, SphereRadius, XP, - YP, Angle2D, ThetaAbs[ICell]); + DetermineSphericalPatchGeometry( + NEdges, edgesOnCell, VerticesOnEdge, XCell, YCell, ZCell, XVertex, + YVertex, ZVertex, CellList, XP, YP, Angle2D, ThetaAbs[ICell]); } else { // On an x-y plane const Array1DI4 edgesOnCell = Kokkos::subview(EdgesOnCell, ICell, Kokkos::ALL); @@ -271,7 +270,6 @@ class SecondDerivativeOnCell { static constexpr R8 Pii = 3.141592653589793_Real; const bool OnSphere; - const Real SphereRadius; ///< radius (m) of the sphere the mesh lies on const I4 NCellsAll; const I4 MaxEdges; Array1DI4 NEdgesOnCell; @@ -341,10 +339,11 @@ class SecondDerivativeOnCell { const Array2DI4 VerticesOnEdge, const Array1DReal XCell, const Array1DReal YCell, const Array1DReal ZCell, const Array1DReal XVertex, const Array1DReal YVertex, - const Array1DReal ZVertex, const Array1DI4 CellList, - const Real sphereRadius, Array1DReal XP, Array1DReal YP, - Array1DReal Angle2D, Real &ThetaAbs) { + const Array1DReal ZVertex, const Array1DI4 CellList, Array1DReal XP, + Array1DReal YP, Array1DReal Angle2D, Real &ThetaAbs) { const Real length_scale = 1._Real; + // the mesh is required to lie on the sphere of radius REarth + const Real sphereRadius = REarth; Real XC[MaxMaxEdges] = {}; Real YC[MaxMaxEdges] = {}; Real ZC[MaxMaxEdges] = {}; diff --git a/components/omega/src/timeStepping/TimeStepper.cpp b/components/omega/src/timeStepping/TimeStepper.cpp index 58e12908208a..a480ba7b561d 100644 --- a/components/omega/src/timeStepping/TimeStepper.cpp +++ b/components/omega/src/timeStepping/TimeStepper.cpp @@ -639,7 +639,6 @@ void TimeStepper::prescribeVelocity(OceanState *State1, int TimeLevel1, OMEGA_SCOPE(LatEdge, Mesh->LatEdge); OMEGA_SCOPE(LonEdge, Mesh->LonEdge); OMEGA_SCOPE(AngleEdge, Mesh->AngleEdge); - const Real SphereRadius = Mesh->SphereRadius; OMEGA_SCOPE(MinLayerEdgeBot, VCoord->MinLayerEdgeBot); OMEGA_SCOPE(MaxLayerEdgeTop, VCoord->MaxLayerEdgeTop); @@ -665,8 +664,8 @@ void TimeStepper::prescribeVelocity(OceanState *State1, int TimeLevel1, 2.0 * Pi * cos(LatEdge(IEdge))); const R8 v = (10.0 / Tau) * sin(2.0 * lon_p) * cos(LatEdge(IEdge)) * cos(Pi * TSim / Tau); - const R8 normalVel = SphereRadius * (u * cos(AngleEdge(IEdge)) + - v * sin(AngleEdge(IEdge))); + const R8 normalVel = REarth * (u * cos(AngleEdge(IEdge)) + + v * sin(AngleEdge(IEdge))); parallelForInner( Team, KRange, INNER_LAMBDA(int KChunk) { @@ -681,7 +680,6 @@ void TimeStepper::prescribeVelocity(OceanState *State1, int TimeLevel1, OMEGA_SCOPE(LatEdge, Mesh->LatEdge); OMEGA_SCOPE(LonEdge, Mesh->LonEdge); OMEGA_SCOPE(AngleEdge, Mesh->AngleEdge); - const Real SphereRadius = Mesh->SphereRadius; OMEGA_SCOPE(MinLayerEdgeBot, VCoord->MinLayerEdgeBot); OMEGA_SCOPE(MaxLayerEdgeTop, VCoord->MaxLayerEdgeTop); @@ -710,8 +708,8 @@ void TimeStepper::prescribeVelocity(OceanState *State1, int TimeLevel1, const R8 v = ((2.5 / Tau) * sin(lon_p) * Kokkos::pow(cos(LatEdge(IEdge)), 3) * cos(Pi * TSim / Tau)); - const R8 normalVel = SphereRadius * (u * cos(AngleEdge(IEdge)) + - v * sin(AngleEdge(IEdge))); + const R8 normalVel = REarth * (u * cos(AngleEdge(IEdge)) + + v * sin(AngleEdge(IEdge))); parallelForInner( Team, KRange, INNER_LAMBDA(int KChunk) { diff --git a/components/omega/test/ocn/AuxiliaryStateTest.cpp b/components/omega/test/ocn/AuxiliaryStateTest.cpp index d8afbb701c93..891f14aa9470 100644 --- a/components/omega/test/ocn/AuxiliaryStateTest.cpp +++ b/components/omega/test/ocn/AuxiliaryStateTest.cpp @@ -27,7 +27,7 @@ using namespace OMEGA; struct TestSetup { // radius of spherical mesh - Real Radius = HorzMesh::getDefault()->SphereRadius; + Real Radius = REarth; KOKKOS_FUNCTION Real pseudoThickness(Real Lon, Real Lat) const { return (2 + std::cos(Lon) * std::pow(std::cos(Lat), 4)); diff --git a/components/omega/test/ocn/AuxiliaryVarsTest.cpp b/components/omega/test/ocn/AuxiliaryVarsTest.cpp index 0f34953a99cc..260ae8cbbf87 100644 --- a/components/omega/test/ocn/AuxiliaryVarsTest.cpp +++ b/components/omega/test/ocn/AuxiliaryVarsTest.cpp @@ -147,7 +147,7 @@ struct TestSetupPlane { struct TestSetupSphere { // radius of spherical mesh - Real Radius = HorzMesh::getDefault()->SphereRadius; + Real Radius = REarth; ErrorMeasures ExpectedKineticEnergyErrors = {0.014358161778181201, 0.0068116910597281113}; diff --git a/components/omega/test/ocn/HorzOperatorsTest.cpp b/components/omega/test/ocn/HorzOperatorsTest.cpp index c29dff8291b1..bab6fc845b6c 100644 --- a/components/omega/test/ocn/HorzOperatorsTest.cpp +++ b/components/omega/test/ocn/HorzOperatorsTest.cpp @@ -80,7 +80,7 @@ struct TestSetupPlane { struct TestSetupSphere1 { // radius of spherical mesh - Real Radius = HorzMesh::getDefault()->SphereRadius; + Real Radius = REarth; ErrorMeasures ExpectedDivErrors = {0.013659556526126423, 0.0036698023569596279}; @@ -131,7 +131,7 @@ struct TestSetupSphere1 { struct TestSetupSphere2 { // radius of spherical mesh - Real Radius = HorzMesh::getDefault()->SphereRadius; + Real Radius = REarth; ErrorMeasures ExpectedDivErrors = {1.37734693033362766e-10, 0.000484370621558727582}; @@ -749,9 +749,8 @@ int testsecondderivativeoncellDetermineSphericalPatchGeometry(Real RTol) { {20, 21, 25, 29, 28, 24}, {23, 24, 28, 32, 31, 27}, {26, 27, 31, 35, 34, 30}, {28, 29, 33, 37, 36, 32}}; const int CtoC = 1039230; - // synthetic patch on its own sphere, independent of the test mesh - const double R = OMEGA::REarth; - const R8 Pii = 3.141592653589793_Real; + const double R = OMEGA::REarth; + const R8 Pii = 3.141592653589793_Real; // Project coordinates to sphere. @@ -898,13 +897,11 @@ int testsecondderivativeoncellDetermineSphericalPatchGeometry(Real RTol) { const Array2DI4 VerticiesOnEdge, const Array1DReal XCell, const Array1DReal YCell, const Array1DReal ZCell, const Array1DReal XVertex, const Array1DReal YVertex, - const Array1DReal ZVertex, const Array1DI4 CellList, - const Real SphereRadius, Array1DReal XP, Array1DReal YP, - Array1DReal Angle2D, Real &ThetaAbs) { + const Array1DReal ZVertex, const Array1DI4 CellList, Array1DReal XP, + Array1DReal YP, Array1DReal Angle2D, Real &ThetaAbs) { SecondDerivativeOnCell::DetermineSphericalPatchGeometry( NEdges, EdgesOnCell, VerticiesOnEdge, XCell, YCell, ZCell, XVertex, - YVertex, ZVertex, CellList, SphereRadius, XP, YP, Angle2D, - ThetaAbs); + YVertex, ZVertex, CellList, XP, YP, Angle2D, ThetaAbs); } }; { @@ -926,7 +923,7 @@ int testsecondderivativeoncellDetermineSphericalPatchGeometry(Real RTol) { Real ThetaAbs = {}; SecondDerivativeOnCellTest::Test( NEdges, edgesOnCell, VerticesOnEdge, XCell, YCell, ZCell, - XVertex, YVertex, ZVertex, CellList, R, XP, YP, Angle2D, + XVertex, YVertex, ZVertex, CellList, XP, YP, Angle2D, ThetaAbs); }); auto XPH = createHostMirrorCopy(XP); diff --git a/components/omega/test/ocn/TendenciesTest.cpp b/components/omega/test/ocn/TendenciesTest.cpp index 1451b12301a4..c20d9cf8a936 100644 --- a/components/omega/test/ocn/TendenciesTest.cpp +++ b/components/omega/test/ocn/TendenciesTest.cpp @@ -32,7 +32,7 @@ using namespace OMEGA; struct TestSetup { // radius of spherical mesh - Real Radius = HorzMesh::getDefault()->SphereRadius; + Real Radius = REarth; KOKKOS_FUNCTION Real pseudoThickness(Real Lon, Real Lat) const { return (2 + std::cos(Lon) * std::pow(std::cos(Lat), 4)); diff --git a/components/omega/test/ocn/TendencyTermsTest.cpp b/components/omega/test/ocn/TendencyTermsTest.cpp index 525540afb25a..0f2c00ce8e6e 100644 --- a/components/omega/test/ocn/TendencyTermsTest.cpp +++ b/components/omega/test/ocn/TendencyTermsTest.cpp @@ -203,7 +203,7 @@ struct TestSetupPlane { struct TestSetupSphere { // radius of spherical mesh - Real Radius = HorzMesh::getDefault()->SphereRadius; + Real Radius = REarth; ErrorMeasures ExpectedDivErrors = {0.013659556526126423, 0.0036698023569596279}; From d8fecb898c1a24304cfcbaed122037212f3d5513 Mon Sep 17 00:00:00 2001 From: Xylar Asay-Davis Date: Tue, 15 Sep 2026 03:54:30 -0500 Subject: [PATCH 09/12] State the smoothness requirement on scalarB without history Co-Authored-By: Claude Opus 5 (1M context) --- components/omega/test/ocn/TendencyTermsTest.cpp | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/components/omega/test/ocn/TendencyTermsTest.cpp b/components/omega/test/ocn/TendencyTermsTest.cpp index 0f2c00ce8e6e..13c5aae1fada 100644 --- a/components/omega/test/ocn/TendencyTermsTest.cpp +++ b/components/omega/test/ocn/TendencyTermsTest.cpp @@ -317,9 +317,7 @@ struct TestSetupSphere { // scalarB is differentiated by the tracer diffusion test, so it must be // smooth on the sphere: cos(Lon)*cos(Lat)*sin(Lat) = x*z/R^2 vanishes at - // the poles, whereas cos(Lon)*sin(Lat) has no limit there and makes the - // exact diffusion at a polar cell depend on the arbitrary longitude of - // that cell's center. + // the poles, where Lon is undefined. KOKKOS_FUNCTION Real scalarB(Real Lon, Real Lat) const { return 2. + std::cos(Lon) * std::cos(Lat) * std::sin(Lat); } From 630d17f93f88ead6a4dee53c36196b523a74848a Mon Sep 17 00:00:00 2001 From: Xylar Asay-Davis Date: Tue, 15 Sep 2026 03:54:59 -0500 Subject: [PATCH 10/12] Express the sphere radius check with OMEGA_REQUIRE Co-Authored-By: Claude Opus 5 (1M context) --- components/omega/src/ocn/HorzMesh.cpp | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/components/omega/src/ocn/HorzMesh.cpp b/components/omega/src/ocn/HorzMesh.cpp index 703be96cbd87..838287269cf3 100644 --- a/components/omega/src/ocn/HorzMesh.cpp +++ b/components/omega/src/ocn/HorzMesh.cpp @@ -192,12 +192,12 @@ HorzMesh::HorzMesh(const std::string &Name, //< [in] Name for new mesh // depend on the build precision; the tolerance leaves room for a mesh // file that stores the radius with fewer significant digits. constexpr R8 RadiusRelTol = 1.e-8; - if (std::abs((pcd::mean_radius - TmpRadius) / pcd::mean_radius) > - RadiusRelTol) - ABORT_ERROR("Input mesh sphere radius {} does not match the Earth " - "radius {} from the Physical Constants Dictionary " - "(relative tolerance {})", - TmpRadius, pcd::mean_radius, RadiusRelTol); + OMEGA_REQUIRE(std::abs((pcd::mean_radius - TmpRadius) / + pcd::mean_radius) <= RadiusRelTol, + "Input mesh sphere radius {} does not match the Earth " + "radius {} from the Physical Constants Dictionary " + "(relative tolerance {})", + TmpRadius, pcd::mean_radius, RadiusRelTol); IsPeriodic = false; XPeriod = 0.0; YPeriod = 0.0; From 676292f15eda222cc16d44d8a3e466f93a677c61 Mon Sep 17 00:00:00 2001 From: Xylar Asay-Davis Date: Wed, 23 Sep 2026 09:59:38 -0500 Subject: [PATCH 11/12] Use the mesh edge lengths in the FCT low-order flux test The FCT low-order flux test built its reference solution from a single hardcoded edge length per geometry, 0.012028130608117 on the plane and 277371 m on the sphere. The spherical value was the mean dvEdge of an earlier icosahedral CTest mesh, so the test measured how far each edge departs from that mean and its expected errors had to be re-derived whenever the mesh changed. Use Mesh->DvEdge instead. The low-order flux of a uniform tracer is then exactly DvEdge times the normal velocity, and the expected errors are zero on both meshes. Co-Authored-By: Claude Opus 5.5 (1M context) --- .../omega/test/ocn/TendencyTermsTest.cpp | 61 +++++++++---------- 1 file changed, 29 insertions(+), 32 deletions(-) diff --git a/components/omega/test/ocn/TendencyTermsTest.cpp b/components/omega/test/ocn/TendencyTermsTest.cpp index 13c5aae1fada..2e01e7b36c55 100644 --- a/components/omega/test/ocn/TendencyTermsTest.cpp +++ b/components/omega/test/ocn/TendencyTermsTest.cpp @@ -64,8 +64,7 @@ struct TestSetupPlane { ErrorMeasures ExpectedFCTHNew = {0.00, 0.00}; ErrorMeasures ExpectedFCT_High = {1.0408340855860843e-17, 1.2836244918860014e-17}; - ErrorMeasures ExpectedFCT_Low = {1.0095554103513882e-15, - 1.0095554103514163e-15}; + ErrorMeasures ExpectedFCT_Low = {0, 0}; ErrorMeasures ExpectedFCTToNonFCT = {9.524542988823467e-16, 5.336197203674528e-16}; ErrorMeasures ExpectedTrDel2Errors = {0.00334357193650093847, @@ -205,32 +204,32 @@ struct TestSetupSphere { // radius of spherical mesh Real Radius = REarth; - ErrorMeasures ExpectedDivErrors = {0.013659556526126423, - 0.0036698023569596279}; - ErrorMeasures ExpectedPVErrors = {0.021937128290005589, - 0.012253081309866708}; - ErrorMeasures ExpectedGradErrors = {0.0018790366180557886, - 0.0014984647360648147}; - ErrorMeasures ExpectedLaplaceErrors = {0.2819223346929331, - 0.27053032782204006}; - ErrorMeasures ExpectedTrHAdvErrors = {0.013259084476801913, - 0.004075236217689366}; - ErrorMeasures ExpectedFCTErrors = {0.00, 0.00}; - ErrorMeasures ExpectedFCTHProv = {3.0542657508680904e-05, - 1.0779233406323438e-06}; - ErrorMeasures ExpectedFCTHInv = {3.0541724683419424e-05, - 1.0779233406323438e-06}; - ErrorMeasures ExpectedFCTHNew = {3.0541724683419424e-05, - 1.0779233406323438e-06}; - ErrorMeasures ExpectedFCT_High = {0.00146484375, 16478.526025524854}; - ErrorMeasures ExpectedFCT_Low = {0.2958192441638326, 0.09980000200239414}; - ErrorMeasures ExpectedFCTToNonFCT = {1.5120472343676013e-15, - 6.215666716489436e-16}; - ErrorMeasures ExpectedTrDel2Errors = {0.015620987792697782, - 0.0032365878896264888}; - ErrorMeasures ExpectedTrDel4Errors = {0.00081985237645818541, - 0.00064699727012728107}; - ErrorMeasures ExpectedSurfTrRestErrors = {0, 0}; + ErrorMeasures ExpectedDivErrors = {0.013659556526126423, + 0.0036698023569596279}; + ErrorMeasures ExpectedPVErrors = {0.021937128290005589, + 0.012253081309866708}; + ErrorMeasures ExpectedGradErrors = {0.0018790366180557886, + 0.0014984647360648147}; + ErrorMeasures ExpectedLaplaceErrors = {0.2819223346929331, + 0.27053032782204006}; + ErrorMeasures ExpectedTrHAdvErrors = {0.013259084476801913, + 0.004075236217689366}; + ErrorMeasures ExpectedFCTErrors = {0.00, 0.00}; + ErrorMeasures ExpectedFCTHProv = {3.0542657508680904e-05, + 1.0779233406323438e-06}; + ErrorMeasures ExpectedFCTHInv = {3.0541724683419424e-05, + 1.0779233406323438e-06}; + ErrorMeasures ExpectedFCTHNew = {3.0541724683419424e-05, + 1.0779233406323438e-06}; + ErrorMeasures ExpectedFCT_High = {0.00146484375, 16478.526025524854}; + ErrorMeasures ExpectedFCT_Low = {0, 0}; + ErrorMeasures ExpectedFCTToNonFCT = {1.5120472343676013e-15, + 6.215666716489436e-16}; + ErrorMeasures ExpectedTrDel2Errors = {0.015620987792697782, + 0.0032365878896264888}; + ErrorMeasures ExpectedTrDel4Errors = {0.00081985237645818541, + 0.00064699727012728107}; + ErrorMeasures ExpectedSurfTrRestErrors = {0, 0}; ErrorMeasures ExpectedSfcStressForcingErrors = {0, 0}; ErrorMeasures ExpectedBottomDragErrors = {0.0015343843060931499, 0.0014861805971789941}; @@ -1481,14 +1480,12 @@ int testFCTTracerHorzAdvOnCell(int NVertLayers, int NTracers, Real RTol) { Kokkos::fence(); Array2DReal LowOrderFlx = TrHorzAdvOnC.GetLowOrderFlx(); + const auto DvEdge = Mesh->DvEdge; parallelFor( {Mesh->NEdgesHaloH(1), NVertLayers}, KOKKOS_LAMBDA(const int IEdge, const int K) { - const double DvEdgeFromMeshFile = - Geom == Geometry::Planar ? 0.012028130608117 : 277371; - FluxSubView(IEdge, K) = - DvEdgeFromMeshFile * NormalVelocity(IEdge, K); + FluxSubView(IEdge, K) = DvEdge(IEdge) * NormalVelocity(IEdge, K); }); const Real ATol = 1.0e-10; From 1a564eba3b31f716ea9ee0a2a27e670e7e231113 Mon Sep 17 00:00:00 2001 From: Xylar Asay-Davis Date: Wed, 23 Sep 2026 10:29:30 -0500 Subject: [PATCH 12/12] Update the spherical FCT provisional thickness expected errors The spherical FCT test advects with a constant lon/lat velocity, so the provisional layer thickness error is dt |div u| = 3 dt tan(lat) / R and scales with 1/R. On the regenerated icosahedral CTest mesh, with the Physical Constants Dictionary radius, the L2 errors of the provisional thickness and its inverse rise by 2.6e-5 relative (3.3e-5 from the radius alone) and the LInf errors by 1e-5 to 7e-5. Record the new values. The FCTHNewInv check still yields inf on every mesh; that is separate from the mesh change and left alone here. Co-Authored-By: Claude Opus 5.5 (1M context) --- components/omega/test/ocn/TendencyTermsTest.cpp | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/components/omega/test/ocn/TendencyTermsTest.cpp b/components/omega/test/ocn/TendencyTermsTest.cpp index 2e01e7b36c55..8e4624fa09ab 100644 --- a/components/omega/test/ocn/TendencyTermsTest.cpp +++ b/components/omega/test/ocn/TendencyTermsTest.cpp @@ -215,10 +215,10 @@ struct TestSetupSphere { ErrorMeasures ExpectedTrHAdvErrors = {0.013259084476801913, 0.004075236217689366}; ErrorMeasures ExpectedFCTErrors = {0.00, 0.00}; - ErrorMeasures ExpectedFCTHProv = {3.0542657508680904e-05, - 1.0779233406323438e-06}; - ErrorMeasures ExpectedFCTHInv = {3.0541724683419424e-05, - 1.0779233406323438e-06}; + ErrorMeasures ExpectedFCTHProv = {3.054304092320326e-05, + 1.077950920692842e-06}; + ErrorMeasures ExpectedFCTHInv = {3.054397382906693e-05, + 1.0779509220589763e-06}; ErrorMeasures ExpectedFCTHNew = {3.0541724683419424e-05, 1.0779233406323438e-06}; ErrorMeasures ExpectedFCT_High = {0.00146484375, 16478.526025524854};