Require the PCD Earth radius for spherical meshes

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 ..
 ```
 

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 <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,

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. |

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;
+      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;

components/omega/src/ocn/HorzOperators.h

@@ -342,6 +342,7 @@ class SecondDerivativeOnCell {
        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]    = {};

components/omega/src/timeStepping/TimeStepper.cpp

@@ -688,7 +688,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(

components/omega/test/ocn/AuxiliaryStateTest.cpp

@@ -26,6 +26,7 @@
 using namespace OMEGA;
 
 struct TestSetup {
+   // radius of spherical mesh
    Real Radius = REarth;
 
    KOKKOS_FUNCTION Real pseudoThickness(Real Lon, Real Lat) const {

components/omega/test/ocn/AuxiliaryVarsTest.cpp

@@ -146,41 +146,42 @@ struct TestSetupPlane {
 
 struct TestSetupSphere {
 
+   // radius of spherical mesh
    Real Radius = REarth;
 
-   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));

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);

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

components/omega/test/ocn/HorzOperatorsTest.cpp

@@ -80,23 +80,22 @@ struct TestSetupPlane {
 
 struct TestSetupSphere1 {
    // radius of spherical mesh
-   // TODO: get this from the mesh
    Real Radius = REarth;
 
-   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);
@@ -132,7 +131,6 @@ struct TestSetupSphere1 {
 
 struct TestSetupSphere2 {
    // radius of spherical mesh
-   // TODO: get this from the mesh
    Real Radius = REarth;
 
    ErrorMeasures ExpectedDivErrors         = {1.37734693033362766e-10,

components/omega/test/ocn/TendenciesTest.cpp

@@ -31,6 +31,7 @@
 using namespace OMEGA;
 
 struct TestSetup {
+   // radius of spherical mesh
    Real Radius = REarth;
 
    KOKKOS_FUNCTION Real pseudoThickness(Real Lon, Real Lat) const {