Traits

[rafaqz]

@gaelforget this is an example of traits for the interface

Theyre not meant to be representative of the actual traits, just an example of the structure that we can use for more appropriate things

# How are the polygons ordered?
abstract type GridOrder end
abstract type OrderedGrid <: GridOrder end

struct NorthSouth <: OrderedGrid end
struct SouthNorth <: OrderedGrid end
struct Unordered <: GridOrder end

gridorder(x) = Unordered()

gridorder(x::MyGrid) = Unordered()
gridorder(x::MyOtherGrid) = NorthSouth()

gridalign(source, dest) = gridalign(gridorder(source), gridorder(dest), source, dest)

function gridalign(::Unordered, ::Unordered, source, dest)
    # algorithm here
end
function gridalign(::OrderedGrid, ::Unordered, source, dest)
    # algorithm here
end
function gridalign(::Unordered, ::OrderedGrid, source, dest)
    # algorithm here
end
function gridalign(::OrderedGrid, ::OrderedGrid, source, dest)
    # algorithm here
end

maybe_flip_north(::NorthSouth, grid) = grid
maybe_flip_north(::SouthNorth, grid) = reverse(grid)

# Are the polygons all quadrilaterals?
struct Quad end
struct NotQuad end

isquad(::MyGrid) = Quad()
isquad(::MyOtherGrid) = NotQuad()

[asinghvi17]

cc @milankl @simone-silvestri

[milankl]

struct NorthSouth <: OrderedGrid end
struct SouthNorth <: OrderedGrid end
struct Unordered <: GridOrder end

This is not a problem to define if helpful for any algorithms here. Oceananigan's grids are "[east:west, south:north]" (in column-major), I believe, whereas SpeedyWeather uses "[east:west, north:south]" so a reverse could do it, but would also swap east:west to west:east. Grids like HEALPix also have a nested ordering e.g. Gorski et al. 2005

SpeedyWeather uses consistently the RING SCHEME (top panel) for all grids. But I reckon the cubed sphere would also likely employ something similar to the NESTED SCHEME outlined here. As far as I understand it, making assumptions about the order could provide performance potential as only a smaller subset of polygons would need to be intersected as you know that all others are far away and therefore wouldn't intersect anymore.

[felixcremer]

@danlooo might be interested in this from the DGGS.jl perspective.

[milankl]

# Are the polygons all quadrilaterals?
struct Quad end
struct NotQuad end

isquad(::MyGrid) = Quad()
isquad(::MyOtherGrid) = NotQuad()

Yes, all grids of SpeedyWeather and Oceananigans consist of only quadrilaterals. But I was wondering whether for the sake of generality we should just do a nvertices(::MyGrid) = 4 so that every grid can define how many vertices they have. We could say we only support grids which have a same number of vertices for all cells. I guess that would allow one to do a lot with matrices instead of vectors of vectors?

[gaelforget]

also maybe https://github.com/eth-cscs/ImplicitGlobalGrid.jl devs might be interested

[asinghvi17]

About nvertices: I was thinking that the getpolygon methods can just return a tuple or static vector of points, then the number of points is statically encoded. If the grid has multiple face types (some pentagons and some hexagons, for example) then the nvertices abstraction breaks down...

Also, matrices probably don't work in the general case. We can always allow high level methods to override for performance optimizations, or create new traits to describe whether the grids are MatrixRepresentable or VectorRepresentable.

[danlooo]

Yes, one can not tessellate the globe using only hexagons, there will always be 12 pentagons, requiring grids of mixed shape (e.g. ISEA4H discrete global grid system or Uber H3).
Furthermore, some grid indices (e.g, DGGRID Q2DI) are represented by not a single matrix but a vector of matrices. Hereby, the global grid is based on a icosahedron in which 2 triangular faces are put together to one matrix.

[asinghvi17]

For the multi-matrix case, I was thinking we could use a flattened vector to store values, with some accelerator index based on which "face set" the series of indices is in. Kind of like spatial chunking :D - Gael has a similar structure I believe, though Simone's grid (ClimaOcean) is simply representable as a matrix.

But GridInterface probably doesn't need to know about that, beyond the implementing methods having some acceleration.