Support Robin boundary conditions for finite difference derivative operators

[valentinaschueller]

Is your feature request related to a problem? Please describe.
The finite difference derivative operators support SetValue and SetGradient boundary operators.
These allow to solve PDE problems with Neumann and Dirichlet boundary conditions (BC), see for example this tutorial.
A third, standard boundary condition for PDEs is the Robin boundary condition; a linear combination of Neumann and Dirichlet BCs given by:
$\left.(a \nabla u \cdot n + b u)\right|_{\partial\Omega} = g$,
with $a,b$ some given constants, $u$ our unknown on $\Omega$, $n$ the normal vector, $\partial\Omega$ the boundary of the domain, and $g$ a function on said boundary.

Right now, as far as I understand it is not possible (or at least not documented and not obvious to me) how one could enforce such a boundary condition.
Robin-type boundary conditions appear in climate modeling, especially in simplified models suitable for testing and analysis.¹

Describe the solution you'd like
It should be possible to provide such a boundary condition for a diffusion-like PDE.
Since it consists of a combination of Dirichlet and Neumann BC, and since Dirichlet BC have to be implemented via SetValue for gradient operators, I think it makes sense to implement this for the gradient operators.
I personally need this for the combination GradientC2F x DivergenceF2C.

With the stencils used, for the lower boundary this will take the form:

a*(u[1] - u[0]) + b u[0] = g
a*u[1] + (b-a) * u[0] = g
=> u[0] = (g - a*u[1]) / (b-a)

and thus can be implemented as

G(1/2) = u[1] - (g - a*u[1]) / (b-a)    (1)

Note that the Dirichlet/Neumann cases are recovered for a=0,b=1, a=1,b=0, respectively.

Describe alternatives you've considered
These are not really alternatives but here are some thoughts I have on this and what make it a bit complicated to get running:

• Implementing this as SetRobin would deviate in spirit from the existing finite difference boundary conditions which use the terms "value" and "gradient" instead of "Dirichlet" and "Neumann". SetLinearCombination might be an alternative but it is not obvious to understand
• SetGradient expects a vector, SetValue a scalar. This is part of why I don't think it is obvious how I could implement a Robin-type BC with the current feature set. Note that Neumann BC and Robin BC typically explicitly set not the gradient but $\nabla u \cdot n$, which, again, is a scalar. In the finite difference framework in ClimaCoupler, as far as I understand from Generalize FiniteDifference operators to work over any axis instead of the hard-coded 3-axis #202, one implicitly assumes to live in the vertical dimension, so $n=(0,0,\pm 1)^T$.
• Note that $b-a=0$ for $b=a$. It might be reasonable to rescale $\tilde{g}=g/a, \tilde{b}=b/a, \tilde{a}=1$ to prevent division by zero in (1).

Footnotes

1. Examples: Bulk conditions at the atmosphere-ocean interface (e.g., Sockwell et al., 2020) or the surface boundary condition for vertical sea ice thermodynamics (e.g., Maykut & Untersteiner, 1971) ↩