Skip to content

Commit

Permalink
let's go
Browse files Browse the repository at this point in the history
  • Loading branch information
@simone-silvestri
simone-silvestri committed Dec 15, 2024
1 parent 5724ee6 commit 2364025
Show file tree
Hide file tree
Showing 3 changed files with 93 additions and 19 deletions.
10 changes: 5 additions & 5 deletions src/Advection/vector_invariant_advection.jl
View file
Original file line number Diff line number Diff line change
@@ -1,6 +1,6 @@
using Oceananigans.Operators
using Oceananigans.Operators: flux_div_xyᶜᶜᶜ, Γᶠᶠᶜ
using Oceananigans.Operators: mask_inactive_points_ℑxyᶠᶜᵃ, mask_inactive_points_ℑxyᶜᶠᵃ
using Oceananigans.Operators: ℑxyMᶠᶜᵃ, ℑxyMᶜᶠᵃ

# These are also used in Coriolis/hydrostatic_spherical_coriolis.jl
struct EnergyConserving{FT} <: AbstractAdvectionScheme{1, FT} end
Expand Down Expand Up @@ -372,7 +372,7 @@ end

Sζ = scheme.vorticity_stencil

@inbounds v̂ = mask_inactive_points_ℑxyᶠᶜᵃ(i, j, k, grid, Δx_qᶜᶠᶜ, v) / Δxᶠᶜᶜ(i, j, k, grid)
@inbounds v̂ = ℑxyMᶠᶜᶜ(i, j, k, grid, Δx_qᶜᶠᶜ, v) / Δxᶠᶜᶜ(i, j, k, grid)
ζᴿ = _biased_interpolate_yᵃᶜᵃ(i, j, k, grid, scheme, scheme.vorticity_scheme, bias(v̂), ζ₃ᶠᶠᶜ, Sζ, u, v)

return - v̂ * ζᴿ
Expand All @@ -382,7 +382,7 @@ end

Sζ = scheme.vorticity_stencil

@inbounds û = mask_inactive_points_ℑxyᶜᶠᵃ(i, j, k, grid, Δy_qᶠᶜᶜ, u) / Δyᶜᶠᶜ(i, j, k, grid)
@inbounds û = ℑxyMᶜᶠᶜ(i, j, k, grid, Δy_qᶠᶜᶜ, u) / Δyᶜᶠᶜ(i, j, k, grid)
ζᴿ = _biased_interpolate_xᶜᵃᵃ(i, j, k, grid, scheme, scheme.vorticity_scheme, bias(û), ζ₃ᶠᶠᶜ, Sζ, u, v)

return + û * ζᴿ
Expand All @@ -394,7 +394,7 @@ end

@inline function U_dot_∇u(i, j, k, grid, advection::AbstractAdvectionScheme, U)

v̂ = ℑxᶠᵃᵃ(i, j, k, grid, ℑyᵃᶜᵃ, Δx_qᶜᶠᶜ, U.v) / Δxᶠᶜᶜ(i, j, k, grid)
v̂ = ℑxyMᶠᶜᶜ(i, j, k, grid, Δx_qᶜᶠᶜ, U.v) / Δxᶠᶜᶜ(i, j, k, grid)
û = @inbounds U.u[i, j, k]

return div_𝐯u(i, j, k, grid, advection, U, U.u) -
Expand All @@ -404,7 +404,7 @@ end

@inline function U_dot_∇v(i, j, k, grid, advection::AbstractAdvectionScheme, U)

û = ℑyᵃᶠᵃ(i, j, k, grid, ℑxᶜᵃᵃ, Δy_qᶠᶜᶜ, U.u) / Δyᶜᶠᶜ(i, j, k, grid)
û = ℑxyMᶜᶠᶜ(i, j, k, grid, Δy_qᶠᶜᶜ, U.u) / Δyᶜᶠᶜ(i, j, k, grid)
v̂ = @inbounds U.v[i, j, k]

return div_𝐯v(i, j, k, grid, advection, U, U.v) +
Expand Down
6 changes: 3 additions & 3 deletions src/Coriolis/hydrostatic_spherical_coriolis.jl
View file
Original file line number Diff line number Diff line change
@@ -1,5 +1,5 @@
using Oceananigans.Grids: LatitudeLongitudeGrid, OrthogonalSphericalShellGrid, peripheral_node, φnode
using Oceananigans.Operators: Δx_qᶜᶠᶜ, Δy_qᶠᶜᶜ, Δxᶠᶜᶜ, Δyᶜᶠᶜ, hack_sind, mask_inactive_points_ℑxyᶠᶜᵃ, mask_inactive_points_ℑxyᶜᶠᵃ
using Oceananigans.Operators: Δx_qᶜᶠᶜ, Δy_qᶠᶜᶜ, Δxᶠᶜᶜ, Δyᶜᶠᶜ, hack_sind, ℑxyMᶠᶜᶜ, ℑxyMᶜᶠᶜ
using Oceananigans.Advection: EnergyConserving, EnstrophyConserving
using Oceananigans.BoundaryConditions
using Oceananigans.Fields
Expand Down Expand Up @@ -69,11 +69,11 @@ const CoriolisActiveCellEnstrophyConserving = HydrostaticSphericalCoriolis{<:Act

@inline x_f_cross_U(i, j, k, grid, coriolis::CoriolisActiveCellEnstrophyConserving, U) =
@inbounds - ℑyᵃᶜᵃ(i, j, k, grid, fᶠᶠᵃ, coriolis) *
mask_inactive_points_ℑxyᶠᶜᵃ(i, j, k, grid, Δx_qᶜᶠᶜ, U[2]) / Δxᶠᶜᶜ(i, j, k, grid)
ℑxyMᶠᶜᶜ(i, j, k, grid, Δx_qᶜᶠᶜ, U[2]) / Δxᶠᶜᶜ(i, j, k, grid)

@inline y_f_cross_U(i, j, k, grid, coriolis::CoriolisActiveCellEnstrophyConserving, U) =
@inbounds + ℑxᶜᵃᵃ(i, j, k, grid, fᶠᶠᵃ, coriolis) *
mask_inactive_points_ℑxyᶜᶠᵃ(i, j, k, grid, Δy_qᶠᶜᶜ, U[1]) / Δyᶜᶠᶜ(i, j, k, grid)
ℑxyMᶜᶠᶜ(i, j, k, grid, Δy_qᶠᶜᶜ, U[1]) / Δyᶜᶠᶜ(i, j, k, grid)

#####
##### Enstrophy-conserving scheme
Expand Down
96 changes: 85 additions & 11 deletions src/Operators/topology_aware_operators.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -66,18 +66,92 @@ const AGYL = AbstractUnderlyingGrid{FT, <:Any, LeftConnected} where FT
@inline ∂xTᶠᶜᶠ(i, j, k, grid, f::Function, args...) = δxTᶠᵃᵃ(i, j, k, grid, f, args...) / Δxᶠᶜᶠ(i, j, k, grid)
@inline ∂yTᶜᶠᶠ(i, j, k, grid, f::Function, args...) = δyTᵃᶠᵃ(i, j, k, grid, f, args...) / Δyᶜᶠᶠ(i, j, k, grid)

# Masking interpolation operators

####
#### Masking interpolation operators
####

# Interpolation operators that mask inactive points.
# They assume that inactive points are zero.

@inline not_peripheral_node(args...) = !peripheral_node(args...)

@inline function mask_inactive_points_ℑxyᶠᶜᵃ(i, j, k, grid, f::Function, args...)
neighboring_active_nodes = ℑxyᶠᶜᵃ(i, j, k, grid, not_peripheral_node, Center(), Face(), Center())
return ifelse(neighboring_active_nodes == 0, zero(grid),
ℑxyᶠᶜᵃ(i, j, k, grid, f, args...) / neighboring_active_nodes)
@inline flip(::Type{Center}) = Face
@inline flip(::Type{Face}) = Center

for LX in (:Center, :Face), LY in (:Center, :Face), LZ in (:Center, :Face)
LXe = @eval $LX
LYe = @eval $LY
LZe = @eval $LZ

LXf = flip(LXe)
LYf = flip(LYe)
LZf = flip(LZe)

ℑxˡᵃᵃ = Symbol(:ℑx, location_code(LXe, nothing, nothing))
ℑyᵃˡᵃ = Symbol(:ℑy, location_code(nothing, LYe, nothing))
ℑzᵃᵃˡ = Symbol(:ℑz, location_code(nothing, nothing, LZe))

ℑxMˡˡˡ = Symbol(:ℑxM, location_code(LXe, LYe, LZe))
ℑyMˡˡˡ = Symbol(:ℑyM, location_code(LXe, LYe, LZe))
ℑzMˡˡˡ = Symbol(:ℑzM, location_code(LXe, LYe, LZe))

@eval begin
@inline function $ℑxMˡˡˡ(i, j, k, grid, f::Function, args...)
neighboring_active_nodes = $ℑxˡᵃᵃ(i, j, k, grid, not_peripheral_node, $LXf(), $LYe(), $LZe())
return ifelse(neighboring_active_nodes == 0, zero(grid),
$ℑxˡᵃᵃ(i, j, k, grid, f, args...) / neighboring_active_nodes)
end

@inline function $ℑyMˡˡˡ(i, j, k, grid, f::Function, args...)
neighboring_active_nodes = $ℑyᵃˡᵃ(i, j, k, grid, not_peripheral_node, $LXe(), $LYf(), $LZe())
return ifelse(neighboring_active_nodes == 0, zero(grid),
$ℑyᵃˡᵃ(i, j, k, grid, f, args...) / neighboring_active_nodes)
end

@inline function $ℑzMˡˡˡ(i, j, k, grid, f::Function, args...)
neighboring_active_nodes = $ℑzᵃᵃˡ(i, j, k, grid, not_peripheral_node, $LXe(), $LYe(), $LZf())
return ifelse(neighboring_active_nodes == 0, zero(grid),
$ℑzᵃᵃˡ(i, j, k, grid, f, args...) / neighboring_active_nodes)
end
end

ℑxyˡˡᵃ = Symbol(:ℑxy, location_code(LXe, LYe, nothing))
ℑyzᵃˡˡ = Symbol(:ℑyz, location_code(nothing, LYe, LZe))
ℑxzˡᵃˡ = Symbol(:ℑxz, location_code(LXe, nothing, LZe))

ℑxyMˡˡˡ = Symbol(:ℑxyM, location_code(LXe, LYe, LZe))
ℑyzMˡˡˡ = Symbol(:ℑyzM, location_code(LXe, LYe, LZe))
ℑxzMˡˡˡ = Symbol(:ℑxzM, location_code(LXe, LYe, LZe))

@eval begin
@inline function $ℑxyMˡˡˡ(i, j, k, grid, f::Function, args...)
neighboring_active_nodes = $ℑxyˡˡᵃ(i, j, k, grid, not_peripheral_node, $LXf(), $LYf(), $LZe())
return ifelse(neighboring_active_nodes == 0, zero(grid),
$ℑxyˡˡᵃ(i, j, k, grid, f, args...) / neighboring_active_nodes)
end

@inline function $ℑyzMˡˡˡ(i, j, k, grid, f::Function, args...)
neighboring_active_nodes = $ℑyzᵃˡˡ(i, j, k, grid, not_peripheral_node, $LXe(), $LYf(), $LZf())
return ifelse(neighboring_active_nodes == 0, zero(grid),
$ℑyzᵃˡˡ(i, j, k, grid, f, args...) / neighboring_active_nodes)
end

@inline function $ℑxzMˡˡˡ(i, j, k, grid, f::Function, args...)
neighboring_active_nodes = $ℑxzˡᵃˡ(i, j, k, grid, not_peripheral_node, $LXf(), $LYe(), $LZf())
return ifelse(neighboring_active_nodes == 0, zero(grid),
$ℑxzˡᵃˡ(i, j, k, grid, f, args...) / neighboring_active_nodes)
end
end

ℑxyzˡˡˡ = Symbol(:ℑxyz, location_code(LXe, LYe, LZe))
ℑxyzMˡˡˡ = Symbol(:ℑxyzM, location_code(LXe, LYe, LZe))

@eval begin
@inline function $ℑxyzMˡˡˡ(i, j, k, grid, f::Function, args...)
neighboring_active_nodes = $ℑxyzˡˡˡ(i, j, k, grid, not_peripheral_node, $LXf(), $LYf(), $LZf())
return ifelse(neighboring_active_nodes == 0, zero(grid),
$ℑxyzˡˡˡ(i, j, k, grid, f, args...) / neighboring_active_nodes)
end
end
end

@inline function mask_inactive_points_ℑxyᶜᶠᵃ(i, j, k, grid, f::Function, args...)
neighboring_active_nodes = @inbounds ℑxyᶜᶠᵃ(i, j, k, grid, not_peripheral_node, Face(), Center(), Center())
return ifelse(neighboring_active_nodes == 0, zero(grid),
ℑxyᶜᶠᵃ(i, j, k, grid, f, args...) / neighboring_active_nodes)
end

0 comments on commit 2364025

Please sign in to comment.