What API should we implement for specifying fluxes across `GridFittedImmersedBoundary`?

[glwagner]

GridFittedImmersedBoundary uses a masking technique to immerse a grid-fitted boundary into some primary grid. For this purpose it elides diffusive fluxes across immersed interfaces (and soon will elide advective fluxes via #1719).

This permits simulations that conserve tracers and momentum. However, users typically want to enforce boundary conditions across immersed boundaries. The simplest boundary condition to support is FluxBoundaryCondition.

But there are a number of complications. Our current boundary condition API really bakes in the assumption of directionally-aligned boundaries, since fluxes are defined in the x, y, z direction, according to the boundary in question. This only makes sense for boundaries that are aligned with x, y, z, respectively. For boundaries that can point in any direction, it probably makes more sense to define fluxes normal to the boundary. But then if we want to have a consistent API for specifying boundary conditions for both immersed boundaries and grid boundaries, we need to change the current convention for specifying grid boundary conditions. Fluxes would no longer point in the positive direction but inwards or outwards (in the former case, this means that negative tracer fluxes at the top boundary would lead to a decrease in tracer, eg cooling in the case of temperature). Not hard to implement, but certainly a major breaking change.

Another question is where we should store the boundary conditions. Right now boundary conditions are stored in Field.boundary_conditions for each field. If we continue with this approach we need to generalize / redesign FieldBoundaryConditions. It might make sense to flatten the structure and have boundary conditions for west, east, south, north, bottom, top[, immersed_boundary_label]. This has the potential to permit us to specify different boundary conditions on different "parts" of the immersed boundary, if we come up with an abstraction for partitioning an immersed boundary.

With the two changes above we might have an API that looks something like

c_bcs = TracerBoundaryConditions(grid, top = NormalFluxBoundaryCondition(args...),
    boundary_1 = NormalFluxBoundaryCondition(args...), boundary_2 = NormalFluxBoundaryCondition(more_args...))

We also may want to change how boundary conditions are implemented. Currently we launch a bunch of 2D kernels to enforce boundary conditions. But enforcing boundary conditions across immersed boundaries either requires a 3D kernel or a 1D kernel that traverses every fluid-solid interface (probably requiring some tedious bookkeeping to generate the iteration). From the software perspective, it's probably simpler to move the enforcement of flux boundary condition inside the 3D interior tendency kernel and evaluate boundary conditions at the same time that interior tendency contributions are evaluated (in other words, replace diffusive fluxes with the specified boundary normal flux when appropriate).

cc @whitleyv this affects your work.

[glwagner]

I guess an even more grave and challenging issue is that we need an abstraction for vectors, since users will likely want to specify fluxes of momentum rather than fluxes of individual components. Specifying a flux of momentum is what we need, for example, to implement drag boundary conditions.

[francispoulin]

This all sounds great and certainly the right way to go, even though it will be some work.

At the moment we specify 3 types of boundary conditions: Value (Dirichlet), Gradient (Neumann) and Flux. The first one does not need to change, I don't think. The second and third need to be generalized in that not only do we specify a value but we need to specify a direction, but the direction is already contained in grid so we get that for free essentially.

I wonder if this would be easier to try in the vector invariant form?

[glwagner]

Ah, it's true that in the continuous equations a Dirichlet condition can be specified without invoking the boundary normal. I'm not sure that holds for us though because we use the weak form of the equations. This means that all the boundary conditions boil down to specification of a Flux. In a sense Value and Gradient are really convenience features that infer what the flux needs to be (assuming second order evaluation of the gradient across the boundary). So I think we need to know about the boundary normal for all boundary conditions...

The name "vector invariant" is confusing but only refers to the advection of momentum, so I don't think it helps.

After thinking about it a little bit I think we can implement a hack that will allow us to proceed for most cases. I think we can prescribe a priori a zero component of the boundary normal component of the viscous momentum flux. This will allow users to specify drag conditions. For example, a drag condition specifies only boundary-tangential viscous fluxes, and has the form flux = C * u | u | for boundaries with both y and z normal vectors (where | u | can be correctly evaluated due to masking).

To specify a no-slip condition we would also have a normal component of the viscous flux, so we'll have to get a little smarter for that case.

[glwagner]

Ah ok @francispoulin I think I appreciate what you're saying. There's no need to worry about the convention for vectors / fluxes when specifying a Value boundary condition. Everything happens under the hood in the case and we don't need to worry about changing the API since users will still specify ValueBoundaryCondition with the same expectations that they have now.

[francispoulin]

Thanks for confirming that I understood that part correctly. I guess Value conditions are easier to enforce and maybe the first thing to check.

[glwagner]

The difficulty with Value conditions is that they depend on the model / turbulence closure being used (in the simplest case, we can use the user-specification to calculate a gradient, and then infer the cross boundary flux with a diffusivity). We can implement this by implementing some standard notation for the turbulence closures (right now there is a function viscosity, for example, and z_viscosity. We need the x and y components as well).

In the grid-aligned case we use halos to enforce Value boundary conditions, but this approach doesn't work with immersed boundaries.

The Flux case is a bit more straightforward since it doesn't depend on the closure, but does require some reasoning about boundary normal.

[whitleyv]

@glwagner Sorry for a late response. The current implementation I have going projects the velocity into tangential and normal directions (with respect to the immersed surface) before interpolating over the boundary. Since the boundary can be pretty variable in made more sense for us to take this approach. So the user would specify what they want the two tangential and 1 normal direction to be instead of (u,v,w) specifically. Obviously with tracers there is no projection because it's not a vector like that. I know this is a little different than the fitted version, because it's a little harder to specify a normal direction in that case, and also different from what you do for the usual boundary condition. My interpolation can handle Neumann and Dirichlet ie. Value and Gradient and the original idea was to infer Flux from Gradient, but right now everything is hard coded in. The normal and tangential though, would not hold with the current boundary condition set up with North or East or Top etc, but I don't think it would be hard to have whatever implementation is decided and just have 'Tang' and 'Normal' or whatever as well?

[glwagner]

No worries about the late response! Here's a couple of comments, might have more later.

The current implementation I have going projects the velocity into tangential and normal directions (with respect to the immersed surface) before interpolating over the boundary. Since the boundary can be pretty variable in made more sense for us to take this approach. So the user would specify what they want the two tangential and 1 normal direction to be instead of (u,v,w) specifically.

Right, this is what I mean by having an "abstraction for vectors" --- awesome! If the momentum equation is treated in vector form then the three components are coupled. User specification is then on the vector momentum equation; users will specify VelocityBoundaryConditions (rather than boundary conditions for each component, and VectorForcing, rather than forcing on each component. It may also make sense to coalesce the kernels that compute tendencies for each velocity component (but I'm less sure about that). Either way this is a major change to the API at the very least, but probably necessary and something we also need for complex domains for GCM simulations, like the cubed sphere.

the original idea was to infer Flux from Gradient

This makes sense for a continuous immersed boundary where it's not possible to discretely specify fluxes. If we can't discretely specify fluxes, we have to rely on a diffusivity extracted from the specified turbulence closure. This is simple for closures that have isotropic diffusivities, but gets more complicated for closures with anisotropic / tensor diffusivities. Even worse is supporting flux boundary conditions for the case that a closure doesn't use a diffusivity at all...

Note also that flux boundary conditions are used almost exclusively except for direct numerical simulation, so this is indeed an important consideration.

I know this is a little different than the fitted version, because it's a little harder to specify a normal direction in that case,

For GridFittedImmersedBoundary the normal direction can be easily and efficiently inferred from the masking function by evaluating it at offset i, j, k. If fluid is adjacent to solid in any of those directions, you are on a boundary with a normal in the x, y, z direction (respectively), and the sign of the normal direction can be inferred from the offset.

[glwagner]

tl;dr

1. Should we change the flux convention in Oceananigans so that FluxBoundaryCondition specifies the flux along the inward pointing normal? This reverses the convention we use on "right" boundaries (west, north, top). For example, a negative buoyancy flux would cause cooling at the surface with this convention. This is a huge change, but I think is necessary to do things like specifying a drag boundary condition on vertical velocity (I can discuss further why this is the case if needed.)

2. How do we feel about wrapping boundary conditions for immersed boundaries in FieldBoundaryConditions (this requires refactoring FieldBoundaryConditions under the hood, but does not change the API).

We can do 2. without 1. In that case we'll be able to support a lot of common cases, including ordinary oceanic cases (where boundaries do not overhang, and we only specify drag on tangential components). It's really LES cases that become tricky because we can't specify drag on vertical velocity components, or on tangential components when the boundary overhangs.

[francispoulin]

1. I think using other outward or inward normal is a good idea, and since you suggeted inward normal I say yes.
2. I feel fine with this but also don't understand the implications of doing this.

[glwagner]

There's some additional discussion here: https://github.com/CliMA/Oceananigans.jl/pull/2275#issuecomment-1050094675

[glwagner]

Note also that since this issue was posted, we do have the property FieldBoundaryConditions.immersed as a placeholder for boundary conditions on immersed boundaries.

[tomchor]

1. Should we change the flux convention in Oceananigans so that FluxBoundaryCondition specifies the flux along the inward pointing normal? This reverses the convention we use on "right" boundaries (west, north, top). For example, a negative buoyancy flux would cause cooling at the surface with this convention. This is a huge change, but I think is necessary to do things like specifying a drag boundary condition on vertical velocity (I can discuss further why this is the case if needed.)

I quite like the current convention for specifying fluxes since (imho) it's more mathematically intuitive, so my vote goes to maintaining that convention if we can (although I agree that it makes it harder to apply it generally for immersed solids...)

One way to do it would be to leave most of the heavy work for the user, which would have to specify fluxes in each of the immersed solid's interfaces separately. The user's script would look like what I did here: https://github.com/CliMA/Oceananigans.jl/blob/faed0c4ac85409cb94811b0bdba2bbb7becf330a/sandbox/drag_test.jl#L25-L41

The downside is that this won't work for non-grid-fitted boundaries and requires more user-written code. The upside is that could (I think) re-use much of the inner-workings of the current BC implementation.

We can do 2. without 1. In that case we'll be able to support a lot of common cases, including ordinary oceanic cases (where boundaries do not overhang, and we only specify drag on tangential components). It's really LES cases that become tricky because we can't specify drag on vertical velocity components, or on tangential components when the boundary overhangs.

Not sure I understand this point. Do you mind clarifying @glwagner? By "overhang" do you mean that there's no fluid-to-the-bottom-immersed-boundary-to-the-top interface?

[tomchor]

1. Should we change the flux convention in Oceananigans so that FluxBoundaryCondition specifies the flux along the inward pointing normal? This reverses the convention we use on "right" boundaries (west, north, top). For example, a negative buoyancy flux would cause cooling at the surface with this convention. This is a huge change, but I think is necessary to do things like specifying a drag boundary condition on vertical velocity (I can discuss further why this is the case if needed.)

I quite like the current convention for specifying fluxes since (imho) it's more mathematically intuitive, so my vote goes to maintaining that convention if we can (although I agree that it makes it harder to apply it generally for immersed solids...)

Also making this change would mean that a bunch of scripts would fail silently. Meaning that, depending on what/where fluxes are prescribed, the fluxes would flip, producing a completely different result, without Oceananigans throwing an error (since I think we wouldn't change the interface). So if we do go that way I think we'd need to be really careful about it

[glwagner]

Well, we have an API now.