Skip to content

Commit 4aa81b5

Browse files
tomchorclaudesimone-silvestri
authored
(0.105.5) Add unit test for zero momentum flux at immersed peripheral nodes (#5402)
Co-authored-by: Claude Sonnet 4.6 <noreply@anthropic.com> Co-authored-by: Simone Silvestri <silvestri.simone0@gmail.com>
1 parent cca93d2 commit 4aa81b5

2 files changed

Lines changed: 127 additions & 1 deletion

File tree

‎src/Advection/immersed_advective_fluxes.jl‎

Lines changed: 29 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -41,6 +41,35 @@ end
4141
@inline _advective_tracer_flux_y(i, j, k, ibg::IBG, ::Nothing, args...) = zero(ibg)
4242
@inline _advective_tracer_flux_z(i, j, k, ibg::IBG, ::Nothing, args...) = zero(ibg)
4343

44+
# Fallback for `nothing` advection
45+
@inline _advective_momentum_flux_Uu(i, j, k, ibg::IBG, ::Nothing, args...) = zero(ibg)
46+
@inline _advective_momentum_flux_Uv(i, j, k, ibg::IBG, ::Nothing, args...) = zero(ibg)
47+
@inline _advective_momentum_flux_Uw(i, j, k, ibg::IBG, ::Nothing, args...) = zero(ibg)
48+
@inline _advective_momentum_flux_Vu(i, j, k, ibg::IBG, ::Nothing, args...) = zero(ibg)
49+
@inline _advective_momentum_flux_Vv(i, j, k, ibg::IBG, ::Nothing, args...) = zero(ibg)
50+
@inline _advective_momentum_flux_Vw(i, j, k, ibg::IBG, ::Nothing, args...) = zero(ibg)
51+
@inline _advective_momentum_flux_Wu(i, j, k, ibg::IBG, ::Nothing, args...) = zero(ibg)
52+
@inline _advective_momentum_flux_Wv(i, j, k, ibg::IBG, ::Nothing, args...) = zero(ibg)
53+
@inline _advective_momentum_flux_Ww(i, j, k, ibg::IBG, ::Nothing, args...) = zero(ibg)
54+
55+
# dx(uu), dy(vu), dz(wu)
56+
# ccc, ffc, fcf
57+
@inline _advective_momentum_flux_Uu(i, j, k, ibg::IBG, scheme, U, u) = conditional_flux_ccc(i, j, k, ibg, zero(ibg), advective_momentum_flux_Uu(i, j, k, ibg, scheme, U, u))
58+
@inline _advective_momentum_flux_Vu(i, j, k, ibg::IBG, scheme, V, u) = conditional_flux_ffc(i, j, k, ibg, zero(ibg), advective_momentum_flux_Vu(i, j, k, ibg, scheme, V, u))
59+
@inline _advective_momentum_flux_Wu(i, j, k, ibg::IBG, scheme, W, u) = conditional_flux_fcf(i, j, k, ibg, zero(ibg), advective_momentum_flux_Wu(i, j, k, ibg, scheme, W, u))
60+
61+
# dx(uv), dy(vv), dz(wv)
62+
# ffc, ccc, cff
63+
@inline _advective_momentum_flux_Uv(i, j, k, ibg::IBG, scheme, U, v) = conditional_flux_ffc(i, j, k, ibg, zero(ibg), advective_momentum_flux_Uv(i, j, k, ibg, scheme, U, v))
64+
@inline _advective_momentum_flux_Vv(i, j, k, ibg::IBG, scheme, V, v) = conditional_flux_ccc(i, j, k, ibg, zero(ibg), advective_momentum_flux_Vv(i, j, k, ibg, scheme, V, v))
65+
@inline _advective_momentum_flux_Wv(i, j, k, ibg::IBG, scheme, W, v) = conditional_flux_cff(i, j, k, ibg, zero(ibg), advective_momentum_flux_Wv(i, j, k, ibg, scheme, W, v))
66+
67+
# dx(uw), dy(vw), dz(ww)
68+
# fcf, cff, ccc
69+
@inline _advective_momentum_flux_Uw(i, j, k, ibg::IBG, scheme, U, w) = conditional_flux_fcf(i, j, k, ibg, zero(ibg), advective_momentum_flux_Uw(i, j, k, ibg, scheme, U, w))
70+
@inline _advective_momentum_flux_Vw(i, j, k, ibg::IBG, scheme, V, w) = conditional_flux_cff(i, j, k, ibg, zero(ibg), advective_momentum_flux_Vw(i, j, k, ibg, scheme, V, w))
71+
@inline _advective_momentum_flux_Ww(i, j, k, ibg::IBG, scheme, W, w) = conditional_flux_ccc(i, j, k, ibg, zero(ibg), advective_momentum_flux_Ww(i, j, k, ibg, scheme, W, w))
72+
4473
# Disambiguation for `FluxForm` momentum fluxes....
4574
@inline _advective_momentum_flux_Uu(i, j, k, ibg::IBG, scheme::FluxFormAdvection, U, u) = _advective_momentum_flux_Uu(i, j, k, ibg, scheme.x, U, u)
4675
@inline _advective_momentum_flux_Vu(i, j, k, ibg::IBG, scheme::FluxFormAdvection, V, u) = _advective_momentum_flux_Vu(i, j, k, ibg, scheme.y, V, u)

‎test/test_forcings.jl‎

Lines changed: 98 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -3,7 +3,7 @@ include("dependencies_for_runtests.jl")
33
using Oceananigans.BoundaryConditions: ImpenetrableBoundaryCondition
44
using Oceananigans.Fields: Field
55
using Oceananigans.Forcings: MultipleForcings
6-
using Oceananigans.ImmersedBoundaries: mask_immersed_field!
6+
using Oceananigans.ImmersedBoundaries: mask_immersed_field!, immersed_peripheral_node, peripheral_node
77

88
""" Take one time step with three forcing arrays on u, v, w. """
99
function time_step_with_forcing_array(arch)
@@ -253,6 +253,96 @@ function seven_forcings(arch)
253253
return true
254254
end
255255

256+
# For each node (i, j, k), store whether the condition
257+
# "if immersed peripheral node then flux == 0" holds.
258+
# That is: !immersed_peripheral_node | (flux == 0).
259+
# @test all(interior(t**)) then verifies zero flux everywhere it is required.
260+
@kernel function _populate_momentum_flux_tests!(tuu, tuv, tuw, tvu, tvv, tvw, twu, twv, tww,
261+
grid, scheme, u, v, w)
262+
i, j, k = @index(Global, NTuple)
263+
264+
Fuu = Oceananigans.Advection._advective_momentum_flux_Uu(i, j, k, grid, scheme, u, u)
265+
Fuv = Oceananigans.Advection._advective_momentum_flux_Uv(i, j, k, grid, scheme, u, v)
266+
Fuw = Oceananigans.Advection._advective_momentum_flux_Uw(i, j, k, grid, scheme, u, w)
267+
Fvu = Oceananigans.Advection._advective_momentum_flux_Vu(i, j, k, grid, scheme, v, u)
268+
Fvv = Oceananigans.Advection._advective_momentum_flux_Vv(i, j, k, grid, scheme, v, v)
269+
Fvw = Oceananigans.Advection._advective_momentum_flux_Vw(i, j, k, grid, scheme, v, w)
270+
Fwu = Oceananigans.Advection._advective_momentum_flux_Wu(i, j, k, grid, scheme, w, u)
271+
Fwv = Oceananigans.Advection._advective_momentum_flux_Wv(i, j, k, grid, scheme, w, v)
272+
Fww = Oceananigans.Advection._advective_momentum_flux_Ww(i, j, k, grid, scheme, w, w)
273+
274+
c = Center()
275+
f = Face()
276+
277+
@inbounds begin
278+
tuu[i, j, k] = !peripheral_node(i, j, k, grid, c, c, c) | (Fuu == 0)
279+
tuv[i, j, k] = !peripheral_node(i, j, k, grid, f, f, c) | (Fuv == 0)
280+
tuw[i, j, k] = !peripheral_node(i, j, k, grid, f, c, f) | (Fuw == 0)
281+
tvu[i, j, k] = !peripheral_node(i, j, k, grid, f, f, c) | (Fvu == 0)
282+
tvv[i, j, k] = !peripheral_node(i, j, k, grid, c, c, c) | (Fvv == 0)
283+
tvw[i, j, k] = !peripheral_node(i, j, k, grid, c, f, f) | (Fvw == 0)
284+
twu[i, j, k] = !peripheral_node(i, j, k, grid, f, c, f) | (Fwu == 0)
285+
twv[i, j, k] = !peripheral_node(i, j, k, grid, c, f, f) | (Fwv == 0)
286+
tww[i, j, k] = !peripheral_node(i, j, k, grid, c, c, c) | (Fww == 0)
287+
end
288+
end
289+
290+
"""
291+
Test that all 9 advective momentum flux kernels are zero at every immersed peripheral node,
292+
for a grid with a random bottom boundary. This verifies that removing the explicit
293+
`conditional_flux` zeroing from `_advective_momentum_flux_*` on `ImmersedBoundaryGrid` is
294+
safe: velocity masking in immersed cells and the no-penetration boundary condition at the
295+
immersed boundary face are sufficient to guarantee zero momentum flux at all peripheral
296+
nodes, independently of the velocity magnitude in the active fluid region.
297+
"""
298+
function test_momentum_flux_zero_at_peripheral_nodes(scheme)
299+
Nx, Ny, Nz = 8, 8, 8
300+
underlying_grid = RectilinearGrid(CPU(), size=(Nx, Ny, Nz), extent=(1, 1, 1), halo=(4, 4, 4))
301+
302+
# Random bottom creates a varied immersed boundary topology, stressing all peripheral
303+
# node configurations (not just a uniform flat slab).
304+
bottom_height = -1 .+ rand(Nx, Ny) .* 0.8 # varies between -1 and -0.2
305+
ibg = ImmersedBoundaryGrid(underlying_grid, GridFittedBottom(bottom_height))
306+
model = NonhydrostaticModel(ibg, advection=scheme)
307+
308+
# Set non-zero velocities everywhere, then zero out immersed cells via masking.
309+
set!(model, u=1, v=1, w=1)
310+
mask_immersed_field!(model.velocities.u)
311+
mask_immersed_field!(model.velocities.v)
312+
mask_immersed_field!(model.velocities.w)
313+
fill_halo_regions!(model.velocities)
314+
315+
u, v, w = model.velocities
316+
317+
# Test fields: store true (1) where the "zero-flux at peripheral node" property holds,
318+
# false (0) where it is violated.
319+
tuu = Field{Center, Center, Center}(ibg)
320+
tuv = Field{Face, Face, Center}(ibg)
321+
tuw = Field{Face, Center, Face}(ibg)
322+
tvu = Field{Face, Face, Center}(ibg)
323+
tvv = Field{Center, Center, Center}(ibg)
324+
tvw = Field{Center, Face, Face}(ibg)
325+
twu = Field{Face, Center, Face}(ibg)
326+
twv = Field{Center, Face, Face}(ibg)
327+
tww = Field{Center, Center, Center}(ibg)
328+
329+
launch!(CPU(), ibg, KernelParameters(0:9, 0:9, 0:9), _populate_momentum_flux_tests!,
330+
tuu, tuv, tuw, tvu, tvv, tvw, twu, twv, tww,
331+
ibg, scheme, u, v, w)
332+
333+
@test all(!iszero, Array(interior(tuu)))
334+
@test all(!iszero, Array(interior(tuv)))
335+
@test all(!iszero, Array(interior(tuw)))
336+
@test all(!iszero, Array(interior(tvu)))
337+
@test all(!iszero, Array(interior(tvv)))
338+
@test all(!iszero, Array(interior(tvw)))
339+
@test all(!iszero, Array(interior(twu)))
340+
@test all(!iszero, Array(interior(twv)))
341+
@test all(!iszero, Array(interior(tww)))
342+
343+
return true
344+
end
345+
256346
function test_settling_tracer_comparison(arch; open_bottom=true)
257347
"""
258348
Test that compares settling tracer simulations on regular vs immersed boundary grids.
@@ -426,6 +516,13 @@ end
426516
@test seven_forcings(arch)
427517
end
428518

519+
@testset "Momentum flux zero at immersed peripheral nodes" begin
520+
@info " Testing momentum flux is zero at immersed peripheral nodes..."
521+
for scheme in (Centered(order=2), UpwindBiased(order=3), WENO(order=5))
522+
@test test_momentum_flux_zero_at_peripheral_nodes(scheme)
523+
end
524+
end
525+
429526
@testset "FieldTimeSeries forcing on [$A]" begin
430527
@info " Testing FieldTimeSeries forcing [$A]..."
431528
@test time_step_with_field_time_series_forcing(arch)

0 commit comments

Comments
 (0)