Better grids, including better RegularCartesianGrid

[glwagner]

I see a few things to improve with the grid right now.

1. We need to add some kind of type information that indicates whether a given dimension is 'flat'. One way to do this could simply be to add the grid sizes in each dimension as parameters in the abstract type Grid, aka:

abstract type Grid{T, Nx, Ny, Nz} end

struct RegularCartesianGrid{T, R, Nx, Ny, Nz} <: Grid{T, Nx, Ny, Nz}
    ...
end

Functions can then dispatch when one of Nx, etc is 1 (including halo-filling functions, which I think may fail when the size of the halo is 0). Another option is to use flags for each dimension rather than the actual size of the grid. In my opinion using the size makes the most sense. Using the actual size could have future advantages; for example, if some optimizations are possible when Nz=2. It is also nice to see the size of the grid from the type signature. A disadvantage is that we then could not have mutable grid types, but I'm not sure we want that.

2. There is a lot of redundant information in the RegularCartesianGrid struct: cell areas, volumes, total number of grid points, etc. I think it would be better --- meaning that our code would be shorter, easier to read, easier to maintain, easier to reason about (since storing them implies they cannot be computed, which is incorrect) and more computationally efficient --- to add functions that compute these quantities on the fly, rather than storing them in memory.

Related: #287.

[ali-ramadhan]

Just wanted to point out that PR #283 addresses your second point.

https://github.com/climate-machine/Oceananigans.jl/blob/2176be86be80fb4b9de2c3051c4ea09c74ed00c9/src/operators/areas_and_volumes.jl#L1-L38

[glwagner]

@ali-ramadhan this looks nice. I think we should also change Grid to AbstractGrid (complies with YASGuide and makes me feel more comfortable). Also functions like

@inline Δx(i, j, k, grid::Grid) = @inbounds grid.Δx[i, j, k] 

that cannot be used probably should not be part of the code.

[ali-ramadhan]

Originally I was against types like AbstractGrid and still am, but it seems to be what the Julia ecosystem uses so would be good to adopt to avoid surprises.

I agree . That PR should only introduce operators that will be used (otherwise they can't be tested so what's the point).

[glwagner]

I have a revision to my original post.

Rather than add the size of the grid as types in the RegularCartesianGrid signature, we perhaps should come up with a Topology type indicator, eg:

abstract type AbstractTopology end

struct Bounded <: AbstractTopology end
struct Flat <: AbstractTopology end
struct Periodic <: AbstractTopology end

AbstractGrid then becomes

abstract type AbstractGrid{FT, Tx, Ty, Tz} end

where (Tx, Ty, Tz) denote the "topology" of the grid in their respect directions.

For Flat grids, differentiation will return 0, and interpolation will be reduced to an identity function. Distinguishing between Periodic and Bounded topologies will be useful for the implementation of higher-order operators that need to limit to a second-order derivative near non-periodic boundaries (this is important for avoiding the dreaded "halo filling explosion").

I think this also corresponds more closely to the mathematical meaning of a "Periodic" and "Flat" dimension.

Another thought is that we should make the constructor more verbose, eg

RegularCartesianGrid(size=(Nx, Ny, Nz), length=(Lx, Ly, Lz), topology=(Tx, Ty, Tz)).