Support unstructured grid (Issue #125)

[gauravity30]

Hey !
I am looked into uxarray and attempted two methods to include support for unstructured grids. The fundamental working of both the methods is same.

Proposed Methodology

• Step 1: Use readnek to import the data
• Step 2: Create an array of nodes, elements, and fields
• Step 3: Create a grid and xarray dataset from the data array
• Step 4: Create a uxarray dataset and return it.

Method 1
Extend the current xarray method to uxarray as uxarray also supports the creation of dataset through load_store(). But unfortunately, it doesn't work. I get the following error.

conflicting sizes for dimension 'y': length 6 on 'xmesh' and length 36 on {'x': 'x', 'y': 'y', 'z': 'z'}.

Method 2
It is manual method to extract data from elements and store them as numpy arrays, use them to form a xarray dataset. I am able to make an xarray dataset which can be used for most post-processing stuff like plotting, POD calculations, etc. It lacks the element connectivity data, so is it required ? If not, this method can be used to handle structured and unstructured arrays both.

extract_elem_data is the implementation of Method 2. We may not require the use of uxarray.

A uxarray Grid can be constructed from the coordinates and connectivity data through uxarray.Grid.from_dataset and used with the available values to create a uxarray dataset using uxarray.UxDataset. Can someone help me in extracting the connectivity, so I can attempt creation of an uxarray dataset.

Thanks !
Gaurav

[ashwinvis]

You got an error with the existing _NekDataStore because of this method:

    def meshgrid_to_dim(self, mesh):
        """Reverse of np.meshgrid. This method extracts one-dimensional
        coordinates from a cubical array format for every direction
        """

I don't see a proper way of building 1D coordinate arrays for an unstructured arrays. We can add the DataStore later, but a working implementation is more important. Good that you chose to go for simpler functions.

[ashwinvis]

You can get the connectivity data at least for the EXODUS mesh like this:

import xarray as xr
mesh = xr.open_dataset("naca(10x4).e")
print(mesh.data_vars)
mesh.connect1

This integer data array has the dimensions connect1 (num_el_in_blk1, num_nod_per_el1) int32 766kB ... and looks like this:

array([[    3,     2,     1, ...,     5,     8,     7],
       [   10,     9,     2, ...,    11,     6,    13],
       [   15,    14,     9, ...,    16,    12,    18],
       ...,
       [72683, 72686, 32640, ..., 72732, 32642, 72731],
       [72686, 72689, 32645, ..., 72733, 32647, 72732],
       [72689, 72692, 32650, ..., 72734, 32652, 72733]], dtype=int32)

I am guessing it denotes one cell per row, and along the column denotes the 8 cells which surround it: 4 on each corner and 4 on each face (you need to double check it). Somebody wrote a short blog post looking at this using netCDF4, but above I used xarray (which also uses netCDF4 / h5netcdf behind the scenes).

https://johnfoster.pge.utexas.edu/blog/posts/extracting-exodus-information-with-netcdf-python/

[ashwinvis]

One thing I was not able to get it working was opening the EXODUS mesh file directly from Uxarray, although they claim it to be supported:

https://uxarray.readthedocs.io/en/latest/user-guide/grid-formats.html#exodus

[gauravity30]

Yeah ! But I think extracting mesh or any mesh data from the exodus file will not provide us with the correct mesh data. The NEK5000 visualization file has the polynomial nodes data along with the nodes of the original mesh. Correct me, if I am wrong but I think this is the case only. So, in that case how to extract the connectivity data from the visualization file.

Thanks, I will try to implement a working function first then we can establish a class like for xarray.

[gauravity30]

What was your error ?

One thing I was not able to get it working was opening the EXODUS mesh file directly from Uxarray, although they claim it to be supported:

https://uxarray.readthedocs.io/en/latest/user-guide/grid-formats.html#exodus

[ashwinvis]

import uxarray as ux
ux.open_grid("naca(10x4).e")

It assumes a 3D grid

---------------------------------------------------------------------------
KeyError                                  Traceback (most recent call last)
~/Sources/foss/pymech/.venv/lib/python3.12/site-packages/xarray/core/dataset.py in ?(self, name)
   1482             variable = self._variables[name]
   1483         except KeyError:
-> 1484             _, name, variable = _get_virtual_variable(self._variables, name, self.sizes)
   1485

KeyError: 'node_z'

During handling of the above exception, another exception occurred:

KeyError                                  Traceback (most recent call last)
~/Sources/foss/pymech/.venv/lib/python3.12/site-packages/xarray/core/dataset.py in ?(self, key)
   1585                 if isinstance(key, tuple):
   1586                     message += f"\nHint: use a list to select multiple variables, for example `ds[{[d for d in key]}]`"
-> 1587                 raise KeyError(message) from e
   1588

~/Sources/foss/pymech/.venv/lib/python3.12/site-packages/xarray/core/dataset.py in ?(self, name)
   1482             variable = self._variables[name]
   1483         except KeyError:
-> 1484             _, name, variable = _get_virtual_variable(self._variables, name, self.sizes)
   1485

~/Sources/foss/pymech/.venv/lib/python3.12/site-packages/xarray/core/dataset.py in ?(variables, key, dim_sizes)
    216     split_key = key.split(".", 1)
    217     if len(split_key) != 2:
--> 218         raise KeyError(key)
    219

KeyError: 'node_z'

The above exception was the direct cause of the following exception:

KeyError                                  Traceback (most recent call last)
<ipython-input-7-18c901787b96> in ?()
----> 1 ux.open_grid("naca(10x4).e")

~/Sources/foss/pymech/.venv/lib/python3.12/site-packages/uxarray/core/api.py in ?(grid_filename_or_obj, latlon, use_dual, **kwargs)
     84         try:
     85             grid_ds = xr.open_dataset(grid_filename_or_obj, **kwargs)
     86
     87             uxgrid = Grid.from_dataset(grid_ds, use_dual=use_dual)
---> 88         except ValueError:
     89             raise ValueError("Inputted grid_filename_or_obj not supported.")
     90
     91     return uxgrid

~/Sources/foss/pymech/.venv/lib/python3.12/site-packages/uxarray/grid/grid.py in ?(cls, dataset, use_dual, **kwargs)
    250             # parse to detect source grid spec
    251
    252             source_grid_spec = _parse_grid_type(dataset)
    253             if source_grid_spec == "Exodus":
--> 254                 grid_ds, source_dims_dict = _read_exodus(dataset)
    255             elif source_grid_spec == "Scrip":
    256                 grid_ds, source_dims_dict = _read_scrip(dataset)
    257             elif source_grid_spec == "UGRID":

~/Sources/foss/pymech/.venv/lib/python3.12/site-packages/uxarray/io/_exodus.py in ?(ext_ds)
    101     )
    102
    103     # populate lon/lat coordinates
    104     lon, lat = _xyz_to_lonlat_deg(
--> 105         ds["node_x"].values, ds["node_y"].values, ds["node_z"].values
    106     )
    107
    108     # populate dataset

~/Sources/foss/pymech/.venv/lib/python3.12/site-packages/xarray/core/dataset.py in ?(self, key)
   1583                 message = f"No variable named {key!r}. Variables on the dataset include {shorten_list_repr(list(self.variables.keys()), max_items=10)}"
   1584                 # If someone attempts `ds['foo' , 'bar']` instead of `ds[['foo', 'bar']]`
   1585                 if isinstance(key, tuple):
   1586                     message += f"\nHint: use a list to select multiple variables, for example `ds[{[d for d in key]}]`"
-> 1587                 raise KeyError(message) from e
   1588
   1589         if utils.iterable_of_hashable(key):
   1590             return self._copy_listed(key)

KeyError: "No variable named 'node_z'. Variables on the dataset include ['node_x', 'node_y', 'face_node_connectivity']"

[ashwinvis]

Yeah ! But I think extracting mesh or any mesh data from the exodus file will not provide us with the correct mesh data. The NEK5000 visualization file has the polynomial nodes data along with the nodes of the original mesh. Correct me, if I am wrong but I think this is the case only. So, in that case how to extract the connectivity data from the visualization file.

That is true. Nek5000's element numbering may be different. In that case we need its connectivity. Do you use gencon command or the genmap? If it is the latter, #40 implemented a reader for the .ma2 files. I didn't get a chance to explore it, but now is the time!

[ashwinvis]

It could also be that .ma2 file does not contain this information. If getting the connectivity out of Nek5000 is too hard, we could tap into uxarray's methods to do that. The test suite TestConnectivity suggest that it should be possible:

https://github.com/UXARRAY/uxarray/blob/bd8901f665d455cf324f9773790af118f2f48ef3/test/test_grid.py#L885-L900

[gauravity30]

Yeah ! But I think extracting mesh or any mesh data from the exodus file will not provide us with the correct mesh data. The NEK5000 visualization file has the polynomial nodes data along with the nodes of the original mesh. Correct me, if I am wrong but I think this is the case only. So, in that case how to extract the connectivity data from the visualization file.

That is true. Nek5000's element numbering may be different. In that case we need its connectivity. Do you use gencon command or the genmap? If it is the latter, #40 implemented a reader for the .ma2 files. I didn't get a chance to explore it, but now is the time!

Yeah, it is essential for simulations, so we build the .ma2. Although, I am also not sure whether it will contain the exact nodes numbering as that of the visualization file. As far as I know, it only depends on the input mesh and doesn't account for the polynomial nodes since it is generated prior to simulations.

We should work on extracting connectivity data from the .f00001 files itself. The uxarray method _build_face_face_connectivity can be helpful but no mention about it in the documentation. I will look into this. Maybe ask in the discussion section of uxarray.

[gauravity30]

@ashwinvis You probably missed my question, I am able to make an xarray dataset which can be used for most post-processing stuff like plotting, POD calculations, etc. It lacks the element connectivity data, so is it required ? If not, this method can be used to handle structured and unstructured arrays both. It is a very simple approach.

[ashwinvis]

I don't know. It depends on what you can do now. Without connectivity you may be able to plot, select a slice.

Sometimes you need to integrate or differentiate (by finite difference to roughly calculate a gradient) or interpolate along a direction. Without connectivity information that will be incredibly inconvenient.

[gauravity30]

You can get the connectivity data at least for the EXODUS mesh like this:

import xarray as xr
mesh = xr.open_dataset("naca(10x4).e")
print(mesh.data_vars)
mesh.connect1

This integer data array has the dimensions connect1 (num_el_in_blk1, num_nod_per_el1) int32 766kB ... and looks like this:

array([[    3,     2,     1, ...,     5,     8,     7],
       [   10,     9,     2, ...,    11,     6,    13],
       [   15,    14,     9, ...,    16,    12,    18],
       ...,
       [72683, 72686, 32640, ..., 72732, 32642, 72731],
       [72686, 72689, 32645, ..., 72733, 32647, 72732],
       [72689, 72692, 32650, ..., 72734, 32652, 72733]], dtype=int32)

I am guessing it denotes one cell per row, and along the column denotes the 8 cells which surround it: 4 on each corner and 4 on each face (you need to double check it). Somebody wrote a short blog post looking at this using netCDF4, but above I used xarray (which also uses netCDF4 / h5netcdf behind the scenes).

https://johnfoster.pge.utexas.edu/blog/posts/extracting-exodus-information-with-netcdf-python/

I have an idea on how to build the connection data. The each face in the exodus file is further divided by polynomial nodes. So, each edge is a line with the polynomial as well as geometric nodes. We can get the slope of the edge and find which polynomial nodes lie on these edge and sort them based on the distance. This way we can find out how the nodes in each element object is form the face. This will work provided we can get connectivity of geometric nodes from exodus file.

Rough diagram of a 2D element (can be extended to 3D)

Red nodes are the polynomial nodes that we get from the .f00001 files and blue nodes are the geometric nodes.

[ashwinvis]

@airwarriorg91 it has been a while. Did you work more on this?

I know it is a hard problem. It will be nice if we can try to complete this or at least document what the next steps are. Just so that we don't forget.

[gauravity30]

@airwarriorg91 it has been a while. Did you work more on this?

I know it is a hard problem. It will be nice if we can try to complete this or at least document what the next steps are. Just so that we don't forget.

Hey !
I was a little busy with the semester end. I need to try the above idea and check its feasibility. I will probably start working on it again after Dec 15 and let you know.

Thanks,
Gaurav

[ashwinvis]

Ok sounds good! It would make for a nice xmas project 🎄

[gauravity30]

Hello Ashwin !
Sorry for being away for a long time, but now I finally have time on weekends to work on these problems. My college is finished and have joined job, all the free time will be dedicated to solving this problem. I will revisit the problem with a fresh mind.

Thanks !

[gauravity30]

Hey !
I did a revision of the problem cause its been a long time, I also checked for update in the uxarray method _build_face_face_connectivity, the PRs are closed but in draft (unmerged), also it attempts to create a grid from points on sphere. It might be not useful for us. But my idea to build the connectivity data from the exodus file is hopeful. Although, instead of using the slope based method, I might just the idea that the polynomial nodes are GLL nodes and arrange them (Less computationally expensive).

This is how I plan to proceed:

• Read the exodus data and plot it.
• Find the polynomial order data automatically from the visualization file.
• Implement a method to estimate the GLL nodes on each edge, find them in the nodes data and arrange them.
• Create a connectivity data from them.
• Create a ugrid from this for further use.

Although I might need help in optimizing the code for speed and memory.

Looking forward to completing this PR.

Thanks,
Gaurav

[ashwinvis]

Awesome @gauravity30 and thank you for finding time for this. I presume you are talking about this:

UXARRAY/uxarray#1195

Let me give you a better base to start from so that the all the current tests pass. I will also try to add a single test suite which based on the NACA airfoil data you send long time back (I think).

[gauravity30]

Most welcome @ashwinvis ! I was talking about this UXARRAY/uxarray#349

Please, let me know after you have finished your part. I can see the checks are failing or cancelled and yes, I had sent the naca airfoil test data a year ago, when we started working on this.

[ashwinvis]

@gauravity30 I have added some basic changes to the functions and a small test. Could you do this in your branch?

git remote add upstream https://github.com/eX-Mech/pymech/
git fetch upstream pr-126-support-unstruct-grid
git merge upstream pr-126-support-unstruct-grid

[gauravity30]

@gauravity30 I have added some basic changes to the functions and a small test. Could you do this in your branch?

git remote add upstream https://github.com/eX-Mech/pymech/
git fetch upstream pr-126-support-unstruct-grid
git merge upstream pr-126-support-unstruct-grid

@ashwinvis I was trying it in my fork of PyMech but it was failling. I saw you created a new branch in this repo. Should I also create a branch with same name and try it ? I have not worked with branches in github before so a little confused.

Also my fork is outdated, need to merge the commits which happened after creation of this PR.

EDIT: I managed to resync my fork with the branch you created. I think i can start the work forward.