Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Set UnstructuredGrid in meshes built by MeshBuilder #150

Closed
fmahebert opened this issue Jun 22, 2023 · 9 comments · Fixed by #152
Closed

Set UnstructuredGrid in meshes built by MeshBuilder #150

fmahebert opened this issue Jun 22, 2023 · 9 comments · Fixed by #152

Comments

@fmahebert
Copy link
Contributor

Is your feature request related to a problem? Please describe.

Some JEDI applications will want to use the MeshBuilder but also to call mesh.grid(). Currently the grid is not set.

Describe the solution you'd like

Setting the grid requires a global list of point coordinates.

I propose this change: MeshBuilder takes an MPI comm as argument, does an allGather of the (owned) point lons/lats, and constructs an UnstructuredGrid from the global list of points.

But there are some alternative options:

  • MeshBuilder takes a global list of lons/lats as argument (instead of currently taking local points) and uses these to construct an UnstructuredGrid. For setting up the connectivities, getting the right subset from the global list will I think require either a new indexing argument or requiring a structure in the global_indices (e.g.: 1-based, continuous indexing) of the points so that they can be reliably indexed.
  • MeshBuilder takes a Grid as argument, putting the responsibility on the user.

I think the propose change is a more straightforward design and easier to think about. The alternative options put more burden on the user, and are a more complicated design (build a global list but then need to find a subset), but do avoid adding an MPI comm to the argument list.

@wdeconinck do you have a preferred design?

Describe alternatives you've considered

No response

Additional context

No response

Organisation

UCAR/JCSDA

@wdeconinck
Copy link
Member

Even though I understand the convenience of mesh.grid(), there's also a downside to it for UnstructuredGrid:
every mesh partition will have a full copy of all the global grid coordinates. On HPC nodes with a lot of MPI ranks per node, and with (extremely) large grids this could put quite some strain on the memory.
So you may not necessarily want to rely on this for unstructured grids.

However your proposed solution (Solution 1) is indeed one of the better ones, also concerning memory if we can keep the grid optional.
I would prefer you add a configuration argument to set the grid.
If the grid is a known named grid, then you don't even need to assemble it. It could e.g. be a structured grid without memory. For unstructured grids you can also set it to type unstructured inferring the instructions to allGather the owned points as you described. If it is not configured, then no grid will be set.
Having an optional strongly typed Grid argument (solution 3) could be an alternative to configuring the grid as well.

The mpi communicator should also be an optional configuration parameter where its value is a string key in the eckit::Comm map. If not passed, the default communicator is used. This is in line with the other issue #131

@fmahebert
Copy link
Contributor Author

fmahebert commented Jun 27, 2023

Thanks for your comments on this, @wdeconinck.

I've been thinking more about our use-case in JEDI, and we should be able to handle the unstructured case separately without constructing the full UnstructuredGrid in the MeshBuilder. I would still want to call mesh.grid() in cases where the Grid exists (perhaps a NodeColumns constructed from something with a StructuredGrid), but handle separately the case where the Grid does not exist.

How do you recommend we set up the MeshBuilder to allow this? Currently the MeshBuilder creates a Mesh whose Grid is not initialized, so any calls to it will be a memory error. I haven't yet found an easy way to flag the Mesh's Grid as being invalid from within the Mesh itself. I thought of initializing the Grid like ~ mesh.setGrid(UnstructuredGrid{}); hoping to enable call something like mesh.grid().valid(), but that didn't work (because of trying to copy the grid.projection() in MeshImpl).

Do you have any advice on how we can more gracefully initialize the mesh's grid, or more gracefully check whether it's safe to use mesh.grid()?

@fmahebert
Copy link
Contributor Author

fmahebert commented Jun 29, 2023

or more gracefully check whether it's safe to use mesh.grid()

I've found it: I can use ObjectHandleBase::operator bool() to just call if (mesh.grid()) to determine if I can query the Grid or not. Lovely!

Using this feature to identify if the Grid is initialized + a JEDI workaround for uninitialized (and implicitly unstructured) grids, I think the develop code suites our needs at present.

Independently of this, would you like to see the MeshBuilder set the Mesh's Grid, @wdeconinck ? Note that I anticipate most use cases for MeshBuilder would have an UnstructuredGrid, because if the Grid were easily constructible from atlas then there are more elegant interfaces to get to a Mesh.

@wdeconinck
Copy link
Member

Indeed sorry that was not obvious before that you could use if (mesh.grid())
You could add a configuration option like:

  grid: 
    type: unstructured

You can do this programmatically with config.set("grid.type","unstructured")
The unstructured grid configuration normally also has a xy member:

grid:
  type: unstructured
  xy : [
    180,0,
    90,0,
    -90,0,
    0,90,
    0,-90,
    0,0,
    ... 
  ]

We could agree that IF the grid type is "unstructured", and IF the "xy" entry is missing, then you have to assemble it via AllGather of owned nodes and global_index yourself within the MeshBuilder.
You will then have to call setGrid internally inside the MeshBuilder. You have to make the MeshBuilder a friend class of the Mesh to access the private setGrid function, similar to the MeshGenerator.
If the "grid" configuration is missing altogether, don't create the Grid at all.

@fmahebert
Copy link
Contributor Author

fmahebert commented Jul 3, 2023

This all makes sense, and looks a lot like what I have on a branch. But how do you want to handle the case where a grid is given including its grid points (perhaps unstructured with "xy" list, or some named structured grid)... in a probabilistic sense, most mesh-related input arguments to the MeshBuilder will not be consistent with a particular named atlas grid, and it might be nice to prevent the user from constructing something ill-posed because the grid and mesh are inconsistent. I see three avenues:

  1. [most permissive] proceed to build the grid from config, the mesh from function arguments, and hope for the best
  2. [most reasonable but most work] check the nodes in the function arguments are consistent with the grid in the config
  3. [most safe] throw an exception warning the user that as the grid can be constructed in atlas, maybe the mesh connectivity should be as well, and avoid the problem

What do you think is the way to proceed? I suggest starting with (3) and writing the code for (2) if/when needed...

@wdeconinck
Copy link
Member

I think 2. is a good approach for safety indeed. Each partition, using its global index, could verify the grid coordinates are as expected.
For 3. you may have used an alternative mesh generation tool within the model that you wish to describe the connectivities, so I would not do this.
So in summary, I think 1. is OK and 2. is a little bit more work, and is more of a debugging option. However there are grids like e.g. the HEALPix grid or the Cubed sphere, where the nodes are not the grid points, rather the cell centres represent the grid points. Any check then is going to fail.

@fmahebert
Copy link
Contributor Author

fmahebert commented Jul 10, 2023

Each partition, using its global index, could verify the grid coordinates are as expected.

@wdeconinck I want to verify my understanding of how this can be implemented...

  • as far as the MeshBuilder is concerned, the global_index can be some arbitrary labeling (possibly non-contiguous) over the nodes, and I could use essentially any choice of global_index as long as the indices are unique (and perhaps not 0-based, for gmesh?).
  • for indexing into the Grid, the global_index (as used in UnstructuredGrid::lonlat(idx_t) or StructuredGrid::index2ij(gidx_t, indx_t, indx_t)) is a 0-based contiguous index in the more traditional sense of the word "index"

But now I want to check consistency between the Mesh and the Grid, perhaps with code like

for (size_t n = 0; n < nb_nodes; ++n) {
  const gidx_t gidx = global_indices[n];  // global_indices passed to MeshBuilder, the arbitrary label
  double lonlat[2];
  // assuming an UnstructuredGrid (can write similar code for Structured):
  grid.lonlat(gidx, lonlat);  // here the grid assumes a 0-based contiguous index
  // check lonlat matches lons[n],lats[n]  // lons, lats passed to MeshBuilder
}

My first question is: does the above understanding seem correct?

If all the above seems correct, then I need some way to connect these two global_indices definitions. Probably by doing one of...

  • requiring that the MeshBuilder's input global_indices has the same structure as the Grid's ordering (perhaps at odds with gmesh?)
  • adding a new argument to the config used to build the Grid; some type of vector<int> to serve as a mapping between the two different global_indices.

Alternatively, I can use a more complex validation algorithm: instead of linearly iterating over the points, do a search to check every point passed to the MeshBuilder is a point in the Grid. A bit computationally demanding perhaps, especially for UnstructuredGrids.

My second question is: do you have an opinion, or an alternative suggestion, to the above?

@wdeconinck
Copy link
Member

For nearly all meshes, the nodes.global_index() will indeed follow the grid ordering with 1-base. This is the assumption that should be made, because it is also how the fields are ordered via MPI Gather. So grid.lonlat(gidx-1, lonlat) in your example.
You have to make sure that the global_index is constructed that way.
It may not always be possible to match mesh points with grid points either. A periodic node for instance has a global_index which is not present in the grid. These points have a global index which is higher than the grid.size(), and could be discarded from the check.

Only for some meshes like HEALPix, or the CubedSphere, the grid-points denote the cell-centres of the mesh. In that case it should be cells.global_index() which follows the grid ordering with 1-base. You should probably disable the validation then, as cells don't have a lonlat field, or at least not a field that is uniquely defined (probably a check if the grid point is contained in the cell could work).

Whether to do any validation check at all could be a configuration option passed to the MeshBuilder constructor.
I would not use more complex validation algorithms.

@fmahebert
Copy link
Contributor Author

Thanks @wdeconinck for the info and opinions above. I agree with sticking to a simple validation that is opt-in via config, assumes the global_index is a 1-based contiguous index, also assumes the mesh is connecting the grid nodes, and just won't work in fully general cases.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging a pull request may close this issue.

2 participants