Skip to content

Add the possibility for an immersed drag #84

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
wants to merge 8 commits into
base: main
Choose a base branch
from
Open

Add the possibility for an immersed drag #84

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
8 commits
Select commit Hold shift + click to select a range
990b0ea
add immersed drag
simone-silvestri Jul 25, 2025
9d3f341
add the ice-advected on coastline to the examples
simone-silvestri Jul 25, 2025
63cfbdb
Update ice_advected_on_coastline.jl
simone-silvestri Jul 25, 2025
584193f
bugfix
simone-silvestri Jul 25, 2025
8f44e16
do this for the moment
simone-silvestri Jul 25, 2025
3998b58
bugfix
simone-silvestri Jul 25, 2025
54bcaea
another bugfix
simone-silvestri Jul 25, 2025
7395fb7
another bugfix
simone-silvestri Jul 25, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
4 changes: 2 additions & 2 deletions docs/make.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -17,7 +17,7 @@ example_scripts = [
"melting_in_spring.jl",
"freezing_of_a_lake.jl",
"ice_advected_by_anticyclone.jl",
# "ice_advected_on_coastline.jl",
"ice_advected_on_coastline.jl",
"arctic_basin_seasonal_cycle.jl"
]

Expand All @@ -31,7 +31,7 @@ example_pages = [
"Melting in Spring" => "literated/melting_in_spring.md",
"Freezing of a Lake" => "literated/freezing_of_a_lake.md",
"Ice advected by anticyclone" => "literated/ice_advected_by_anticyclone.md",
# "Ice advected on coastline" => "literated/ice_advected_on_coastline.md",
"Ice advected on coastline" => "literated/ice_advected_on_coastline.md",
"Arctic basin seasonal cycle" => "literated/arctic_basin_seasonal_cycle.md"
]

Expand Down
24 changes: 20 additions & 4 deletions examples/ice_advected_on_coastline.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -64,15 +64,31 @@ dynamics = SeaIceMomentumEquation(grid;
rheology = ElastoViscoPlasticRheology(),
solver = SplitExplicitSolver(substeps=150))

u_bcs = FieldBoundaryConditions(top = nothing, bottom = nothing,
@inline immersed_u_drag(i, j, k, grid, clock, fields, D) = @inbounds - D * fields.u[i, j, k]
@inline immersed_v_drag(i, j, k, grid, clock, fields, D) = @inbounds - D * fields.v[i, j, k]

immersed_u_bc = FluxBoundaryCondition(immersed_u_drag, discrete_form=true, parameters=3e-1)
immersed_v_bc = FluxBoundaryCondition(immersed_v_drag, discrete_form=true, parameters=3e-1)

immersed_u_bc = ImmersedBoundaryCondition(top=nothing, bottom=nothing, west=nothing, east=nothing,
south=immersed_u_bc, north=immersed_u_bc)

immersed_v_bc = ImmersedBoundaryCondition(top=nothing, bottom=nothing, south=nothing, north=nothing,
west=immersed_v_bc, east=immersed_v_bc)

u_bcs = FieldBoundaryConditions(grid, (Face, Center, Nothing);
north = ValueBoundaryCondition(0),
south = ValueBoundaryCondition(0))
south = ValueBoundaryCondition(0),
immersed = immersed_u_bc)

v_bcs = FieldBoundaryConditions(grid, (Center, Face, Nothing);
immersed = immersed_v_bc)

#Define the model!
model = SeaIceModel(grid;
advection = WENO(order=7),
dynamics = dynamics,
boundary_conditions = (; u=u_bcs),
boundary_conditions = (; u=u_bcs, v=v_bcs),
ice_thermodynamics = nothing)

# We start with a concentration of ℵ = 1 everywhere
Expand Down Expand Up @@ -159,4 +175,4 @@ CairoMakie.record(fig, "sea_ice_advected_on_coastline.mp4", 1:Nt, framerate = 8)
end
nothing #hide

# ![](sea_ice_advected_on_coastline.mp4)
# ![](sea_ice_advected_on_coastline.mp4)
61 changes: 38 additions & 23 deletions src/Rheologies/ice_stress_divergence.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,6 +1,7 @@
using Oceananigans.ImmersedBoundaries: ImmersedBoundaryGrid,
ImmersedBoundaryCondition

using Oceananigans.BoundaryConditions: FBC, getbc
using Oceananigans.Advection:conditional_flux_ccc, conditional_flux_ffc
using Oceananigans.Operators: index_left, index_right

Expand Down Expand Up @@ -52,44 +53,58 @@ end

@inline function immersed_∂ⱼ_σ₁ⱼ(i, j, k, ibg::IBG, u_bc::IBC, rheology, clock, fields)
# Fetch fluxes across immersed boundary
q̃ᵂ = ib_ice_stress_ux(i, j, k, ibg, u_bc.west, rheology, clock, fields)
q̃ᴱ = ib_ice_stress_ux(i+1, j, k, ibg, u_bc.east, rheology, clock, fields)
q̃ˢ = ib_ice_stress_uy(i, j-1, k, ibg, u_bc.south, rheology, clock, fields)
q̃ᴺ = ib_ice_stress_uy(i, j, k, ibg, u_bc.north, rheology, clock, fields)
q̃ᵂ = ib_ice_stress_ux_west(i, j, k, ibg, u_bc.west, rheology, clock, fields)
q̃ᴱ = ib_ice_stress_ux_east(i+1, j, k, ibg, u_bc.east, rheology, clock, fields)
q̃ˢ = ib_ice_stress_uy_south(i, j-1, k, ibg, u_bc.south, rheology, clock, fields)
q̃ᴺ = ib_ice_stress_uy_north(i, j, k, ibg, u_bc.north, rheology, clock, fields)

iᵂ, jˢ, _ = map(index_left, (i, j, k), (c, c, c)) # Broadcast instead of map causes inference failure
iᴱ, jᴺ, _ = map(index_right, (i, j, k), (f, f, c))

# Impose i) immersed fluxes if we're on an immersed boundary or ii) zero otherwise.
qᵂ = conditional_flux_ccc(iᵂ, j, k, ibg, q̃ᵂ, zero(ibg)) * Axᶜᶜᶜ(iᵂ, j, k, grid)
qᴱ = conditional_flux_ccc(iᴱ, j, k, ibg, q̃ᴱ, zero(ibg)) * Axᶜᶜᶜ(iᴱ, j, k, grid)
qˢ = conditional_flux_ffc(i, jˢ, k, ibg, q̃ˢ, zero(ibg)) * Ayᶠᶠᶜ(i, jˢ, k, grid)
qᴺ = conditional_flux_ffc(i, jᴺ, k, ibg, q̃ᴺ, zero(ibg)) * Ayᶠᶠᶜ(i, jᴺ, k, grid)
# Impose i) immersed fluxes if we're on an immersed boundary or ii) zero oibgwise.
qᵂ = conditional_flux_ccc(iᵂ, j, k, ibg, q̃ᵂ, zero(ibg)) * Axᶜᶜᶜ(iᵂ, j, k, ibg)
qᴱ = conditional_flux_ccc(iᴱ, j, k, ibg, q̃ᴱ, zero(ibg)) * Axᶜᶜᶜ(iᴱ, j, k, ibg)
qˢ = conditional_flux_ffc(i, jˢ, k, ibg, q̃ˢ, zero(ibg)) * Ayᶠᶠᶜ(i, jˢ, k, ibg)
qᴺ = conditional_flux_ffc(i, jᴺ, k, ibg, q̃ᴺ, zero(ibg)) * Ayᶠᶠᶜ(i, jᴺ, k, ibg)

return (qᴱ - qᵂ + qᴺ - qˢ) / Vᶠᶜᶜ(i, j, k, grid)
return (qᴱ - qᵂ + qᴺ - qˢ) / Vᶠᶜᶜ(i, j, k, ibg)
end

@inline function immersed_∂ⱼ_σ₂ⱼ(i, j, k, ibg::IBG, v_bc::IBC, rheology, clock, fields)
# Fetch fluxes across immersed boundary
q̃ᵂ = ib_ice_stress_vx(i-1, j, k, ibg, v_bc.west, rheology, clock, fields)
q̃ᴱ = ib_ice_stress_vx(i, j, k, ibg, v_bc.east, rheology, clock, fields)
q̃ˢ = ib_ice_stress_vy(i, j, k, ibg, v_bc.south, rheology, clock, fields)
q̃ᴺ = ib_ice_stress_vy(i, j+1, k, ibg, v_bc.north, rheology, clock, fields)
q̃ᵂ = ib_ice_stress_vx_west(i-1, j, k, ibg, v_bc.west, rheology, clock, fields)
q̃ᴱ = ib_ice_stress_vx_east(i, j, k, ibg, v_bc.east, rheology, clock, fields)
q̃ˢ = ib_ice_stress_vy_south(i, j-1, k, ibg, v_bc.south, rheology, clock, fields)
q̃ᴺ = ib_ice_stress_vy_north(i, j, k, ibg, v_bc.north, rheology, clock, fields)

iᵂ, jˢ, _ = map(index_left, (i, j, k), (f, f, c)) # Broadcast instead of map causes inference failure
iᴱ, jᴺ, _ = map(index_right, (i, j, k), (c, c, c))

# Impose i) immersed fluxes if we're on an immersed boundary or ii) zero otherwise.
qᵂ = conditional_flux_ffc(iᵂ, j, k, ibg, q̃ᵂ, zero(ibg)) * Axᶠᶠᶜ(iᵂ, j, k, grid)
qᴱ = conditional_flux_ffc(iᴱ, j, k, ibg, q̃ᴱ, zero(ibg)) * Axᶠᶠᶜ(iᴱ, j, k, grid)
qˢ = conditional_flux_ccc(i, jˢ, k, ibg, q̃ˢ, zero(ibg)) * Ayᶜᶜᶜ(i, jˢ, k, grid)
qᴺ = conditional_flux_ccc(i, jᴺ, k, ibg, q̃ᴺ, zero(ibg)) * Ayᶜᶜᶜ(i, jᴺ, k, grid)
qᵂ = conditional_flux_ffc(iᵂ, j, k, ibg, q̃ᵂ, zero(ibg)) * Axᶠᶠᶜ(iᵂ, j, k, ibg)
qᴱ = conditional_flux_ffc(iᴱ, j, k, ibg, q̃ᴱ, zero(ibg)) * Axᶠᶠᶜ(iᴱ, j, k, ibg)
qˢ = conditional_flux_ccc(i, jˢ, k, ibg, q̃ˢ, zero(ibg)) * Ayᶜᶜᶜ(i, jˢ, k, ibg)
qᴺ = conditional_flux_ccc(i, jᴺ, k, ibg, q̃ᴺ, zero(ibg)) * Ayᶜᶜᶜ(i, jᴺ, k, ibg)

return (qᴱ - qᵂ + qᴺ - qˢ) / Vᶜᶠᶜ(i, j, k, grid)
return (qᴱ - qᵂ + qᴺ - qˢ) / Vᶜᶠᶜ(i, j, k, ibg)
end

# TODO: Implement immersed fluxes (0 for the moment)
@inline ib_ice_stress_ux(i, j, k, grid, args...) = zero(grid)
@inline ib_ice_stress_vx(i, j, k, grid, args...) = zero(grid)
@inline ib_ice_stress_uy(i, j, k, grid, args...) = zero(grid)
@inline ib_ice_stress_vy(i, j, k, grid, args...) = zero(grid)
@inline ib_ice_stress_ux_west(i, j, k, ibg, args...) = zero(ibg)
@inline ib_ice_stress_ux_east(i, j, k, ibg, args...) = zero(ibg)
@inline ib_ice_stress_vx_west(i, j, k, ibg, args...) = zero(ibg)
@inline ib_ice_stress_vx_east(i, j, k, ibg, args...) = zero(ibg)
@inline ib_ice_stress_uy_south(i, j, k, ibg, args...) = zero(ibg)
@inline ib_ice_stress_uy_north(i, j, k, ibg, args...) = zero(ibg)
@inline ib_ice_stress_vy_south(i, j, k, ibg, args...) = zero(ibg)
@inline ib_ice_stress_vy_north(i, j, k, ibg, args...) = zero(ibg)

# Only supporting FluxBoundaryConditions for now
@inline ib_ice_stress_ux_west(i, j, k, ibg, bc::FBC, rheology, args...) = + getbc(bc, i, j, k, ibg, args...)
@inline ib_ice_stress_ux_east(i, j, k, ibg, bc::FBC, rheology, args...) = - getbc(bc, i, j, k, ibg, args...)
@inline ib_ice_stress_vx_west(i, j, k, ibg, bc::FBC, rheology, args...) = + getbc(bc, i, j, k, ibg, args...)
@inline ib_ice_stress_vx_east(i, j, k, ibg, bc::FBC, rheology, args...) = - getbc(bc, i, j, k, ibg, args...)
@inline ib_ice_stress_uy_south(i, j, k, ibg, bc::FBC, rheology, args...) = + getbc(bc, i, j, k, ibg, args...)
@inline ib_ice_stress_uy_north(i, j, k, ibg, bc::FBC, rheology, args...) = - getbc(bc, i, j, k, ibg, args...)
@inline ib_ice_stress_vy_south(i, j, k, ibg, bc::FBC, rheology, args...) = + getbc(bc, i, j, k, ibg, args...)
@inline ib_ice_stress_vy_north(i, j, k, ibg, bc::FBC, rheology, args...) = - getbc(bc, i, j, k, ibg, args...)
Loading