From d9529a211cf122476869c7e98cff723ae537ec19 Mon Sep 17 00:00:00 2001 From: simone-silvestri Date: Mon, 5 Aug 2024 16:36:31 -0400 Subject: [PATCH 01/12] fix a bug in the adapt --- src/ImmersedBoundaries/partial_cell_bottom.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/ImmersedBoundaries/partial_cell_bottom.jl b/src/ImmersedBoundaries/partial_cell_bottom.jl index cce15e6e06..fef90d6c4c 100644 --- a/src/ImmersedBoundaries/partial_cell_bottom.jl +++ b/src/ImmersedBoundaries/partial_cell_bottom.jl @@ -77,10 +77,10 @@ function on_architecture(arch, ib::PartialCellBottom{<:Field}) return PartialCellBottom(new_bottom_height, ib.minimum_fractional_cell_height) end -Adapt.adapt_structure(to, ib::PartialCellBottom) = PartialCellBottom(adapt(to, ib.bottom_height.data), +Adapt.adapt_structure(to, ib::PartialCellBottom) = PartialCellBottom(adapt(to, ib.bottom_height), ib.minimum_fractional_cell_height) -on_architecture(to, ib::PartialCellBottom) = PartialCellBottom(on_architecture(to, ib.bottom_height.data), +on_architecture(to, ib::PartialCellBottom) = PartialCellBottom(on_architecture(to, ib.bottom_height), on_architecture(to, ib.minimum_fractional_cell_height)) """ From 6e5b3077eee2bfc5c6def73ee8d6026cf54693df Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Wed, 7 Aug 2024 20:09:45 -0400 Subject: [PATCH 02/12] first fix --- src/ImmersedBoundaries/partial_cell_bottom.jl | 11 +++++++---- 1 file changed, 7 insertions(+), 4 deletions(-) diff --git a/src/ImmersedBoundaries/partial_cell_bottom.jl b/src/ImmersedBoundaries/partial_cell_bottom.jl index fef90d6c4c..080db5c92c 100644 --- a/src/ImmersedBoundaries/partial_cell_bottom.jl +++ b/src/ImmersedBoundaries/partial_cell_bottom.jl @@ -95,10 +95,13 @@ Criterion is h ≥ z - ϵ Δz """ @inline function _immersed_cell(i, j, k, underlying_grid, ib::PartialCellBottom) - # Face node above current cell - z = znode(i, j, k+1, underlying_grid, c, c, f) - h = @inbounds ib.bottom_height[i, j, 1] - return z ≤ h + # Face node below current cell + z = znode(i, j, k, underlying_grid, c, c, f) + h = @inbounds ib.bottom_height[i, j, 1] + ϵ = ib.minimum_fractional_cell_height + # z + Δz is equal to the face above the current cell + Δz = Δzᶜᶜᶜ(i, j, k, ibg.underlying_grid) + return z + Δz * ϵ ≤ h end @inline bottom_cell(i, j, k, ibg::PCBIBG) = !immersed_cell(i, j, k, ibg.underlying_grid, ibg.immersed_boundary) & From ed6d2520bf64c7f9287df354ddefa19574d09390 Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Thu, 29 Aug 2024 11:00:34 -0400 Subject: [PATCH 03/12] bugfix --- src/ImmersedBoundaries/partial_cell_bottom.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/ImmersedBoundaries/partial_cell_bottom.jl b/src/ImmersedBoundaries/partial_cell_bottom.jl index 080db5c92c..8e6c7b852d 100644 --- a/src/ImmersedBoundaries/partial_cell_bottom.jl +++ b/src/ImmersedBoundaries/partial_cell_bottom.jl @@ -100,7 +100,7 @@ Criterion is h ≥ z - ϵ Δz h = @inbounds ib.bottom_height[i, j, 1] ϵ = ib.minimum_fractional_cell_height # z + Δz is equal to the face above the current cell - Δz = Δzᶜᶜᶜ(i, j, k, ibg.underlying_grid) + Δz = Δzᶜᶜᶜ(i, j, k, underlying_grid) return z + Δz * ϵ ≤ h end From ffca52b365c762d530c002947ba33e75fe3313e4 Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Thu, 29 Aug 2024 12:16:31 -0400 Subject: [PATCH 04/12] these should be all the bugs, still need to add the tests --- src/ImmersedBoundaries/partial_cell_bottom.jl | 14 +++++++++++++- test/test_immersed_boundary_grid.jl | 16 ++++++++++++++++ 2 files changed, 29 insertions(+), 1 deletion(-) create mode 100644 test/test_immersed_boundary_grid.jl diff --git a/src/ImmersedBoundaries/partial_cell_bottom.jl b/src/ImmersedBoundaries/partial_cell_bottom.jl index 8e6c7b852d..e8909a26c8 100644 --- a/src/ImmersedBoundaries/partial_cell_bottom.jl +++ b/src/ImmersedBoundaries/partial_cell_bottom.jl @@ -11,7 +11,7 @@ struct PartialCellBottom{H, E} <: AbstractGridFittedBottom{H} minimum_fractional_cell_height :: E end -const PCBIBG = ImmersedBoundaryGrid{<:Any, <:Any, <:Any, <:Any, <:Any, <:PartialCellBottom} +const PCBIBG{FT, TX, TY, TZ} = ImmersedBoundaryGrid{FT, TX, TY, TZ, <:Any, <:PartialCellBottom} where {FT, TX, TY, TZ} function Base.summary(ib::PartialCellBottom) hmax = maximum(parent(ib.bottom_height)) @@ -140,6 +140,10 @@ end return Δz end +# Make sure Δz works for Flat topologies. We assume PCBIBG requires a Bounded Z-topology. +XFlatPCBIBG = ImmersedBoundaryGrid{<:Any, <:Flat, <:Any, <:Any, <:Any, <:PartialCellBottom} +YFlatPCBIBG = ImmersedBoundaryGrid{<:Any, <:Any, <:Flat, <:Any, <:Any, <:PartialCellBottom} + @inline Δzᶠᶜᶜ(i, j, k, ibg::PCBIBG) = min(Δzᶜᶜᶜ(i-1, j, k, ibg), Δzᶜᶜᶜ(i, j, k, ibg)) @inline Δzᶜᶠᶜ(i, j, k, ibg::PCBIBG) = min(Δzᶜᶜᶜ(i, j-1, k, ibg), Δzᶜᶜᶜ(i, j, k, ibg)) @inline Δzᶠᶠᶜ(i, j, k, ibg::PCBIBG) = min(Δzᶠᶜᶜ(i, j-1, k, ibg), Δzᶠᶜᶜ(i, j, k, ibg)) @@ -148,4 +152,12 @@ end @inline Δzᶜᶠᶠ(i, j, k, ibg::PCBIBG) = min(Δzᶜᶜᶠ(i, j-1, k, ibg), Δzᶜᶜᶠ(i, j, k, ibg)) @inline Δzᶠᶠᶠ(i, j, k, ibg::PCBIBG) = min(Δzᶠᶜᶠ(i, j-1, k, ibg), Δzᶠᶜᶠ(i, j, k, ibg)) +@inline Δzᶠᶜᶜ(i, j, k, ibg::XFlatPCBIBG) = Δzᶜᶜᶜ(i, j, k, ibg) +@inline Δzᶠᶜᶠ(i, j, k, ibg::XFlatPCBIBG) = Δzᶜᶜᶠ(i, j, k, ibg) +@inline Δzᶜᶠᶜ(i, j, k, ibg::YFlatPCBIBG) = Δzᶜᶜᶜ(i, j, k, ibg) + +@inline Δzᶜᶠᶠ(i, j, k, ibg::YFlatPCBIBG) = Δzᶜᶜᶠ(i, j, k, ibg) +@inline Δzᶠᶠᶜ(i, j, k, ibg::XFlatPCBIBG) = Δzᶜᶠᶜ(i, j, k, ibg) +@inline Δzᶠᶠᶜ(i, j, k, ibg::YFlatPCBIBG) = Δzᶠᶜᶜ(i, j, k, ibg) + @inline z_bottom(i, j, ibg::PCBIBG) = @inbounds ibg.immersed_boundary.bottom_height[i, j, 1] diff --git a/test/test_immersed_boundary_grid.jl b/test/test_immersed_boundary_grid.jl new file mode 100644 index 0000000000..5be526b3be --- /dev/null +++ b/test/test_immersed_boundary_grid.jl @@ -0,0 +1,16 @@ +include("dependencies_for_runtests.jl") + +grid = RectilinearGrid(; size=(2, 2, 2), 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) + ibg = ImmersedBoundaryGrid(grid, GridFittedBottom(bottom)) + + # Unit test (bottom is at the right position) + + @info "Testing stably stratified initial conditions..." + +end From 0fd398613170ebb0c0b82b4e55b5bb6659373712 Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Thu, 29 Aug 2024 12:24:08 -0400 Subject: [PATCH 05/12] another bugfix --- examples/internal_tide.jl | 13 +++++++------ src/ImmersedBoundaries/partial_cell_bottom.jl | 2 +- 2 files changed, 8 insertions(+), 7 deletions(-) diff --git a/examples/internal_tide.jl b/examples/internal_tide.jl index 7ef6322a57..48ccb6f956 100644 --- a/examples/internal_tide.jl +++ b/examples/internal_tide.jl @@ -14,6 +14,7 @@ using Oceananigans using Oceananigans.Units +using Oceananigans.ImmersedBoundaries: PartialCellBottom # ## Grid @@ -46,7 +47,7 @@ width = 20kilometers hill(x) = h₀ * exp(-x^2 / 2width^2) bottom(x) = - H + hill(x) -grid = ImmersedBoundaryGrid(underlying_grid, GridFittedBottom(bottom)) +grid = ImmersedBoundaryGrid(underlying_grid, PartialCellBottom(bottom)) # Let's see how the domain with the bathymetry is. @@ -54,7 +55,7 @@ x = xnodes(grid, Center()) bottom_boundary = interior(grid.immersed_boundary.bottom_height, :, 1, 1) top_boundary = 0 * x -using CairoMakie +using GLMakie fig = Figure(size = (700, 200)) ax = Axis(fig[1, 1], @@ -132,7 +133,7 @@ bᵢ(x, z) = Nᵢ² * z set!(model, u=uᵢ, b=bᵢ) -# Now let's built a `Simulation`. +# Now let's build a `Simulation`. Δt = 5minutes stop_time = 4days @@ -234,9 +235,9 @@ n = Observable(1) title = @lift @sprintf("t = %1.2f days = %1.2f T₂", round(times[$n] / day, digits=2) , round(times[$n] / T₂, digits=2)) -u′ₙ = @lift u′_t[$n] - wₙ = @lift w_t[$n] -N²ₙ = @lift N²_t[$n] +u′ₙ = @lift interior(u′_t[$n], :, 1, :) + wₙ = @lift interior( w_t[$n], :, 1, :) +N²ₙ = @lift interior(N²_t[$n], :, 1, :) axis_kwargs = (xlabel = "x [km]", ylabel = "z [km]", diff --git a/src/ImmersedBoundaries/partial_cell_bottom.jl b/src/ImmersedBoundaries/partial_cell_bottom.jl index e8909a26c8..c6983f18ee 100644 --- a/src/ImmersedBoundaries/partial_cell_bottom.jl +++ b/src/ImmersedBoundaries/partial_cell_bottom.jl @@ -101,7 +101,7 @@ Criterion is h ≥ z - ϵ Δz ϵ = ib.minimum_fractional_cell_height # z + Δz is equal to the face above the current cell Δz = Δzᶜᶜᶜ(i, j, k, underlying_grid) - return z + Δz * ϵ ≤ h + return z + Δz * (1 - ϵ) ≤ h end @inline bottom_cell(i, j, k, ibg::PCBIBG) = !immersed_cell(i, j, k, ibg.underlying_grid, ibg.immersed_boundary) & From 0ed035bab572ae9bdb78ed74b0ae1465d96f02f9 Mon Sep 17 00:00:00 2001 From: Gregory Wagner Date: Tue, 24 Sep 2024 19:03:46 -0600 Subject: [PATCH 06/12] Clamp bottom_height to top and bottom of grid --- src/ImmersedBoundaries/grid_fitted_bottom.jl | 12 ++++++++++++ 1 file changed, 12 insertions(+) diff --git a/src/ImmersedBoundaries/grid_fitted_bottom.jl b/src/ImmersedBoundaries/grid_fitted_bottom.jl index f456476f61..4b542e1071 100644 --- a/src/ImmersedBoundaries/grid_fitted_bottom.jl +++ b/src/ImmersedBoundaries/grid_fitted_bottom.jl @@ -87,12 +87,24 @@ Computes `ib.bottom_height` and wraps it in a Field. function ImmersedBoundaryGrid(grid, ib::GridFittedBottom) bottom_field = Field{Center, Center, Nothing}(grid) set!(bottom_field, ib.bottom_height) + launch!(architecture(grid), grid, :xy, clip_bottom_height!, ib.bottom_height, grid) fill_halo_regions!(bottom_field) new_ib = GridFittedBottom(bottom_field, ib.immersed_condition) TX, TY, TZ = topology(grid) return ImmersedBoundaryGrid{TX, TY, TZ}(grid, new_ib) end +const c = Center() +const f = Face() + +@kernel function clip_bottom_height!(z, grid) + i, j, k, = @index(Global, NTuple) + Nz = size(grid, 3) + z_min = znode(i, j, 1, grid, c, c, f) + z_max = znode(i, j, Nz+1, grid, c, c, f) + @inbounds z[i, j, 1] = clamp(z[i, j, 1], z_min, z_max) +end + @inline function _immersed_cell(i, j, k, underlying_grid, ib::GridFittedBottom{<:Any, <:InterfaceImmersedCondition}) z = znode(i, j, k+1, underlying_grid, c, c, f) h = @inbounds ib.bottom_height[i, j, 1] From 00ffa579bc4f4a5e979b478c64301cac3d97f05a Mon Sep 17 00:00:00 2001 From: Gregory Wagner Date: Tue, 24 Sep 2024 19:23:45 -0600 Subject: [PATCH 07/12] Bugfix --- src/ImmersedBoundaries/grid_fitted_bottom.jl | 22 ++++++++++---------- 1 file changed, 11 insertions(+), 11 deletions(-) diff --git a/src/ImmersedBoundaries/grid_fitted_bottom.jl b/src/ImmersedBoundaries/grid_fitted_bottom.jl index 4b542e1071..1c2787c371 100644 --- a/src/ImmersedBoundaries/grid_fitted_bottom.jl +++ b/src/ImmersedBoundaries/grid_fitted_bottom.jl @@ -54,15 +54,15 @@ Keyword Arguments GridFittedBottom(bottom_height) = GridFittedBottom(bottom_height, CenterImmersedCondition()) function Base.summary(ib::GridFittedBottom) - hmax = maximum(ib.bottom_height) - hmin = minimum(ib.bottom_height) - hmean = mean(ib.bottom_height) + zmax = maximum(ib.bottom_height) + zmin = minimum(ib.bottom_height) + zmean = mean(ib.bottom_height) summary1 = "GridFittedBottom(" - summary2 = string("mean(z)=", prettysummary(hmean), - ", min(z)=", prettysummary(hmin), - ", max(z)=", prettysummary(hmax)) + summary2 = string("mean(z)=", prettysummary(zmean), + ", min(z)=", prettysummary(zmin), + ", max(z)=", prettysummary(zmax)) summary3 = ")" @@ -87,7 +87,7 @@ Computes `ib.bottom_height` and wraps it in a Field. function ImmersedBoundaryGrid(grid, ib::GridFittedBottom) bottom_field = Field{Center, Center, Nothing}(grid) set!(bottom_field, ib.bottom_height) - launch!(architecture(grid), grid, :xy, clip_bottom_height!, ib.bottom_height, grid) + launch!(architecture(grid), grid, :xy, clip_bottom_height!, bottom_field, grid) fill_halo_regions!(bottom_field) new_ib = GridFittedBottom(bottom_field, ib.immersed_condition) TX, TY, TZ = topology(grid) @@ -98,11 +98,11 @@ const c = Center() const f = Face() @kernel function clip_bottom_height!(z, grid) - i, j, k, = @index(Global, NTuple) + i, j = @index(Global, NTuple) Nz = size(grid, 3) - z_min = znode(i, j, 1, grid, c, c, f) - z_max = znode(i, j, Nz+1, grid, c, c, f) - @inbounds z[i, j, 1] = clamp(z[i, j, 1], z_min, z_max) + zmin = znode(i, j, 1, grid, c, c, f) + zmax = znode(i, j, Nz+1, grid, c, c, f) + @inbounds z[i, j, 1] = clamp(z[i, j, 1], zmin, zmax) end @inline function _immersed_cell(i, j, k, underlying_grid, ib::GridFittedBottom{<:Any, <:InterfaceImmersedCondition}) From 2ed8f0140d5e363263f27df855b100e2050b96ce Mon Sep 17 00:00:00 2001 From: Gregory Wagner Date: Wed, 25 Sep 2024 06:12:53 -0600 Subject: [PATCH 08/12] Implement clamping for both GF and PC bottom --- .../abstract_grid_fitted_boundary.jl | 12 ++++ src/ImmersedBoundaries/grid_fitted_bottom.jl | 13 +--- src/ImmersedBoundaries/partial_cell_bottom.jl | 59 ++++++++++++------- 3 files changed, 51 insertions(+), 33 deletions(-) diff --git a/src/ImmersedBoundaries/abstract_grid_fitted_boundary.jl b/src/ImmersedBoundaries/abstract_grid_fitted_boundary.jl index 219d41e2cb..a4d113b32a 100644 --- a/src/ImmersedBoundaries/abstract_grid_fitted_boundary.jl +++ b/src/ImmersedBoundaries/abstract_grid_fitted_boundary.jl @@ -53,3 +53,15 @@ const AGFB = AbstractGridFittedBoundary @inline immersed_cell(i, j, k, grid::AbstractGrid{<:Any, Flat, Flat, Flat}, ib::AGFB) = _immersed_cell(1, 1, 1, grid, ib) end +const c = Center() +const f = Face() + +@kernel function clamp_bottom_height!(z, grid) + i, j = @index(Global, NTuple) + Nz = size(grid, 3) + zmin = znode(i, j, 1, grid, c, c, f) + zmax = znode(i, j, Nz+1, grid, c, c, f) + @inbounds z[i, j, 1] = clamp(z[i, j, 1], zmin, zmax) +end + + diff --git a/src/ImmersedBoundaries/grid_fitted_bottom.jl b/src/ImmersedBoundaries/grid_fitted_bottom.jl index 1c2787c371..cee9ef5424 100644 --- a/src/ImmersedBoundaries/grid_fitted_bottom.jl +++ b/src/ImmersedBoundaries/grid_fitted_bottom.jl @@ -87,24 +87,13 @@ Computes `ib.bottom_height` and wraps it in a Field. function ImmersedBoundaryGrid(grid, ib::GridFittedBottom) bottom_field = Field{Center, Center, Nothing}(grid) set!(bottom_field, ib.bottom_height) - launch!(architecture(grid), grid, :xy, clip_bottom_height!, bottom_field, grid) + launch!(architecture(grid), grid, :xy, clamp_bottom_height!, bottom_field, grid) fill_halo_regions!(bottom_field) new_ib = GridFittedBottom(bottom_field, ib.immersed_condition) TX, TY, TZ = topology(grid) return ImmersedBoundaryGrid{TX, TY, TZ}(grid, new_ib) end -const c = Center() -const f = Face() - -@kernel function clip_bottom_height!(z, grid) - i, j = @index(Global, NTuple) - Nz = size(grid, 3) - zmin = znode(i, j, 1, grid, c, c, f) - zmax = znode(i, j, Nz+1, grid, c, c, f) - @inbounds z[i, j, 1] = clamp(z[i, j, 1], zmin, zmax) -end - @inline function _immersed_cell(i, j, k, underlying_grid, ib::GridFittedBottom{<:Any, <:InterfaceImmersedCondition}) z = znode(i, j, k+1, underlying_grid, c, c, f) h = @inbounds ib.bottom_height[i, j, 1] diff --git a/src/ImmersedBoundaries/partial_cell_bottom.jl b/src/ImmersedBoundaries/partial_cell_bottom.jl index c6983f18ee..51b2f75552 100644 --- a/src/ImmersedBoundaries/partial_cell_bottom.jl +++ b/src/ImmersedBoundaries/partial_cell_bottom.jl @@ -14,15 +14,15 @@ end const PCBIBG{FT, TX, TY, TZ} = ImmersedBoundaryGrid{FT, TX, TY, TZ, <:Any, <:PartialCellBottom} where {FT, TX, TY, TZ} function Base.summary(ib::PartialCellBottom) - hmax = maximum(parent(ib.bottom_height)) - hmin = minimum(parent(ib.bottom_height)) - hmean = mean(parent(ib.bottom_height)) + zmax = maximum(parent(ib.bottom_height)) + zmin = minimum(parent(ib.bottom_height)) + zmean = mean(parent(ib.bottom_height)) summary1 = "PartialCellBottom(" - summary2 = string("mean(z)=", prettysummary(hmean), - ", min(z)=", prettysummary(hmin), - ", max(z)=", prettysummary(hmax), + summary2 = string("mean(zb)=", prettysummary(zmean), + ", min(zb)=", prettysummary(zmin), + ", max(zb)=", prettysummary(zmax), ", ϵ=", prettysummary(ib.minimum_fractional_cell_height)) summary3 = ")" @@ -63,6 +63,7 @@ end function ImmersedBoundaryGrid(grid, ib::PartialCellBottom) bottom_field = Field{Center, Center, Nothing}(grid) set!(bottom_field, ib.bottom_height) + launch!(architecture(grid), grid, :xy, clamp_bottom_height!, bottom_field, grid) fill_halo_regions!(bottom_field) new_ib = PartialCellBottom(bottom_field, ib.minimum_fractional_cell_height) TX, TY, TZ = topology(grid) @@ -84,28 +85,42 @@ on_architecture(to, ib::PartialCellBottom) = PartialCellBottom(on_architecture(t on_architecture(to, ib.minimum_fractional_cell_height)) """ - - --x-- - ∘ k+1 - k+1 --x-- ↑ <- node z - ∘ k | Δz - k --x-- ↓ + immersed underlying + + --x-- --x-- + + + ∘ ↑ ∘ k+1 + | + | + k+1 --x-- | k+1 --x-- ↑ <- node z + ∘ ↓ | + zb ⋅⋅x⋅⋅ | + | + ∘ k | Δz + | + | + k --x-- ↓ -Criterion is h ≥ z - ϵ Δz +Criterion is zb ≥ z - ϵ Δz """ @inline function _immersed_cell(i, j, k, underlying_grid, ib::PartialCellBottom) # Face node below current cell z = znode(i, j, k, underlying_grid, c, c, f) - h = @inbounds ib.bottom_height[i, j, 1] + zb = @inbounds ib.bottom_height[i, j, 1] ϵ = ib.minimum_fractional_cell_height # z + Δz is equal to the face above the current cell Δz = Δzᶜᶜᶜ(i, j, k, underlying_grid) - return z + Δz * (1 - ϵ) ≤ h + return (z + Δz * (1 - ϵ)) ≤ zb end -@inline bottom_cell(i, j, k, ibg::PCBIBG) = !immersed_cell(i, j, k, ibg.underlying_grid, ibg.immersed_boundary) & - immersed_cell(i, j, k-1, ibg.underlying_grid, ibg.immersed_boundary) +@inline function bottom_cell(i, j, k, ibg::PCBIBG) + grid = ibg.underlying_grid + ib = ibg.immersed_boundary + # This one's not immersed, but the next one down is + return !immersed_cell(i, j, k, grid, ib) & immersed_cell(i, j, k-1, grid, ib) +end @inline function Δzᶜᶜᶜ(i, j, k, ibg::PCBIBG) underlying_grid = ibg.underlying_grid @@ -140,10 +155,6 @@ end return Δz end -# Make sure Δz works for Flat topologies. We assume PCBIBG requires a Bounded Z-topology. -XFlatPCBIBG = ImmersedBoundaryGrid{<:Any, <:Flat, <:Any, <:Any, <:Any, <:PartialCellBottom} -YFlatPCBIBG = ImmersedBoundaryGrid{<:Any, <:Any, <:Flat, <:Any, <:Any, <:PartialCellBottom} - @inline Δzᶠᶜᶜ(i, j, k, ibg::PCBIBG) = min(Δzᶜᶜᶜ(i-1, j, k, ibg), Δzᶜᶜᶜ(i, j, k, ibg)) @inline Δzᶜᶠᶜ(i, j, k, ibg::PCBIBG) = min(Δzᶜᶜᶜ(i, j-1, k, ibg), Δzᶜᶜᶜ(i, j, k, ibg)) @inline Δzᶠᶠᶜ(i, j, k, ibg::PCBIBG) = min(Δzᶠᶜᶜ(i, j-1, k, ibg), Δzᶠᶜᶜ(i, j, k, ibg)) @@ -152,6 +163,11 @@ YFlatPCBIBG = ImmersedBoundaryGrid{<:Any, <:Any, <:Flat, <:Any, <:Any, <:Partial @inline Δzᶜᶠᶠ(i, j, k, ibg::PCBIBG) = min(Δzᶜᶜᶠ(i, j-1, k, ibg), Δzᶜᶜᶠ(i, j, k, ibg)) @inline Δzᶠᶠᶠ(i, j, k, ibg::PCBIBG) = min(Δzᶠᶜᶠ(i, j-1, k, ibg), Δzᶠᶜᶠ(i, j, k, ibg)) +# Make sure Δz works for horizontally-Flat topologies. +# (There's no point in using z-Flat with PartialCellBottom). +XFlatPCBIBG = ImmersedBoundaryGrid{<:Any, <:Flat, <:Any, <:Any, <:Any, <:PartialCellBottom} +YFlatPCBIBG = ImmersedBoundaryGrid{<:Any, <:Any, <:Flat, <:Any, <:Any, <:PartialCellBottom} + @inline Δzᶠᶜᶜ(i, j, k, ibg::XFlatPCBIBG) = Δzᶜᶜᶜ(i, j, k, ibg) @inline Δzᶠᶜᶠ(i, j, k, ibg::XFlatPCBIBG) = Δzᶜᶜᶠ(i, j, k, ibg) @inline Δzᶜᶠᶜ(i, j, k, ibg::YFlatPCBIBG) = Δzᶜᶜᶜ(i, j, k, ibg) @@ -161,3 +177,4 @@ YFlatPCBIBG = ImmersedBoundaryGrid{<:Any, <:Any, <:Flat, <:Any, <:Any, <:Partial @inline Δzᶠᶠᶜ(i, j, k, ibg::YFlatPCBIBG) = Δzᶠᶜᶜ(i, j, k, ibg) @inline z_bottom(i, j, ibg::PCBIBG) = @inbounds ibg.immersed_boundary.bottom_height[i, j, 1] + From a200fd0599df4e57a8f04dd37faaf22d317d0bd4 Mon Sep 17 00:00:00 2001 From: Gregory Wagner Date: Wed, 25 Sep 2024 06:13:13 -0600 Subject: [PATCH 09/12] Convert internal tide example to use extension --- examples/internal_tide.jl | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/examples/internal_tide.jl b/examples/internal_tide.jl index 48ccb6f956..7fa4f64c68 100644 --- a/examples/internal_tide.jl +++ b/examples/internal_tide.jl @@ -235,9 +235,9 @@ n = Observable(1) title = @lift @sprintf("t = %1.2f days = %1.2f T₂", round(times[$n] / day, digits=2) , round(times[$n] / T₂, digits=2)) -u′ₙ = @lift interior(u′_t[$n], :, 1, :) - wₙ = @lift interior( w_t[$n], :, 1, :) -N²ₙ = @lift interior(N²_t[$n], :, 1, :) +u′n = @lift u′_t[$n] + wn = @lift w_t[$n] +N²n = @lift N²_t[$n] axis_kwargs = (xlabel = "x [km]", ylabel = "z [km]", @@ -249,15 +249,15 @@ fig = Figure(size = (700, 900)) fig[1, :] = Label(fig, title, fontsize=24, tellwidth=false) ax_u = Axis(fig[2, 1]; title = "u'-velocity", axis_kwargs...) -hm_u = heatmap!(ax_u, xu, zu, u′ₙ; colorrange = (-umax, umax), colormap = :balance) +hm_u = heatmap!(ax_u, xu, zu, u′n; nan_color=:gray, colorrange=(-umax, umax), colormap=:balance) Colorbar(fig[2, 2], hm_u, label = "m s⁻¹") ax_w = Axis(fig[3, 1]; title = "w-velocity", axis_kwargs...) -hm_w = heatmap!(ax_w, xw, zw, wₙ; colorrange = (-wmax, wmax), colormap = :balance) +hm_w = heatmap!(ax_w, xw, zw, wn; nan_color=:gray, colorrange=(-wmax, wmax), colormap=:balance) Colorbar(fig[3, 2], hm_w, label = "m s⁻¹") ax_N² = Axis(fig[4, 1]; title = "stratification N²", axis_kwargs...) -hm_N² = heatmap!(ax_N², xN², zN², N²ₙ; colorrange = (0.9Nᵢ², 1.1Nᵢ²), colormap = :thermal) +hm_N² = heatmap!(ax_N², xN², zN², N²n; nan_color=:gray, colorrange=(0.9Nᵢ², 1.1Nᵢ²), colormap=:magma) Colorbar(fig[4, 2], hm_N², label = "s⁻²") fig From 47eeff0db47d62cc3713c6a5fc0f0468dcaf7ca8 Mon Sep 17 00:00:00 2001 From: Gregory Wagner Date: Wed, 25 Sep 2024 07:10:31 -0600 Subject: [PATCH 10/12] apply_regionally *sigh* --- src/ImmersedBoundaries/abstract_grid_fitted_boundary.jl | 5 +++++ src/ImmersedBoundaries/grid_fitted_bottom.jl | 2 +- src/ImmersedBoundaries/partial_cell_bottom.jl | 2 +- 3 files changed, 7 insertions(+), 2 deletions(-) diff --git a/src/ImmersedBoundaries/abstract_grid_fitted_boundary.jl b/src/ImmersedBoundaries/abstract_grid_fitted_boundary.jl index a4d113b32a..f3fb9bfd44 100644 --- a/src/ImmersedBoundaries/abstract_grid_fitted_boundary.jl +++ b/src/ImmersedBoundaries/abstract_grid_fitted_boundary.jl @@ -53,6 +53,11 @@ const AGFB = AbstractGridFittedBoundary @inline immersed_cell(i, j, k, grid::AbstractGrid{<:Any, Flat, Flat, Flat}, ib::AGFB) = _immersed_cell(1, 1, 1, grid, ib) end +function clamp_bottom_height!(bottom_field, grid) + launch!(architecture(grid), grid, :xy, clamp_bottom_height!, bottom_field, grid) + return nothing +end + const c = Center() const f = Face() diff --git a/src/ImmersedBoundaries/grid_fitted_bottom.jl b/src/ImmersedBoundaries/grid_fitted_bottom.jl index cee9ef5424..e007f461a1 100644 --- a/src/ImmersedBoundaries/grid_fitted_bottom.jl +++ b/src/ImmersedBoundaries/grid_fitted_bottom.jl @@ -87,7 +87,7 @@ Computes `ib.bottom_height` and wraps it in a Field. function ImmersedBoundaryGrid(grid, ib::GridFittedBottom) bottom_field = Field{Center, Center, Nothing}(grid) set!(bottom_field, ib.bottom_height) - launch!(architecture(grid), grid, :xy, clamp_bottom_height!, bottom_field, grid) + @apply_regionally clamp_bottom_height!(bottom_field, grid) fill_halo_regions!(bottom_field) new_ib = GridFittedBottom(bottom_field, ib.immersed_condition) TX, TY, TZ = topology(grid) diff --git a/src/ImmersedBoundaries/partial_cell_bottom.jl b/src/ImmersedBoundaries/partial_cell_bottom.jl index 51b2f75552..d08f6db815 100644 --- a/src/ImmersedBoundaries/partial_cell_bottom.jl +++ b/src/ImmersedBoundaries/partial_cell_bottom.jl @@ -63,7 +63,7 @@ end function ImmersedBoundaryGrid(grid, ib::PartialCellBottom) bottom_field = Field{Center, Center, Nothing}(grid) set!(bottom_field, ib.bottom_height) - launch!(architecture(grid), grid, :xy, clamp_bottom_height!, bottom_field, grid) + @apply_regionally clamp_bottom_height!(bottom_field, grid) fill_halo_regions!(bottom_field) new_ib = PartialCellBottom(bottom_field, ib.minimum_fractional_cell_height) TX, TY, TZ = topology(grid) From 9b1244a514c12a3d27bdbe564f68b760e7f9671c Mon Sep 17 00:00:00 2001 From: Gregory Wagner Date: Wed, 25 Sep 2024 09:18:35 -0600 Subject: [PATCH 11/12] Disambiguate --- src/ImmersedBoundaries/abstract_grid_fitted_boundary.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/ImmersedBoundaries/abstract_grid_fitted_boundary.jl b/src/ImmersedBoundaries/abstract_grid_fitted_boundary.jl index f3fb9bfd44..dcc442b6ac 100644 --- a/src/ImmersedBoundaries/abstract_grid_fitted_boundary.jl +++ b/src/ImmersedBoundaries/abstract_grid_fitted_boundary.jl @@ -54,14 +54,14 @@ const AGFB = AbstractGridFittedBoundary end function clamp_bottom_height!(bottom_field, grid) - launch!(architecture(grid), grid, :xy, clamp_bottom_height!, bottom_field, grid) + launch!(architecture(grid), grid, :xy, _clamp_bottom_height!, bottom_field, grid) return nothing end const c = Center() const f = Face() -@kernel function clamp_bottom_height!(z, grid) +@kernel function _clamp_bottom_height!(z, grid) i, j = @index(Global, NTuple) Nz = size(grid, 3) zmin = znode(i, j, 1, grid, c, c, f) From 9807bd06b29d5fb5413e355475270beb8226fec2 Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Wed, 25 Sep 2024 12:35:00 -0400 Subject: [PATCH 12/12] using CairoMakie --- examples/internal_tide.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/internal_tide.jl b/examples/internal_tide.jl index 7fa4f64c68..5b31f18418 100644 --- a/examples/internal_tide.jl +++ b/examples/internal_tide.jl @@ -55,7 +55,7 @@ x = xnodes(grid, Center()) bottom_boundary = interior(grid.immersed_boundary.bottom_height, :, 1, 1) top_boundary = 0 * x -using GLMakie +using CairoMakie fig = Figure(size = (700, 200)) ax = Axis(fig[1, 1],