Be careful of using end in forcing functions and boundary conditions

[ali-ramadhan]

The docstrings for forcing functions and boundary functions with i, j, k signatures should tell users not to use end as it'll reach into the halo regions when they probably meant grid.Nz, etc.

cc @sandreza

[glwagner]

I think we are moving towards hiding the "discrete form" of boundary condition functions and forcing functions beneath an abstraction layer. When this change is complete, we can write some documentation on how to use the "discrete form"s if needed that includes this information.

[glwagner]

A few additional thoughts:

This problem could be solved if we ever figure out how to adapt fields to work on the GPU (#746). We can then define lastindex for fields properly in a special way.

A second solution is to define a thin wrapper FieldData around OffsetArray with a special lastindex method that covers the same purpose, but does not contain pointers to grid or boundary_conditions (therefore being more likely to work on the GPU, unlike Field).

[ali-ramadhan]

That would be nice, especially #746, but on the topic of defining lastindex: is it possible that someone might want a forcing function to access halo information? Using grid.H{x,y,z} would be better than end but it would be weird if end < grid.Nx + 2grid.Hx when working with halos.

I can't think of a possible use case but might still be good to add a recommendation against using end to preserve flexibility.

[glwagner]

Yea agree that it might be considered an "abuse of abstraction" to define end to mean "the grid point corresponding to the physical edge of the domain". That's what @sandreza expected, however? #UsersComeFirst

Maybe what we need is a section in the docs about working with underlying discrete data structures that includes a practical discussion about halos and the staggered grid, with examples. And we'll recommend using continuous interfaces as much as possible to avoid mistakes.

[glwagner]

it would be weird if end < grid.Nx + 2grid.Hx when working with halos.

I agree with this now. This is better fixed with documentation. It should not be surprising that end does not refer to the final physical grid point. When working with discrete data, a knowledge of the discrete grid is required. With a solid interface for working with continuous function input, we can create the impression that discrete input requires care and careful consideration, which will hopefully mitigate the misunderstanding.

[glwagner]

An alternative solution is to supply our own syntax for working with the boundary of fields, similar to

https://github.com/JuliaArrays/EndpointRanges.jl/blob/master/src/EndpointRanges.jl

For example, we can define east, west, north, south, top, bottom, so that one can write

top_flux(i, j, grid, clock, fields) = fields.c[i, j, top]

This will work best if we generalize adapt_structure for fields so that field location information is preserved, which will allow us to accurately interpret these indices when applied to fields at cell interfaces (currently we throw away the entire field wrapper and keep only data.)

[glwagner]

We could solve this by unwrapping field data into a view rather than exposing the entire OffsetArray. Not sure we want to do that though.

[glwagner]

This is documented.