Bugfix in compute bottom height (#3940)

* add the tests

* fix docs

docs/src/grids.md

@@ -174,7 +174,7 @@ mountain_grid = ImmersedBoundaryGrid(grid, GridFittedBottom(mountain))
 
 # output
 20×20×20 ImmersedBoundaryGrid{Float64, Bounded, Bounded, Bounded} on CPU with 3×3×3 halo:
-├── immersed_boundary: GridFittedBottom(mean(z)=6.28318, min(z)=1.58939e-8, max(z)=93.9413)
+├── immersed_boundary: GridFittedBottom(mean(z)=4.5, min(z)=0.0, max(z)=100.0)
 ├── underlying_grid: 20×20×20 RectilinearGrid{Float64, Bounded, Bounded, Bounded} on CPU with 3×3×3 halo
 ├── Bounded  x ∈ [-5000.0, 5000.0] regularly spaced with Δx=500.0
 ├── Bounded  y ∈ [-5000.0, 5000.0] regularly spaced with Δy=500.0

src/ImmersedBoundaries/grid_fitted_bottom.jl

@@ -118,7 +118,7 @@ compute_numerical_bottom_height!(bottom_field, grid, ib) =
         z⁺ = znode(i, j, k+1, grid, c, c, f)
         z  = znode(i, j, k,   grid, c, c, c)
         bottom_cell = ifelse(condition isa CenterImmersedCondition, z ≤ zb, z⁺ ≤ zb)
-        @inbounds bottom_field[i, j, 1] = ifelse(bottom_cell, z⁺, zb)
+        @inbounds bottom_field[i, j, 1] = ifelse(bottom_cell, z⁺, bottom_field[i, j, 1])
     end
 end
 

test/runtests.jl

@@ -46,6 +46,7 @@ CUDA.allowscalar() do
     if group == :unit || group == :all
         @testset "Unit tests" begin
             include("test_grids.jl")
+            include("test_immersed_boundary_grid.jl")
             include("test_operators.jl")
             include("test_vector_rotation_operators.jl")
             include("test_boundary_conditions.jl")

test/test_immersed_boundary_grid.jl

@@ -1,16 +1,20 @@
 include("dependencies_for_runtests.jl")
 
-grid = RectilinearGrid(; size=(2, 2, 2), extent = (1, 1, 1))
+grid = RectilinearGrid(; size=(100, 100, 100), extent = (1, 1, 1))
 
 @testset "Testing Immersed Boundaries" begin
 
     @info "Testing the immersed boundary construction..."
-
-    bottom(x, y) = -1 + 0.5 * exp(-x^2 - y^2)
+    bottom(x, y) = rand()
+
     ibg = ImmersedBoundaryGrid(grid, GridFittedBottom(bottom))
 
-    # Unit test (bottom is at the right position)
+    # Test that the bottom is at the same position
+    bottom_height = interior(ibg.immersed_boundary.bottom_height)
 
-    @info "Testing stably stratified initial conditions..."
+    zfaces = znodes(ibg, Face())
 
+    for i in 1:size(ibg, 1), j in 1:size(ibg, 2)
+        @test bottom_height[i, j, 1] ∈ zfaces
+    end
 end