Outputting/Evaluating fields at the bottom of an immersed boundary grid

[francispoulin]

Sorry that I don't know this but when using an immersed boundary method, what's the best way to evaluate a field at the bottom, right at the bathymetry?

[simone-silvestri]

I am not sure I understood correctly, but do you want a 2D field that in i, j has the value that (in 3D) corresponds to the cell just above bathymetry? If this is the case, we don't have a ready-to-go functionality, but you can use conditional operations

The cell just above bathymetry is the only cell in the column for which immersed_peripheral_node(i, j, k-1, grid, LX, LY, LZ) == true (note this holds only for 2D immersed map, not for a 3D mask)

So you can do

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

julia> ibg = ImmersedBoundaryGrid(grid, GridFittedBottom((x, y) -> - x))
10×10×10 ImmersedBoundaryGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo:
├── immersed_boundary: GridFittedBottom(min(h)=-9.50e-01, max(h)=-5.00e-02)
├── underlying_grid: 10×10×10 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo
├── Periodic x ∈ [0.0, 1.0)  regularly spaced with Δx=0.1
├── Periodic y ∈ [0.0, 1.0)  regularly spaced with Δy=0.1
└── Bounded  z ∈ [-1.0, 0.0] regularly spaced with Δz=0.1

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

julia> set!(field, (x, y, z) -> x)
10×10×10 Field{Center, Center, Center} on ImmersedBoundaryGrid on CPU
├── grid: 10×10×10 ImmersedBoundaryGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo
├── boundary conditions: FieldBoundaryConditions
│   └── west: Periodic, east: Periodic, south: Periodic, north: Periodic, bottom: ZeroFlux, top: ZeroFlux, immersed: ZeroFlux
└── data: 16×16×16 OffsetArray(::Array{Float64, 3}, -2:13, -2:13, -2:13) with eltype Float64 with indices -2:13×-2:13×-2:13
    └── max=0.95, min=0.15, mean=0.683333

julia> using Oceananigans.ImmersedBoundaries: immersed_peripheral_node

julia> @inline bottom_condition(i, j, k, grid) = immersed_peripheral_node(i, j, k-1, grid, Center(), Center(), Center())
bottom_condition (generic function with 1 method)

julia> using Oceananigans.Fields: condition_operand

julia> bottom_height = sum(condition_operand(identity, field, bottom_condition, 0.0), dims = 3)
10×10×1 Field{Center, Center, Nothing} reduced over dims = (3,) on ImmersedBoundaryGrid on CPU
├── grid: 10×10×10 ImmersedBoundaryGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo
├── boundary conditions: FieldBoundaryConditions
│   └── west: Periodic, east: Periodic, south: Periodic, north: Periodic, bottom: Nothing, top: Nothing, immersed: ZeroFlux
└── data: 16×16×1 OffsetArray(::Array{Float64, 3}, -2:13, -2:13, 1:1) with eltype Float64 with indices -2:13×-2:13×1:1
    └── max=0.95, min=0.15, mean=0.55

julia> interior(bottom_height)
10×10×1 view(::Array{Float64, 3}, 4:13, 4:13, 1:1) with eltype Float64:
[:, :, 1] =
 0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0
 0.15  0.15  0.15  0.15  0.15  0.15  0.15  0.15  0.15  0.15
 0.25  0.25  0.25  0.25  0.25  0.25  0.25  0.25  0.25  0.25
 0.35  0.35  0.35  0.35  0.35  0.35  0.35  0.35  0.35  0.35
 0.45  0.45  0.45  0.45  0.45  0.45  0.45  0.45  0.45  0.45
 0.55  0.55  0.55  0.55  0.55  0.55  0.55  0.55  0.55  0.55
 0.65  0.65  0.65  0.65  0.65  0.65  0.65  0.65  0.65  0.65
 0.75  0.75  0.75  0.75  0.75  0.75  0.75  0.75  0.75  0.75
 0.85  0.85  0.85  0.85  0.85  0.85  0.85  0.85  0.85  0.85
 0.95  0.95  0.95  0.95  0.95  0.95  0.95  0.95  0.95  0.95

[francispoulin]

Your points are well taken. I am all for clarity and efficiency.

I feel like I have seen a function that finds the index right above the topography before, or maybe even wrote it once? I think it was brute force, which is not very expensive but since you only do it once at the beginning, why not?

I like this idea. Thanks.

[glwagner]

Above, @simone-silvestri says that

immersed_peripheral_node(i, j, k-1, grid, LX, LY, LZ) == true

implies that i, j, k is right above topography, so we can use that.

[navidcy]

Just to make a note, whether is the sum with conditional operations or anything else we settle on, we should add some user functionality so that users can output bottom fields without need to know the algorithm behind it? So let's find the most efficient way to do it!

(I agree it took me a bit to understand what the conditional sum was achieving here.)

Why don't we compute the map with the indices that correspond to the bottom values can be computed when we construct the immersed boundary grid and then be part of the IBG properties? Then this would be part of all fields that live on the IBG. Wouldn't that be helpful?

[francispoulin]

I think having a function that outputs the fields along the immersed boundary is a great idea. Then a second function, that can sum over the physical domain, would also be very helpful.

I am happy to work on this with someone else if there is interest.

[glwagner]

Generally our philosophy is only to compute things as they are needed and store information locally. So I'm not sure we want to precompute the index map "just in case" a user wants it. The other problem is where to put it. GridFittedBottom is in principle agnostic to the grid (ie the bottom height could be a function of x, y). On the other hand ImmersedBoundaryGrid is more general and can be used for 3D immersed boundaries (holes, caves, bridges), not just "2.5D" bathymetry.

I like the name MappedIndexOperation

[navidcy]

OK, fair point. We do not need to precompute anything. Just create the functionality for the user to do it just by prescribing a kwarg or something.