`OnSphere` is currently not being read from the initial condition.

[hyungyukang]

OnSphere, which is used in SecondDerivativeOnCell in HorzOperators.h, is not currently configured from the initial condition. Instead, it is hard-coded to true in HorzOperators.cpp:

 30 SecondDerivativeOnCell::SecondDerivativeOnCell(HorzMesh const *Mesh)
 31     : OnSphere(true), NCellsAll(Mesh->NCellsAll), MaxEdges(1 + Mesh->MaxEdges),

As a result, SecondDerivativeOnCell always precomputes the second derivative coefficients for tracer advection on a sphere, even when the geometry is planar.

234       if (OnSphere) {
235          const Array1DI4 edgesOnCell =
236              Kokkos::subview(EdgesOnCell, ICell, Kokkos::ALL);
237          DetermineSphericalPatchGeometry(
238              NEdges, edgesOnCell, VerticesOnEdge, XCell, YCell, ZCell, XVertex,
239              YVertex, ZVertex, CellList, XP, YP, Angle2D, ThetaAbs[ICell]);
240       } else { // On an x-y plane
241          const Array1DI4 edgesOnCell =
242              Kokkos::subview(EdgesOnCell, ICell, Kokkos::ALL);
243          DeterminePlanerPatchGeometry(ICell, NEdges, edgesOnCell, CellsOnEdge,
244                                       AngleEdge, DcEdge, XP, YP, Angle2D);
245       }

This eventually leads to NaN values in AdvCoefs and AdvCoefs3rd when high-order horizontal tracer advection is used in the baroclinic channel and the overflow test.

IEdge I AdvCoefs(I, IEdge) AdvCoefs3rd(I, IEdge)
    2 1                nan                   nan
    2 2                nan                   nan
    2 3                nan                   nan
    2 4                nan                   nan
    2 5                nan                   nan
    2 6                nan                   nan
    2 7                nan                   nan

When OnSphere(false) is used for the tests on a plane domain, I confirmed that AdvCoefs and AdvCoefs3rd don't have NaN values and that the test completes successfully.

IEdge I AdvCoefs(I, IEdge) AdvCoefs3rd(I, IEdge)
    2 1            481.125                96.225
    2 2              769.8               -192.45
    2 3            -96.225                96.225
    2 4        3.66054e-15          -3.66054e-15
    2 5        2.06117e-14          -2.06117e-14
    2 6       -6.59004e-15           6.59004e-15
    2 7         7.0732e-15           -7.0732e-15

Therefore, OnSphere should be read from the initial condition file, as in MPAS-O, using the global attribute on_a_sphere, or handled through a simpler alternative if available.

[philipwjones]

@hyungyukang Thanks for catching this - there are a few global attributes we need to add to HorzMesh and read in from the input mesh file. That includes OnSphere, SphereRadius, IsPeriodic, XPeriod, YPeriod. I think I can add this to #382

[xylar]

@philipwjones, that sounds great!

[hyungyukang]

@philipwjones , sounds great. Thanks a lot!

[philipwjones]

@hyungyukang I've added these metadata to both the mesh structure and to the mesh file read routines so they are now initialized correctly once #382 is merged. I also passed the proper argument into the constructor.

While doing that, I noticed that the variable sphereRadius is being set to REarth in the HorzOperators. Should I change that to the mesh variable SphereRadius now?

[xylar]

While doing that, I noticed that the variable sphereRadius is being set to REarth in the HorzOperators. Should I change that to the mesh variable SphereRadius now?

It feels to me like we should raise an error if SphereRadius and REarth aren't the same (probably to a tolerance?). We don't want to use SphereRadius over REarth. We want Physical Constants Dictionary to be the source of truth. But it would be good to know if there's been a screw up on the Polaris side and the SphereRadius attribute isn't what it should be.

[hyungyukang]

@philipwjones , Thanks a lot for working on this! The change looks good.

I did a quick check and found a discrepancy between REarth from PCD and sphere_radius in some mesh files. REarth from PCD is 6371008.7714, while sphere_radius in ocean.QU.240km.151209.nc and oRRS18to6v3.171116.nc is 6371229..

I ran QU240 tests using both values, and the total KE at day 30 is different:

REarth_from_PCD       : 4.358571259295178e+17
SphereRadius_from_mesh: 4.363688674501280e+17

@xylar, then we keep the current code using REarth from PCD for now, and probably sphere_radius in the mesh files can be updated later for Omega?

[xylar]

@hyungyukang, yes, those are not files we would ideally be using for omega testing because of the discrepancy in Earth radius. It's a good point, though. @katsmith133 is working on a conversion script for MPAS-O to Omega and it could make sense to scale the Earth radius there, too. @katsmith133 , what do you think?

[katsmith133]

@xylar yup, I can add in a change to update sphere_radius in both the Omega init file and the new MPAS-O file if they are not consistent with REarth.

[xylar]

I think it's more complicated than just changing the metadata. I think we need to scale xCell, xVertex, etc. I'll comment on the Polaris PR.

[xylar]

@hyungyukang, @philipwjones and anyone else with an opinion, do you think we should have the capability to rescale the coordinates in Omega if sphere_radius doesn't match REarth or should it just be an error? I lean toward the latter but obviously it would be limiting in that existing initial conditions and meshes would need to be rescaled before Omega could use them.

[philipwjones]

While I suspect the errors are small enough that we could scale various mesh quantities and spread the error around, it's probably safer to regenerate a mesh with the proper radius so that we make sure we retain all the other mesh properties. So guess I have a slight preference for throwing an error and forcing a new mesh file.