implement a procedure to carry out a regridding from high to low resolution (binning)

[mo-lormi]

Description

This PR is intended to extend the group of the interpolation operations, in particular to:

• implement a procedure to evaluate the areas of the grid cells (grid type: CS-LFR)
• implement a procedure to carry out a regridding from high to low resolution (binning)
• add new tests to the test battery
  
  
  

Rigridding from high to low resolution via binning ...

method: binning
binning matrix: B = N A^T W
area weights matrix: W
interpolation matrix: A
normalization factor matrix: N

Acceptance Criteria (Definition of Done)

All the tests pass.

[mo-lormi]

@wdeconinck , @odlomax This PR covers the implementation of the regridding procedure (from high to low resolution).
At this stage, I implemented a simple test associated with the procedure to generate data in Gmsh format for data visualization. Further tests will be added.

[odlomax]

Thanks, @mo-lormi

Could you add some of the outputs to this PR so we can easily have a look?

[mo-lormi]

The following plots refer to two fields:

1. field: field_01_s ; regridding type: high to low; source grid: CS-LFR-100
2. field: field_01_t ; regridding type: high to low; target grid: CS-LFR-50

[twsearle]

This looks very nice! I am curious as a bystander why you haven't used the interpolation elements. Perhaps if you had it would have worked for any grid? https://github.com/ecmwf/atlas/blob/develop/src/atlas/interpolation/element/Quad2D.h

Maybe out of scope though!

[mo-lormi]

I carried out further testing. The procedure also works for a different type of grid (function space: StructuredColumns).

The following plots refer to two fields:

1. field: field_02_s ; regridding type: high to low; source grid: O32
2. field: field_02_t ; regridding type: high to low; target grid: O16

(note -- It seems that a set of contiguous cells is missing in both plots. Thus, further investigation is needed)

[wdeconinck]

@mo-lormi concerning the "missing" part in the Gmsh plot, is nothing to worry about. That's because the StructuredColumns fields are unlike NodeColumns and don't have extra duplicated nodes at 360 degrees.

[mo-lormi]

@mo-lormi concerning the "missing" part in the Gmsh plot, is nothing to worry about. That's because the StructuredColumns fields are unlike NodeColumns and don't have extra duplicated nodes at 360 degrees.

@wdeconinck thanks for this info ... To take into account this point, I will need to update the procedure.

[wdeconinck]

Thank you and sorry for the delayed review. Please see my in-code comments.

[wdeconinck]

In order to encapsulate as much as possible the grid from objects like Field, FunctionSpace, Mesh,
I would avoid this function.

Possibly you could create a different function

double CubedSphere::gridCellArea(idx_t index);

The grid itself should have all possible information needed to compute just a single cell value.

Then where needed fill a container such as std::vector or atlas::Field or ...
But in this case here, I think you could just compute the cellArea within the loop where you compute the weights, avoiding the entire allocation of a container.

[mo-lormi]

@wdeconinck ... good point. However, there is reason why the method CubedSphere::gridCellArea() has been designed and implemented in this way.

So far, at the MO, we have identified at least two workflow steps where this function is needed: one is the evaluation of the area weights matrix within the interpolation/binning procedure (this PR); the other one is the evaluation of the Jc term, a penalty term associated with the gravity wave constraints for the cost function.
Thus, we thought that having a function signature for CubedSphere::gridCellArea() based on the FunctionSpace data type could make it general enough to be used in different contexts.

[wdeconinck]

I managed to boil this down to this standalone function, independent of functionspace

// Area around nodes of cubed sphere mesh
Field CubedSphereNodalArea(const Mesh& mesh) {

  constexpr auto deg2rad = M_PI / 180.;
  const auto& proj = mesh.projection();
  auto lonlat = array::make_view<double, 2>(mesh.nodes().lonlat());
  auto gcell_area_field = Field("grid_cell_areas", make_datatype<double>(), array::make_shape(mesh.nodes().size()));  
  auto gcell_area_fview = array::make_view<double, 1>(gcell_area_field);

  double gcell_area_cs = [&] {
      /******* CUBED_SPHERE SPECIFIC *********/
      // (grid_res * grid_res) = no. of cells on a tile
      ATLAS_ASSERT(CubedSphereGrid(mesh.grid()));
      auto grid_res = CubedSphereGrid(mesh.grid()).N();
      // area of a grid cell (cubed-sphere coord. system)
      return M_PI/(2*grid_res) * M_PI/(2*grid_res);
   }();

  for (size_t i = 0; i < gcell_area_fview.size(); i++) {
    PointLonLat loc = PointLonLat(lonlat(i, atlas::LON), lonlat(i, atlas::LAT));
    double cos_lat = std::cos(deg2rad * loc.lat());
    double grid_jac_det = 1/proj.jacobian(loc).determinant();
    // area of a grid cell (geographic coord. system)
    gcell_area_fview(i) = grid_jac_det * gcell_area_cs * cos_lat ;   /****** USE gcell_area_cs here ******/
  }

  return gcell_area_field;
}

Possibly with more thought this could be made even more independent of CubedSphere.

What happens to the Binning procedure if the combination does not happen to be "grid == CubedSphere && functionspace == NodeColumns" ?

[mo-lormi]

@wdeconinck ... thanks for this.
Where do you think this function should be stored? In which file, class,...?
If you tell me where I should define the function, I will do a test with it. Thanks ...

[wdeconinck]

It could be a mesh-action, similar to https://github.com/ecmwf/atlas/blob/develop/src/atlas/mesh/actions/BuildCellCentres.h
So putting it parallel to this file, and following its design is a good idea.
You can then "register" this field within the mesh.nodes() and could be reused later.

What happens to the Binning procedure if the combination does not happen to be "grid == CubedSphere && functionspace == NodeColumns" ?

[mo-lormi]

functionspace == NodeColumns

The binning procedure works also for "functionspace == StructuredColumns" (grid type: "O").
I implemented a test to check this scenario, see here.

[mo-lormi]

@wdeconinck thanks again for the comments and suggestions.
I re-implemented the procedure to evaluate the areas of the grid cells based on your proposal.
Thanks also to @odlomax for helping to identify the cause of a segmentation fault failure.

[wdeconinck]

Thanks so much for this contribution and addressing all my suggestions. 🙏

[mo-lormi]

Thanks so much for this contribution and addressing all my suggestions. 🙏

@wdeconinck thanks again.

[wdeconinck]

This will be released in 0.38.0 some time this week.