Derivatives of a reduced field over an Immersed boundary

[simone-silvestri]

From the discussion #3587 @scott-conn found out that our definition of immersed derivatives ~~does not work very well~~ is wrong when using reduced fields.

Consider, for example, a field reduced in z. If we try performing a derivative in the x-direction it will check the immersed condition https://github.com/CliMA/Oceananigans.jl/blob/7a4b3f04e402be70d45fcb775a4dedef087f3bb0/src/ImmersedBoundaries/conditional_differences.jl#L21 on k == 1. This might be ~~very~~ wrong if the field is a sum (or a mean) over the column.
For example, for integrals we want to check if the whole column is immersed. If we are dealing with a GridFittedBottom, instead, we just need to check the upmost grid cell (k == Nz)

We could procede in a couple of ways from here:

• use simple derivatives for AsbtractOperations that do not account for immersed boundaries.
• We can augment operations on Immersed Reduced fields by attaching to them a condition like in conditional operations

[glwagner]

Can you try to describe the problem? I'm confused by saying "this might be very wrong". Is it wrong, or not? I don't understand.

Why can't we use different conditions for different reductions? Isn't the problem that we shouldn't use immersed_inactive_node for reductions? If it's a reduction then k should not be used at all. The fact that k matters is the root of the issue.

[simone-silvestri]

The problem is not with the reduction but with derivatives that act on reduced fields. The reductions are not an issue, because by performing a reduction we know what operation leads us to a reduced field, so we can perform the reduction accordingly (for example we exclude immersed cells from reductions).

When performing a derivative we use all three indices regardless of the field being reduced or not. In this case we get a funky result where we are trying to evaluate a derivative at k == 1 for a reduced field that does not necessarily live at k == 1.

In other words, the assumption that the field lives at k == 1 is wrong for a reduced field; the right solution, on the other hand, is not so clear cut:

If we assume that reduced fields "lives" nowhere in the reduced direction, then we can remove the k index and just perform an "immersed-boundary-unaware" derivative.

On the other hand, if the reduced field lives on the whole reduced column (like for example an integral) then we need to be aware of the "immersed column", because if the whole column is immersed then the derivative should return a zero.

[glwagner]

Isn't correct behavior to return zero when the whole dimension is immersed? Otherwise, return the derivative.

[glwagner]

We could procede in a couple of ways from here:

Let me provide a few other options:

• Compute the 2D boolean masks for reduced operations in x, y, and z when forming ImmersedBoundaryGrid. Then, use those masks for conditional differencing of reduced fields, rather than calling immersed_inactive_node.

• For GridFittedBottom and PartialCellBottom, as a stop gap, define these methods:

@inline conditional_δx_f(ℓy, ℓz, i, j, k, ibg::GFBIBG, δx, r::ZIRF, args...) = ifelse(immersed_inactive_node(i,   j, ibg.Nz, ibg, c, ℓy, ℓz) |
                                                                                      immersed_inactive_node(i-1, j, ibg.Nz, ibg, c, ℓy, ℓz),
                                                                                      zero(ibg),
                                                                                      δx(i, j, k, ibg.underlying_grid, r, args...))

That use a condition based on whether k == Nz is immersed. This will fix fields that are reduced in z without increasing memory storage or doing a computation in the immersed boundary grid constructor. And that's the most common case anyways.

• When building a ReducedField, compute the mask that has to be applied to abstract operations. Then, extend the _derivative constructor for the combination of a ReducedField argument + ImmersedBoundaryGrid using conditional operation. This has the advantage of avoiding the mask computation in ImmersedBoundaryGrid (since its only needed to do operations on reduced fields). The disadvantage is that different reduced fields have to redo the computation. Also, this only fixes abstract operations and does not fix the internal operators. We also have to update the conditional operators to throw away the immersed boundary grid for reduced fields, or throw away the immersed boundary grid inside the abstract operation.

• A variant on the above approach is to compute the mask when forming Derivative. But then a new mask is computed for every operation.

There's probably a lot of other options. Keep the brainstorming coming.

[glwagner]

If we assume that reduced fields "lives" nowhere in the reduced direction, then we can remove the k index and just perform an "immersed-boundary-unaware" derivative.

On the other hand, if the reduced field lives on the whole reduced column (like for example an integral) then we need to be aware of the "immersed column", because if the whole column is immersed then the derivative should return a zero.

What is the difference between "living nowhere" and "living on the whole column"?

I think the point of reduced fields is that they are derived from a reduction over one or more dimensions. Since they are derived from a reduction, they invoke values from every element in the column / reduced direction. I don't understand what it means to "live nowhere".