Grid Struct Template start for MOLE 2.0

[jbrzensk]

For MOLE 2.0, we are exploring implementing structs to store key information about the grid and boundary conditions. This simplifies calls to the operators, but also simplifies request by the user.

Many of the grid objects are generated automatically, based on an initial call to the grid. The grid struct stores all given and derived information for the operators. Below is a sample of what the grid struct could be. It is a work in progress.

---

MOLE Grid Struct Reference (MATLAB/Octave)

---

1. Overview

The grid struct is the central data object passed to all MOLE mimetic operators
(grad, div, lap, curl, nodal, addScalarBC, interpol). It bundles
everything an operator needs to know about the physical domain: how many cells
exist in each direction, how large each cell is, where nodes/faces/cell-centers
sit in physical space, and what type of boundary treatment is expected. (It is expected that interpol will fill in fields in the grid struct, like facesToNode, NodesToFaces, etc.).

Two public functions build and validate the struct:

grid = makeGrid('m', 64, 'n', 64, 'dx', 1/64, 'dy', 1/64, 'bc', bc)
grid = validateGrid(grid)   % explicit re-validation / enrichment (fill-in-missing)

makeGrid accepts either name-value pairs or a partially filled struct and
always calls validateGrid before returning. validateGrid infers missing
metadata fields, normalises boundary coefficients, and populates all coordinate
arrays. Callers should treat the returned struct as immutable: operators read
from it but never write to it.

---

2. Core Fields

2.1 dim — dimensionality

Value	Meaning
1	1-D problem (x only)
2	2-D problem (x, y)
3	3-D problem (x, y, z)

Inferred automatically from whichever size fields are present:

• o or dz present → dim = 3
• n or dy present → dim = 2
• otherwise → dim = 1

2.2 type — grid topology

Value	Meaning
'uniform'	Constant spacing (dx/dy/dz scalars); no periodic axes
'periodic'	Constant spacing; one or more axes are periodic (set by bc?)
'curvilinear'	Physical coordinates supplied by the caller as node arrays

Inferred automatically:

• X, Y, or Z fields present at struct root → 'curvilinear'
• x, y, or z fields (lowercase) present → 'nonuniform' (reserved, not yet active)
• otherwise → 'uniform' or 'periodic' depending on boundary coefficients

2.3 Cell counts: m, n, o

Field	Axis	Required for
m	x	1-D, 2-D, 3-D
n	y	2-D, 3-D
o	z	3-D

These are cell counts, not node counts. A grid with m cells in x has
m+1 nodes along x.

For curvilinear grids, m/n/o must equal the node-array dimension minus one:
nodes.X must be (m+1) × (n+1) for 2-D, (m+1) × (n+1) × (o+1) for 3-D.

2.4 Cell spacings: dx, dy, dz

Scalar values give the uniform physical spacing between adjacent nodes along
each axis. Infered if passed x,y,z grid values. Required for uniform and periodic grids; not used for
curvilinear grids (physical distances are encoded in the node coordinate
arrays instead).

2.5 bc — boundary condition sub-struct

The bc sub-struct encodes Robin boundary coefficients and periodicity flags.
See the companion Boundary Conditions Reference document for full details.
The fields relevant to the grid struct itself are:

• bc.isPeriodic — logical vector (1×1 for 1-D, 2×1 for 2-D, 3×1 for 3-D).
  Each element is true if the corresponding axis wraps around. Used by
  validateGrid to set grid.type = 'periodic'.

---

3. Coordinate Arrays

After validateGrid, the struct contains three sub-structs holding physical
coordinates for every point class in the staggered mimetic grid:

Sub-struct	Point class	Role
nodes	Grid corners	Where physical boundaries and node values live
faces	Cell-face midpoints	Where flux (vector field) unknowns live
centers	Cell interiors	Where scalar (potential) unknowns live

All arrays use ndgrid layout: the first index varies in x, the second in y,
the third in z. This is the opposite of MATLAB's default meshgrid convention.

3.1 One-dimensional grid

Field	Size	Description
nodes.X	(m+1) × 1	0, dx, 2·dx, …, m·dx
faces.X	(m+1) × 1	Identical to nodes.X (in 1-D, face midpoints coincide with nodes)
centers.X	(m+2) × 1	0, 0.5·dx, 1.5·dx, …, (m−0.5)·dx, m·dx

The centers array has two ghost entries: the first element (0) and the
last element (m·dx) mirror the domain boundary. This is required so that
the mimetic boundary condition operator can augment the interior matrix rows
with boundary data without resizing.

3.2 Two-dimensional grid

Face arrays are split by normal direction — u faces are perpendicular to x
(x-normal), v faces are perpendicular to y (y-normal).

Field	Size	Description
nodes.X, nodes.Y	(m+1) × (n+1)	ndgrid of node positions
centers.X, centers.Y	(m+2) × (n+2)	ndgrid with ghost padding on all sides
faces.u.X, faces.u.Y	(m+1) × n	x-normal face midpoints
faces.v.X, faces.v.Y	m × (n+1)	y-normal face midpoints

u-face coordinate ranges (uniform):

faces.u.X = 0, dx, 2·dx, …, m·dx          (m+1 values — same as node x-positions)
faces.u.Y = 0.5·dy, 1.5·dy, …, (n−0.5)·dy (n values — face centers between node rows)

v-face coordinate ranges (uniform):

faces.v.X = 0.5·dx, 1.5·dx, …, (m−0.5)·dx (m values)
faces.v.Y = 0, dy, 2·dy, …, n·dy           (n+1 values — same as node y-positions)

3.3 Three-dimensional grid

Face arrays split into three components: u (x-normal), v (y-normal),
w (z-normal).

Field	Size	Description
nodes.X/Y/Z	(m+1) × (n+1) × (o+1)	ndgrid of node positions
centers.X/Y/Z	(m+2) × (n+2) × (o+2)	ghost-padded ndgrid
faces.u.X/Y/Z	(m+1) × n × o	x-normal faces
faces.v.X/Y/Z	m × (n+1) × o	y-normal faces
faces.w.X/Y/Z	m × n × (o+1)	z-normal faces

The pattern is consistent: the face component named after a direction (u→x,
v→y, w→z) has a full node count in that direction and cell-center count in
the other two directions.

---

4. Curvilinear Grids

For non-uniform physical domains the caller constructs a partially filled
struct, populates the node coordinate arrays, then passes it to validateGrid:

g = struct();
g.m = m - 1;   % cell count = node count - 1
g.n = n - 1;
g.type = 'curvilinear';
g.nodes.X = X';   % (m+1) × (n+1) node array, transposed to ndgrid layout
g.nodes.Y = Y';
g = validateGrid(g);

4.1 What the caller must supply

Dimension	Required fields
2-D	m, n, type='curvilinear', nodes.X, nodes.Y
3-D	m, n, o, type='curvilinear', nodes.X, nodes.Y, nodes.Z

Node arrays must be exact size (m+1) × (n+1) (2-D) or
(m+1) × (n+1) × (o+1) (3-D) in ndgrid layout (x-index first).
validateGrid throws validateGrid:SizeMismatch if the sizes are wrong, and
validateGrid:CurvilinearMissingNodes if the arrays are absent.

4.2 What validateGrid derives

All coordinate arrays are computed by averaging the caller-supplied node
positions. No dx/dy/dz is computed or stored; physical distances are
implicit in the node positions themselves.

u-faces (x-normal): average of adjacent nodes along the n-direction (column pairs):

faces.u.X(:, j) = 0.5 * (nodes.X(:, j) + nodes.X(:, j+1))   for j = 1…n
faces.u.Y(:, j) = 0.5 * (nodes.Y(:, j) + nodes.Y(:, j+1))

Size: (m+1) × n.

v-faces (y-normal): average of adjacent nodes along the m-direction (row pairs):

faces.v.X(i, :) = 0.5 * (nodes.X(i, :) + nodes.X(i+1, :))   for i = 1…m
faces.v.Y(i, :) = 0.5 * (nodes.Y(i, :) + nodes.Y(i+1, :))

Size: m × (n+1).

Cell centers: bilinear average of the four surrounding corner nodes:

centers.X(i, j) = 0.25 * (nodes.X(i,j) + nodes.X(i+1,j) + nodes.X(i,j+1) + nodes.X(i+1,j+1))
centers.Y(i, j) = 0.25 * (nodes.Y(i,j) + nodes.Y(i+1,j) + nodes.Y(i,j+1) + nodes.Y(i+1,j+1))

Size: m × n.

For 3-D curvilinear, the same edge-midpoint and trilinear-average rules
apply in all three axis directions. The node arrays are
(m+1) × (n+1) × (o+1) and the derived face/center arrays follow the same
naming and sizing conventions as the uniform 3-D case above.

4.3 Layout note for operator internals

Curvilinear operators (gradCurv_impl, divCurv_impl, nodalCurv_impl)
internally need meshgrid layout for the Jacobian computation. They
transpose or permute the ndgrid arrays on entry:

X = grid.nodes.X';            % 2-D: (n+1)×(m+1) meshgrid
X = permute(grid.nodes.X, [2,1,3]);  % 3-D: (n+1)×(m+1)×(o+1)

This conversion is the operator's responsibility, not the grid struct's.
The struct always stores coordinates in ndgrid layout.

---

5. Quick-Reference Table

Field	1-D	2-D	3-D	Curvilinear	Notes
dim	1	2	3	2 or 3	Auto-inferred
type	str	str	str	'curvilinear'	'uniform'/'periodic'/'curvilinear'
m	scalar	scalar	scalar	scalar	Cell count in x
n	—	scalar	scalar	scalar	Cell count in y
o	—	—	scalar	scalar	Cell count in z
dx	scalar	scalar	scalar	—	x cell spacing
dy	—	scalar	scalar	—	y cell spacing
dz	—	—	scalar	—	z cell spacing
nodes.X	(m+1)×1	(m+1)×(n+1)	(m+1)×(n+1)×(o+1)	same (user-supplied)	ndgrid layout
nodes.Y	—	(m+1)×(n+1)	(m+1)×(n+1)×(o+1)	same (user-supplied)
nodes.Z	—	—	(m+1)×(n+1)×(o+1)	same (user-supplied)
centers.X	(m+2)×1	(m+2)×(n+2)	(m+2)×(n+2)×(o+2)	m×n or m×n×o	Staggered grid for MOLE scalars
centers.Y	—	(m+2)×(n+2)	(m+2)×(n+2)×(o+2)	m×n or m×n×o
centers.Z	—	—	(m+2)×(n+2)×(o+2)	m×n×o
faces.X	(m+1)×1	—	—	—	1-D only; equals nodes.X
faces.u.X	—	(m+1)×n	(m+1)×n×o	(m+1)×n	x-normal faces
faces.u.Y	—	(m+1)×n	(m+1)×n×o	(m+1)×n
faces.u.Z	—	—	(m+1)×n×o	—	3-D only
faces.v.X	—	m×(n+1)	m×(n+1)×o	m×(n+1)	y-normal faces
faces.v.Y	—	m×(n+1)	m×(n+1)×o	m×(n+1)
faces.v.Z	—	—	m×(n+1)×o	—	3-D only
faces.w.X	—	—	m×n×(o+1)	—	z-normal faces
faces.w.Y	—	—	m×n×(o+1)	—
faces.w.Z	—	—	m×n×(o+1)	—
bc.isPeriodic	1×1	2×1	3×1	—	One flag per axis

Keep-open request

• I am requesting maintainer review for keep-open.

[jbrzensk]

As an example of what I have working, here is the code for the grid on elliptic1DPeriodicBC:

k = 2;
m = 20; 
dx = 1/m;
xn = (0:dx:1)';

grid = struct();
grid.m = m;
grid.dx = dx;
grid.nodes.X = xn;
grid.bc = struct('dc', [0;0], 'nc', [0;0]);  % periodic ends
grid = validateGrid(grid,true);

Here the boundary conditions are coded as a struct for the grid, so the grid knows it is periodic. validategrid.m checks the grid for consistency, and generates any missing information. Note here I only gave the nodal grid.

>> disp(grid)
           m: 20
          dx: 0.0500
       nodes: [1×1 struct]
          bc: [1×1 struct]
         dim: 1
        type: 'periodic'
           x: [21×1 double]
       faces: [1×1 struct]
     centers: [1×1 struct]
    periodic: [1×1 struct]

The function has filled in faces, centers, dim, type, and periodic. The structs are structs to hold the coordinates of the dimensions. Fo 1D they are 1x1, but for 2D they would be 1x2 and 3D 1x3.

>> disp(grid.nodes)
    X: [21×1 double]

>> disp(grid.faces)
    X: [21×1 double]

>> disp(grid.centers)
    X: [22×1 double]

And the periodic is there because those grids, currently, are different, and you need these to calculate the b for the RHS.

>> disp(grid.periodic)
      nodes: [1×1 struct]
      faces: [1×1 struct]
    centers: [1×1 struct]

>> disp(grid.periodic.nodes)
    X: [20×1 double]

>> disp(grid.periodic.faces)
    X: [20×1 double]

>> disp(grid.periodic.centers)
    X: [20×1 double]

[mdumett]

My two cents on this:

I hope to convey here why square operators are necessary for PDEs with periodic boundary conditions. You will find several argument about this in what follows. The most important aspect of these square operators is not the name. It is that they are able to properly solve problems that the non-square operators cannot. These square operators not only include gradients and divergences but also interpolations and quadratures that make them mimetic.

One take out of what follows is that the new basic gradient, and divergence should be represented by square matrices (i.e., for a domain without boundaries or an infinite domain; periodic boundary conditions imply an infinite domain, but one represents only one period in the discrete domain). In case there are finite boundaries, one should extend the square gradients and divergences to non-square operators in order to include boundaries, creating additional rows or columns to deal with the boundary conditions for finite domains.

A secondary topic is what would be a convenient name these square divergences and gradients. Probably "basic operators" or "base operators".

Having "basic operators" which are modified internally to include finite boundaries will hide the square and non-square versions from the user. The user should know about the number of cells in each direction (m,n,o), and not to know that there are (m+2,n+2,o+2) centers and (m+1,n+1,o+1) faces. Nevertheless, when an "advanced" user is debugging (which will occur), he/she will discover the different internal representations of divergences and gradients.

The following 1D examples exhibit that non-square gradients and divergences are not adequate but square ones are:

Consider the equation du/dx = f with periodic boundary conditions.
a) If one uses Gu = f, with rectangular G, then the first and last component of u have to be the same. How to enforce that without duplicating that calculation (one for the left vertex and one for the right vertex) and destroying the calculation of the derivatives at the boundaries?

b) If one uses Du = f, with rectangular D, then the boundary first and last rows in D are set to zero. Enforcing the boundary condition will duplicate the calculation (one for the left vertex and one for the right vertex).

Consider the PDE d3u/dx3 = f with periodic boundary conditions.
In this case, periodic boundary conditions implies continuity of u, of du/dn, and of a certain second-derivative. These are three conditions and there are only two rows to enforce them. How could this be done?

Consider the PDE d4u/dx =f with periodic boundary conditionsIn this case, there are four conditions related to the continuity of u and its normal derivative. There are only two rows to enforce them. There exists a trick to convert this into a problem where two Poisson equations are solved (with appropriate boundary conditions for each them), but the trick will not work for a general elliptic operator of fourth-order.

Square operators can handle Ben's example (see the figure Ben added before in this conversation) and rectangular ones cannot.

All these issues are solved if one chop off one of the input vertices. Moreover,
a) the expected self-adjoint properties for the interior of D and G are automatically satisfied by the operators as a whole (D = -G'),
b) the non-mimetic usual property of symmetry of L is attained (L' = L), and
c) the symmetry between the interpolation from centers to faces (Icf) and to the interpolation from faces to centers (Ifc) holds, i.e., Icf = Ifc'.

These operators have been used successfully for solving PDEs with anti-symmetric structure and with 2D and 3D curls in the context of periodic boundary conditions.

Enforcing periodic boundary conditions to the rectangular operators has its limitations.

Some history:
a) I spent months thinking about this while implementing the addScalerBC feature and finally figuring out the convenience of having square matrices in the presence of boundary condition (see months of conversations with Angel and Jared before approving the addScalarBC PR).
b) Initially, I tried to impose periodicity on the boundary of the non-square gradient and divergence but I discovered that there were situations where that approach would not work.
c) I ended up changing the non-square operators for the case of periodic boundary conditions for somme of the reasons exposed above.

[Tony-Drummond]

Hi @jbrzensk ,
Thanks for sharing this design, it looks almost complete and it has sufficient information for the language implementations of the operators.
It actually also provides an answer to @mdumett's above, since the matricial representation of the operators are going to be determined by the grid.type, and when the grid.type is periodic, we will be building the correct square operators.

Your example building a grid on elliptic1DPeriodicBC, was very useful to understand the abstractions that you described in the previous sections. However, we have a few questions for you:

1. When building the example, you did not use the makeGrid but instead used the equivalent of a constructor-like for the Grid struct, my understanding is that this constructor is the one that determines the grid.type property (?) Can a user also set this attribute?

2. we are also operating under the assumption that 'uniform' = nonperiodic+uniform, correct? and also that 'periodic' = periodic+uniform, correct?

3. The nonuniform grids will be an extension that will be implemented in the future (?), are these a completely different grid struct (e.g. NUGrid struct)? or an extension of this proposed struct?

Thanks,

[Tony-Drummond]

One more question, in our current abstract design of the operators, we had considering these calls:
G = grad(k, grid, bc), G gradient
D = div(k, grid, bc), D divergence
L = laplacian(k, grid, bc), Laplacian
C = curl(k, grid, bc), C curl (location (i.e., center, edges) defined by values in grid
[J, T] = jacobian(k, grid), J = Jacobian, T is a set of tensors

Should we remove bc from the argument lists since bc is now a substruct of grid?

[mdumett]

I would say grid.bctype instead of grid.bc, since dc and nc are for choosing the bc type.

When considering a general non-periodic boundary condition, we will have something like this: a(x) u + b(x) gradient u = g(x), where x is a boundary point, a(x) is a vector containing the values of a(x) at each boundary point, b(x) is a full 1D, 2D, or 3D square matrix (possibly a multiple of the diagonal, a diagonal, a symmetric positive definite matrix, or a general matrix at each boundary point), and g(x) is a function that varies point-to-point

[mdumett]

We could remove the bc from the arguments since grid will have a bc type field.
However, curl operators will have a different way of handling bc, depending one how are they constructed (in terms of divergence for example), and where the input and out are.

I was thinking probably each vector calculus operator should have an input grid (where the input field data is) and an output grid (where are the calculation being done) as arguments. This way the user just provides where he wants the input and output. This could also allow the user to have their own gradient and divergence implementations.

The best way to facilitate that might be a fully-staggered implementation of all vector calculus operators where all vector calculus operators use data from all meshes (centers and faces) and their output is everywhere too (center and faces). If we internally calculate everywhere then the output could be restricted at the locations where the user wants it.