Using chmy.jl for global ocean grids

[briochemc]

Saw the JuliaCon Local Paris 2025 talk! What a cool package! Congrats! I have two questions:

1. If I have some grid files from an ocean model (that define all the geometries of a staggered grid, including all nodes and wet/dry cells), can I "load" that grid in Chmy.jl to build all the finite-difference operators on it? Or is there some fundamental difference between how Chmy.jl does its finite differencing and how ocean models work?

2. Also if I heard correctly there was some mention of sparsity patterns in the talk, but I'm not sure what was meant. Could you elaborate?

---

Context of my interest: I build what is called "Transport Matrices" from ocean model outputs (https://github.com/TMIP-code/OceanTransportMatrixBuilder.jl), which allow for approximate but fast offline simulations. To build these matrices, I do the spatial finite differencing I need "by hand", potentially with some mistakes, so I would love to rely on some more robust Julia package to do that for me. And Chmy.jl looks like the thing I need. But the question is: Is it ready for my use case?

[utkinis]

Hi @briochemc ! Thank you for the interest in the package. Please note that the v0.2 is in the heavy development, no PR even opened yet 😅 The current doc isn't representative of what features will be in the overhauled package. I plan to open a PR in the coming weeks. To answer your questions:

1. This depends on how the grid is arranged. I'm not familiar with ocean models, but we'll support standard Arakawa staggered grids, and it should be easy to construct custom staggered arrangement. More subtle issue is the topology of the grid. As of now, we only support topologically rectangular grids, so if your grid is something like a "cubed sphere", it might take some time to implement this functionality, and I don't have a clear plan on how to do this yet. Other than that, if your grid is topologically a cube, it should be easy to load this in Chmy.jl.

2. Constructing finite difference operators is also supported, with the only limitation that they should be local, i.e. implementing certain compact/mimetic FD schemes could be impossible. Otherwise, you can create a 1D derivative, and then generate nD partial derivative from it, higher-order derivatives by repeated application, flux reconstructions etc.

3. Chmy v0.2 can infer the stencil from the finite-difference operator. From that stencil it can construct a sparsity pattern of the Jacobian of this operator. One could then use graph colouring algorithms and automatic differentiation (AD) to compute sparse Jacobians from FD operators to solve e.g. Navier-Stokes equations. The colouring and AD are not part of Chmy.jl, but there several packages in Julia ecosystem to perform this steps (see SparseMatrixColorings.jl, Enzyme.jl).

In short, please stay tuned for the upcoming release 😉

[briochemc]

On ocean model grids: Most use standard Arakawa grids (although a few are unstructured, like MPAS-O). In terms of topology, most are regular (logical rectangles), although for global grids they need to wrap around the globe (so periodic in the longitude dimension; related #42). In addition, to avoid singularities (vanishing cells) in the Arctic, most either displace the North Pole onto land (over Greenland), or "fold" in the Artic with 2 displaced poles (so there is a "seam" on the "logical north" boundary). In both cases the grid is then curvilinear (but still logically rectangular) Maybe this image is clearer than my words:

Figure copied from https://gmd.copernicus.org/articles/17/4603/2024/. Caption: Examples of two curvilinear grids typically used for ocean modeling. (a) A tripolar grid with two northern poles (one over Eurasia and one over North America). (b) A displaced dipolar grid with one northern pole displaced over Greenland.

And then of course land points are out of the domain. Hierarchical grids, like cubed spheres are not that common as far as I know (e.g., MITgcm uses that).

About the operator: Could I use Chmy.jl to construct the full operator itself rather than its sparsity pattern? That's exactly what I do BTW: I use face-fluxes/velocities to construct the full operator of the flux divergence (which is what I call the "Transport Matrix").