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. | diff --git a/components/omega/src/ocn/HorzMesh.cpp b/components/omega/src/ocn/HorzMesh.cpp index 78f041584e1d..838287269cf3 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; + 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; diff --git a/components/omega/src/ocn/HorzOperators.h b/components/omega/src/ocn/HorzOperators.h index bcb96b12dcab..a78eba11fe5d 100644 --- a/components/omega/src/ocn/HorzOperators.h +++ b/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] = {}; diff --git a/components/omega/src/timeStepping/TimeStepper.cpp b/components/omega/src/timeStepping/TimeStepper.cpp index 759fdaae3198..a480ba7b561d 100644 --- a/components/omega/src/timeStepping/TimeStepper.cpp +++ b/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( diff --git a/components/omega/test/ocn/AuxiliaryStateTest.cpp b/components/omega/test/ocn/AuxiliaryStateTest.cpp index f9bd3633bfae..891f14aa9470 100644 --- a/components/omega/test/ocn/AuxiliaryStateTest.cpp +++ b/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 { diff --git a/components/omega/test/ocn/AuxiliaryVarsTest.cpp b/components/omega/test/ocn/AuxiliaryVarsTest.cpp index 4f916d998c7d..260ae8cbbf87 100644 --- a/components/omega/test/ocn/AuxiliaryVarsTest.cpp +++ b/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)); 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/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 diff --git a/components/omega/test/ocn/HorzOperatorsTest.cpp b/components/omega/test/ocn/HorzOperatorsTest.cpp index f4af88fd369c..bab6fc845b6c 100644 --- a/components/omega/test/ocn/HorzOperatorsTest.cpp +++ b/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, diff --git a/components/omega/test/ocn/TendenciesTest.cpp b/components/omega/test/ocn/TendenciesTest.cpp index cf53a858d2eb..c20d9cf8a936 100644 --- a/components/omega/test/ocn/TendenciesTest.cpp +++ b/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 { diff --git a/components/omega/test/ocn/TendencyTermsTest.cpp b/components/omega/test/ocn/TendencyTermsTest.cpp index aa4698b87b5f..8e4624fa09ab 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, @@ -203,42 +202,41 @@ struct TestSetupPlane { struct TestSetupSphere { // radius of spherical mesh - // TODO: get this from the mesh Real Radius = REarth; - 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 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.04865718541236144, - 0.005105510870642706}; - ErrorMeasures ExpectedTrDel4Errors = {0.0008646345116716073, - 0.0007118574326665881}; - 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.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}; + 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.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; } @@ -316,17 +314,27 @@ 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, where Lon is undefined. 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; } @@ -1472,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;