A new boundary condition API

[glwagner]

This PR implements a new boundary condition API.

It's a bit of a work in progress. Before merging we need:

• tests (simple 1D test with temperature to see if budgets are correct?)
• export some appropriate objects / functions to user
• some API development to make boundary condition specification easier
• documentation.

Criticism welcome!

Summary of implementation: the BoundaryCondition

The core type of this PR is BoundaryCondition. The BoundaryCondition has one field, calc, and is made callable via

(bc::BoundaryCondition)(args...) = bc.calc(args...)

A BoundaryCondition has one parameter: BCType. The BCType can be Default, Flux, or Value. We include Default because the specification of boundary conditions is baked in to the algorithm. In the case of 'default' boundary conditions, no calculations are made.

nboundary_conditions = nfields x ndims x 2.

In principle, the user can specify 5 x 3 x 2 boundary conditions: on each of the 5 fields (which we currently have --- unfortunately, we may want more in the future), the user can specify a condition on any of the 6 boundaries. The model is initialized with all 30 boundary conditions set to Default.

We allow for all possible boundaries to have conditions by introducing three structs:

• CoordinateBoundaryConditions with fields left and right
• FieldBoundaryConditions with fields x, y, and z
• BoundaryConditions with fields u, v, w, T, and S.

Specifying boundary conditions and user API

A boundary condition is now specified by defining a function that calculates it within a kernel that loops over the boundary in question. The arguments of a flux boundary conditions on a top or bottom boundary (the only situation supported so far) are

u, v, w, T, S, t, step, Nx, Ny, Nz, Δx, Δy, Δz, i, j

This will probably change in the future.

The user API is rather threadbare at the moment. To define a constant flux boundary condition, for example, the user might write (to specify the boundary condition as a function)

model = Model(...)

const constant_flux = # something
surface_temperature_flux(args...) = constant_flux
flux_bc = BoundaryCondition{Flux}(surface_temperature_flux)
model.boundary_conditions.T.z.right = flux_bc

to specify the boundary condition as a constant,

model = Model(...)

flux_bc = BoundaryCondition(Flux, constant_flux)
model.boundary_conditions.T.z.right = flux_bc

There are probably some sugary things that we can implement to smooth this specification.

Caveats

• We support only flux boundary conditions at top and bottom at the moment
• We can only handle fluxes in terms of the constant vertical viscosities and diffusivities
• If the viscosities and diffusivities are allowed to vary in the vertical at some point in the future, we will have to re-configure the implementation
• We cannot intercept the specification of a flux if we have, for example, a parameterization like KPP that takes care of the surface fluxes with a closure (rather than as a boundary condition).

A side note

I think that, at some point, it would be preferable to move the implementation of boundary conditions and equations to a new file / module so that new physics can be easily implemented my modifying that module / file.

[codecov]

Codecov Report

Merging #118 into master will increase coverage by 4.58%.
The diff coverage is 72.97%.

@@            Coverage Diff             @@
##           master     #118      +/-   ##
==========================================
+ Coverage   52.71%   57.29%   +4.58%     
==========================================
  Files          19       19              
  Lines         645      644       -1     
==========================================
+ Hits          340      369      +29     
+ Misses        305      275      -30

Impacted Files	Coverage Δ
src/Oceananigans.jl	100% <ø> (ø)	⬆️
src/equation_of_state.jl	100% <ø> (+90.9%)	⬆️
src/models.jl	84% <50%> (-3.5%)	⬇️
src/time_steppers.jl	70.34% <69.09%> (+2.38%)	⬆️
src/boundary_conditions.jl	88.23% <88.23%> (-11.77%)	⬇️
src/utils.jl	0% <0%> (-72.73%)	⬇️
src/poisson_solvers.jl	43.08% <0%> (-3.37%)	⬇️
src/diagnostics.jl	0% <0%> (ø)	⬆️
src/output_writers.jl	0% <0%> (ø)	⬆️
... and 1 more

---

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 3307031...b5d3983. Read the comment docs.

[glwagner]

@ali-ramadhan @jm-c I realized that I may have specified the flux boundary condition wrong --- I think that I am double counting the flux coming in/out of the cell from the interior (this part is calculated during the main time-stepping loop. I guess the sanity check is that it shouldn't do anything if flux is zero.

I'll add you as collaborators on my repo in case you want to change anything (I probably won't have time till Monday). Also looks like Travis is failing; not sure why that is.

Feel free to criticize my design. This is just a first pass as I see it.

[glwagner]

I propose changing the struct BoundaryConditions to ModelBoundaryConditions for clarity.

[ali-ramadhan]

Thanks for embarking on this @glwagner!

It's a bit of a work in progress. Before merging we need:

• tests (simple 1D test with temperature to see if budgets are correct?)

Just added "WIP" (work in progress) to the PR title.

We can work on the tests together. For testing, we should be able to extend the 1D column model example: https://github.com/climate-machine/Oceananigans.jl/blob/master/examples/column_model.jl

Also looks like Travis is failing; not sure why that is.

Looks like it's just the one test in the time stepping section that's failing. Maybe something to do with how the boundary conditions are called during the time stepping? The Model test passes so it's not erroring during model initialization/construction.

I propose changing the struct BoundaryConditions to ModelBoundaryConditions for clarity.

I'll always vote for clarity!

Just one initial question: I might be misunderstanding the purpose of bc.calc but why not bc.impose(args...) instead of bc.calc(args...) as we usually say that we impose boundary conditions?

Have to run out but will have a more detailed look later today.

Otherwise, looks really neat! On a more practical note, might end up being easier to discuss this PR in person.

[ali-ramadhan]

Just another thought. If it's implemented in this way with

• CoordinateBoundaryConditions on the left and right parts of the domain.
• FieldBoundaryConditions along the x, y, and z dimensions.
• BoundaryConditions on the u, v, w, T, and S fields.

with the final boundary condition being a function, would there be 30+ nested parameterized types for ModelBoundaryConditions. Is/should that be a concern?

[johncmarshall54]

It's not obvious to me what a field boundary condition is, vs a coordinate boundary condition.
Fundamentally we are applying boundary conditions to u, v, w, T, and S fields, and that's it. It all seems a bit too complicated.

[glwagner]

@johncmarshall54 each field (u, v, w, T, S) can have boundary conditions in each direction (x, y, z).

We could restrict ourselves to specific combinations of boundary conditions instead, which would reduce the number of possibilities. For example, we might have just doubly periodic in (x, y) plus flux in z on all fields, or singly-periodic in x and flux in (y, z) on all fields.

The current design is meant to give the user maximum flexibility, but I'm open to simplifying it. That would restrict what users can do, but perhaps that's ok?

@jm-c feedback welcome!

[glwagner]

would there be 30+ nested parameterized types for ModelBoundaryConditions. Is/should that be a concern?

@ali-ramadhan do you mean with regards to performance? I'm not sure. With multiple dispatch being core to julia it seems this scenario is not uncommon (30+ may not be very large).

Just one initial question: I might be misunderstanding the purpose of bc.calc but why not bc.impose(args...) instead of bc.calc(args...) as we usually say that we impose boundary conditions?

The function calc does not actually impose a boundary condition --- the imposition of boundary condition depends on, for example, the viscosity and diffusivity, and is a property of the equation (or turbulent closure) being implemented. Again for example, the K-Profile-Parameterization includes a modification of how a flux boundary condition is implemented. In other words, the "specification of flux" is separate from the "imposition of a boundary condition". The former is determined by the user. The latter is determined by the model/governing equation.

[glwagner]

It should be noted that this PR also makes progress on addressing #100.

[ali-ramadhan]

We could restrict ourselves to specific combinations of boundary conditions instead, which would reduce the number of possibilities. For example, we might have just doubly periodic in (x, y) plus flux in z on all fields, or singly-periodic in x and flux in (y, z) on all fields.

Agree that CoordinateBoundaryConditions might be a weird name but yeah, maximum flexibility would be very powerful. Maybe the common use case isn't to impose each of the 30 boundary conditions one-by-one but we can just have nice helper functions/abstractions like

model.boundary conditions += HorizontallyPeriodic()

@ali-ramadhan do you mean with regards to performance? I'm not sure. With multiple dispatch being core to julia it seems this scenario is not uncommon (30+ may not be very large).

I'm still pretty new to Julia so yeah don't know if that will be an issue, especially on the GPU. Only way to find out is to try and benchmark! Maybe you're right and 30+ isn't a lot.

@vchuravy any idea on whether 30+ parameterized types for a struct is too many? Would this hurt performance on the GPU?

The function calc does not actually impose a boundary condition --- the imposition of boundary condition depends on, for example, the viscosity and diffusivity, and is a property of the equation (or turbulent closure) being implemented. Again for example, the K-Profile-Parameterization includes a modification of how a flux boundary condition is implemented. In other words, the "specification of flux" is separate from the "imposition of a boundary condition". The former is determined by the user. The latter is determined by the model/governing equation.

I see. So if no parameterizations are being used, are the boundary conditions actually being imposed then? Even with KPP, isn't the boundary condition still being imposed only to later be modified by KPP?

I still feel like bc.calc() feels obscure, I'm not sure why a boundary condition should have to calculated. Perhaps it's just semantics but it would be nice to see boundary conditions be imposed in the code and see something like bc.impose(...).

[glwagner]

So if no parameterizations are being used, are the boundary conditions actually being imposed then?

I believe yes: if the viscosity/diffusivity is 0, then flux or no-slip boundary conditions cannot be imposed (mathematically, the order of the PDE is reduced in that case). We then can only impose no-penetration or periodic conditions.

I still feel like bc.calc() feels obscure

My primitive logic: for a flux boundary condition, bc.calc() "calculates" the flux at the given grid point and time-step. For a "value" boundary condition, bc.calc() "calculates" the value of the boundary condition at the given grid point and time-step. But I agree it is a weird name. What is a better name?

[glwagner]

Other names:

• bc.get
• bc.specification
• ... ?

Agree that CoordinateBoundaryConditions might be a weird name but yeah, maximum flexibility would be very powerful. Maybe the common use case isn't to impose each of the 30 boundary conditions one-by-one but we can just have nice helper functions/abstractions

precisely, and I like the syntax you propose. We can do lots of stuff here to make our user's lives easy.

There are two issues: the backend, and the user interface. Maybe the title of this PR is confusing, because I think it's primarily about the backend.

[ali-ramadhan]

My primitive logic: for a flux boundary condition, bc.calc() "calculates" the flux at the given grid point and time-step. For a "value" boundary condition, bc.calc() "calculates" the value of the boundary condition at the given grid point and time-step. But I agree it is a weird name. What is a better name?

It may not stand up to mathematical rigor but I still like bc.impose(). In your two examples, I feel like a flux is being imposed and a value is being imposed.

If bc.calc() calculates a number then it should be used as

something = bc.calc(args...)

while bc.impose should be used like

bc.impose!(args...)

But now we're just arguing semantics instead of what's important.

There are two issues: the backend, and the user interface. Maybe the title of this PR is confusing, because I think it's primarily about the backend.

I think so too. API suggests more front-end.

I also think discussing these dense and complicated issues (e.g. this PR and #120) among multiple busy people is difficult on GitHub. Little far ahead but maybe the Monday CliMA meetings are a good place to get high-level feedback?

[glwagner]

Sure, we can talk about code design at the CliMa meeting.

Let me summarize what we've discussed so far:

The current design

A "boundary condition" is a flux, gradient, or value of a field (u, v, w, T, S) at the end (left, right) in a given direction (x, y, z).

Pros

• Generality: any boundary condition can be specified
• Similarity to the mathematics: this design corresponds to how boundary conditions are specified in PDEs
• Extensibility: the design accommodates changes in types of boundary conditions or changes in dimensionality. (Adding new solution variables requires code modification, but that is true throughout the code)

Concerns

• It's more complicated than allowing only 3-4 possibilities to users; as @johncmarshall54 says:

It's not obvious to me what a field boundary condition is, vs a coordinate boundary condition.
Fundamentally we are applying boundary conditions to u, v, w, T, and S fields, and that's it. It all seems a bit too complicated.

• The hierarchical dependence of this design on parameterized types creates a lot of possibilities for different boundary conditions to be specified, and this is a concern for performance (@ali-ramadhan)

• We aren't sure what to call the function (currently bc.calc) whose purpose is to return the value of the boundary condition (either the field value, flux, gradient, or whatever) at a given grid point and simulation time

Other miscellaneous thoughts

• This PR implements a backend / abstraction system for specifying boundary conditions. For users we can add as much sugar on top as we want.

• As we have discussed, the way we implement boundary conditions is intimately connected to the way that we specify equations. This can be seen in the current code (though it is commented out) in that the lines that specify a no-slip condition depend on the "vertical viscosity". Adding other diffusive terms (hyperdiffusivity, Leith diffusivity, various turbulent diffusivity parameterizations/closures) can affect how the boundary condition is implemented. We are not ready to develop an equation abstraction system yet, but we must keep in mind the difference between implementing a boundary condition in the context of a set of equations, and specifying a boundary condition. This PR attempts to solve the latter problem.

• It would be nice to know how MITgcm implements different boundary conditions.

[glwagner]

@SandreOuza has requested a Gradient BCType.

Also here are some proposed tests:

0. Constant initial conditions + zero flux = average value of quantity doesn't change
1. Linear temperature profile + prescribed identical fluxes at top and bottom = linear temperature profile
2. Diffusive solution to the free-convection example

[christophernhill]

That sounds like a potentially nice abstraction. Can we get together tomorrow to sync up on details? Chris

…

On Thu, Mar 14, 2019 at 20:27 Gregory L. Wagner ***@***.***> wrote: @SandreOuza <https://github.com/SandreOuza> has requested a Gradient BCType. Also here are some proposed tests: 1. Constant initial conditions + zero flux = average value of quantity doesn't change 2. Linear temperature profile + prescribed identical fluxes at top and bottom = linear temperature profile 3. Diffusive solution to the free-convection example — You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHub <#118 (comment)>, or mute the thread <https://github.com/notifications/unsubscribe-auth/ADXx4GPeeGVVSIRI-T7YAVmhrd4tnTj9ks5vWujkgaJpZM4bmbWT> .

[glwagner]

Ok, sounds good! I’m going to try to write some of those tests this morning.

[ali-ramadhan]

Typo here? Should it be apply_z_bcs! instead of apply_x_bcs!?