Skip to content

Commit

Permalink
Different advection schemes in different directions (#3732)
Browse files Browse the repository at this point in the history
* should be done

* add the dicussion

* disambiguation

* remove repeat code

* precopmiles

* keep these changes for later

* correct the outside buffer

* fix tests

* change the show method

* back to previous buffer scheme

* allow for changes in the advection

* switch position of functions

* better comment

* better info

* bugfix

* better info

* make sure N==1 -> nothing

* Update src/Advection/momentum_advection_operators.jl

Co-authored-by: Gregory L. Wagner <[email protected]>

* add topologies

* add topology argument

* better adaptation

* formatting

* whoops

* remove topology and less metaprogramming

* add advection

* dont adapt nothing

* retain same behavior as before

* fix multi region advection

* bugfix scalar diffusivities parameter types

* bugfix

* enzyme error

* make sure constructors are type-stable...

* force required_halo_size to be an integer

* add better comment

* add a couple of unit tests

* test biharmonic diffusivity

* fix enzyme tests

* correct adaptbetter explanation

* add a test for adapt_advection

* whoops wrong test

* revert changes

* add couple more tests

* fix tests

* without the extra space

* fix enzyme test

* fix enzyme tests

* fix nonhydrostatic tests

* rmove one newline

* simplify test

* Update src/Models/HydrostaticFreeSurfaceModels/hydrostatic_free_surface_model.jl

Co-authored-by: Gregory L. Wagner <[email protected]>

* remove extra dispatches

* relaunch the tests

* change comment

* switch to B

* new syntax for KernelParameters

* bugfix

* another bugfix

* make sure the range is correct

* Update src/Advection/adapt_advection_order.jl

Co-authored-by: Gregory L. Wagner <[email protected]>

* Update src/Advection/adapt_advection_order.jl

Co-authored-by: Gregory L. Wagner <[email protected]>

* Update src/Advection/adapt_advection_order.jl

Co-authored-by: Gregory L. Wagner <[email protected]>

* remove comments

* this should be type stable for the moment

* bugfix

* bugfix

* better comment

* fix the docs

* whoops

---------

Co-authored-by: Gregory L. Wagner <[email protected]>
  • Loading branch information
simone-silvestri and glwagner authored Oct 11, 2024
1 parent 989108a commit 0210e66
Show file tree
Hide file tree
Showing 24 changed files with 487 additions and 219 deletions.
11 changes: 8 additions & 3 deletions src/Advection/Advection.jl
Original file line number Diff line number Diff line change
Expand Up @@ -21,7 +21,7 @@ export
UpwindBiased, UpwindBiasedFirstOrder, UpwindBiasedThirdOrder, UpwindBiasedFifthOrder,
WENO, WENOThirdOrder, WENOFifthOrder,
VectorInvariant, WENOVectorInvariant,
TracerAdvection,
FluxFormAdvection,
EnergyConserving,
EnstrophyConserving

Expand All @@ -38,7 +38,7 @@ using Oceananigans.Architectures: architecture, CPU
using Oceananigans.Operators

import Base: show, summary
import Oceananigans.Grids: required_halo_size
import Oceananigans.Grids: required_halo_size_x, required_halo_size_y, required_halo_size_z
import Oceananigans.Architectures: on_architecture

abstract type AbstractAdvectionScheme{B, FT} end
Expand All @@ -55,9 +55,12 @@ abstract type AbstractUpwindBiasedAdvectionScheme{B, FT} <: AbstractAdvectionSch
# Note that it is not possible to compile schemes for `advection_buffer = 41` or higher.
const advection_buffers = [1, 2, 3, 4, 5, 6]

@inline required_halo_size(::AbstractAdvectionScheme{B}) where B = B
@inline Base.eltype(::AbstractAdvectionScheme{<:Any, FT}) where FT = FT

@inline required_halo_size_x(::AbstractAdvectionScheme{B}) where B = B
@inline required_halo_size_y(::AbstractAdvectionScheme{B}) where B = B
@inline required_halo_size_z(::AbstractAdvectionScheme{B}) where B = B

include("centered_advective_fluxes.jl")
include("upwind_biased_advective_fluxes.jl")

Expand All @@ -72,12 +75,14 @@ include("vector_invariant_upwinding.jl")
include("vector_invariant_advection.jl")
include("vector_invariant_self_upwinding.jl")
include("vector_invariant_cross_upwinding.jl")
include("flux_form_advection.jl")

include("flat_advective_fluxes.jl")
include("topologically_conditional_interpolation.jl")
include("momentum_advection_operators.jl")
include("tracer_advection_operators.jl")
include("positivity_preserving_tracer_advection_operators.jl")
include("cell_advection_timescale.jl")
include("adapt_advection_order.jl")

end # module
98 changes: 98 additions & 0 deletions src/Advection/adapt_advection_order.jl
Original file line number Diff line number Diff line change
@@ -0,0 +1,98 @@
using Oceananigans.Grids: topology

"""
adapt_advection_order(advection, grid::AbstractGrid)
Adapts the advection operator `advection` based on the grid `grid` by adjusting the order of advection in each direction.
For example, if the grid has only one point in the x-direction, the advection operator in the x-direction is set to first order
upwind or 2nd order centered scheme, depending on the original user-specified advection scheme. A high order advection sheme
is reduced to a lower order advection scheme if the grid has fewer points in that direction.
# Arguments
- `advection`: The original advection scheme.
- `grid::AbstractGrid`: The grid on which the advection scheme is applied.
If the order of advection is changed in at least one direction, the adapted advection scheme with adjusted advection order returned
by this function is a `FluxFormAdvection`.
"""
function adapt_advection_order(advection, grid::AbstractGrid)
advection_x = x_advection(advection)
advection_y = y_advection(advection)
advection_z = z_advection(advection)

new_advection_x = adapt_advection_order(advection_x, size(grid, 1), grid)
new_advection_y = adapt_advection_order(advection_y, size(grid, 2), grid)
new_advection_z = adapt_advection_order(advection_z, size(grid, 3), grid)

# Check that we indeed changed the advection operator
changed_x = new_advection_x != advection_x
changed_y = new_advection_y != advection_y
changed_z = new_advection_z != advection_z

new_advection = FluxFormAdvection(new_advection_x, new_advection_y, new_advection_z)
changed_advection = any((changed_x, changed_y, changed_z))

if changed_x
@info "Using the advection scheme $(summary(new_advection.x)) in the x-direction because size(grid, 1) = $(size(grid, 1))"
end
if changed_y
@info "Using the advection scheme $(summary(new_advection.y)) in the y-direction because size(grid, 2) = $(size(grid, 2))"
end
if changed_z
@info "Using the advection scheme $(summary(new_advection.z)) in the x-direction because size(grid, 3) = $(size(grid, 3))"
end

return ifelse(changed_advection, new_advection, advection)
end


x_advection(flux_form::FluxFormAdvection) = flux_form.x
y_advection(flux_form::FluxFormAdvection) = flux_form.y
z_advection(flux_form::FluxFormAdvection) = flux_form.z

x_advection(advection) = advection
y_advection(advection) = advection
z_advection(advection) = advection

# For the moment, we do not adapt the advection order for the VectorInvariant advection scheme
adapt_advection_order(advection::VectorInvariant, grid::AbstractGrid) = advection
adapt_advection_order(advection::Nothing, grid::AbstractGrid) = nothing
adapt_advection_order(advection::Nothing, N::Int, grid::AbstractGrid) = nothing

#####
##### Directional adapt advection order
#####

function adapt_advection_order(advection::Centered{B}, N::Int, grid::AbstractGrid) where B
if N >= B
return advection
else
return Centered(; order = 2N)
end
end

function adapt_advection_order(advection::UpwindBiased{B}, N::Int, grid::AbstractGrid) where B
if N >= B
return advection
else
return UpwindBiased(; order = 2N - 1)
end
end

"""
new_weno_scheme(grid, order, bounds, XT, YT, ZT)
Constructs a new WENO scheme based on the given parameters. `XT`, `YT`, and `ZT` is the type of the precomputed weno coefficients in the
x-direction, y-direction and z-direction. A _non-stretched_ WENO scheme has `T` equal to `Nothing` everywhere. In case of a non-stretched WENO scheme,
we rebuild the advection without passing the grid information, otherwise we use the grid to account for stretched directions.
"""
new_weno_scheme(::WENO, grid, order, bounds, ::Type{Nothing}, ::Type{Nothing}, ::Type{Nothing},) = WENO(; order, bounds)
new_weno_scheme(::WENO, grid, order, bounds, XT, YT, ZT) = WENO(grid; order, bounds)

function adapt_advection_order(advection::WENO{B, FT, XT, YT, ZT}, N::Int, grid::AbstractGrid) where {B, FT, XT, YT, ZT}
if N >= B
return advection
else
return new_weno_scheme(advection, grid, 2N - 1, advection.bounds, XT, YT, ZT)
end
end
58 changes: 58 additions & 0 deletions src/Advection/flux_form_advection.jl
Original file line number Diff line number Diff line change
@@ -0,0 +1,58 @@
using Oceananigans.Operators: Vᶜᶜᶜ
using Oceananigans.Fields: ZeroField

struct FluxFormAdvection{N, FT, A, B, C} <: AbstractAdvectionScheme{N, FT}
x :: A
y :: B
z :: C

FluxFormAdvection{N, FT}(x::A, y::B, z::C) where {N, FT, A, B, C} = new{N, FT, A, B, C}(x, y, z)
end

"""
function FluxFormAdvection(x, y, z)
Builds a `FluxFormAdvection` type with reconstructions schemes `x`, `y`, and `z` to be applied in
the x, y, and z direction, respectively.
"""
function FluxFormAdvection(x_advection, y_advection, z_advection)
Hx = required_halo_size_x(x_advection)
Hy = required_halo_size_y(y_advection)
Hz = required_halo_size_z(z_advection)

FT = eltype(x_advection)
H = max(Hx, Hy, Hz)

return FluxFormAdvection{H, FT}(x_advection, y_advection, z_advection)
end

Base.show(io::IO, scheme::FluxFormAdvection) =
print(io, "FluxFormAdvection with reconstructions: ", " \n",
" ├── x: ", summary(scheme.x), "\n",
" ├── y: ", summary(scheme.y), "\n",
" └── z: ", summary(scheme.z))

@inline required_halo_size_x(scheme::FluxFormAdvection) = required_halo_size_x(scheme.x)
@inline required_halo_size_y(scheme::FluxFormAdvection) = required_halo_size_y(scheme.y)
@inline required_halo_size_z(scheme::FluxFormAdvection) = required_halo_size_z(scheme.z)

Adapt.adapt_structure(to, scheme::FluxFormAdvection{N, FT}) where {N, FT} =
FluxFormAdvection{N, FT}(Adapt.adapt(to, scheme.x),
Adapt.adapt(to, scheme.y),
Adapt.adapt(to, scheme.z))

@inline _advective_tracer_flux_x(i, j, k, grid, advection::FluxFormAdvection, args...) = _advective_tracer_flux_x(i, j, k, grid, advection.x, args...)
@inline _advective_tracer_flux_y(i, j, k, grid, advection::FluxFormAdvection, args...) = _advective_tracer_flux_y(i, j, k, grid, advection.y, args...)
@inline _advective_tracer_flux_z(i, j, k, grid, advection::FluxFormAdvection, args...) = _advective_tracer_flux_z(i, j, k, grid, advection.z, args...)

@inline _advective_momentum_flux_Uu(i, j, k, grid, advection::FluxFormAdvection, args...) = _advective_momentum_flux_Uu(i, j, k, grid, advection.x, args...)
@inline _advective_momentum_flux_Vu(i, j, k, grid, advection::FluxFormAdvection, args...) = _advective_momentum_flux_Vu(i, j, k, grid, advection.y, args...)
@inline _advective_momentum_flux_Wu(i, j, k, grid, advection::FluxFormAdvection, args...) = _advective_momentum_flux_Wu(i, j, k, grid, advection.z, args...)

@inline _advective_momentum_flux_Uv(i, j, k, grid, advection::FluxFormAdvection, args...) = _advective_momentum_flux_Uv(i, j, k, grid, advection.x, args...)
@inline _advective_momentum_flux_Vv(i, j, k, grid, advection::FluxFormAdvection, args...) = _advective_momentum_flux_Vv(i, j, k, grid, advection.y, args...)
@inline _advective_momentum_flux_Wv(i, j, k, grid, advection::FluxFormAdvection, args...) = _advective_momentum_flux_Wv(i, j, k, grid, advection.z, args...)

@inline _advective_momentum_flux_Uw(i, j, k, grid, advection::FluxFormAdvection, args...) = _advective_momentum_flux_Uw(i, j, k, grid, advection.x, args...)
@inline _advective_momentum_flux_Vw(i, j, k, grid, advection::FluxFormAdvection, args...) = _advective_momentum_flux_Vw(i, j, k, grid, advection.y, args...)
@inline _advective_momentum_flux_Ww(i, j, k, grid, advection::FluxFormAdvection, args...) = _advective_momentum_flux_Ww(i, j, k, grid, advection.z, args...)
31 changes: 20 additions & 11 deletions src/Advection/momentum_advection_operators.jl
Original file line number Diff line number Diff line change
Expand Up @@ -28,17 +28,9 @@ const ZeroU = NamedTuple{(:u, :v, :w), Tuple{ZeroField, ZeroField, ZeroField}}
@inline div_𝐯v(i, j, k, grid, advection, U, ::ZeroField) = zero(grid)
@inline div_𝐯w(i, j, k, grid, advection, U, ::ZeroField) = zero(grid)

@inline div_𝐯u(i, j, k, grid, ::Nothing, U, u) = zero(grid)
@inline div_𝐯v(i, j, k, grid, ::Nothing, U, v) = zero(grid)
@inline div_𝐯w(i, j, k, grid, ::Nothing, U, w) = zero(grid)

@inline div_𝐯u(i, j, k, grid, ::Nothing, ::ZeroU, u) = zero(grid)
@inline div_𝐯v(i, j, k, grid, ::Nothing, ::ZeroU, v) = zero(grid)
@inline div_𝐯w(i, j, k, grid, ::Nothing, ::ZeroU, w) = zero(grid)

@inline div_𝐯u(i, j, k, grid, ::Nothing, U, ::ZeroField) = zero(grid)
@inline div_𝐯v(i, j, k, grid, ::Nothing, U, ::ZeroField) = zero(grid)
@inline div_𝐯w(i, j, k, grid, ::Nothing, U, ::ZeroField) = zero(grid)
@inline div_𝐯u(i, j, k, grid, advection, ::ZeroU, ::ZeroField) = zero(grid)
@inline div_𝐯v(i, j, k, grid, advection, ::ZeroU, ::ZeroField) = zero(grid)
@inline div_𝐯w(i, j, k, grid, advection, ::ZeroU, ::ZeroField) = zero(grid)

"""
div_𝐯u(i, j, k, grid, advection, U, u)
Expand Down Expand Up @@ -89,3 +81,20 @@ which ends up at the location `ccf`.
δyᵃᶜᵃ(i, j, k, grid, _advective_momentum_flux_Vw, advection, U[2], w) +
δzᵃᵃᶠ(i, j, k, grid, _advective_momentum_flux_Ww, advection, U[3], w))
end

#####
##### Fallback advection fluxes!
#####

# Fallback for `nothing` advection
@inline _advective_momentum_flux_Uu(i, j, k, grid, ::Nothing, args...) = zero(grid)
@inline _advective_momentum_flux_Uv(i, j, k, grid, ::Nothing, args...) = zero(grid)
@inline _advective_momentum_flux_Uw(i, j, k, grid, ::Nothing, args...) = zero(grid)

@inline _advective_momentum_flux_Vu(i, j, k, grid, ::Nothing, args...) = zero(grid)
@inline _advective_momentum_flux_Vv(i, j, k, grid, ::Nothing, args...) = zero(grid)
@inline _advective_momentum_flux_Vw(i, j, k, grid, ::Nothing, args...) = zero(grid)

@inline _advective_momentum_flux_Wu(i, j, k, grid, ::Nothing, args...) = zero(grid)
@inline _advective_momentum_flux_Wv(i, j, k, grid, ::Nothing, args...) = zero(grid)
@inline _advective_momentum_flux_Ww(i, j, k, grid, ::Nothing, args...) = zero(grid)
27 changes: 18 additions & 9 deletions src/Advection/topologically_conditional_interpolation.jl
Original file line number Diff line number Diff line change
Expand Up @@ -26,14 +26,23 @@ const AUGXYZ = AUG{<:Any, <:Bounded, <:Bounded, <:Bounded}
# Left-biased buffers are smaller by one grid point on the right side; vice versa for right-biased buffers
# Center interpolation stencil look at i + 1 (i.e., require one less point on the left)

@inline outside_symmetric_haloᶠ(i, N, adv) = (i >= required_halo_size(adv) + 1) & (i <= N + 1 - required_halo_size(adv))
@inline outside_symmetric_haloᶜ(i, N, adv) = (i >= required_halo_size(adv)) & (i <= N + 1 - required_halo_size(adv))

@inline outside_biased_haloᶠ(i, N, adv) = (i >= required_halo_size(adv) + 1) & (i <= N + 1 - (required_halo_size(adv) - 1)) & # Left bias
(i >= required_halo_size(adv)) & (i <= N + 1 - required_halo_size(adv)) # Right bias
@inline outside_biased_haloᶜ(i, N, adv) = (i >= required_halo_size(adv)) & (i <= N + 1 - (required_halo_size(adv) - 1)) & # Left bias
(i >= required_halo_size(adv) - 1) & (i <= N + 1 - required_halo_size(adv)) # Right bias

for dir in (:x, :y, :z)
outside_symmetric_haloᶠ = Symbol(:outside_symmetric_halo_, dir, :ᶠ)
outside_symmetric_haloᶜ = Symbol(:outside_symmetric_halo_, dir, :ᶜ)
outside_biased_haloᶠ = Symbol(:outside_biased_halo_, dir, :ᶠ)
outside_biased_haloᶜ = Symbol(:outside_biased_halo_, dir, :ᶜ)
required_halo_size = Symbol(:required_halo_size_, dir)

@eval begin
@inline $outside_symmetric_haloᶠ(i, N, adv) = (i >= $required_halo_size(adv) + 1) & (i <= N + 1 - $required_halo_size(adv))
@inline $outside_symmetric_haloᶜ(i, N, adv) = (i >= $required_halo_size(adv)) & (i <= N + 1 - $required_halo_size(adv))

@inline $outside_biased_haloᶠ(i, N, adv) = (i >= $required_halo_size(adv) + 1) & (i <= N + 1 - ($required_halo_size(adv) - 1)) & # Left bias
(i >= $required_halo_size(adv)) & (i <= N + 1 - $required_halo_size(adv)) # Right bias
@inline $outside_biased_haloᶜ(i, N, adv) = (i >= $required_halo_size(adv)) & (i <= N + 1 - ($required_halo_size(adv) - 1)) & # Left bias
(i >= $required_halo_size(adv) - 1) & (i <= N + 1 - $required_halo_size(adv)) # Right bias
end
end
# Separate High order advection from low order advection
const HOADV = Union{WENO,
Tuple(Centered{N} for N in advection_buffers[2:end])...,
Expand Down Expand Up @@ -61,7 +70,7 @@ for bias in (:symmetric, :biased)
@eval @inline $alt1_interp(i, j, k, grid::$GridType, scheme::LOADV, args...) = $interp(i, j, k, grid, scheme, args...)
end

outside_buffer = Symbol(:outside_, bias, :_halo, loc)
outside_buffer = Symbol(:outside_, bias, :_halo_, ξ, loc)

# Conditional high-order interpolation in Bounded directions
if ξ == :x
Expand Down
72 changes: 5 additions & 67 deletions src/Advection/tracer_advection_operators.jl
Original file line number Diff line number Diff line change
@@ -1,81 +1,17 @@
using Oceananigans.Operators: Vᶜᶜᶜ
using Oceananigans.Fields: ZeroField

struct TracerAdvection{N, FT, A, B, C} <: AbstractAdvectionScheme{N, FT}
x :: A
y :: B
z :: C

TracerAdvection{N, FT}(x::A, y::B, z::C) where {N, FT, A, B, C} = new{N, FT, A, B, C}(x, y, z)
end

"""
function TracerAdvection(x, y, z)
Builds a `TracerAdvection` type with reconstructions schemes `x`, `y`, and `z` to be applied in
the x, y, and z direction, respectively.
"""
function TracerAdvection(x_advection, y_advection, z_advection)
Hx = required_halo_size(x_advection)
Hy = required_halo_size(y_advection)
Hz = required_halo_size(z_advection)

FT = eltype(x_advection)
H = max(Hx, Hy, Hz)

return TracerAdvection{H, FT}(x_advection, y_advection, z_advection)
end

Adapt.adapt_structure(to, scheme::TracerAdvection{N, FT}) where {N, FT} =
TracerAdvection{N, FT}(Adapt.adapt(to, scheme.x),
Adapt.adapt(to, scheme.y),
Adapt.adapt(to, scheme.z))

@inline _advective_tracer_flux_x(args...) = advective_tracer_flux_x(args...)
@inline _advective_tracer_flux_y(args...) = advective_tracer_flux_y(args...)
@inline _advective_tracer_flux_z(args...) = advective_tracer_flux_z(args...)

@inline _advective_tracer_flux_x(i, j, k, grid, advection::TracerAdvection, args...) =
_advective_tracer_flux_x(i, j, k, grid, advection.x, args...)

@inline _advective_tracer_flux_y(i, j, k, grid, advection::TracerAdvection, args...) =
_advective_tracer_flux_y(i, j, k, grid, advection.y, args...)

@inline _advective_tracer_flux_z(i, j, k, grid, advection::TracerAdvection, args...) =
_advective_tracer_flux_z(i, j, k, grid, advection.z, args...)
#####
##### Fallback tracer fluxes!
#####

# Fallback for `nothing` advection
@inline _advective_tracer_flux_x(i, j, k, grid, ::Nothing, args...) = zero(grid)
@inline _advective_tracer_flux_y(i, j, k, grid, ::Nothing, args...) = zero(grid)
@inline _advective_tracer_flux_z(i, j, k, grid, ::Nothing, args...) = zero(grid)

# Fallback for `nothing` advection and `ZeroField` tracers and velocities
@inline _advective_tracer_flux_x(i, j, k, grid, ::Nothing, ::ZeroField, ::ZeroField) = zero(grid)
@inline _advective_tracer_flux_y(i, j, k, grid, ::Nothing, ::ZeroField, ::ZeroField) = zero(grid)
@inline _advective_tracer_flux_z(i, j, k, grid, ::Nothing, ::ZeroField, ::ZeroField) = zero(grid)

@inline _advective_tracer_flux_x(i, j, k, grid, ::Nothing, U, ::ZeroField) = zero(grid)
@inline _advective_tracer_flux_y(i, j, k, grid, ::Nothing, V, ::ZeroField) = zero(grid)
@inline _advective_tracer_flux_z(i, j, k, grid, ::Nothing, W, ::ZeroField) = zero(grid)
@inline _advective_tracer_flux_x(i, j, k, grid, ::Nothing, ::ZeroField, c) = zero(grid)
@inline _advective_tracer_flux_y(i, j, k, grid, ::Nothing, ::ZeroField, c) = zero(grid)
@inline _advective_tracer_flux_z(i, j, k, grid, ::Nothing, ::ZeroField, c) = zero(grid)

# Fallback for `ZeroField` tracers and velocities
@inline _advective_tracer_flux_x(i, j, k, grid, scheme, ::ZeroField, ::ZeroField) = zero(grid)
@inline _advective_tracer_flux_y(i, j, k, grid, scheme, ::ZeroField, ::ZeroField) = zero(grid)
@inline _advective_tracer_flux_z(i, j, k, grid, scheme, ::ZeroField, ::ZeroField) = zero(grid)

# Fallback for `ZeroField` tracers
@inline _advective_tracer_flux_x(i, j, k, grid, scheme, U, ::ZeroField) = zero(grid)
@inline _advective_tracer_flux_y(i, j, k, grid, scheme, V, ::ZeroField) = zero(grid)
@inline _advective_tracer_flux_z(i, j, k, grid, scheme, W, ::ZeroField) = zero(grid)

# Fallback for `ZeroField` velocities
@inline _advective_tracer_flux_x(i, j, k, grid, scheme, ::ZeroField, c) = zero(grid)
@inline _advective_tracer_flux_y(i, j, k, grid, scheme, ::ZeroField, c) = zero(grid)
@inline _advective_tracer_flux_z(i, j, k, grid, scheme, ::ZeroField, c) = zero(grid)

#####
##### Tracer advection operator
#####
Expand All @@ -100,7 +36,9 @@ end
# Fallbacks for zero velocities, zero tracer and `nothing` advection
@inline div_Uc(i, j, k, grid, advection, ::ZeroU, c) = zero(grid)
@inline div_Uc(i, j, k, grid, advection, U, ::ZeroField) = zero(grid)
@inline div_Uc(i, j, k, grid, advection, ::ZeroU, ::ZeroField) = zero(grid)

@inline div_Uc(i, j, k, grid, ::Nothing, U, c) = zero(grid)
@inline div_Uc(i, j, k, grid, ::Nothing, ::ZeroU, c) = zero(grid)
@inline div_Uc(i, j, k, grid, ::Nothing, U, ::ZeroField) = zero(grid)
@inline div_Uc(i, j, k, grid, ::Nothing, ::ZeroU, ::ZeroField) = zero(grid)
Loading

0 comments on commit 0210e66

Please sign in to comment.