Feature request: stretched grid file attributes and run-time error checking

[lizziel]

The MAPL stretched grid capability is now available for use by the general public in GCHP and we are starting to receive feedback. This week a user asked a question about the model not crashing when using a restart on the wrong grid. The grid size was correct but the restart was not generated using the same grid stretching as the simulation. The question was would this output be scientifically valid. The answer was no, there would be a silent bug.

To avoid this problem there could be (1) stretched grid information (target lat, lon, stretch factor) stored in the checkpoint file, and (2) grid validation of the restart file upon simulation start beyond just size. Having this information in the checkpoint files for easy referencing would help users stay organized and avoid mistakes. Having this information in History files would also be beneficial for the same reasons.

Pinging @LiamBindle, @bena-nasa, @atrayano, @tclune, @sdeastham for discussion.

[tclune]

I think this issue came up previously. Certainly such errors should be trapped, but there might be reasons that this is nontrivial.

[bena-nasa]

I'll give this a crack. It should not be that hard other than adding more complexity to what is already way too long of a code in MAPL_IO.F90. Perhaps I could add another function to the grid factories, some sort of check and this might make sense to do a level higher than MAPL_IO so avoid adding more to that, perhaps in the call in MAPL_Generic. Turns out the stretch parameters are already being set in the grid so let me first confirm I can retrieve them when needed.

[LiamBindle]

Just an idea, but an alternative method to enforcing this via attributes would be inspecting one of the corner coordinates of the grid. Specifically, one can calculate the stretch-factor, target lat, and target lon from a single point. If the calculated stretch-factor is 1.0, the target lat is -90, and the target lon is 170 then it is a normal cubed-sphere grid. The advantage would be it doesn't require new attributes.

Sorry for the late comment, but before I didn't realize it was pretty easily doable. I'm happy to write a routine that does this if it would be helpful.

[tclune]

A bit of a kludge, but at first glance it does look like it would be reasonably robust. I've not looked at the stretch formulation recently, so I don't know how easily invertible it is. Or do we just need to check that a few points are identical in which case the stretch parameters are identical. In theory we need to be a bit careful as there should be some degeneracies, but I doubt those are an issue in practical cases.

I'll await @bena-nasa 's analysis.

[sdeastham]

I've not quite thought it through, but could you use the fact that you have multiple workers to make this non-degenerate? If you always have >= 6 workers so that there is always one worker whose [1,1] grid cell is at a cube corner, I think that if you have every worker check its corner then you'll always check all 6 corners (at minimum) and fully locate the grid?

[LiamBindle]

The stretching procedure is

1. apply schmidt transform [can be solved for the stretch-factor]
2. rotate about (0,1,0) [reversible]
3. rotate about (0,0,1) [reversible]

so I think it's doable. The rotations can be determined by inspecting the center of face 6; the center of face 6 is well-defined (south pole) for the GMAO/ESMF cubed-sphere so for any grid we can figure out how much it was rotated. Once we know how much the grid was rotated we can apply the reverse rotations to a single point and then plug it in to the Schmidt transform solved for the stretch-factor.

[LiamBindle]

Here's a proof-of-concept for getting the stretching parameters according to the coordinates of two points (center of face 6 and the coordinate of some grid-box on face 6): get_stretching_params_demo.zip. Btw, this POC doesn't handle the special case where the center of face 6 is at one of the poles, but it could be added.

[tclune]

ok - thanks. We'll look into it. Note that even determining on any given process, whether or not it owns the points in question is a bit nontrivial. And presumably we'll allow some numerical tolerance so that we don't reject things where the numerical inverse is not exact.

[lizziel]

I'm still in favor of having stretch information in the file metadata for reference. It seems like basic important information that should be easily accessible without having to do math.

[LiamBindle]

I agree. That way you could ncdump the file to see the grid info which would be useful.

Sorry, to clarify, I just meant to share that the run-time check could be implemented with only the grid's coordinates, in case it's useful.

[bena-nasa]

I've been working on this. The strategy I'm working on to

1. Add the stretch coordinates to the restarts when writing.
2. Modify the cubed-sphere grid factory to when it tests equality it takes these into account
3. When we read the restart make a factory from the file and check it it matches the application grid

This is a little more involved that I thought it would be for technical reasons.

[bena-nasa]

This has been implemented, turned out to be harder than I thought and unfortunately only works for cubed-sphere (although as you pointed out that is all GCHP uses so you should be good). Note that once you take this, users who are currently running stretched grid simulations and have restarts already will need to add these attributes to the restart so it doesn't fail. Not a problem if you are bootstrapping the restart but if they already have them, then they would need to do this.

[bena-nasa]

@lizziel I'm going to close this issue as it is implemented and is in the 2.8.2 release

[lizziel]

Thanks @bena-nasa. Heads up on this @LiamBindle since it means the attributes needs to be added to the stretched CS regrid tool you developed.

[LiamBindle]

@lizziel Thanks for the reminder! Do you know what GCHP version this will come in?