Grid representations

[wpbonelli]

Some notes on representing grids.

DISU

A DISU grid is a directed graph. There are several options to represent these.

Currently MODFLOW 6 uses an adjacency list, where ja specifies adjacent cells and ihc indicates direction.

Flopy currently follows MF6 in using the adjacency list representation. It computes neighbors lazily and caches them to speed up subsequent searches.

@deltamarnix has looked into interop with networkx in #58, which provides an adjacency matrix representation — storing these in dense numpy arrays and sparse scipy arrays, as well as in xarray. It seems we could manually store a sparse array in an xarray dataset in CSR format, or use a sparse array since xarray supports these and they can be converted to/from scipy sparse arrays.

[wpbonelli]

reproducing some things said by @deltamarnix in #49 about #58

I've researched some of the packages that seem to support UGRID. It's hard to find a library that supports 3D unstructured grids. I found the following so far:

• UGridPy, a Deltares product. Only supports 1D and 2D.
• gridded is a package that only supports 2D grids (lat and lon). And it seems to store everything in memory.
• XUGrid, a Deltares product by Huite. Only supports 2D, and can extend with adding an extra XArray dimension for layer in order to support DISV.
• UXArray, currently being worked on by a masters student at NCAR. Supports 2D with XArray, so would likely also be able to add an extra dimension layer for DISV.

I haven't found any DISU compatible libraries. And I was thinking that it might be wise to keep the FloPy design as close to the MODFLOW6 input format as possible, like it is right now. So perhaps it's best as @mjreno suggested, to use numpy-like arrays so we can support dask where possible as input formats. Preferable as lazy arrays. And the moment that we need to write the MF6 files, we compute() the values and write them directly to the file format that is needed.

I will try to make a little prototype that shows how that works.

---

I'm still playing in my head with the fact that DISU is actually just a graph, and should be treated that way. The graph itself might fit into memory, but the data on the nodes and edges might not.

There area some core packages in python that represent graphs: python graphblas, networkx, scipy. And on top of that there is a package that can translate graphs from one to another, or call algorithms from one on another type: metagraph (metagraph was last released in 2021, so not as interesting anymore).

So far we want to stay as close as possible with the packages to the input of modflow, but we might want to use these packages in generator helper functions networkx -> gwfdisu. Or we might accept a networkx-graph as input if it is set up correctly and convert it to our own representation internally.

For the output we can do something similar. By default we return a single list of data points based on the cell index, but if there are use cases for traversal, we can add a translation layer on top to return a networkx or graphblas object with the data linked to the correct nodes and edges.

---

I'm diving a little deeper into what XArray supports for representing a graph (DISU) in its data format. [XArray](https://docs.xarray.dev/en/latest/user-guide/duckarrays.html) supports [numpy sparse arrays](https://sparse.pydata.org/en/stable/).

NetworkX has the option to transform its representation to numpy arrays to_numpy_array. This information can be directly stored in an XArray object. And the other way around, a network representation in XArray can be transformed back to a NetworkX object with from_numpy_array.

---

I made a little prototype in #58 to show how we could store networkx graph via a sparse scipy matrix into an XArray Dataset, and get the structure back into networkx too.

[wpbonelli]

Huite recommends we represent DISU with flat arrays in the same way as mf6 for now. We can add cell id as an xarray index for some convenience. Over the coming months the fate of xugrid vs uxarray will become clearer, at which point we can decide which to adopt for DISV, and at that point maybe reconsider options for DISU.

[wpbonelli]

prior art: grid classes in xgcm, xugrid, uxarray

• https://xgcm.readthedocs.io/en/latest/grids.html
• https://uxarray.readthedocs.io/en/latest/user-guide/data-structures.html#grid
• https://deltares.github.io/xugrid/api/xugrid.Ugrid2d.html

[wpbonelli]

Some conversions we might want to do with a grid class:

• Create a package from a grid. When defining models programmatically, we like the existing pattern/recommendation to make a grid first, then use it to create package/model/simulation.
• Create an xarray dataset from a grid.
• Create a grid from a package or xarray dataset. This lets us e.g. load a DIS package from a netCDF file and infer a grid from it.

[wpbonelli]

Component classes would presumably continue to have a .modelgrid or .grid property. Datasets could also have something like that, either via subclassing (as uxarray and xugrid do) or accessor registration. While you can save state on accessors, it seems like we might want to keep the state (i.e. grid) in an index and consume it from an accessor. This also lets us build a grid-aware index. This is the path Huite seems to want to eventually take with xugrid.

[wpbonelli]

In the near term, maybe we can start with a simple index which stores a flopy.discretization.grid.Grid and does some basic things like k/i/j indexing (#108). In the mid term, we could prototype CRS-awareness and CF (or other) conventions. Whenever we know how to proceed with xugrid/uxarray (perhaps they will adopt the custom index approach) we can consider what a full solution might look like.

[wpbonelli]

This project is a good example of how to use the accessor extension pattern.

[wpbonelli]

I played around with the index mechanisms and got some guidance from an xarray developer pydata/xarray#10257

[wpbonelli]

an xarray development to watch pydata/xarray#1475.. cell boundaries as coordinates