GOcean1.0 stencil metadata

[rupertford]

For the GOcean api to support OmpSs and MPI distributed memory we need to add stencil metadata to the GOcean1.0 API.

I propose this is done in two stages

1) define appropriate names to describe stencils and the metadata changes 2) add parsing support for the new metadata (+ tests, examples etc).

[rupertford]

In terms of finite difference stencil descriptions we have ...

2D single depth 5-point and 9-point stencils.

forward and backward finite difference.

We could describe (depth 1) stencils in terms of N, S, E, W, NE, SE, NW, SW

We could describe stencils as 1's and 0's e.g. 010,111,010 for a 5 point stencil.

[hiker]

I wonder if it would make sense to declare the stencil as a string? This way we would not require additional types in GOcean, and it would be intuitive names without polluting the namespace with "N", "SE". Example:

          arg(WRITE,      CT, "N-S")   ! North/South Halo

Other examples: "N-S-W-E" (5 point), "NN-SS-WW-EE" (9 point) or "NN-SS-W-E" (not sure if this makes sense in reality, but it would be a (strange) 7 point stencil).

Advantage: easier to understand than the binary flags; avoids name pollution. Disadvantage: incorrect stencils will only be detected when psyclone runs (unless we modify kernels, in which case it would be guaranteed that psyclone runs before the compiler anyway ;) ).

[rupertford]

I asked @thomasmelvin for his thoughts by email and he said

"Finite difference notation is a bit ad-hoc so I don't think there are really many standard names for what you are looking for. For stencils I would also include a 16 point stencil (this would be centred on a vertex) as this is widely used for cubic interpolation (such as in ENDGame's semi-lagrangian scheme). I would also include central differences to your finite difference list.

Where there is more standard nomenclature is in the grid staggering, i.e. A-grid, C-grid, D-grid etc, in the horizontal and co-located, Lorenz, Charney-Philips in the vertical but I don't know if this is of use to you"

[TeranIvy]

To add to Tom's comment and naming confusion, FD stencils of the same order representing spatial derivatives are usually called "centred" (symmetric) and "upwind" or "downwind" (asymmetric). "forward" and "backward" more often refer to the FD discretisation of time derivatives.

There are widely used names for some FD schemes which mainly arose from solving initial value problems (hence referring to temporal discretisations). For instance, the first order Euler method can be "forward" (explicit) or "backward" (implicit) in time or "upwind" or "downwind" in space. The symmetric second-order scheme in time is called "leapfrog" and in space usually just "centred". I can dig out a table of the commonly used names if that helps.

Odd ordered schemes (even number of gridpoints/cells) are asymmetric whereas the even ordered (odd number of gridpoints/cells) can be made centred or asymmetric (for a high enough stencil extent) depending on which stencil point is chosen to derive the scheme around (however they are usually centred).

In space, the schemes are also called by the degree of polynomial used to derive them ("linear", "parabolic", "cubic", etc.), with the prefix "bi" or "tri" for 2D or 3D applications. Perhaps here they could be denoted by (polynomial degree, centred/upwind/downwind) for spatial and (polynomial degree, centred/forward/backward) for temporal discretisations? For instance, "linear upwind" (1D 2 point stencil), "biparabolic centred" (2D 5 point centred stencil), "biquadratic centred" (2D 9 point centred stencil), etc.

[rupertford]

Many thanks for the input @TeranIvy. One of the issues I see is that we need to be able to support users writing stencil accesses from domains beyond just atmospheric and ocean dynamics so we might be limiting ourselves by only supporting certain names. We might also get cases where "strange" stencils are required.

I certainly like the idea of supporting names that make sense to the scientist but I think we might (also?) need a general way of specifying a stencil, as @hiker was suggesting. My current view is that the simplest way to do this is using coordinates. This is verbose but has the advantage of being able to specify any stencil shape in any number of dimensions and with any depth.

For example ...

((-1,0),(1,0)) for linear_centred? and ((-1,0),(1,0),(0,-1),(0,1)) for biparabolic_centred?

However, what would ((0,-1),(0,1)) be? linear_centred too? We would need to distinguish.

So, I'm currently thinking we could support both forms (explicit shape specification and numerical notation), if you think you can flesh out a set of terms and keywords. This would mean we would have more that one way of specifying the same thing but I don't think that is too big a deal.

Perhaps we could allow users to specify names for different stencil shapes so the terminology could be tailored to the problem and domain, rather than supporting a fixed vocabulary?

Obviously we would also benefit from input from @arporter for how ocean models fit in.

[thomasmelvin]

Just to add to the 'strange' stencils list, I have added a 2d region stencil in https://code.metoffice.gov.uk/trac/lfric/ticket/919 and this has the 'nice' property that it can be a different size at different points on the grid (e.g. around corners of a cuber sphere).

An extra comment on the naming, I think calling them something like linear_centred or biparabolic_centred is confusing what the the stencils might be used for (e.g. linear approximation) with what the stencil is providing, your examples seem to me to be: ((-1,0),(1,0)) = depth1_left_right ((0,-1),(0,1)) = depth1_up_down

Secondly, do you ever plan to support non-quadrilateral based grids? in which case the (x,y) stencil id doesn't make as much sense sine the grid is not logically rectangular.

Maybe you could base the stencil id's on how many neighbours a cell has, i.e. for a quad with four neighbours your examples would be: (1,0,1,0) : linear centred (in x-direction) (0,1,0,1) : linear centred (in y-direction) (1,1,1,1) : biparabolic_centred and for a triangle you might have (1,1,1) : biparabolic_centred

Some more thought would have to go into how to do larger 2d stencils...

[rupertford]

Thanks @thomasmelvin. The GOcean API is specialised to logically rectangular grids. I presume we would have a separate API for triangular (but share infrastructure and optimisations where possible).

I guess it boils down to "what is the domain?" of our DSL (what we confusingly call an API). In our case the thinking is that we keep the domain (API) relatively specialised so that we can take advantage of structure for performance. Another way of thinking about it is that we support coded kernels which (when written natively) have an assumed data layout, so we need to specialise down to a domain which naturally uses and allows us to take advantage of this structure for performance.

[TeranIvy]

Speaking of a general way of specifying stencils, multidimensional stencils can have different depths in each dimension as per @hiker's example "NN-SS-W-E". Looking at @rupertford's and @thomasmelvin's suggestions so far, it seems to me that there are two issues in play here.

• There are four "states" that a stencil cell/point can be in: 1a) Focus cell participating in the stencil, 1b) Focus cell not participating in the stencil, 2a) Neighbour cell participating in the stencil, 2b) Neighbour cell not participating in the stencil.

• There can be different combination of cells (schemes) in each dimension, so each dimension would need to be represented separately.

By "focus cell" I mean the one for which the result is calculated (and around which an approximation polynomial is built).

2a) and 2b) can be denoted with 1 and 0 respectively. If the focus cell always participates in the stencil and has equal number of neighbour cells in each direction then this information is enough to represent it. For instance, the 7-point "NN-SS-W-E" stencil can be decomposed as "W-F-E" and "NN-F-SS" (where "F" is the focus cell) and represented as (1,1) tuple in WE (increasing x) direction and (1,1,1,1) in SN (increasing y) direction where it is implicitly assumed than "F" is also 1 so it is not included in the description.

If the focus cell does not participate in the stencil (such as in "leapfrog" scheme) or it has different number of neighbouring cells on each side, then it would need to have be denoted differently from neighbouring cells. Perhaps 1a) and b) can be denoted with 3 and 2 (or some characters)? For instance, the 7-point "NN-SS-W-E" stencil could be represented as (1,3,1) tuple in WE and (1,1,3,1,1) in SN direction. Leapfrog in 1 D could be represented as (1,2,1) and forward Euler (upwind) as (1,3). Alternatively, tuples could be condensed into "number strings" in each dimension, such as (131, 11311) for "NN-SS-W-E" stencil or (131, 131) for 2D biparabolic centred.

[rupertford]

A few comments and questions.

1) I don't think it matters to PSyclone whether the focus cell is, or is not, participating in the stencil as an access to this cell does not result in any communication (which is what we care about), although there is no harm including it.

2) I'm afraid I don't understand what NN and SS actually mean - do they mean N, then N again i.e. (0,1) and (0,2) in coordinate notation, or perhaps just coordinate (0,2), or something else?

3) I don't see how the NN, SS etc notation and your proposed number notation above works for the combinations i.e. NE, SE, NW and SW. i.e. I don't think you can treat EW separately from NS.

4) I think we can treat stencils as 2D (horizontal) + 1D (vertical) as we are limiting our API to this domain. Therefore I think the vertical dimension can be treated separately from the horizontal, just giving us up and down (rather than a full 3D stencil) - at least I hope that is correct.

5) My proposal for coordinates would be difficult to get working with a compiler so I'd be happy to drop it

6) I think that PSyclone just needs to know whether we need to communicate halo information to a neighbour or not and if so, what depth. For a 2D stencil there will be 8 types of neighbour relationships N,S,E,W,NE,NW,SE,SW. For each of these we need to know what depth of halo to exchange (1..n). There are a number of ways to capture this information generically e.g.

stencil(type="N,S,E,W",depth=1111) or stencil(type=(010,111,010) depth=1111), or have depth as part of the type stencil(type=(010,111,010))

In the above case I think naming it e.g. stencil(type=fd_cross, depth=1) would be nicest

For a strange case such as an X with NE accessing depth 2 we could have

stencil(type="NE,SE,SW,NW",depth=2111) or stencil(type=(101,010,101),depth=1211), or have depth as part of the type ... stencil(type=(102,010,101)).

[hiker]

My intention was that each N indicates one (more) halo row on the north: so N would indicate one row to the North, NN two rows etc.. I am not sure if we intend to support diagonals (e.g. a 9 point stencil with all 8 neighbours) - if so, we would need to use: N-NE-E-SE-S-SW-W-NW (N-S-W-E would not exchange the diagonal halo). And it gets messy if we have larger diagonal stencils: N, NENE, ... 😢

BUT: while I really like using N etc, we would need to enforce a certain orientation of the mesh (e.g. N is towards decreasing J, S would be towards increasing J). I have no idea, but I doubt that all applications would agree that N is towards decreasing J coordinates etc. Of course, we could add this information to the grid: N is decreasing J, ... But that seems unnecessary complicated (why does PSyclone need to know where North is) and complicate things by providing one more parameter that either needs to be tested for consistency (all grids must have N towards decreasing j).

It took me a while to get the triplets ( (0,1,0), (1,1,1), (0,1,0) ) - that would only be self-documenting when written in 2d in the code:

stencil(type=(101,  &
              111,  &
              010), depth=...

If we are using this, I would recommend to include the depth (as suggested by Rupert: (102) etc.).

That kind of leaves: stencil(type="N,S,E,W",depth=1111) But it has the same issue as just using the number of 'N's for depths - we enforce an orientation on the grid.

Or use pairs to indicate this: (1,0) -one column at the right (increasing i), and (0,-2) two rows at the top (decreasing J). So the above 5 point stencil would become: ( (-1,0), (0,1), (1,0), (0,-1) ) (I don't think we need to include the 'focus' point (0,0) in this.

Otherwise I admit I don't like using the long/mathematical names ... too much 'natural science' for me :) Not only would I get confused, it makes it harder to support different stencils: Most other options mentioned so far would not require a code change if we suddenly have a stencil with 3 rows top, 1 row bottom, 2 left, and 4 right or so - they would just work. And I don't like supporting two options either - makes things unnecessary complicated, harder to enforce coding style to users.

TBH, I am ok with anything except for the last paragraph.

[TeranIvy]

Good point about the representation of "diagonals", I completely forgot about this! Regarding long/mathematical names, perhaps we could have a "mapping" between internal PSyclone representation and the commonly used names for each API/science branch the API represents in the documentation?

P.S. In most of the atmospheric models based on some version of lat-lon grid I have seen, the positive y axis goes from South to North so N is the direction of increasing J. Plus, when working with non-quad cells (or quad cells with non-standard orientation) the terms like N and S do not really make sense - important things are which edges (and/or other entity?) to inspect for neighbours, and whether the contribution from a neigbouring cells is added (+1) or subtracted (-1).

[rupertford]

@TeranIvy, I agree about supporting mappings. This would allow @hiker and like minded "un-natural scientists" to not use names and "natural scientists" to use names. My personal preference would be to use names in general and only fall back to a description where needed but people can choose. If fear my un-naturalness is weakening.

@TeranIvy and @hiker. The idea of using N,S,E,W etc is only meant to be logical not physical. What we actually care about is which dimension of an array we are talking about and whether we are indexing +1 or -1. If using N,S,E,W confuses matters then we should not use it. I would therefore lean toward triplets (which would also need a definition for which dimension is which, but that is easy to define), or coordinates (which would follow the array indexing ordering). I think the most compact, practical and easy to use in Fortran is triplets so I vote for that. I also agree with @hiker that we can simply capture depth information here by using numbers e.g. stencil(020,212,020) for a depth 2 'cross' stencil, so do not need an explicit depth. If we go to 3d=2D+1D we can treat the 1D separately e.g. stencil(h=(010,111,010),v=111). Triplet notation does not necessarily exactly capture a complex stencil with depth greater than 1 but I believe it captures all we require for PSyclone.

[arporter]

Like @hiker, it took me a while to get the triplet notation too. I'm also confused about what the central "1" in the following:

stencil(020, &
        212, &
        020)

means? Presumably that location is the 'owned'/current grid point and therefore is always accessed but depth is meaningless and thus it must always be "1"? If so, I propose maybe using "x" instead to show that it is a special point.

Having said all of that, it seems to me that the triplet notation is the simplest way to capture everything we need and all scientists (whether un-natural or not) will understand it. If there are any agreed names out there (e.g. "five-point") then I think a map would work well.

I also agree with keeping the vertical separate as we are aiming at a domain where the MPI decomposition is in the horizontal.

[rupertford]

@arporter, yes the 1 in the middle is the location of the current point i.e the point at the centre of the stencil. @TeranIvy called this the focus cell and pointed out that it might, or might not, be accessed - e.g. in a leapfrog scheme it is not apparently. I was using 1 to mean it is accessed and would expect 0 to mean it is not. The focus cell is counted in the *-point stencil terminology e.g. 5-point is (010,111,010). However I can't think of a situation where this information would be useful to PSyclone. We could put an x but wouldn't that mean using strings rather than integers? If so, I think the neatest thing would be to limit the value to 0 or 1.

[rupertford]

Created a Stencil_Metadata branch which will add metadata processing support for stencil metadata in the Gocean1.0 api.

[arporter]

169 has been merged to master. Closing issue.