Cartesian grid generation

[thabbott]

Purpose

Implement grid generation for doubly-periodic Cartesian domains.

Remove the sections below which do not apply.

Code changes:

• util/pace/util/grid/generation.py: Add fields to MetricTerms class to store grid spacing and latitude for Cartesian grids, and add an _init_cartesian method to compute Cartesian grid properties. Also fix a bug in vertical grid calculations that was preventing grid translate tests from passing.
• fv3core/tests/savepoint/translate/translate_grid.py: Pass grid type and Cartesian grid spacing and latitude to MetricTerms constructor in InitGrid and InitGridUtils tests.

Requirements changes:

• constraints.txt: Update h5py and netcdf versions. (make build failed with versions in unchanged constraints.txt.)

Checklist

Before submitting this PR, please make sure:

• You have followed the coding standards guidelines established at Code Review Checklist.
• Docstrings and type hints are added to new and updated routines, as appropriate
• All relevant documentation has been updated or added (e.g. README, CONTRIBUTING docs)
• [TODO] Unit tests are added or updated for non-stencil code changes

Additionally, if this PR contains code authored by new contributors:

• The names of all the new contributors have been added to CONTRIBUTORS.md

[FlorianDeconinck]

+1

Can you pass a "fill_all_fields" option to the __init__ default to False that's flipped for test?

[FlorianDeconinck]

LGTM.

Look at putting non-required floats on a bool for test only

[oelbert]

Commenting for now with some small tweaks that could make the code a bit neater and more optimized. The init_cartesian method does work under a different paradigm than the non-cartesian initialization, so I want to think a bit about that

[oelbert]

Since dx, dxc, and dxa are all the same on an orthogonal grid I think you could consolidate these into one _set_dxdy_cartesian method

[thabbott]

I have three separate methods because _dx, _dx_agrid, and _dx_center (for example) sit at different locations on the grid. Could still work on reducing code duplication if you think that's worthwhile, e.g. with a single method like

def _compute_dxdy_cartesian(self, dx_dims, dy_dims):
    dx_64 = self.quantity_factory.zeros(dx_dims, ...)

that can handle all three cases.

[oelbert]

Similar to the dx/dy computation above, these areas are also the same so you could have one _set_area_cartesian method that can do either/both of these instead of two separate methods

[thabbott]

Same as above: areas sit at different locations, but could still reduce code duplication.

[thabbott]

Commenting for now with some small tweaks that could make the code a bit neater and more optimized. The init_cartesian method does work under a different paradigm than the non-cartesian initialization, so I want to think a bit about that

One alternative: could modify _init_cartesian so it only precomputes things that are also precomputed in _init_dgrid and _init_agrid, and add an _is_cartesian field to MetricTerms that's used when computing cached properties:

@property
def dx(self) -> util.Quantity:
    """
    the distance between grid corners along the x-direction
    """
    if self._dx is None:
        self._dx, self._dy = self._compute_dxdy_cartesian() if self._is_cartesian else self._compute_dxdy()
    return self._dx

[oelbert]

One alternative: could modify _init_cartesian so it only precomputes things that are also precomputed in _init_dgrid and _init_agrid, and add an _is_cartesian field to MetricTerms that's used when computing cached properties:

On reflection I think this would be my preference. If the metric terms generation works the same way for each grid type that makes things much simpler to deal with if we need to revisit or modify it going forward.

[thabbott]

One alternative: could modify _init_cartesian so it only precomputes things that are also precomputed in _init_dgrid and _init_agrid, and add an _is_cartesian field to MetricTerms that's used when computing cached properties:

On reflection I think this would be my preference. If the metric terms generation works the same way for each grid type that makes things much simpler to deal with if we need to revisit or modify it going forward.

I just pushed a commit with this change. Let me know what you think. A couple of thoughts:

1. There's more potential for this change to introduce bugs in cubed sphere grid generation. Are there some translate tests I should run to check that this hasn't happened?
2. The grid generation code now has a lot of additional if statements, and some are duplicated. For example:

self._dx, self._dy = (
    self._compute_dxdy_cartesian()
    if self._is_cartesian
    else self._compute_dxdy()
)

appears in dx, dy, and _calculate_divg_del6. Upon further reflection, I think a better choice would be structure things so the above code can just be

self._dx, self._dy = self._compute_dxdy(),

which I could do by changing _compute_dxdy to

def _compute_dxdy(self):
   if self._is_cartesian:
      return self._compute_dxdy_cartesian()
  ... # cubed-sphere calculation 

This seems better to me since if self._is_cartesian only shows up in once place. Does this also seem better to you?

[oelbert]

@FlorianDeconinck @bensonr Thoughts?

[FlorianDeconinck]

There's a way to have the best of both world. The key here is to keep to the "one idea one code" mantra. Call to dx triggering a lazy compute is one concept, needs for different internal numerics is another. So my advice would be to keep configuration in init. E.g.

class MetricTerms:

 def __init__(...):
   ...
   # Configuration block for internal numerics that differ. Requires functions to share signature
   if is_cartesian:
       self._compute_dx_dy = self._compute_dx_dy_cartesian
       ....
   else:
       self._compute_dx_dy = self._compute_dx_dy_cube_sphere
       ...    
   ...

# Lazy function remains the same and will pull on the proper numerics
def dx(self,...):
   ...
   dx = self._compute_dx_dy(...)

[FlorianDeconinck]

If we need to implement other grid in the future, we will break out the numerics in it's object and use OO and ABC to make it a true configuration pattern. With two it's fine to keep it as a if/else

[oelbert]

There's a way to have the best of both world. The key here is to keep to the "one idea one code" mantra. Call to dx triggering a lazy compute is one concept, needs for different internal numerics is another. So my advice would be to keep configuration in init. E.g.

class MetricTerms:

 def __init__(...):
   ...
   # Configuration block for internal numerics that differ. Requires functions to share signature
   if is_cartesian:
       self._compute_dx_dy = self._compute_dx_dy_cartesian
       ....
   else:
       self._compute_dx_dy = self._compute_dx_dy_cube_sphere
       ...    
   ...

# Lazy function remains the same and will pull on the proper numerics
def dx(self,...):
   ...
   dx = self._compute_dx_dy(...)

This looks like a good way to go to me

[thabbott]

@FlorianDeconinck Thanks for the suggestion! I just pushed a commit that implements it, and adds a unit test that checks some basic properties of Cartesian grids. These changes pass the InitGrid and InitGridUtils translate tests using test data from either Cartesian or cube sphere grids---happy to document this in some form if you'd like.

[oelbert]

Thanks for all the revisions, this is great!

[FlorianDeconinck]

Looking good, thanks for the utest.

All this will need a good documentation pass, but it needs to be done globally so it's out-of-scope here and we will probably attack this post framework/model refactor.

Great to see new contributor and clean code going in!