RobinBoundaryCondition to complement Gradient and Value boundary conditions

[glwagner]

We need to implement RobinBoundaryCondition, which is a boundary condition of the form

$$ \partial_n c = a c + b $$

for two inputs $a$ and $b$. On a staggered grid and a "right boundary" (eg $$i=N+1$$) this condition looks like

$$ \frac{c_{N+1} - c_N}{\Delta} = a \frac{c_{N+1} + c_N}{2} + b $$

which implies a rule for the halo point,

$$ c_{N+1} = \frac{c_N \left ( 1 + \frac{a \Delta}{2} \right ) + b \Delta}{1 + \frac{a \Delta}{2}} $$

@joernc, @shriyafruitwala

Right now we have two boundary conditions which invoke rules for filling halos: Gradient and Value. The way it's currently implemented, these boundary conditions imply rules for determining the across-boundary gradient, eg:

Oceananigans.jl/src/BoundaryConditions/fill_halo_regions_value_gradient.jl

Line 9 in 2dc85f0

@inline left_gradient(bc::GBC, c¹, Δ, i, j, args...) = getbc(bc, i, j, args...)

for GradientBoundaryCondition, where getbc returns g for user input GradientBoundaryCondition(g). Another example is ValueBoundaryCondition(v),

Oceananigans.jl/src/BoundaryConditions/fill_halo_regions_value_gradient.jl

Line 12 in 2dc85f0

@inline left_gradient(bc::VBC, c¹, Δ, i, j, args...) = ( c¹ - getbc(bc, i, j, args...) ) / (Δ/2)

This implementation has the effect that reconstruction/interpolation across the boundary returns v.

We can rewrite the algebra for RobinBoundaryCondition to fit this pattern, which (I think) yields

$$ \frac{c_{N+1} - c_N}{\Delta} = \frac{b}{1 - \frac{a \Delta}{2}} $$

So we can either refactor how boundary conditions are implemented (first formula) or re-use the infrastructure by just adding a line (above formula).

If we want extra credit, we can allow the getbc infrastructure to work for the two Robin values $a$ and $b$. I think for that we actually want to refactor getbc first --- it's weird that getbc takes in the BoundaryCondition, when it looks like a simpler implementation would invoke

getcondition(bc.condition, i, j, grid, args...)

so the first step might be to do this refactor, then implement RobinBoundaryCondition.