Storing boundary conditions inside fields

[ali-ramadhan]

The issue of whether fields should store their boundary conditions was brought up in #601. We've discussed this idea before but I'm opening this issue to see if this is something we want to pursue.

Note: Moving boundary conditions to fields should happen alongside or after the introduce of grid topologies (#489) which would decouple the grid topology from the boundary conditions (although some boundary conditions will be determined by the grid topology).

I think it's a good idea that will pay off in code clarity and ease of use.

1. Functions such as fill_halo_regions! will simplify as we no longer need to pass in both a field and boundary conditions, just the field itself.
2. No need to prepare boundary condition tuples to pass to functions, for example, in complete_pressure_correction_step!.
3. Oceananigans has have a pretty deep hierarchy of boundary conditions: BC -> Coordinate BCs -> field BCs -> Solution BCs -> Model BCs. Moving BCs to fields will eliminate the need for solution BCs and model BCs, which will create a simpler hierarchy of just BC -> coordinate BCs -> field BCs. This will make it easier to developers and users to access and interact with boundary conditions, especially as ModelBoundaryConditions can only grow in complexity.
4. Storing fields in one struct and boundary conditions in another struct feels antithetical to orthogonal design.

Another reason to consider having fields carry their boundary conditions around is to avoid having to build more boundary condition machinery, i.e. more solution_and_model_boundary_conditions.jl when new fields are added or new models are added (e.g. #605).

User API may have to change as well, although I don't think we should get rid of the current boundary_conditions Model kwarg unless we have a better interface.

[ali-ramadhan]

Copy pasting a comment by @glwagner from #601 so it doesn't get lost in a closed PR:

It's an interesting philosophical question whether fields should carry around their boundary conditions or not. @kburns might have something to say. I think there are some good arguments in favor of this view.

There are some oddities: for example, the boundary conditions on pressure are not zero gradient; however we use these "effective" boundary conditions for convenience and (at least eventually) will include inhomogeneous parts of the pressure boundary condition on the RHS of the pressure Poisson equation. So at least for pressure the "boundary conditions" we would specify are not actually the true boundary conditions.

Once we have topological information embedded in our AbstractGrid, adding default boundary conditions for fields may not be too difficult. But it may not just be a refactor. It seems it might change the user API, too.

Originally posted by @glwagner in #601 (comment)

[glwagner]

Both the grid specification and the boundary condition API will change when we have the grid topology, because there won't be a need for monstrosities like HorizontallyPeriodicSolutionBCs. So, maybe this would be a nice change to make immediately following the upgrade to both grids. I agree that it would be nice to de-complexify the deep hierarchical structure of the current boundary_conditions container.

I agree we should probably keep the boundary_conditions keyword argument to Model --- but note that boundary_conditions would no longer be a Model field.

The Model constructor will get a bit more complicated, because we can no longer initialize fields within the function signature of Model. On the other hand, I have noticed that doing everything inside function signatures as keyword arguments has the downside of making errors that occur there much more difficult to locate and diagnose. So maybe this isn't such a bad thing. We will have to instead initialize field tuples as nothing, and then build the field if it isn't explicitly specified. I'm ok with this --- I suspect we will be working on smoothing out kinks in the user API for some time.

We will also need a function that does

new_field = with_boundary_conditions(old_field, new_field_boundary_conditions)

that creates a new Field wrapper around the array in old_field, but with specified rather than default boundary conditions.

The correct implementation of Flux boundary conditions (eg determining gradients / ghost point values when flux boundary conditions are specified) is still something of a puzzle, but the situation for fields-with-boundary-conditions would be no worse than it is with the current design.

[glwagner]

Another change that we may want to wrap into the changes discussed in this issue (because it also affects fill_halo_regions! and field initialization) is to the convention we use for face-centered directions. Currently, the left halo of a face-centered field lies outside the domain, while the right "halo point" of a face-centered field is actually located on the boundary.

We need to change this definition so that all halo points are located outside the domain for face-centered fields. In other words, parent arrays need to have N + 1 + 2H points along face-centered dimensions, and N + 2H points along cell-centered directions. Because we use OffsetArrays for indexing I don't think this change will affect any code except field initialization and fill_halo_regions! (which makes assumptions about indexing due to the need to reference the parent array directly).

This change is necessary for the correct evaluation of the pressure Poisson equation in bounded domains and for general equations / boundary conditions. There will also need to be concomitant changes in fill_halo_regions! for no-penetration boundary conditions.

[kburns]

Just a quick thought: mathematically it sometimes makes sense to use boundary conditions to reduce the size of the function space a variable can live in, i.e. considering a certain field to be in the space of functions that have zero value (or zero normal derivative) at a boundary.

Physically, though, I'm not sure this is so useful, since many systems have non-trivial boundary conditions, including nonlinear conditions and conditions involving multiple fields. Then the notion of a certain BC even "belonging" to a certain field is a little less clear... you may just have a set of fields that need to collectively satisfy a set of conditions on the boundary.

E.g. thermal radiation conditions (heat flux depends nonlinearly on temperature)
E.g. nontrivial stress BCs (maybe where viscosity is temperature-dependent, like ice)