trying to mask with a new immersed boundary method `PartialCellBottom`

[francispoulin]

Following Adroft et al. (1997), we are developing a partial cell method in the following branch: fjp/partial_cell_immersed_boundaries

Presently, we can build the mesh without any complaints, but nothing else.

I thought a good first step would be to try and use mask_immersed_field! on a scalar field. When I try with the current set up, it doesn't complain, but it also doesn't modify the entries either. Clearly more work needs to be done.

I don't understand all of mask_immersed_field.jl but I wonder if we need to build a new scalar_mask. At present, it only works for grids that are AbstractGridFittedBoundary, but the new grid is not of that type.

Any advice on whether we need to create a new scalar_mask for dispatch? I tried something pretty simple but that had no effect so I suspect it's not being used,

@inline function scalar_mask(i, j, k, grid, ::PartialCellBottom, LX, LY, LZ, value, field)
    return @inbounds ifelse(solid_interface(LX, LY, LZ, i, j, k, grid),
                            value,
                            field[i, j, k])
end

Or maybe the problem should be solved by going into mask_immersed_field.jl and modifying solid_interface?

Any suggestions @simone-silvestri ?

[simone-silvestri]

Good idea to do a partial cell Immersed boundary!

weird, from what I see here the mask_immersed_field! should work properly

have you tried to do

struct PartialCellBottom{B, E} <: AbstractGridFittedBoundary
    bottom_height :: B
    minimum_fractional_partial_Δz :: E
end

Like that it seems to be working for me...

julia> grid = RectilinearGrid(size=(1, 1, 10), extent=(1, 1, 1))
1×1×10 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 1×1×1 halo
├── Periodic x ∈ [0.0, 1.0)           regularly spaced with Δx=1.0
├── Periodic y ∈ [0.0, 1.0)           regularly spaced with Δy=1.0
└── Bounded  z ∈ [-1.0, -1.29526e-16] regularly spaced with Δz=0.1

julia> ibg = ImmersedBoundaryGrid(grid, PartialCellBottom((x, y)-> -0.5, 1))
┌ Warning: ImmersedBoundaryGrid is unvalidated and may produce incorrect results. Help validate ImmersedBoundaryGrid by reporting any bugs or unexpected behavior to https://github.com/CliMA/Oceananigans.jl/issues.
└ @ Oceananigans.ImmersedBoundaries ~/Oceananigans.jl/src/ImmersedBoundaries/ImmersedBoundaries.jl:110
ImmersedBoundaryGrid on: 
    architecture: CPU()
            grid: 1×1×10 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 1×1×1 halo
   with immersed: PartialCellBottom{var"#1#2", Int64}

julia> c = CenterField(ibg)
1×1×10 Field{Center, Center, Center} on ImmersedBoundaryGrid on CPU
├── grid: 1×1×10 ImmersedBoundaryGrid{Float64, Periodic, Periodic, Bounded} on CPU with 1×1×1 halo
├── boundary conditions: west=Periodic, east=Periodic, south=Periodic, north=Periodic, bottom=ZeroFlux, top=ZeroFlux, immersed=ZeroFlux
└── data: 3×3×12 OffsetArray(::Array{Float64, 3}, 0:2, 0:2, 0:11) with eltype Float64 with indices 0:2×0:2×0:11
    └── max=0.0, min=0.0, mean=0.0

julia> fill!(c, 1)
3×3×12 Array{Float64, 3}:
[:, :, 1] =
 1.0  1.0  1.0
 1.0  1.0  1.0
 1.0  1.0  1.0

[:, :, 2] =
 1.0  1.0  1.0
 1.0  1.0  1.0
 1.0  1.0  1.0

[:, :, 3] =
 1.0  1.0  1.0
 1.0  1.0  1.0
 1.0  1.0  1.0

;;; … 

[:, :, 10] =
 1.0  1.0  1.0
 1.0  1.0  1.0
 1.0  1.0  1.0

[:, :, 11] =
 1.0  1.0  1.0
 1.0  1.0  1.0
 1.0  1.0  1.0

[:, :, 12] =
 1.0  1.0  1.0
 1.0  1.0  1.0
 1.0  1.0  1.0

julia> event = mask_immersed_field!(c)
KernelAbstractions.CPUEvent(Task (done) @0x000000010f388a20)

julia> wait(event)

julia> interior(c)
1×1×10 view(::Array{Float64, 3}, 2:2, 2:2, 2:11) with eltype Float64:
[:, :, 1] =
 0.0

[:, :, 2] =
 0.0

[:, :, 3] =
 0.0

[:, :, 4] =
 0.0

[:, :, 5] =
 0.0

[:, :, 6] =
 1.0

[:, :, 7] =
 1.0

[:, :, 8] =
 1.0

[:, :, 9] =
 1.0

[:, :, 10] =
 1.0

Ah sorry, i saw now that you don't want the new immersed boundary to be of type AbstractGridFittedBoundary (why is that?).

Then just dispatching mask_immersed_field! on ::PartialCellBottom should do the trick

But at this point I would suggest a new abstract type, something like AbstractPartialCellBoundary and dispatch on that so that you can have more concrete types such as PartialCellBottom and PartialCellBoundary

[glwagner]

I guess we have a decision here, but I think PartialCellBottom probably deserves to subtype AbstractGridFittedBoundary. The masking algorithm is the same... the only additional feature that PartialCellBottom brings is a change in the calculation of the vertical cell spacing in the bottommost cell. For example I think the masking of the velocity field is the same for both partial cells and full cells.

[francispoulin]

Thanks @simone-silvestri. I can confirm that does solve the problem for me as well. Great!

Yes, @glwagner , raises a good point. PartialCellBottom is not the same as GridFittedBoundary, and so the name we give to both doesn't really work.

Would we want to find a name that works for both GridFitted and PartialCells?

[simone-silvestri]

But the solid_interface function would have to be different right?

I ll give an example:

     Solid    Fluid
    -------.................
   |     ∘     |     ∘     |
   f     c     f     c     f
  i-1   i-1    i     i    i+1

In this case z(c[i-1]) < immersed_boundary (i.e. solid_interface(f, i, j, k, grid) = true) but here I think we do not want to mask the normal velocity which resides on f[i]. Actually, what are we going to do with that velocity?

I think usually people interpolate using f[i+1] and assuming that u is zero at the actual immersed boundary (inbetween c[i-1] and f[i])....

[francispoulin]

Very good questions. I don't have very good answers, but happy to share my thoughts.

I don't know how things will change for velocity but it should be clearer after we figure out how to deal with the tracers.

As for how the velocity is set, we always set the normal velocity to be zero at the interface. The only difference is the ParitalCell has a bottom that can be higher and we need to take this into account when computing the fluxes in and out of this region. We haven't gotten there yet but it's part of the plan.

There aren't any essential differences between the current method, what we call GridFitted and Adcroft et al. call a FullCell method, and the PartialCell method. The former is the limit when we choose the height of the topography in a cell interface to be zero, whereas the latter can pick it to be between a tolerance (currently 0.2) and 1 of the full depth. So the GridFitted boundary must have the boundary between fluid and solid at a face, whereas PartialCell can have it anywhere in the interior. That means that the interface is really at i + alpha where 0 \le alpha < 1, or something like that. I can imagine that sometime in the future they are part of the same function, but much testing must happen before.

I don't believe we need to change how solid_interface works but I need to look over how it's defined to make sure.

[simone-silvestri]

Ok, so from what I understand, the immersed_boundary would be the same in case of GridFitted and PartialCell but the latter would "add" an additional height on top of that? (which modifies fluxes in the cell)

In that case, then we wouldn't have to change the masking procedure...

Maybe another doubt, what do we do if the α is larger than half (or even 3/4) of a cell in terms of masking? we still do not mask a cell which is (mostly) solid?

[glwagner]

I think we will use logic such that i-1 is not immersed / solid. Only i-2 is solid. (It'd be better to use k here because I think we are talking about the vertical index.)

[glwagner]

So the logic is slightly different for is_immersed with PartialCellBottom. is_immersed is false for cells that are "partially immersed". It's only true for cells that are "fully immersed".

[glwagner]

Maybe another doubt, what do we do if the α is larger than half (or even 3/4) of a cell in terms of masking? we still do not mask a cell which is (mostly) solid?

We don't mask the partial cells at all. We only mask cells that are fully immersed.

There is also a parameter that controls the "minimum fractional height" of a cell. MITgcm uses 0.2 for this parameter, which means that the partial cell has to be at least 20% fluid; otherwise, it becomes immersed.

[glwagner]

Ok, so from what I understand, the immersed_boundary would be the same in case of GridFitted and PartialCell but the latter would "add" an additional height on top of that? (which modifies fluxes in the cell)

The two changes are:

1. A cell is_immersed if the interface above the cell center is below the bottom (or nearly so). Actually the precise condition requires using minimum_fractional_Δz as noted above; I think we should have something like

function is_immersed(i, j, k, ibg)
    ϵ = Δzᶜᶜᶜ(i, j, k, ibg) * ibg.immersed_boundary.minimum_fractional_Δz
    z_above = znode(c, c, f, i, j, k+1, ibg.grid)
    return z_above - ϵ < get_bottom_height(i, j, k, ibg.grid, ibg.immersed_boundary.bottom_height)
end

function get_bottom_height(i, j, k, grid, bottom_height)
    x, y, z = nodes(c, c, c, i, j, k, grid)
    return bottom_height(x, y)
end

get_bottom_height(i, j, k, grid, bottom_height::AbstractMatrix) = @inbounds bottom_height[i, j]

Should be inspected closely, might be an error...

2. The vertical grid spacings Δzᶜᶜᶜ are modified for the cells that are just above the bottom (ie, is_immersed(k) is false, but is_immersed(k-1) is true).

Note that the interface spacings Δzᶠᶜᶜ have to take the minimum between the Δzᶜᶜᶜ in adjacent cells.

[glwagner]

I don't believe we need to change how solid_interface works but I need to look over how it's defined to make sure.

I think we only need an accurate definition of is_immersed.

[jm-c]

Just a minor comment here, default minimum fraction might be 0.2 but we more often use a smaller value (e.g., 0.1 or smaller).

On Wed, Feb 16, 2022 at 06:18:05PM -0800, Gregory L. Wagner wrote:

Maybe another doubt, what do we do if the ?? is larger than half (or even 3/4) of a cell in terms of masking? we still do not mask a cell which is (mostly) solid?

We don't mask the partial cells at all. We only mask cells that are fully immersed.

There is also a parameter that controls the "minimum fractional height" of a cell. MITgcm uses 0.2 for this parameter, which means that the partial cell has to be at least 20% fluid; otherwise, it becomes immersed.

-- Reply to this email directly or view it on GitHub: https://github.com/CliMA/Oceananigans.jl/issues/2251#issuecomment-1042500683 You are receiving this because you are subscribed to this thread.

Message ID: @.***>

[simone-silvestri]

function is_immersed(i, j, k, ibg) ϵ = Δzᶜᶜᶜ(i, j, k, ibg) * ibg.immersed_boundary.minimum_fractional_Δz z_above = znode(c, c, f, i, j, k+1, ibg.grid) return z_above - ϵ < get_bottom_height(i, j, k, ibg.grid, ibg.immersed_boundary.bottom_height) end

But in this way wouldn't c[i-1] (from my example above) return is_immersed = true and be treated as an immersed cell?

I think maybe is_immersed should stay the same (i.e. does not "know" about the partial component, only about the bottom as in GridFittedBottom).

We just need a way to modify vertical spacings in case is_immersed(k-1)==true && is_immersed(k)==false which can also be expressed with solid_interface(face, k) && !solid_node(center, k)

[glwagner]

function is_immersed(i, j, k, ibg) ϵ = Δzᶜᶜᶜ(i, j, k, ibg) * ibg.immersed_boundary.minimum_fractional_Δz z_above = znode(c, c, f, i, j, k+1, ibg.grid) return z_above - ϵ < get_bottom_height(i, j, k, ibg.grid, ibg.immersed_boundary.bottom_height) end

But in this way wouldn't c[i-1] (from my example above) return is_immersed = true and be treated as an immersed cell?

I think maybe is_immersed should stay the same (i.e. does not "know" about the partial component, only about the bottom as in GridFittedBottom).

We just need a way to modify vertical spacings in case is_immersed(k-1)==true && is_immersed(k)==false which can also be expressed with solid_interface(face, k) && !solid_node(center, k)

The cell i-1 is a partial cell right?

If the bottom is above the center of the cell, then it is immersed for GridFittedBottom, but is not immersed for PartialCellBottom. I might be misinterpreting the figure...

The point is perhaps what "being immersed" means. I think if a cell is immersed it is not part of the domain. The partial cells are definitely part of the domain; the values of the field there contribute to the volume mean, we evaluate tendency there, there are fluxes through the sides. The thing that changes is the cell height (and thus areas and volumes).

[francispoulin]

I am copying this out from immersed_grid_metrics.jl since it seems a little ofset above, and switching i's for k's, since it is the vertical direction.

     Immersed      Fluid
    ----------- ...........
   |     ∘     |     ∘
   f     c     f     c
  k-1   k-1    k     k

My understanding was that in GridFittedBottom we had a transition from solid to fluid at f_k. If that's the case then I thought the partial cell would be the one above with c_k in the itnerior. That's why I thought c_{k-1} is immersed and c_k is at the interface. Agreed?

The grid fitted problem as a special case where the height in the immersed cell is 0, where as in partial cells it can be anything above zero and up to the top of that cell (or within a tolerance). I hope that we can use the same functions for both, otherwise, much more confusion can arise and I am glad we are having this discussion.

Thank you @jm-c for the comment. I will now change the default to 0.1, since that's what is currently used. I presume this means we don't don't want to have partial cells that are in the bottom 20% of the cell or the top 20% of the cell as well?

[glwagner]

I believe the only issue is with cells that are too thin. Reducing the cell height by 10% (90% of original height) is ok, but reducing by 90% (to 10% of original height) is not.

[glwagner]

I am copying this out from immersed_grid_metrics.jl since it seems a little ofset above, and switching i's for k's, since it is the vertical direction.

     Immersed      Fluid

    ----------- ...........

   |     ∘     |     ∘

   f     c     f     c

  k-1   k-1    k     k

My understanding was that in GridFittedBottom we had a transition from solid to fluid at f_k. If that's the case then I thought the partial cell would be the one above with c_k in the itnerior. That's why I thought c_{k-1} is immersed and c_k is at the interface. Agreed?

The question is where bottom_height lies in these diagrams.

GridFittedBottom is not a limiting case of PartialCellBottom in the way we have coded it, because GridFittedBottom uses the cell center to determine whether the cell is immersed or not.

[francispoulin]

Thanks @glwagner . I think that PartialCellBottom should also use the cell center to determine if a cell is immersed or not.

I agree that GrifFittedBottom is not a limiting case right now but it would be nice if we get this as a limit of PartialCellBottom. We could do this by instead of using the depth in the centre, we use the depth in the center truncated down to the lower cell. Something I will try and keep in mind as we move forward.

[glwagner]

Thanks @glwagner . I think that PartialCellBottom should also use the cell center to determine if a cell is immersed or not.

I agree that GrifFittedBottom is not a limiting case right now but it would be nice if we get this as a limit of PartialCellBottom. We could do this by instead of using the depth in the centre, we use the depth in the center truncated down to the lower cell. Something I will try and keep in mind as we move forward.

If we use

bottom_height(x, y) > znode(c, c, f, i, j, k, grid)

As criteria, then GridFittedBottom is a limiting case when the minimum fractional height is 1.0.

But this doesn't seem like a sensible formulation for GridFittedBottom. Why do we want them to be limited cases of each other?

[francispoulin]

It seems to me that both methods are trying to do the same thing.

They both identify a cell that the topopgraphy passes through and picks a height to represent the topogrpahy, and it should be the same cell. In the partial cell method we find the height to be something between the bottom and the top of that cell. The grid fitted boundary used the bottom (or the top, if I'm mistaken).

If they are both doing the same thing and the only differences is the particular height we assign, then do we want to have two completely different methods? Instead we could have one general method that allows for the two (or more) options.

We might want to have different options on how to pick the height of the topography.

1. Bottom of the cell (GridFittedBoundary)
2. Top of the cell
3. Midpoint rule (2nd order approximation, and what I am currently playing with)
4. Simpsons rule (4th order approxiamtion, and not that much harder to do).

If it's not that much more effort then why wouldn't we want to give the user a choice in how to represent the topogrpahy?

Also, one method would be less to maintain, I would think, when we decide what that method should be and code it up.

This is just a suggestion and I am of course happy to follow whatever people decide is best moving forward.

[glwagner]

It seems to me that both methods are trying to do the same thing.

They both identify a cell that the topopgraphy passes through and picks a height to represent the topogrpahy, and it should be the same cell. In the partial cell method we find the height to be something between the bottom and the top of that cell. The grid fitted boundary used the bottom (or the top, if I'm mistaken).

If they are both doing the same thing and the only differences is the particular height we assign, then do we want to have two completely different methods? Instead we could have one general method that allows for the two (or more) options.

We might want to have different options on how to pick the height of the topography.

1. Bottom of the cell (GridFittedBoundary)
2. Top of the cell
3. Midpoint rule (2nd order approximation, and what I am currently playing with)
4. Simpsons rule (4th order approxiamtion, and not that much harder to do).

If it's not that much more effort then why wouldn't we want to give the user a choice in how to represent the topogrpahy?

Also, one method would be less to maintain, I would think, when we decide what that method should be and code it up.

This is just a suggestion and I am of course happy to follow whatever people decide is best moving forward.

I think generalizing GridFittedBottom so that we have different methods available to identify the "height of a cell" is a nice idea.

Using bottom_height evaluated at cell centers is a low-order method in a sense. So we can have higher-order methods available. You can add a property to GridFittedBottom that allows the user to control that behavior, something like height_model (although I hate the word "model" so maybe there's a better name).

[francispoulin]

We can run our attempts at a partial cell method without any errors. :)

Today we plan to look at the results and see if the are reasoanble. Thinking about some tests to add would probably be a good idea.

As you know, the two methods are virtually identical. I am happy to go the way of generalizing GridFittedBottom and add in a tolerance and a way of specifying the height. Shouldn't be much work, it's just a matter of deciding exactly what we want this to look like.

I'll share our results later on and maybe that will help.

[glwagner]

As you know, the two methods are virtually identical. I am happy to go the way of generalizing GridFittedBottom and add in a tolerance and a way of specifying the height. Shouldn't be much work, it's just a matter of deciding exactly what we want this to look like.

The only change to existing code you need is a new property in GridFittedBottom.

Then the current methods are the fallback (default), and we extend behavior by defining new functions for specific cases.

[simone-silvestri]

We can run our attempts at a partial cell method without any errors. :)

Today we plan to look at the results and see if the are reasoanble. Thinking about some tests to add would probably be a good idea.

As you know, the two methods are virtually identical. I am happy to go the way of generalizing GridFittedBottom and add in a tolerance and a way of specifying the height. Shouldn't be much work, it's just a matter of deciding exactly what we want this to look like.

I'll share our results later on and maybe that will help.

Great development! I'm excited to see the results! It's great if we can have a better way to represent bathymetry, as deep cells are usually very coarse

[francispoulin]

Thanks @glwagner . Will work on that after we're happy with the results.

Thanks @simone-silvestri . Partial cells are very easy but certain a little better from the current approach. We have hopes of doing shaved cells too, but that will certainly be more work. Will need to think about that a lot more, but I think it would be a great addition.

[francispoulin]

With help from @glwagner , I believe we have a partial cell method in this branch here.

We are going to run some tests to make sure that everything is behaving as normal but wanted to give you an udpate. Hopefully have something to review soon.