(v0.85.0) Switching cells where fluxes are applied in ImmersedBoundaryCondition to match FieldBoundaryConditions' kwargs

[simone-silvestri]

see #3141
closes #3141

The results from validation/immersed_boundaries/immersed_couette_flow.jl from this branch

immersed_couette_flow.mp4

[glwagner]

Wasn't immersed Couette flow doing the right thing previously? Which one is old and which is new?

[simone-silvestri]

Those are c and u from the couette flow simulation in this branch. Couette flow was doing the right thing before also because immersed boundaries is implemented there as

c_immersed_bc = ValueBoundaryCondition(1)
c_top_bc = ValueBoundaryCondition(-1)
c_bcs = FieldBoundaryConditions(immersed=c_immersed_bc, top=c_top_bc)

Effectively putting an immersed boundary on both the immersed top and the immersed bottom.
it would have failed if

c_imm_bc = ValueBoundaryCondition(1)
c_top_bc = ValueBoundaryCondition(-1)
c_immersed_bc = ImmersedBoundaryCondition(bottom = c_imm_bc)
c_bcs = FieldBoundaryConditions(immersed=c_immersed_bc, top=c_top_bc)

[glwagner]

Those are c and u from the couette flow simulation in this branch. Couette flow was doing the right thing before also because immersed boundaries is implemented there as

c_immersed_bc = ValueBoundaryCondition(1)
c_top_bc = ValueBoundaryCondition(-1)
c_bcs = FieldBoundaryConditions(immersed=c_immersed_bc, top=c_top_bc)

Effectively putting an immersed boundary on both the immersed top and the immersed bottom. it would have failed if

c_imm_bc = ValueBoundaryCondition(1)
c_top_bc = ValueBoundaryCondition(-1)
c_immersed_bc = ImmersedBoundaryCondition(bottom = c_imm_bc)
c_bcs = FieldBoundaryConditions(immersed=c_immersed_bc, top=c_top_bc)

I'm confused. How does it get the right result in the first case?

[simone-silvestri]

In main:

since we write the immersed boundary condition as

c_immersed_bc = ValueBoundaryCondition(1)
c_top_bc = ValueBoundaryCondition(-1)
c_bcs = FieldBoundaryConditions(immersed=c_immersed_bc, top=c_top_bc)

immersed boundary conditions are defined on all facets, and in particular, what is important in this case, is that bc are assigned both at the top and at the bottom. Now, in this case bottom assigns a flux to the cell below the immersed boundary and, therefore, it does nothing (there is no immersed ceiling).

I am here referring to the internals of the code (not the API) where we have a divergence of immersed fluxes defined as

@inline Base.div(i, j, k, grid::AbstractGrid, loc, q_west, q_east, q_south, q_north, q_bottom, q_top) =
    1 / volume(i, j, k, grid, loc...) * (δx_Ax_q(i, j, k, grid, loc, q_west, q_east) + 
                                         δy_Ay_q(i, j, k, grid, loc, q_south, q_north) + 
                                         δz_Az_q(i, j, k, grid, loc, q_bottom, q_top))

and by top and bottom I refer at q_bottom and q_top (in this case q_bottom == 0 and q_top = something)

The heavy lifting is done by the top boundary condition, which assigns a flux to the cell above the immersed boundary.
Therefore, if instead of that we would have written:

c_imm_bc = ValueBoundaryCondition(1)
c_top_bc = ValueBoundaryCondition(-1)
c_immersed_bc = ImmersedBoundaryCondition(bottom = c_imm_bc)
c_bcs = FieldBoundaryConditions(immersed=c_immersed_bc, top=c_top_bc)

The simulation would not have had the intended outcome since in this case

julia> c_immersed_bc = ImmersedBoundaryCondition(bottom = c_imm_bc)
ImmersedBoundaryCondition:
├── west: Nothing
├── east: Nothing
├── south: Nothing
├── north: Nothing
├── bottom: ValueBoundaryCondition: 1
└── top: Nothing

with the bottom field being inactive (q_bottom == 0 because it looks at cells below immersed boundaries)

[simone-silvestri]

on main in case of

c_imm_bc = ValueBoundaryCondition(1)
c_top_bc = ValueBoundaryCondition(-1)
c_immersed_bc = ImmersedBoundaryCondition(bottom = c_imm_bc)
c_bcs = FieldBoundaryConditions(immersed=c_immersed_bc, top=c_top_bc)

and

u_imm_bc = ValueBoundaryCondition(-1)
u_top_bc = ValueBoundaryCondition(1)
u_immersed_bc = ImmersedBoundaryCondition(bottom = u_imm_bc)
u_bcs = FieldBoundaryConditions(immersed=u_immersed_bc, top=u_top_bc)

this is the result

immersed_couette_flow.mp4

while

c_imm_bc = ValueBoundaryCondition(1)
c_top_bc = ValueBoundaryCondition(-1)
c_immersed_bc = ImmersedBoundaryCondition(top = c_imm_bc)
c_bcs = FieldBoundaryConditions(immersed=c_immersed_bc, top=c_top_bc)

and

u_imm_bc = ValueBoundaryCondition(-1)
u_top_bc = ValueBoundaryCondition(1)
u_immersed_bc = ImmersedBoundaryCondition(top = u_imm_bc)
u_bcs = FieldBoundaryConditions(immersed=u_immersed_bc, top=u_top_bc)

gives

immersed_couette_flow.mp4

[tomchor]

Nice catch @Yixiao-Zhang and @simone-silvestri!

IMO this is kind of a serious bug. It affected several of my simulations, for example, and I'm just now finding out. So can we release a patch version with this fix soon please?

[glwagner]

Should we discourage using ImmersedBoundaryCondition?

[simone-silvestri]

I think it is good to use it in hydrostatic simulations where the aspect ratio is very high (for example drag on the sides should not be equal to drag on the bottom)

[glwagner]

Yes, I see now that it is basically a convenience feature (ie friendly), which allows users to build models that satisfy global constraints without having to, for example, compute the bottom surface area of their discrete geometry. More discussion on #3141 .

My only misgiving is that I think it obscures the physics of a problem and can lead to confusion or modeling choices that don't make sense. That said, it's very difficult to achieve convenience without also allowing wrong modeling choices.

[glwagner]

@simone-silvestri @tomchor @Yixiao-Zhang I updated the docs on ImmersedBoundaryGrid. See what you think.

[glwagner]

Probably makes sense to add a test to make sure this doesn't get reverted in the future too

[tomchor]

I think it is good to use it in hydrostatic simulations where the aspect ratio is very high (for example drag on the sides should not be equal to drag on the bottom)

Also if we use a Monin-Obukhov-based drag law the drag coefficients will depend on the grid-spacing in that direction, meaning they'd need to be applied separately if the spacing is different in different directions.

[glwagner]

Also if we use a Monin-Obukhov-based drag law the drag coefficients will depend on the grid-spacing in that direction, meaning they'd need to be applied separately if the spacing is different in different directions.

That's an interesting example. This wouldn't generalize to cut cells, right? Perhaps we need an abstraction that represents the distances to the boundary, for use within boundary conditions. That way we can ensure the distance to the boundary is computed correctly regardless of how the boundary is represented.

As an aside, hasn't it been shown that MOST is not valid in complex terrain? (That wouldn't prevent it from still somehow improving model fidelity, at least in principle.)

[tomchor]

Also if we use a Monin-Obukhov-based drag law the drag coefficients will depend on the grid-spacing in that direction, meaning they'd need to be applied separately if the spacing is different in different directions.

That's an interesting example. This wouldn't generalize to cut cells, right? Perhaps we need an abstraction that represents the distances to the boundary, for use within boundary conditions. That way we can ensure the distance to the boundary is computed correctly regardless of how the boundary is represented.

Agreed. Although I think that's for another PR, right?

As an aside, hasn't it been shown that MOST is not valid in complex terrain? (That wouldn't prevent it from still somehow improving model fidelity, at least in principle.)

Yes, although the discrepancy is less severe for neutral stability, which I suspect encompasses most of the simulations people have been using Oceananigans for. Interestingly, there's been a very recent attempt at generalization that looks pretty promising!

[glwagner]

Agreed. Although I think that's for another PR, right?

For sure but if that case actually requires another abstraction that would mean we can safely recommend not using ImmersedBoundaryCondition, then we shouldn't put too much effort into it.

[tomchor]

@simone-silvestri bump patch version and merge?

[tomchor]

Actually now that I think about it, this is a breaking release. Idk how many people use this feature, but I wonder if we should do something similar to #2990 and throw a warning whenever ImmersedBoundaryCondition is used so that previously working user scripts don't break silently.

[glwagner]

I'd be in favor of always emitting a warning for ImmersedBoundaryCondition that it's experimental (no published papers, few validation tests)