Upgrade fields to store boundary conditions

[ali-ramadhan]

This PR (part 3/3) upgrades the field abstraction so fields store their own boundary conditions. This simplifies the model boundary condition hierarchy and generalizes the field and boundary conditions abstractions so they can be used for a compressible model (and any other model we come up with in the future). All future fields will have boundary conditions so PRs like #601 won't be necessary again.

The only change in user interface is that you pass a named tuple to the model constructor now instead of an instance of SolutionBoundaryConditions.

This also works for LES diffusivities so the amount of convoluted scripting gymnastics is much reduced (see test from #601). Setting a diffusivity BC is now almost as easy as a tracer BC.

grid = RegularCartesianGrid(FT, size=(16, 16, 16), length=(1, 1, 1))

buoyancy_bcs = TracerBoundaryConditions(grid, bottom=BoundaryCondition(Gradient, bz))
κₑ_bcs = DiffusivityBoundaryConditions(grid, bottom=BoundaryCondition(Value, κ₀))
model_bcs = (b=buoyancy_bcs, κₑ=(b=κₑ_bcs,))

model = IncompressibleModel(
    grid=grid, architecture=arch, float_type=FT, tracers=:b, buoyancy=BuoyancyTracer(),
    closure=AnisotropicMinimumDissipation(), boundary_conditions=model_bcs
)

Internally: No surprise, this change ended up being pretty invasive. But note that we now have a more flexible and easier to use package with fewer lines of code!

I'm happy to discuss and iterate over the choices that were made in this PR. But glad that I was able to make these changes. Development of the compressible model can continue based on this branch.

Changes:

1. Fields has a new property: field.boundary_conditions.

2. Better pretty printing for fields:

Field located at (Cell, Cell, Cell)
├── data: OffsetArrays.OffsetArray{Float64,3,Array{Float64,3}}, size: (18, 18, 18)
├── grid: RegularCartesianGrid{Float64, Periodic, Periodic, Bounded}(Nx=16, Ny=16, Nz=16)
└── boundary conditions: x=(west=Periodic, east=Periodic), y=(south=Periodic, north=Periodic), z=(bottom=NoFlux, top=NoFlux)

3. Field tuple constructors have two versions now. One that accepts boundary conditions (useful in model constructors) and another that accepts full fields (useful for checkpoint restoration and other custom functionality).

4. You can pass anything into VelocityFields and TracerFields. Perhaps we should check whether everything passed into these functions is a Field or AbstractField. But for now that's extra code and it's not an exported function so I'm okay leaving it as is (we might want the flexibility in the future?). A check was added in Ability to construct field tuples with non-zero data #627.

5. Model no longer has the boundary_conditions property.

6. There is no need for SolutionBoundaryConditions and ModelBoundaryConditions anymore: code is simpler 🎉

7. TurbulentDiffusivities has been renamed to DiffusivityFields for consistency with VelocityFields and TracerFields, etc.

8. For each turbulence closure, DiffusivityFields gets the same two versions (one if you want control over the full field, another if you just want to set boundary conditions).

9. Halo filling and parts of the time stepping have been simplified as there is no need to juggle around and bundle up boundary conditions.

Some comments:

1. It's not clear whether fields for abstract operations should have boundary conditions as they don't have halos and we don't impose boundary conditions on them. Do we want to generally set their boundary conditions to nothing? I did this for the abstract operations tests and they all passed.
2. We could further simplify the time stepping code if we adapt the Field abstraction to be GPU compatible (The future of the Field abstraction #298).
3. This is probably an edge case but I don't think setting diffusivity BCs when using a tuple of closures will actually work. Or rather, I'm not sure how to specify two different diffusivity BCs, one for each closure.
4. Before merging I still need to update the checkpointer and update the documentation.

Resolves #606

Note: This PR branches off from #628.

[glwagner]

We could further simplify the time stepping code if we adapt the Field abstraction to be GPU compatible (#298).

Just a note: I attempted this while working on abstract operations and was unable to solve the problem. But we are using new versions of julia + cuda tools all the time so it is / will be worth revisiting

[glwagner]

This is probably an edge case but I don't think setting diffusivity BCs when using a tuple of closures will actually work. Or rather, I'm not sure how to specify two different diffusivity BCs, one for each closure.

For tuples of closures, everything related to the closure is tupled, including the DiffusivityFields. If there are different diffusivity fields for different closures, then it is possible to set them with different diffusivities?

[glwagner]

For each turbulence closure, DiffusivityFields gets the same two versions (one if you want control over the full field, another if you just want to set boundary conditions).

I don't understand this point, can it be elaborated?

[ali-ramadhan]

We should always use default boundary conditions for grid and the location, eg Face or Cell along each dimension.

I'm okay with doing this but should we be worried that not all {Cell,Cell,Cell} fields are tracer fields and not all {Face,Cell,Cell} fields are u-velocity fields?

Just a note: I attempted this while working on abstract operations and was unable to solve the problem. But we are using new versions of julia + cuda tools all the time so it is / will be worth revisiting

Hmmm, do you remember what you got stuck on? I'd be in favor of reopening the issue and trying to adapt Field to work in GPU kernels.

For tuples of closures, everything related to the closure is tupled, including the DiffusivityFields. If there are different diffusivity fields for different closures, then it is possible to set them with different diffusivities?

Ah sorry I meant that passing a named tuple like

boundary_conditions = (b=buoyancy_bcs, κₑ=(b=κₑ_bcs,))

to the model constructor will apply diffusivity boundary conditions as expected. But if you have two closures and you want to impose different diffusivity boundary conditions on each one then you can't do so by passing a named tuple to the model constructor. You can of course use DiffusivityFields to accomplish this.

We could upgrade the logic for interpreting boundary condition named tuples so that you can do something like

κ1_bcs = (T=some_diffusivity_bcs, S=some_diffusivity_bcs)
κ2_bcs = (T=some_other_diffusivity_bcs, S=some_other_diffusivity_bcs)
boundary_conditions = (T=T_bcs, κₑ=(κ1_bcs, κ2_bcs))

but it felt like a pretty edge case (very very few users would probably use it) so I didn't do it. But we can do it, it's within scope for this PR.

[ali-ramadhan]

For each turbulence closure, DiffusivityFields gets the same two versions (one if you want control over the full field, another if you just want to set boundary conditions).

I don't understand this point, can it be elaborated?

Ah the first one allows you to pass pre-constructed fields as keyword arguments. Users will probably not use this, it's mostly useful for custom diffusivity field tuple constructions. I had restore_from_checkpoint in mind when I wrote this constructor. If you pass it a field for x, it will use it, otherwise it constructs a field with default boundary conditions.

https://github.com/climate-machine/Oceananigans.jl/blob/000158447877ddfcd0f53bff7a3c44fabd29fe88/src/TurbulenceClosures/diffusivity_fields.jl#L34-L41

The second one allows you to pass non-default boundary conditions as a named tuple. Here bcs will be the named tuple users pass to the model constructor. Users will probably not use constructor either. This constructor is meant to be used by model constructors. If bcs contains boundary conditions for x they will be used, otherwise it constructs a field with default boundary conditions.

https://github.com/climate-machine/Oceananigans.jl/blob/000158447877ddfcd0f53bff7a3c44fabd29fe88/src/TurbulenceClosures/diffusivity_fields.jl#L43-L50

It might be possible to merge the two constructors (which have different uses) into one but I couldn't figure out how to do it. VelocityFields and TracerFields has two constructors for the same reason. TendencyFields only has one because users cannot specify tendency boundary conditions by passing a named tuple of FieldBoundaryConditions to a model constructor.

I can add docstrings with these explanations.

[ali-ramadhan]

Another thing to potentially discuss if whether we want to switch from using Face and Cell to Face and Center. I think we agreed to do this at some point, this PR might be a good place to do so.

EDIT: Well, actually this PR does quite a lot already. Better leave this for a future PR.

[glwagner]

I'm okay with doing this but should we be worried that not all {Cell,Cell,Cell} fields are tracer fields and not all {Face,Cell,Cell} fields are u-velocity fields?

Defaults are not about correctness, but convenience --- the concept of a "default" implies that it is not always correct (otherwise we would not use the word "default", because it would never have to change). I think defaults can make user's lives simpler. I think there may be a lot of cases in field definition where boundary conditions are unused.

Hmmm, do you remember what you got stuck on?

Was just getting the dreaded "dynamic function invocation". It's worth a shot and documenting our progress. It seems unlikely there's a fundamental issue. We just have to figure it out.

But if you have two closures and you want to impose different diffusivity boundary conditions on each one then you can't do so by passing a named tuple to the model constructor. You can of course use DiffusivityFields to accomplish this.

Okay thank you I understand. The issue here is the ability to use defaults. For that, we would have to require users to name closures when used in a tuple.

I am happy to leave this functionality as "unsupported". That said, I suspect we will have to confront this issue in the future if we want to do sophisticated GCM stuff, which can have complicated sets of turbulence closures.

[glwagner]

but it felt like a pretty edge case (very very few users would probably use it) so I didn't do it.

It may become common if users do realistic GCM stuff, because there will be many turbulence closures, but only one of them would require explicit boundary-diffusivity modeling to control field gradients and diffusivities along boundaries (for example, prescribing some interesting lateral diffusivity parameterization on top of a vertical mixing parameterization).

EDIT: we don't need to confront these complexities now of course, but it's worth keeping in mind. One brute force solution is just to define convenience constructors for models themselves (OceanCirculationModel) that sets up an IncompressibleModel in "MITgcm" mode (for example) 😏. Then we can do whatever needs to be done there, provided it's at least possible.

[ali-ramadhan]

It may become common if users do realistic GCM stuff, because there will be many turbulence closures, but only one of them would require explicit boundary-diffusivity modeling to control field gradients and diffusivities along boundaries (for example, prescribing some interesting lateral diffusivity parameterization on top of a vertical mixing parameterization).

Ah interesting I didn't think of such cases. I believe we can set up arbitrarily complicated boundary conditions so anything is possible. It's just hard and probably unfeasible/undesirable to extend the user interface to handle all these cases.

[glwagner]

Pre-emptive approval for planned changes.

Thought: we provide the sugary syntax .top and .bottom for boundary conditions in z. Should we also provide east, west, south, and north, and avoid using .left and .right in the code for full clarity?

[ali-ramadhan]

Thought: we provide the sugary syntax .top and .bottom for boundary conditions in z. Should we also provide east, west, south, and north, and avoid using .left and .right in the code for full clarity?

Sounds like a good idea. I think setbc! and getbc aren't fully tested so also worth adding some simple tests there.

I'll open an issue to document this.

[codecov]

Codecov Report

❗ No coverage uploaded for pull request base (master@0001584). Click here to learn what that means.
The diff coverage is 87.03%.

@@            Coverage Diff            @@
##             master     #631   +/-   ##
=========================================
  Coverage          ?   78.11%           
=========================================
  Files             ?      119           
  Lines             ?     2326           
  Branches          ?        0           
=========================================
  Hits              ?     1817           
  Misses            ?      509           
  Partials          ?        0

Impacted Files	Coverage Δ
src/Models/incompressible_model.jl	85.71% <ø> (ø)
src/BoundaryConditions/BoundaryConditions.jl	100% <ø> (ø)
src/Fields/Fields.jl	100% <ø> (ø)
src/Oceananigans.jl	100% <100%> (ø)
src/Fields/set!.jl	67.64% <100%> (ø)
src/Fields/field.jl	76.19% <62.5%> (ø)
src/OutputWriters/checkpointer.jl	91.52% <87.5%> (ø)
src/OutputWriters/output_writer_utils.jl	62.16% <90%> (ø)
...rc/BoundaryConditions/field_boundary_conditions.jl	92% <93.75%> (ø)

---

Continue to review full report at Codecov.

Legend - Click here to learn more
Δ = absolute <relative> (impact), ø = not affected, ? = missing data
Powered by Codecov. Last update 0001584...871ca40. Read the comment docs.