From ea309b72dbf23404c9fd05a6072edd5d7f9d2d8b Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Sun, 25 Aug 2024 20:39:03 -0400 Subject: [PATCH 01/68] should be done --- src/Advection/Advection.jl | 10 ++- src/Advection/flux_form_advection.jl | 52 ++++++++++++++ src/Advection/momentum_advection_operators.jl | 17 +++++ ...topologically_conditional_interpolation.jl | 23 ++++--- src/Advection/tracer_advection_operators.jl | 41 ----------- src/Advection/vector_invariant_advection.jl | 22 +++++- src/Grids/automatic_halo_sizing.jl | 68 ++++++++++++++++--- src/ImmersedBoundaries/conditional_fluxes.jl | 8 +-- ...static_free_surface_boundary_tendencies.jl | 2 - ...pute_nonhydrostatic_boundary_tendencies.jl | 13 ++-- src/Oceananigans.jl | 2 +- src/TurbulenceClosures/TurbulenceClosures.jl | 6 +- src/TurbulenceClosures/closure_tuples.jl | 4 +- .../scalar_biharmonic_diffusivity.jl | 18 ++++- .../scalar_diffusivity.jl | 16 ++++- 15 files changed, 216 insertions(+), 86 deletions(-) create mode 100644 src/Advection/flux_form_advection.jl diff --git a/src/Advection/Advection.jl b/src/Advection/Advection.jl index c54b310475..3fed3242de 100644 --- a/src/Advection/Advection.jl +++ b/src/Advection/Advection.jl @@ -21,7 +21,7 @@ export UpwindBiased, UpwindBiasedFirstOrder, UpwindBiasedThirdOrder, UpwindBiasedFifthOrder, WENO, WENOThirdOrder, WENOFifthOrder, VectorInvariant, WENOVectorInvariant, - TracerAdvection, + FluxFormAdvection, EnergyConserving, EnstrophyConserving @@ -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 @@ -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") @@ -72,6 +75,7 @@ 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") diff --git a/src/Advection/flux_form_advection.jl b/src/Advection/flux_form_advection.jl new file mode 100644 index 0000000000..58369fc19f --- /dev/null +++ b/src/Advection/flux_form_advection.jl @@ -0,0 +1,52 @@ +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 + +@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...) \ No newline at end of file diff --git a/src/Advection/momentum_advection_operators.jl b/src/Advection/momentum_advection_operators.jl index a2aa5c29ed..78b24dea8c 100644 --- a/src/Advection/momentum_advection_operators.jl +++ b/src/Advection/momentum_advection_operators.jl @@ -89,3 +89,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! +##### + +for flux_dir in (:Uu, :Vu, :Wu, :Uv, :Vv, :Wv, :Uw, :Vw, :Ww) + advective_momentum_flux = Symbol(:_advective_momentum_flux_, flux_dir) + + @eval begin + @inline $advective_momentum_flux(i, j, k, grid, ::Nothing, args...) = zero(grid) + @inline $advective_momentum_flux(i, j, k, grid, ::Nothing, ::ZeroField, ::ZeroField) = zero(grid) + @inline $advective_momentum_flux(i, j, k, grid, ::Nothing, U, ::ZeroField) = zero(grid) + @inline $advective_momentum_flux(i, j, k, grid, scheme, ::ZeroField, ::ZeroField) = zero(grid) + @inline $advective_momentum_flux(i, j, k, grid, scheme, U, ::ZeroField) = zero(grid) + @inline $advective_momentum_flux(i, j, k, grid, scheme, ::ZeroField, u) = zero(grid) + end +end \ No newline at end of file diff --git a/src/Advection/topologically_conditional_interpolation.jl b/src/Advection/topologically_conditional_interpolation.jl index 094f5a7253..37f98a76c7 100644 --- a/src/Advection/topologically_conditional_interpolation.jl +++ b/src/Advection/topologically_conditional_interpolation.jl @@ -26,14 +26,21 @@ 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 direction in (:x, :y, :z) + outside_symmetric_haloᶠ = Symbol(:outside_symmetric_halo_, dir, :ᶠ) + outside_symmetric_haloᶜ = Symbol(:outside_symmetric_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])..., diff --git a/src/Advection/tracer_advection_operators.jl b/src/Advection/tracer_advection_operators.jl index fcae1923ce..4808aa60e9 100644 --- a/src/Advection/tracer_advection_operators.jl +++ b/src/Advection/tracer_advection_operators.jl @@ -1,49 +1,8 @@ -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 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) diff --git a/src/Advection/vector_invariant_advection.jl b/src/Advection/vector_invariant_advection.jl index a9a89b2623..e16cd7438a 100644 --- a/src/Advection/vector_invariant_advection.jl +++ b/src/Advection/vector_invariant_advection.jl @@ -116,7 +116,14 @@ function VectorInvariant(; vorticity_scheme = EnstrophyConserving(), upwinding = OnlySelfUpwinding(; cross_scheme = divergence_scheme), multi_dimensional_stencil = false) - N = required_halo_size(vorticity_scheme) + N = max(required_halo_size_x(vorticity_scheme), + required_halo_size_y(vorticity_scheme), + required_halo_size_x(divergence_scheme), + required_halo_size_y(divergence_scheme), + required_halo_size_x(kinetic_energy_gradient_scheme), + required_halo_size_y(kinetic_energy_gradient_scheme), + required_halo_size_z(vertical_scheme)) + FT = eltype(vorticity_scheme) return VectorInvariant{N, FT, multi_dimensional_stencil}(vorticity_scheme, @@ -216,11 +223,20 @@ function WENOVectorInvariant(FT::DataType = Float64; upwinding) end - # Since vorticity itself requires one halo, if we use an upwinding scheme (N > 1) we require one additional # halo for vector invariant advection -required_halo_size(scheme::VectorInvariant{N}) where N = N == 1 ? N : N + 1 +@inline function required_halo_size_x(scheme::VectorInvariant) + Hx₁ = required_halo_size_x(scheme.vorticity_scheme) + Hx₂ = required_halo_size_x(scheme.divergence_scheme) + Hx₃ = required_halo_size_x(scheme.kinetic_energy_gradient_scheme) + + Hx = max(Hx₁, Hx₂, Hx₃) + return Hx == 1 ? Hx : Hx + 1 +end +@inline required_halo_size_y(scheme::VectorInvariant) = required_halo_size_x(scheme) +@inline required_halo_size_z(scheme::VectorInvariant) = required_halo_size_z(scheme.vertical_scheme) + Adapt.adapt_structure(to, scheme::VectorInvariant{N, FT, M}) where {N, FT, M} = VectorInvariant{N, FT, M}(Adapt.adapt(to, scheme.vorticity_scheme), Adapt.adapt(to, scheme.vorticity_stencil), diff --git a/src/Grids/automatic_halo_sizing.jl b/src/Grids/automatic_halo_sizing.jl index a4fd19b526..262728efa1 100644 --- a/src/Grids/automatic_halo_sizing.jl +++ b/src/Grids/automatic_halo_sizing.jl @@ -1,7 +1,7 @@ """ - required_halo_size(tendency_term) + required_halo_size_x(tendency_term) -Return the required size of halos for a term appearing +Return the required size of halos in the x direction for a term appearing in the tendency for a velocity field or tracer field. Example @@ -9,17 +9,61 @@ Example ```jldoctest using Oceananigans.Advection: CenteredFourthOrder -using Oceananigans.Grids: required_halo_size +using Oceananigans.Grids: required_halo_size_x -required_halo_size(CenteredFourthOrder()) +required_halo_size_x(CenteredFourthOrder()) # output 2 """ -function required_halo_size end +function required_halo_size_x end -required_halo_size(tendency_term) = 1 -required_halo_size(::Nothing) = 0 +""" + required_halo_size_y(tendency_term) + +Return the required size of halos in the y direction for a term appearing +in the tendency for a velocity field or tracer field. + +Example +======= + +```jldoctest +using Oceananigans.Advection: CenteredFourthOrder +using Oceananigans.Grids: required_halo_size_y + +required_halo_size_y(CenteredFourthOrder()) + +# output +2 +""" +function required_halo_size_y end + +""" + required_halo_size_z(tendency_term) + +Return the required size of halos in the y direction for a term appearing +in the tendency for a velocity field or tracer field. + +Example +======= + +```jldoctest +using Oceananigans.Advection: CenteredFourthOrder +using Oceananigans.Grids: required_halo_size_z + +required_halo_size_z(CenteredFourthOrder()) + +# output +2 +""" +function required_halo_size_z end + +required_halo_size_x(tendency_term) = 1 +required_halo_size_x(::Nothing) = 0 +required_halo_size_y(tendency_term) = 1 +required_halo_size_y(::Nothing) = 0 +required_halo_size_z(tendency_term) = 1 +required_halo_size_z(::Nothing) = 0 inflate_halo_size_one_dimension(req_H, old_H, _, grid) = max(req_H, old_H) inflate_halo_size_one_dimension(req_H, old_H, ::Type{Flat}, grid) = 0 @@ -27,10 +71,12 @@ inflate_halo_size_one_dimension(req_H, old_H, ::Type{Flat}, grid) = 0 function inflate_halo_size(Hx, Hy, Hz, grid, tendency_terms...) topo = topology(grid) for term in tendency_terms - H = required_halo_size(term) - Hx = inflate_halo_size_one_dimension(H, Hx, topo[1], grid) - Hy = inflate_halo_size_one_dimension(H, Hy, topo[2], grid) - Hz = inflate_halo_size_one_dimension(H, Hz, topo[3], grid) + Hx_required = required_halo_size_x(term) + Hy_required = required_halo_size_y(term) + Hz_required = required_halo_size_z(term) + Hx = inflate_halo_size_one_dimension(Hx_required, Hx, topo[1], grid) + Hy = inflate_halo_size_one_dimension(Hy_required, Hy, topo[2], grid) + Hz = inflate_halo_size_one_dimension(Hz_required, Hz, topo[3], grid) end return Hx, Hy, Hz diff --git a/src/ImmersedBoundaries/conditional_fluxes.jl b/src/ImmersedBoundaries/conditional_fluxes.jl index 260e397c75..02d55dd674 100644 --- a/src/ImmersedBoundaries/conditional_fluxes.jl +++ b/src/ImmersedBoundaries/conditional_fluxes.jl @@ -1,7 +1,7 @@ using Oceananigans.Advection: AbstractAdvectionScheme, advection_buffers using Oceananigans.Operators: ℑxᶠᵃᵃ, ℑxᶜᵃᵃ, ℑyᵃᶠᵃ, ℑyᵃᶜᵃ, ℑzᵃᵃᶠ, ℑzᵃᵃᶜ using Oceananigans.TurbulenceClosures: AbstractTurbulenceClosure, AbstractTimeDiscretization -using Oceananigans.Advection: LOADV, HOADV, WENO, TracerAdvection +using Oceananigans.Advection: LOADV, HOADV, WENO, FluxFormAdvection using Oceananigans.Fields: ZeroField const ATC = AbstractTurbulenceClosure @@ -83,13 +83,13 @@ end @inline _advective_tracer_flux_z(i, j, k, ibg::IBG, args...) = conditional_flux_ccf(i, j, k, ibg, zero(ibg), advective_tracer_flux_z(i, j, k, ibg, args...)) # Disambiguation for tracer fluxes.... -@inline _advective_tracer_flux_x(i, j, k, ibg::IBG, advection::TracerAdvection, args...) = +@inline _advective_tracer_flux_x(i, j, k, ibg::IBG, advection::FluxFormAdvection, args...) = _advective_tracer_flux_x(i, j, k, ibg, advection.x, args...) -@inline _advective_tracer_flux_y(i, j, k, ibg::IBG, advection::TracerAdvection, args...) = +@inline _advective_tracer_flux_y(i, j, k, ibg::IBG, advection::FluxFormAdvection, args...) = _advective_tracer_flux_y(i, j, k, ibg, advection.y, args...) -@inline _advective_tracer_flux_z(i, j, k, ibg::IBG, advection::TracerAdvection, args...) = +@inline _advective_tracer_flux_z(i, j, k, ibg::IBG, advection::FluxFormAdvection, args...) = _advective_tracer_flux_z(i, j, k, ibg, advection.z, args...) # Fallback for `nothing` advection diff --git a/src/Models/HydrostaticFreeSurfaceModels/compute_hydrostatic_free_surface_boundary_tendencies.jl b/src/Models/HydrostaticFreeSurfaceModels/compute_hydrostatic_free_surface_boundary_tendencies.jl index ca7214e82c..4096cf0972 100644 --- a/src/Models/HydrostaticFreeSurfaceModels/compute_hydrostatic_free_surface_boundary_tendencies.jl +++ b/src/Models/HydrostaticFreeSurfaceModels/compute_hydrostatic_free_surface_boundary_tendencies.jl @@ -8,8 +8,6 @@ using Oceananigans.Models.NonhydrostaticModels: boundary_tendency_kernel_paramet boundary_κ_kernel_parameters, boundary_parameters -using Oceananigans.TurbulenceClosures: required_halo_size - # We assume here that top/bottom BC are always synchronized (no partitioning in z) function compute_boundary_tendencies!(model::HydrostaticFreeSurfaceModel) grid = model.grid diff --git a/src/Models/NonhydrostaticModels/compute_nonhydrostatic_boundary_tendencies.jl b/src/Models/NonhydrostaticModels/compute_nonhydrostatic_boundary_tendencies.jl index 933f732d8d..9fb4e44c00 100644 --- a/src/Models/NonhydrostaticModels/compute_nonhydrostatic_boundary_tendencies.jl +++ b/src/Models/NonhydrostaticModels/compute_nonhydrostatic_boundary_tendencies.jl @@ -1,6 +1,6 @@ import Oceananigans.Models: compute_boundary_tendencies! -using Oceananigans.TurbulenceClosures: required_halo_size +using Oceananigans.TurbulenceClosures: equired_halo_size_x, required_halo_size_y using Oceananigans.Grids: XFlatGrid, YFlatGrid # TODO: the code in this file is difficult to understand. @@ -64,15 +64,16 @@ end function boundary_κ_kernel_parameters(grid, closure, arch) Nx, Ny, Nz = size(grid) - B = required_halo_size(closure) + Bx = required_halo_size_x(closure) + By = required_halo_size_y(closure) - Sx = (B+1, Ny, Nz) - Sy = (Nx, B+1, Nz) + Sx = (Bx+1, Ny, Nz) + Sy = (Nx, By+1, Nz) Oxᴸ = (-1, 0, 0) Oyᴸ = (0, -1, 0) - Oxᴿ = (Nx-B, 0, 0) - Oyᴿ = (0, Ny-B, 0) + Oxᴿ = (Nx-Bx, 0, 0) + Oyᴿ = (0, Ny-By, 0) sizes = (Sx, Sy, Sx, Sy) offs = (Oxᴸ, Oyᴸ, Oxᴿ, Oyᴿ) diff --git a/src/Oceananigans.jl b/src/Oceananigans.jl index b91c83b7fd..b1b9f3d73a 100644 --- a/src/Oceananigans.jl +++ b/src/Oceananigans.jl @@ -34,7 +34,7 @@ export UpwindBiased, UpwindBiasedFirstOrder, UpwindBiasedThirdOrder, UpwindBiasedFifthOrder, WENO, WENOThirdOrder, WENOFifthOrder, VectorInvariant, WENOVectorInvariant, EnergyConserving, EnstrophyConserving, - TracerAdvection, + FluxFormAdvection, # Boundary conditions BoundaryCondition, diff --git a/src/TurbulenceClosures/TurbulenceClosures.jl b/src/TurbulenceClosures/TurbulenceClosures.jl index f810ce5e2f..701a7b1927 100644 --- a/src/TurbulenceClosures/TurbulenceClosures.jl +++ b/src/TurbulenceClosures/TurbulenceClosures.jl @@ -50,7 +50,7 @@ using Oceananigans.Utils using Oceananigans.Architectures: AbstractArchitecture, device using Oceananigans.Fields: FunctionField -import Oceananigans.Advection: required_halo_size +import Oceananigans.Advection: required_halo_size_x, required_halo_size_y, required_halo_size_z import Oceananigans.Architectures: on_architecture const VerticallyBoundedGrid{FT} = AbstractGrid{FT, <:Any, <:Any, <:Bounded} @@ -77,7 +77,9 @@ compute_diffusivities!(K, closure::AbstractTurbulenceClosure, args...; kwargs... # point at each side to calculate viscous fluxes at the edge of the domain. # If diffusivity itself requires one halo to be computed (e.g. κ = ℑxᶠᵃᵃ(i, j, k, grid, ℑxᶜᵃᵃ, T), # or `AnisotropicMinimumDissipation` and `SmagorinskyLilly`) then B = 2 -@inline required_halo_size(::AbstractTurbulenceClosure{TD, B}) where {TD, B} = B +@inline required_halo_size_x(::AbstractTurbulenceClosure{TD, B}) where {TD, B} = B +@inline required_halo_size_y(::AbstractTurbulenceClosure{TD, B}) where {TD, B} = B +@inline required_halo_size_z(::AbstractTurbulenceClosure{TD, B}) where {TD, B} = B const ClosureKinda = Union{Nothing, AbstractTurbulenceClosure, AbstractArray{<:AbstractTurbulenceClosure}} add_closure_specific_boundary_conditions(closure::ClosureKinda, bcs, args...) = bcs diff --git a/src/TurbulenceClosures/closure_tuples.jl b/src/TurbulenceClosures/closure_tuples.jl index 6321800b14..b552331bb0 100644 --- a/src/TurbulenceClosures/closure_tuples.jl +++ b/src/TurbulenceClosures/closure_tuples.jl @@ -89,7 +89,9 @@ function add_closure_specific_boundary_conditions(closure_tuple::Tuple, bcs, arg return bcs end -required_halo_size(closure_tuple::Tuple) = maximum(map(required_halo_size, closure_tuple)) +required_halo_size_x(closure_tuple::Tuple) = maximum(map(required_halo_size_x, closure_tuple)) +required_halo_size_y(closure_tuple::Tuple) = maximum(map(required_halo_size_y, closure_tuple)) +required_halo_size_z(closure_tuple::Tuple) = maximum(map(required_halo_size_z, closure_tuple)) ##### ##### Compiler-inferrable time_discretization for tuples diff --git a/src/TurbulenceClosures/turbulence_closure_implementations/scalar_biharmonic_diffusivity.jl b/src/TurbulenceClosures/turbulence_closure_implementations/scalar_biharmonic_diffusivity.jl index 5d5056461c..070c30fb5b 100644 --- a/src/TurbulenceClosures/turbulence_closure_implementations/scalar_biharmonic_diffusivity.jl +++ b/src/TurbulenceClosures/turbulence_closure_implementations/scalar_biharmonic_diffusivity.jl @@ -1,4 +1,4 @@ -import Oceananigans.Grids: required_halo_size +import Oceananigans.Grids: required_halo_size_x, required_halo_size_y, required_halo_size_z using Oceananigans.Utils: prettysummary """ @@ -117,3 +117,19 @@ function on_architecture(to, closure::ScalarBiharmonicDiffusivity{F, <:Any, <:An κ = on_architecture(to, closure.κ) return ScalarBiharmonicDiffusivity{F, N}(ν, κ) end + +required_halo_size_x(::ScalarBiharmonicDiffusivity{<:Any, ThreeDimensionalFormulation, N}) where N = N +required_halo_size_y(::ScalarBiharmonicDiffusivity{<:Any, ThreeDimensionalFormulation, N}) where N = N +required_halo_size_z(::ScalarBiharmonicDiffusivity{<:Any, ThreeDimensionalFormulation, N}) where N = N + +required_halo_size_x(::ScalarBiharmonicDiffusivity{<:Any, VerticalFormulation, N}) where N = 0 +required_halo_size_y(::ScalarBiharmonicDiffusivity{<:Any, VerticalFormulation, N}) where N = 0 +required_halo_size_z(::ScalarBiharmonicDiffusivity{<:Any, VerticalFormulation, N}) where N = N + +required_halo_size_x(::ScalarBiharmonicDiffusivity{<:Any, HorizontalFormulation, N}) where N = N +required_halo_size_y(::ScalarBiharmonicDiffusivity{<:Any, HorizontalFormulation, N}) where N = N +required_halo_size_z(::ScalarBiharmonicDiffusivity{<:Any, HorizontalFormulation, N}) where N = 0 + +required_halo_size_x(::ScalarBiharmonicDiffusivity{<:Any, HorizontalDivergenceFormulation, N}) where N = N +required_halo_size_y(::ScalarBiharmonicDiffusivity{<:Any, HorizontalDivergenceFormulation, N}) where N = N +required_halo_size_z(::ScalarBiharmonicDiffusivity{<:Any, HorizontalDivergenceFormulation, N}) where N = 0 diff --git a/src/TurbulenceClosures/turbulence_closure_implementations/scalar_diffusivity.jl b/src/TurbulenceClosures/turbulence_closure_implementations/scalar_diffusivity.jl index d17f1eb72e..8fdd8be46f 100644 --- a/src/TurbulenceClosures/turbulence_closure_implementations/scalar_diffusivity.jl +++ b/src/TurbulenceClosures/turbulence_closure_implementations/scalar_diffusivity.jl @@ -1,7 +1,7 @@ using Oceananigans.Utils: prettysummary import Adapt -import Oceananigans.Grids: required_halo_size +import Oceananigans.Grids: required_halo_size_x, required_halo_size_y, required_halo_size_z struct ScalarDiffusivity{TD, F, V, K, N} <: AbstractScalarDiffusivity{TD, F, N} ν :: V @@ -173,8 +173,18 @@ Shorthand for a `ScalarDiffusivity` with `HorizontalDivergenceFormulation()`. Se HorizontalScalarDiffusivity(FT::DataType; kwargs...) = ScalarDiffusivity(ExplicitTimeDiscretization(), HorizontalFormulation(), FT; kwargs...) HorizontalDivergenceScalarDiffusivity(FT::DataType; kwargs...) = ScalarDiffusivity(ExplicitTimeDiscretization(), HorizontalDivergenceFormulation(), FT; kwargs...) -required_halo_size(closure::ScalarDiffusivity) = 1 - +required_halo_size_x(::ScalarDiffusivity{<:Any, ThreeDimensionalFormulation, N}) where N = N +required_halo_size_y(::ScalarDiffusivity{<:Any, ThreeDimensionalFormulation, N}) where N = N +required_halo_size_z(::ScalarDiffusivity{<:Any, ThreeDimensionalFormulation, N}) where N = N + +required_halo_size_x(::ScalarDiffusivity{<:Any, VerticalFormulation, N}) where N = 0 +required_halo_size_y(::ScalarDiffusivity{<:Any, VerticalFormulation, N}) where N = 0 +required_halo_size_z(::ScalarDiffusivity{<:Any, VerticalFormulation, N}) where N = N + +required_halo_size_x(::ScalarDiffusivity{<:Any, HorizontalFormulation, N}) where N = N +required_halo_size_y(::ScalarDiffusivity{<:Any, HorizontalFormulation, N}) where N = N +required_halo_size_z(::ScalarDiffusivity{<:Any, HorizontalFormulation, N}) where N = 0 + @inline function with_tracers(tracers, closure::ScalarDiffusivity{TD, F, N}) where {TD, F, N} κ = tracer_diffusivities(tracers, closure.κ) return ScalarDiffusivity{TD, F, N}(closure.ν, κ) From 83f870e11100980b95fb3652b0590bb07f4a783a Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Sun, 25 Aug 2024 21:06:29 -0400 Subject: [PATCH 02/68] add the dicussion --- .../scalar_biharmonic_diffusivity.jl | 6 +++++ .../scalar_diffusivity.jl | 22 +++++++++++++++++++ 2 files changed, 28 insertions(+) diff --git a/src/TurbulenceClosures/turbulence_closure_implementations/scalar_biharmonic_diffusivity.jl b/src/TurbulenceClosures/turbulence_closure_implementations/scalar_biharmonic_diffusivity.jl index 070c30fb5b..677b4ebbc3 100644 --- a/src/TurbulenceClosures/turbulence_closure_implementations/scalar_biharmonic_diffusivity.jl +++ b/src/TurbulenceClosures/turbulence_closure_implementations/scalar_biharmonic_diffusivity.jl @@ -118,6 +118,11 @@ function on_architecture(to, closure::ScalarBiharmonicDiffusivity{F, <:Any, <:An return ScalarBiharmonicDiffusivity{F, N}(ν, κ) end +# We do not have at the moment the ability to distinguish between the halos required in the x, y and z +# direction for a general closure since the closure could contain any user-defined function + +# In the end we probably want something like this: +#= required_halo_size_x(::ScalarBiharmonicDiffusivity{<:Any, ThreeDimensionalFormulation, N}) where N = N required_halo_size_y(::ScalarBiharmonicDiffusivity{<:Any, ThreeDimensionalFormulation, N}) where N = N required_halo_size_z(::ScalarBiharmonicDiffusivity{<:Any, ThreeDimensionalFormulation, N}) where N = N @@ -133,3 +138,4 @@ required_halo_size_z(::ScalarBiharmonicDiffusivity{<:Any, HorizontalFormulation, required_halo_size_x(::ScalarBiharmonicDiffusivity{<:Any, HorizontalDivergenceFormulation, N}) where N = N required_halo_size_y(::ScalarBiharmonicDiffusivity{<:Any, HorizontalDivergenceFormulation, N}) where N = N required_halo_size_z(::ScalarBiharmonicDiffusivity{<:Any, HorizontalDivergenceFormulation, N}) where N = 0 +=# \ No newline at end of file diff --git a/src/TurbulenceClosures/turbulence_closure_implementations/scalar_diffusivity.jl b/src/TurbulenceClosures/turbulence_closure_implementations/scalar_diffusivity.jl index 8fdd8be46f..00473f28fa 100644 --- a/src/TurbulenceClosures/turbulence_closure_implementations/scalar_diffusivity.jl +++ b/src/TurbulenceClosures/turbulence_closure_implementations/scalar_diffusivity.jl @@ -229,3 +229,25 @@ function on_architecture(to, closure::ScalarDiffusivity{TD, F, <:Any, <:Any, N}) κ = on_architecture(to, closure.κ) return ScalarDiffusivity{TD, F, N}(ν, κ) end + +# We do not have at the moment the ability to distinguish between the halos required in the x, y and z +# direction for a general closure since the closure could contain any user-defined function + +# In the end we probably want something like this: +#= +required_halo_size_x(::ScalarDiffusivity{<:Any, ThreeDimensionalFormulation, N}) where N = N +required_halo_size_y(::ScalarDiffusivity{<:Any, ThreeDimensionalFormulation, N}) where N = N +required_halo_size_z(::ScalarDiffusivity{<:Any, ThreeDimensionalFormulation, N}) where N = N + +required_halo_size_x(::ScalarDiffusivity{<:Any, VerticalFormulation, N}) where N = 0 +required_halo_size_y(::ScalarDiffusivity{<:Any, VerticalFormulation, N}) where N = 0 +required_halo_size_z(::ScalarDiffusivity{<:Any, VerticalFormulation, N}) where N = N + +required_halo_size_x(::ScalarDiffusivity{<:Any, HorizontalFormulation, N}) where N = N +required_halo_size_y(::ScalarDiffusivity{<:Any, HorizontalFormulation, N}) where N = N +required_halo_size_z(::ScalarDiffusivity{<:Any, HorizontalFormulation, N}) where N = 0 + +required_halo_size_x(::ScalarDiffusivity{<:Any, HorizontalDivergenceFormulation, N}) where N = N +required_halo_size_y(::ScalarDiffusivity{<:Any, HorizontalDivergenceFormulation, N}) where N = N +required_halo_size_z(::ScalarDiffusivity{<:Any, HorizontalDivergenceFormulation, N}) where N = 0 +=# \ No newline at end of file From ad1ae414c4e76b0f01f9dcb577af392462e81ab5 Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Sun, 25 Aug 2024 21:10:37 -0400 Subject: [PATCH 03/68] disambiguation --- src/Advection/tracer_advection_operators.jl | 56 ++++++------- src/ImmersedBoundaries/conditional_fluxes.jl | 88 +++++++++++++------- 2 files changed, 84 insertions(+), 60 deletions(-) diff --git a/src/Advection/tracer_advection_operators.jl b/src/Advection/tracer_advection_operators.jl index 4808aa60e9..44719bfbc9 100644 --- a/src/Advection/tracer_advection_operators.jl +++ b/src/Advection/tracer_advection_operators.jl @@ -3,37 +3,35 @@ @inline _advective_tracer_flux_y(args...) = advective_tracer_flux_y(args...) @inline _advective_tracer_flux_z(args...) = advective_tracer_flux_z(args...) -# 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 tracer fluxes! +##### -# 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) +for flux_dir in (:x, :y, :z) + advective_tracer_flux = Symbol(:_advective_tracer_flux_, flux_dir) + + @eval begin + @inline $advective_tracer_flux(i, j, k, grid, ::Nothing, args...) = zero(grid) + @inline $advective_tracer_flux(i, j, k, grid, ::Nothing, ::ZeroField, ::ZeroField) = zero(grid) + @inline $advective_tracer_flux(i, j, k, grid, ::Nothing, U, ::ZeroField) = zero(grid) + @inline $advective_tracer_flux(i, j, k, grid, scheme, ::ZeroField, ::ZeroField) = zero(grid) + @inline $advective_tracer_flux(i, j, k, grid, scheme, U, ::ZeroField) = zero(grid) + @inline $advective_tracer_flux(i, j, k, grid, scheme, ::ZeroField, c) = zero(grid) + end +end -# 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) +for flux_dir in (:Uu, :Vu, :Wu, :Uv, :Vv, :Wv, :Uw, :Vw, :Ww) + advective_momentum_flux = Symbol(:_advective_momentum_flux_, flux_dir) + + @eval begin + @inline $advective_momentum_flux(i, j, k, grid, ::Nothing, args...) = zero(grid) + @inline $advective_momentum_flux(i, j, k, grid, ::Nothing, ::ZeroField, ::ZeroField) = zero(grid) + @inline $advective_momentum_flux(i, j, k, grid, ::Nothing, U, ::ZeroField) = zero(grid) + @inline $advective_momentum_flux(i, j, k, grid, scheme, ::ZeroField, ::ZeroField) = zero(grid) + @inline $advective_momentum_flux(i, j, k, grid, scheme, U, ::ZeroField) = zero(grid) + @inline $advective_momentum_flux(i, j, k, grid, scheme, ::ZeroField, u) = zero(grid) + end +end ##### ##### Tracer advection operator diff --git a/src/ImmersedBoundaries/conditional_fluxes.jl b/src/ImmersedBoundaries/conditional_fluxes.jl index 02d55dd674..03be5fb0f9 100644 --- a/src/ImmersedBoundaries/conditional_fluxes.jl +++ b/src/ImmersedBoundaries/conditional_fluxes.jl @@ -92,37 +92,63 @@ end @inline _advective_tracer_flux_z(i, j, k, ibg::IBG, advection::FluxFormAdvection, args...) = _advective_tracer_flux_z(i, j, k, ibg, advection.z, args...) -# Fallback for `nothing` advection -@inline _advective_tracer_flux_x(i, j, k, ibg::IBG, ::Nothing, args...) = zero(ibg) -@inline _advective_tracer_flux_y(i, j, k, ibg::IBG, ::Nothing, args...) = zero(ibg) -@inline _advective_tracer_flux_z(i, j, k, ibg::IBG, ::Nothing, args...) = zero(ibg) - -# Fallback for `nothing` advection and `ZeroField` tracers and velocities -@inline _advective_tracer_flux_x(i, j, k, ibg::IBG, ::Nothing, ::ZeroField, ::ZeroField) = zero(ibg) -@inline _advective_tracer_flux_y(i, j, k, ibg::IBG, ::Nothing, ::ZeroField, ::ZeroField) = zero(ibg) -@inline _advective_tracer_flux_z(i, j, k, ibg::IBG, ::Nothing, ::ZeroField, ::ZeroField) = zero(ibg) - -@inline _advective_tracer_flux_x(i, j, k, ibg::IBG, ::Nothing, U, ::ZeroField) = zero(ibg) -@inline _advective_tracer_flux_y(i, j, k, ibg::IBG, ::Nothing, V, ::ZeroField) = zero(ibg) -@inline _advective_tracer_flux_z(i, j, k, ibg::IBG, ::Nothing, W, ::ZeroField) = zero(ibg) -@inline _advective_tracer_flux_x(i, j, k, ibg::IBG, ::Nothing, ::ZeroField, c) = zero(ibg) -@inline _advective_tracer_flux_y(i, j, k, ibg::IBG, ::Nothing, ::ZeroField, c) = zero(ibg) -@inline _advective_tracer_flux_z(i, j, k, ibg::IBG, ::Nothing, ::ZeroField, c) = zero(ibg) - -# Fallback for `ZeroField` tracers and velocities -@inline _advective_tracer_flux_x(i, j, k, ibg::IBG, scheme, ::ZeroField, ::ZeroField) = zero(ibg) -@inline _advective_tracer_flux_y(i, j, k, ibg::IBG, scheme, ::ZeroField, ::ZeroField) = zero(ibg) -@inline _advective_tracer_flux_z(i, j, k, ibg::IBG, scheme, ::ZeroField, ::ZeroField) = zero(ibg) - -# Fallback for `ZeroField` tracers -@inline _advective_tracer_flux_x(i, j, k, ibg::IBG, scheme, U, ::ZeroField) = zero(ibg) -@inline _advective_tracer_flux_y(i, j, k, ibg::IBG, scheme, V, ::ZeroField) = zero(ibg) -@inline _advective_tracer_flux_z(i, j, k, ibg::IBG, scheme, W, ::ZeroField) = zero(ibg) - -# Fallback for `ZeroField` velocities -@inline _advective_tracer_flux_x(i, j, k, ibg::IBG, scheme, ::ZeroField, c) = zero(ibg) -@inline _advective_tracer_flux_y(i, j, k, ibg::IBG, scheme, ::ZeroField, c) = zero(ibg) -@inline _advective_tracer_flux_z(i, j, k, ibg::IBG, scheme, ::ZeroField, c) = zero(ibg) +# Disambiguation for momentum fluxes in x.... +@inline _advective_momentum_flux_Uu(i, j, k, ibg::IBG, advection::FluxFormAdvection, args...) = + _advective_momentum_flux_Uu(i, j, k, ibg, advection.x, args...) + +@inline _advective_momentum_flux_Uv(i, j, k, ibg::IBG, advection::FluxFormAdvection, args...) = + _advective_momentum_flux_Uv(i, j, k, ibg, advection.x, args...) + +@inline _advective_momentum_flux_Uw(i, j, k, ibg::IBG, advection::FluxFormAdvection, args...) = + _advective_momentum_flux_Uw(i, j, k, ibg, advection.x, args...) + +# Disambiguation for momentum fluxes in y.... +@inline _advective_momentum_flux_Vu(i, j, k, ibg::IBG, advection::FluxFormAdvection, args...) = + _advective_momentum_flux_Vu(i, j, k, ibg, advection.y, args...) + +@inline _advective_momentum_flux_Vv(i, j, k, ibg::IBG, advection::FluxFormAdvection, args...) = + _advective_momentum_flux_Vv(i, j, k, ibg, advection.y, args...) + +@inline _advective_momentum_flux_Vw(i, j, k, ibg::IBG, advection::FluxFormAdvection, args...) = + _advective_momentum_flux_Vw(i, j, k, ibg, advection.y, args...) + +# Disambiguation for momentum fluxes in z.... +@inline _advective_momentum_flux_Wu(i, j, k, ibg::IBG, advection::FluxFormAdvection, args...) = + _advective_momentum_flux_Wu(i, j, k, ibg, advection.z, args...) + +@inline _advective_momentum_flux_Wv(i, j, k, ibg::IBG, advection::FluxFormAdvection, args...) = + _advective_momentum_flux_Wv(i, j, k, ibg, advection.z, args...) + +@inline _advective_momentum_flux_Ww(i, j, k, ibg::IBG, advection::FluxFormAdvection, args...) = + _advective_momentum_flux_Ww(i, j, k, ibg, advection.z, args...) + +# Fallback for `nothing` Advection + +for flux_dir in (:x, :y, :z) + advective_tracer_flux = Symbol(:_advective_tracer_flux_, flux_dir) + + @eval begin + @inline $advective_tracer_flux(i, j, k, ibg::IBG, ::Nothing, args...) = zero(ibg) + @inline $advective_tracer_flux(i, j, k, ibg::IBG, ::Nothing, ::ZeroField, ::ZeroField) = zero(ibg) + @inline $advective_tracer_flux(i, j, k, ibg::IBG, ::Nothing, U, ::ZeroField) = zero(ibg) + @inline $advective_tracer_flux(i, j, k, ibg::IBG, scheme, ::ZeroField, ::ZeroField) = zero(ibg) + @inline $advective_tracer_flux(i, j, k, ibg::IBG, scheme, U, ::ZeroField) = zero(ibg) + @inline $advective_tracer_flux(i, j, k, ibg::IBG, scheme, ::ZeroField, c) = zero(ibg) + end +end + +for flux_dir in (:Uu, :Vu, :Wu, :Uv, :Vv, :Wv, :Uw, :Vw, :Ww) + advective_momentum_flux = Symbol(:_advective_momentum_flux_, flux_dir) + + @eval begin + @inline $advective_momentum_flux(i, j, k, ibg::IBG, ::Nothing, args...) = zero(ibg) + @inline $advective_momentum_flux(i, j, k, ibg::IBG, ::Nothing, ::ZeroField, ::ZeroField) = zero(ibg) + @inline $advective_momentum_flux(i, j, k, ibg::IBG, ::Nothing, U, ::ZeroField) = zero(ibg) + @inline $advective_momentum_flux(i, j, k, ibg::IBG, scheme, ::ZeroField, ::ZeroField) = zero(ibg) + @inline $advective_momentum_flux(i, j, k, ibg::IBG, scheme, U, ::ZeroField) = zero(ibg) + @inline $advective_momentum_flux(i, j, k, ibg::IBG, scheme, ::ZeroField, u) = zero(ibg) + end +end ##### ##### "Boundary-aware" reconstruct From 7843e0d7a70f74b9612e47e3ab63024cf53072ce Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Mon, 26 Aug 2024 07:51:29 -0400 Subject: [PATCH 04/68] remove repeat code --- src/Advection/tracer_advection_operators.jl | 13 ------------- 1 file changed, 13 deletions(-) diff --git a/src/Advection/tracer_advection_operators.jl b/src/Advection/tracer_advection_operators.jl index 44719bfbc9..239628b481 100644 --- a/src/Advection/tracer_advection_operators.jl +++ b/src/Advection/tracer_advection_operators.jl @@ -20,19 +20,6 @@ for flux_dir in (:x, :y, :z) end end -for flux_dir in (:Uu, :Vu, :Wu, :Uv, :Vv, :Wv, :Uw, :Vw, :Ww) - advective_momentum_flux = Symbol(:_advective_momentum_flux_, flux_dir) - - @eval begin - @inline $advective_momentum_flux(i, j, k, grid, ::Nothing, args...) = zero(grid) - @inline $advective_momentum_flux(i, j, k, grid, ::Nothing, ::ZeroField, ::ZeroField) = zero(grid) - @inline $advective_momentum_flux(i, j, k, grid, ::Nothing, U, ::ZeroField) = zero(grid) - @inline $advective_momentum_flux(i, j, k, grid, scheme, ::ZeroField, ::ZeroField) = zero(grid) - @inline $advective_momentum_flux(i, j, k, grid, scheme, U, ::ZeroField) = zero(grid) - @inline $advective_momentum_flux(i, j, k, grid, scheme, ::ZeroField, u) = zero(grid) - end -end - ##### ##### Tracer advection operator ##### From abecb2819aa6e2804ac5c0892bab2954015975f9 Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Mon, 26 Aug 2024 08:03:48 -0400 Subject: [PATCH 05/68] precopmiles --- src/Advection/topologically_conditional_interpolation.jl | 4 +++- .../compute_nonhydrostatic_boundary_tendencies.jl | 2 +- 2 files changed, 4 insertions(+), 2 deletions(-) diff --git a/src/Advection/topologically_conditional_interpolation.jl b/src/Advection/topologically_conditional_interpolation.jl index 37f98a76c7..2eb526eef1 100644 --- a/src/Advection/topologically_conditional_interpolation.jl +++ b/src/Advection/topologically_conditional_interpolation.jl @@ -26,9 +26,11 @@ 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) -for direction in (:x, :y, :z) +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 diff --git a/src/Models/NonhydrostaticModels/compute_nonhydrostatic_boundary_tendencies.jl b/src/Models/NonhydrostaticModels/compute_nonhydrostatic_boundary_tendencies.jl index 9fb4e44c00..f6a5cbf1a9 100644 --- a/src/Models/NonhydrostaticModels/compute_nonhydrostatic_boundary_tendencies.jl +++ b/src/Models/NonhydrostaticModels/compute_nonhydrostatic_boundary_tendencies.jl @@ -1,6 +1,6 @@ import Oceananigans.Models: compute_boundary_tendencies! -using Oceananigans.TurbulenceClosures: equired_halo_size_x, required_halo_size_y +using Oceananigans.TurbulenceClosures: required_halo_size_x, required_halo_size_y using Oceananigans.Grids: XFlatGrid, YFlatGrid # TODO: the code in this file is difficult to understand. From 9d87f01cb9b51ced800ffbbe47340560b66843fc Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Mon, 26 Aug 2024 08:23:51 -0400 Subject: [PATCH 06/68] keep these changes for later --- src/TurbulenceClosures/TurbulenceClosures.jl | 2 +- .../scalar_diffusivity.jl | 12 ------------ 2 files changed, 1 insertion(+), 13 deletions(-) diff --git a/src/TurbulenceClosures/TurbulenceClosures.jl b/src/TurbulenceClosures/TurbulenceClosures.jl index 701a7b1927..a99d3481f7 100644 --- a/src/TurbulenceClosures/TurbulenceClosures.jl +++ b/src/TurbulenceClosures/TurbulenceClosures.jl @@ -50,7 +50,7 @@ using Oceananigans.Utils using Oceananigans.Architectures: AbstractArchitecture, device using Oceananigans.Fields: FunctionField -import Oceananigans.Advection: required_halo_size_x, required_halo_size_y, required_halo_size_z +import Oceananigans.Grids: required_halo_size_x, required_halo_size_y, required_halo_size_z import Oceananigans.Architectures: on_architecture const VerticallyBoundedGrid{FT} = AbstractGrid{FT, <:Any, <:Any, <:Bounded} diff --git a/src/TurbulenceClosures/turbulence_closure_implementations/scalar_diffusivity.jl b/src/TurbulenceClosures/turbulence_closure_implementations/scalar_diffusivity.jl index 00473f28fa..c5f4a71de2 100644 --- a/src/TurbulenceClosures/turbulence_closure_implementations/scalar_diffusivity.jl +++ b/src/TurbulenceClosures/turbulence_closure_implementations/scalar_diffusivity.jl @@ -173,18 +173,6 @@ Shorthand for a `ScalarDiffusivity` with `HorizontalDivergenceFormulation()`. Se HorizontalScalarDiffusivity(FT::DataType; kwargs...) = ScalarDiffusivity(ExplicitTimeDiscretization(), HorizontalFormulation(), FT; kwargs...) HorizontalDivergenceScalarDiffusivity(FT::DataType; kwargs...) = ScalarDiffusivity(ExplicitTimeDiscretization(), HorizontalDivergenceFormulation(), FT; kwargs...) -required_halo_size_x(::ScalarDiffusivity{<:Any, ThreeDimensionalFormulation, N}) where N = N -required_halo_size_y(::ScalarDiffusivity{<:Any, ThreeDimensionalFormulation, N}) where N = N -required_halo_size_z(::ScalarDiffusivity{<:Any, ThreeDimensionalFormulation, N}) where N = N - -required_halo_size_x(::ScalarDiffusivity{<:Any, VerticalFormulation, N}) where N = 0 -required_halo_size_y(::ScalarDiffusivity{<:Any, VerticalFormulation, N}) where N = 0 -required_halo_size_z(::ScalarDiffusivity{<:Any, VerticalFormulation, N}) where N = N - -required_halo_size_x(::ScalarDiffusivity{<:Any, HorizontalFormulation, N}) where N = N -required_halo_size_y(::ScalarDiffusivity{<:Any, HorizontalFormulation, N}) where N = N -required_halo_size_z(::ScalarDiffusivity{<:Any, HorizontalFormulation, N}) where N = 0 - @inline function with_tracers(tracers, closure::ScalarDiffusivity{TD, F, N}) where {TD, F, N} κ = tracer_diffusivities(tracers, closure.κ) return ScalarDiffusivity{TD, F, N}(closure.ν, κ) From 03b0d2825137b1ea18e31e6414cf22b58451a704 Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Mon, 26 Aug 2024 09:14:56 -0400 Subject: [PATCH 07/68] correct the outside buffer --- src/Advection/topologically_conditional_interpolation.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/Advection/topologically_conditional_interpolation.jl b/src/Advection/topologically_conditional_interpolation.jl index 2eb526eef1..9c623c6b57 100644 --- a/src/Advection/topologically_conditional_interpolation.jl +++ b/src/Advection/topologically_conditional_interpolation.jl @@ -70,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 From 70d647514852404da657ec2340ac1e232e6c3fdc Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Mon, 26 Aug 2024 09:56:36 -0400 Subject: [PATCH 08/68] fix tests --- src/Advection/weno_reconstruction.jl | 5 +++-- test/test_immersed_advection.jl | 4 ++-- 2 files changed, 5 insertions(+), 4 deletions(-) diff --git a/src/Advection/weno_reconstruction.jl b/src/Advection/weno_reconstruction.jl index 5ed5a16adc..bc0720e1e3 100644 --- a/src/Advection/weno_reconstruction.jl +++ b/src/Advection/weno_reconstruction.jl @@ -96,7 +96,8 @@ WENO reconstruction order 7 function WENO(FT::DataType=Float64; order = 5, grid = nothing, - bounds = nothing) + bounds = nothing, + buffer_scheme = nothing) if !(grid isa Nothing) FT = eltype(grid) @@ -111,8 +112,8 @@ function WENO(FT::DataType=Float64; N = Int((order + 1) ÷ 2) weno_coefficients = compute_reconstruction_coefficients(grid, FT, :WENO; order = N) - buffer_scheme = WENO(FT; grid, order = order - 2, bounds) advecting_velocity_scheme = Centered(FT; grid, order = order - 1) + buffer_scheme = isnothing(buffer_scheme) ? WENO(FT; grid, order = order - 2, bounds) : buffer_scheme end return WENO{N, FT}(weno_coefficients..., bounds, buffer_scheme, advecting_velocity_scheme) diff --git a/test/test_immersed_advection.jl b/test/test_immersed_advection.jl index aaae16a899..dc36776b6a 100644 --- a/test/test_immersed_advection.jl +++ b/test/test_immersed_advection.jl @@ -10,7 +10,7 @@ using Oceananigans.Advection: _biased_interpolate_xᶠᵃᵃ, _biased_interpolate_yᵃᶜᵃ, _biased_interpolate_yᵃᶠᵃ, - TracerAdvection + FluxFormAdvection advection_schemes = [Centered, UpwindBiased, WENO] @@ -128,7 +128,7 @@ for arch in archs for adv in advection_schemes, buffer in [1, 2, 3, 4, 5] directional_scheme = adv(order = advective_order(buffer, adv)) - scheme = TracerAdvection(directional_scheme, directional_scheme, directional_scheme) + scheme = FluxFormAdvection(directional_scheme, directional_scheme, directional_scheme) for g in [grid, ibg] @info " Testing immersed tracer conservation [$(typeof(arch)), $(summary(scheme)), $(typeof(g).name.wrapper)]" run_tracer_conservation_test(g, scheme) From 42b1f3a8e50f59c33230c7fcc730970046edac70 Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Mon, 26 Aug 2024 12:31:51 -0400 Subject: [PATCH 09/68] change the show method --- src/Advection/flux_form_advection.jl | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/src/Advection/flux_form_advection.jl b/src/Advection/flux_form_advection.jl index 58369fc19f..aa1931a53d 100644 --- a/src/Advection/flux_form_advection.jl +++ b/src/Advection/flux_form_advection.jl @@ -26,6 +26,12 @@ function FluxFormAdvection(x_advection, y_advection, z_advection) 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) From 6181e01f19ef66a4d7360d3c9262d7c6f82683a3 Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Mon, 26 Aug 2024 12:55:45 -0400 Subject: [PATCH 10/68] back to previous buffer scheme --- src/Advection/weno_reconstruction.jl | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/src/Advection/weno_reconstruction.jl b/src/Advection/weno_reconstruction.jl index bc0720e1e3..f91c633557 100644 --- a/src/Advection/weno_reconstruction.jl +++ b/src/Advection/weno_reconstruction.jl @@ -96,8 +96,7 @@ WENO reconstruction order 7 function WENO(FT::DataType=Float64; order = 5, grid = nothing, - bounds = nothing, - buffer_scheme = nothing) + bounds = nothing) if !(grid isa Nothing) FT = eltype(grid) @@ -113,7 +112,7 @@ function WENO(FT::DataType=Float64; weno_coefficients = compute_reconstruction_coefficients(grid, FT, :WENO; order = N) advecting_velocity_scheme = Centered(FT; grid, order = order - 1) - buffer_scheme = isnothing(buffer_scheme) ? WENO(FT; grid, order = order - 2, bounds) : buffer_scheme + buffer_scheme = WENO(FT; grid, order = order - 2, bounds) end return WENO{N, FT}(weno_coefficients..., bounds, buffer_scheme, advecting_velocity_scheme) From 77873b3aeb86d8b3bdedf1d77c3191d2a130be68 Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Mon, 26 Aug 2024 13:37:33 -0400 Subject: [PATCH 11/68] allow for changes in the advection --- src/Advection/Advection.jl | 1 + src/Advection/adapt_advection_order.jl | 120 ++++++++++++++++++ .../hydrostatic_free_surface_model.jl | 8 +- .../nonhydrostatic_model.jl | 7 +- 4 files changed, 134 insertions(+), 2 deletions(-) create mode 100644 src/Advection/adapt_advection_order.jl diff --git a/src/Advection/Advection.jl b/src/Advection/Advection.jl index 3fed3242de..27b5ea6e27 100644 --- a/src/Advection/Advection.jl +++ b/src/Advection/Advection.jl @@ -83,5 +83,6 @@ 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 diff --git a/src/Advection/adapt_advection_order.jl b/src/Advection/adapt_advection_order.jl new file mode 100644 index 0000000000..ebecd3035a --- /dev/null +++ b/src/Advection/adapt_advection_order.jl @@ -0,0 +1,120 @@ + + +function adapt_advection_order_x(advection::Centered{H}, grid::AbstractGrid, N::Int) where H + if N == 1 + return nothing + elseif N >= H + return advection + else + return Centered(; order = N) + end +end + +function adapt_advection_order_x(advection::UpwindBiased{H}, grid::AbstractGrid, N::Int) where H + if N == 1 + return nothing + elseif N >= H + return advection + else + return UpwindBiased(; order = N) + end +end + +adapt_advection_order_y(advection, grid::AbstractGrid, N::Int) = adapt_advection_order_x(advection, grid, N) +adapt_advection_order_z(advection, grid::AbstractGrid, N::Int) = adapt_advection_order_x(advection, grid, N) + +""" + new_weno_scheme(grid, order, bounds, T) + +Constructs a new WENO scheme based on the given parameters. `T` is the type of the weno coefficients. +A _non-stretched_ WENO scheme has `T` equal to `Nothing`. 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(grid, order, bounds, T) = ifelse(T == Nothing, WENO(; order, bounds), WENO(grid; order, bounds)) + +function adapt_advection_order_x(advection::WENO{H, FT, XT, YT, ZT}, grid::AbstractGrid) where {H, FT, XT, YT, ZT} + + if N == 1 + return nothing + elseif N >= H + return advection + else + return new_weno_scheme(grid, N, advection.bounds, XT) + end +end + +function adapt_advection_order_y(advection::WENO{H, FT, XT, YT, ZT}, grid::AbstractGrid, N::Int) where {H, FT, XT, YT, ZT} + + if N > H + return advection + else + return new_weno_scheme(grid, N, advection.bounds, YT) + end +end + +function adapt_advection_order_z(advection::WENO{H, FT, XT, YT, ZT}, grid::AbstractGrid, N::Int) where {H, FT, XT, YT, ZT} + + if N == 1 + return nothing + elseif N >= H + return advection + else + return new_weno_scheme(grid, N, advection.bounds, XT) + end +end + +""" + 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 `nothing`. +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. + +The adapted advection scheme with adjusted advection order returned by this function is a `FluxFormAdvection`. +""" +function adapt_advection_order(advection, grid::AbstractGrid) + advection_x = adapt_advection_order_x(advection, grid, grid.Nx) + advection_y = adapt_advection_order_y(advection, grid, grid.Ny) + advection_z = adapt_advection_order_z(advection, grid, grid.Nz) + + # Check that we indeed changed the advection operator + changed_x = advection_x != advection + changed_y = advection_y != advection + changed_z = advection_z != advection + + new_advection = FluxFormAdvection(advection_x, advection_y, advection_z) + changed_advection = any((changed_x, changed_y, changed_z)) + + if changed_advection + @info "User-defined advection scheme $(advection) reduced to $(new_advection) to comply with grid-size limitations." + end + + return ifelse(changed_advection, new_advection, advection) +end + +function adapt_advection_order(advection::FluxFormAdvection, grid::AbstractGrid) + advection_x = adapt_advection_order_x(advection.x, grid, grid.Nx) + advection_y = adapt_advection_order_y(advection.y, grid, grid.Ny) + advection_z = adapt_advection_order_z(advection.z, grid, grid.Nz) + + # Check that we indeed changed the advection operator + changed_x = advection_x != advection.x + changed_y = advection_y != advection.y + changed_z = advection_z != advection.z + + new_advection = FluxFormAdvection(advection_x, advection_y, advection_z) + changed_advection = any((changed_x, changed_y, changed_z)) + + if changed_advection + @info "User-defined advection scheme $(advection) reduced to $(new_advection) to comply with grid-size limitations." + end + + return ifelse(changed_advection, new_advection, advection) +end + +# For the moment, we do not adapt the advection order for the VectorInvariant advection scheme +adapt_advection_order(advection::VectorInvariant, grid::AbstractGrid) = advection \ No newline at end of file diff --git a/src/Models/HydrostaticFreeSurfaceModels/hydrostatic_free_surface_model.jl b/src/Models/HydrostaticFreeSurfaceModels/hydrostatic_free_surface_model.jl index 588235882e..3c2b501041 100644 --- a/src/Models/HydrostaticFreeSurfaceModels/hydrostatic_free_surface_model.jl +++ b/src/Models/HydrostaticFreeSurfaceModels/hydrostatic_free_surface_model.jl @@ -3,7 +3,7 @@ using OrderedCollections: OrderedDict using Oceananigans.DistributedComputations using Oceananigans.Architectures: AbstractArchitecture -using Oceananigans.Advection: AbstractAdvectionScheme, CenteredSecondOrder, VectorInvariant +using Oceananigans.Advection: AbstractAdvectionScheme, CenteredSecondOrder, VectorInvariant, adapt_advection_order using Oceananigans.BuoyancyModels: validate_buoyancy, regularize_buoyancy, SeawaterBuoyancy, g_Earth using Oceananigans.BoundaryConditions: regularize_field_boundary_conditions using Oceananigans.Biogeochemistry: validate_biogeochemistry, AbstractBiogeochemistry, biogeochemical_auxiliary_fields @@ -130,6 +130,12 @@ function HydrostaticFreeSurfaceModel(; grid, tracers = tupleit(tracers) # supports tracers=:c keyword argument (for example) + # Adjust advection scheme to be valid on a particular grid size. i.e. if the grid size + # is smaller than the advection order, reduce the order of the advection in that particular + # direction + momentum_advection = adapt_advection_order(momentum_advection, grid) + tracer_advection = adapt_advection_order(tracer_advection, grid) + tracers, auxiliary_fields = validate_biogeochemistry(tracers, merge(auxiliary_fields, biogeochemical_auxiliary_fields(biogeochemistry)), biogeochemistry, grid, clock) validate_buoyancy(buoyancy, tracernames(tracers)) buoyancy = regularize_buoyancy(buoyancy) diff --git a/src/Models/NonhydrostaticModels/nonhydrostatic_model.jl b/src/Models/NonhydrostaticModels/nonhydrostatic_model.jl index 688f834f3a..03ff08a517 100644 --- a/src/Models/NonhydrostaticModels/nonhydrostatic_model.jl +++ b/src/Models/NonhydrostaticModels/nonhydrostatic_model.jl @@ -3,7 +3,7 @@ using OrderedCollections: OrderedDict using Oceananigans.Architectures: AbstractArchitecture using Oceananigans.DistributedComputations: Distributed -using Oceananigans.Advection: CenteredSecondOrder +using Oceananigans.Advection: CenteredSecondOrder, adapt_advection_order using Oceananigans.BuoyancyModels: validate_buoyancy, regularize_buoyancy, SeawaterBuoyancy using Oceananigans.Biogeochemistry: validate_biogeochemistry, AbstractBiogeochemistry, biogeochemical_auxiliary_fields using Oceananigans.BoundaryConditions: regularize_field_boundary_conditions @@ -175,6 +175,11 @@ function NonhydrostaticModel(; grid, validate_buoyancy(buoyancy, tracernames(tracers)) buoyancy = regularize_buoyancy(buoyancy) + # Adjust advection scheme to be valid on a particular grid size. i.e. if the grid size + # is smaller than the advection order, reduce the order of the advection in that particular + # direction + advection = adapt_advection_order(advection, grid) + # Adjust halos when the advection scheme or turbulence closure requires it. # Note that halos are isotropic by default; however we respect user-input here # by adjusting each (x, y, z) halo individually. From 9c7ff379784325a7542575f4cc4ff9326d123515 Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Mon, 26 Aug 2024 13:39:46 -0400 Subject: [PATCH 12/68] switch position of functions --- src/Advection/adapt_advection_order.jl | 114 +++++++++++++------------ 1 file changed, 58 insertions(+), 56 deletions(-) diff --git a/src/Advection/adapt_advection_order.jl b/src/Advection/adapt_advection_order.jl index ebecd3035a..d9fbce1b02 100644 --- a/src/Advection/adapt_advection_order.jl +++ b/src/Advection/adapt_advection_order.jl @@ -1,4 +1,62 @@ +""" + 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 `nothing`. +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. + +The adapted advection scheme with adjusted advection order returned by this function is a `FluxFormAdvection`. +""" +function adapt_advection_order(advection, grid::AbstractGrid) + advection_x = adapt_advection_order_x(advection, grid, grid.Nx) + advection_y = adapt_advection_order_y(advection, grid, grid.Ny) + advection_z = adapt_advection_order_z(advection, grid, grid.Nz) + + # Check that we indeed changed the advection operator + changed_x = advection_x != advection + changed_y = advection_y != advection + changed_z = advection_z != advection + + new_advection = FluxFormAdvection(advection_x, advection_y, advection_z) + changed_advection = any((changed_x, changed_y, changed_z)) + if changed_advection + @info "User-defined advection scheme $(advection) reduced to $(new_advection) to comply with grid-size limitations." + end + + return ifelse(changed_advection, new_advection, advection) +end + +function adapt_advection_order(advection::FluxFormAdvection, grid::AbstractGrid) + advection_x = adapt_advection_order_x(advection.x, grid, grid.Nx) + advection_y = adapt_advection_order_y(advection.y, grid, grid.Ny) + advection_z = adapt_advection_order_z(advection.z, grid, grid.Nz) + + # Check that we indeed changed the advection operator + changed_x = advection_x != advection.x + changed_y = advection_y != advection.y + changed_z = advection_z != advection.z + + new_advection = FluxFormAdvection(advection_x, advection_y, advection_z) + changed_advection = any((changed_x, changed_y, changed_z)) + + if changed_advection + @info "User-defined advection scheme $(advection) reduced to $(new_advection) to comply with grid-size limitations." + end + + return ifelse(changed_advection, new_advection, advection) +end + +# For the moment, we do not adapt the advection order for the VectorInvariant advection scheme +adapt_advection_order(advection::VectorInvariant, grid::AbstractGrid) = advection + +##### +##### Directional adapt advection order +##### function adapt_advection_order_x(advection::Centered{H}, grid::AbstractGrid, N::Int) where H if N == 1 @@ -62,59 +120,3 @@ function adapt_advection_order_z(advection::WENO{H, FT, XT, YT, ZT}, grid::Abstr return new_weno_scheme(grid, N, advection.bounds, XT) end end - -""" - 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 `nothing`. -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. - -The adapted advection scheme with adjusted advection order returned by this function is a `FluxFormAdvection`. -""" -function adapt_advection_order(advection, grid::AbstractGrid) - advection_x = adapt_advection_order_x(advection, grid, grid.Nx) - advection_y = adapt_advection_order_y(advection, grid, grid.Ny) - advection_z = adapt_advection_order_z(advection, grid, grid.Nz) - - # Check that we indeed changed the advection operator - changed_x = advection_x != advection - changed_y = advection_y != advection - changed_z = advection_z != advection - - new_advection = FluxFormAdvection(advection_x, advection_y, advection_z) - changed_advection = any((changed_x, changed_y, changed_z)) - - if changed_advection - @info "User-defined advection scheme $(advection) reduced to $(new_advection) to comply with grid-size limitations." - end - - return ifelse(changed_advection, new_advection, advection) -end - -function adapt_advection_order(advection::FluxFormAdvection, grid::AbstractGrid) - advection_x = adapt_advection_order_x(advection.x, grid, grid.Nx) - advection_y = adapt_advection_order_y(advection.y, grid, grid.Ny) - advection_z = adapt_advection_order_z(advection.z, grid, grid.Nz) - - # Check that we indeed changed the advection operator - changed_x = advection_x != advection.x - changed_y = advection_y != advection.y - changed_z = advection_z != advection.z - - new_advection = FluxFormAdvection(advection_x, advection_y, advection_z) - changed_advection = any((changed_x, changed_y, changed_z)) - - if changed_advection - @info "User-defined advection scheme $(advection) reduced to $(new_advection) to comply with grid-size limitations." - end - - return ifelse(changed_advection, new_advection, advection) -end - -# For the moment, we do not adapt the advection order for the VectorInvariant advection scheme -adapt_advection_order(advection::VectorInvariant, grid::AbstractGrid) = advection \ No newline at end of file From 9ca656eee302c1b1317ee5e519a35677f08dbf90 Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Mon, 26 Aug 2024 13:42:11 -0400 Subject: [PATCH 13/68] better comment --- src/Advection/adapt_advection_order.jl | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/Advection/adapt_advection_order.jl b/src/Advection/adapt_advection_order.jl index d9fbce1b02..d4e57b79bb 100644 --- a/src/Advection/adapt_advection_order.jl +++ b/src/Advection/adapt_advection_order.jl @@ -9,7 +9,8 @@ A high order advection sheme is reduced to a lower order advection scheme if the - `advection`: The original advection scheme. - `grid::AbstractGrid`: The grid on which the advection scheme is applied. -The adapted advection scheme with adjusted advection order returned by this function is a `FluxFormAdvection`. +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 = adapt_advection_order_x(advection, grid, grid.Nx) From 4aabf62c7abf1e93d700c8672c31941a12d8aca0 Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Mon, 26 Aug 2024 13:42:36 -0400 Subject: [PATCH 14/68] better info --- src/Advection/adapt_advection_order.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/Advection/adapt_advection_order.jl b/src/Advection/adapt_advection_order.jl index d4e57b79bb..23c2a2f910 100644 --- a/src/Advection/adapt_advection_order.jl +++ b/src/Advection/adapt_advection_order.jl @@ -26,7 +26,7 @@ function adapt_advection_order(advection, grid::AbstractGrid) changed_advection = any((changed_x, changed_y, changed_z)) if changed_advection - @info "User-defined advection scheme $(advection) reduced to $(new_advection) to comply with grid-size limitations." + @info "The user-defined advection scheme $(advection) has been reduced to $(new_advection) to comply with grid-size limitations." end return ifelse(changed_advection, new_advection, advection) From 39a4ef86e6a609a6ddcfef0e06201b830b484f54 Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Mon, 26 Aug 2024 13:44:47 -0400 Subject: [PATCH 15/68] bugfix --- src/Advection/adapt_advection_order.jl | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/src/Advection/adapt_advection_order.jl b/src/Advection/adapt_advection_order.jl index 23c2a2f910..249577c911 100644 --- a/src/Advection/adapt_advection_order.jl +++ b/src/Advection/adapt_advection_order.jl @@ -65,7 +65,7 @@ function adapt_advection_order_x(advection::Centered{H}, grid::AbstractGrid, N:: elseif N >= H return advection else - return Centered(; order = N) + return Centered(; order = N * 2) end end @@ -75,7 +75,7 @@ function adapt_advection_order_x(advection::UpwindBiased{H}, grid::AbstractGrid, elseif N >= H return advection else - return UpwindBiased(; order = N) + return UpwindBiased(; order = N * 2 - 1) end end @@ -91,14 +91,14 @@ we rebuild the advection without passing the grid information, otherwise we use """ new_weno_scheme(grid, order, bounds, T) = ifelse(T == Nothing, WENO(; order, bounds), WENO(grid; order, bounds)) -function adapt_advection_order_x(advection::WENO{H, FT, XT, YT, ZT}, grid::AbstractGrid) where {H, FT, XT, YT, ZT} +function adapt_advection_order_x(advection::WENO{H, FT, XT, YT, ZT}, grid::AbstractGrid, N::Int) where {H, FT, XT, YT, ZT} if N == 1 return nothing elseif N >= H return advection else - return new_weno_scheme(grid, N, advection.bounds, XT) + return new_weno_scheme(grid, N * 2 - 1, advection.bounds, XT) end end @@ -107,7 +107,7 @@ function adapt_advection_order_y(advection::WENO{H, FT, XT, YT, ZT}, grid::Abstr if N > H return advection else - return new_weno_scheme(grid, N, advection.bounds, YT) + return new_weno_scheme(grid, N * 2 - 1, advection.bounds, YT) end end @@ -118,6 +118,6 @@ function adapt_advection_order_z(advection::WENO{H, FT, XT, YT, ZT}, grid::Abstr elseif N >= H return advection else - return new_weno_scheme(grid, N, advection.bounds, XT) + return new_weno_scheme(grid, N * 2 - 1, advection.bounds, XT) end end From 0ba24d6172c4cfa39fc0a2f4f513f01fa736ee05 Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Mon, 26 Aug 2024 13:46:50 -0400 Subject: [PATCH 16/68] better info --- src/Advection/adapt_advection_order.jl | 20 ++++++++++++++++---- 1 file changed, 16 insertions(+), 4 deletions(-) diff --git a/src/Advection/adapt_advection_order.jl b/src/Advection/adapt_advection_order.jl index 249577c911..a273c27f96 100644 --- a/src/Advection/adapt_advection_order.jl +++ b/src/Advection/adapt_advection_order.jl @@ -25,8 +25,14 @@ function adapt_advection_order(advection, grid::AbstractGrid) new_advection = FluxFormAdvection(advection_x, advection_y, advection_z) changed_advection = any((changed_x, changed_y, changed_z)) - if changed_advection - @info "The user-defined advection scheme $(advection) has been reduced to $(new_advection) to comply with grid-size limitations." + if changed_x + @info "User-defined advection scheme $(summary(advection)) reduced to $(summary(new_advection.x)) in the x-direction to comply with grid-size limitations." + end + if changed_y + @info "User-defined advection scheme $(summary(advection)) reduced to $(summary(new_advection.y)) in the y-direction to comply with grid-size limitations." + end + if changed_z + @info "User-defined advection scheme $(summary(advection)) reduced to $(summary(new_advection.z)) in the z-direction to comply with grid-size limitations." end return ifelse(changed_advection, new_advection, advection) @@ -45,8 +51,14 @@ function adapt_advection_order(advection::FluxFormAdvection, grid::AbstractGrid) new_advection = FluxFormAdvection(advection_x, advection_y, advection_z) changed_advection = any((changed_x, changed_y, changed_z)) - if changed_advection - @info "User-defined advection scheme $(advection) reduced to $(new_advection) to comply with grid-size limitations." + if changed_x + @info "User-defined advection scheme $(summary(advection)) reduced to $(summary(new_advection.x)) in the x-direction to comply with grid-size limitations." + end + if changed_y + @info "User-defined advection scheme $(summary(advection)) reduced to $(summary(new_advection.y)) in the y-direction to comply with grid-size limitations." + end + if changed_z + @info "User-defined advection scheme $(summary(advection)) reduced to $(summary(new_advection.z)) in the z-direction to comply with grid-size limitations." end return ifelse(changed_advection, new_advection, advection) From 71dbde1ab526a83c845a38811344618b33150924 Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Mon, 26 Aug 2024 13:48:47 -0400 Subject: [PATCH 17/68] make sure N==1 -> nothing --- src/Advection/adapt_advection_order.jl | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/src/Advection/adapt_advection_order.jl b/src/Advection/adapt_advection_order.jl index a273c27f96..6767463a3a 100644 --- a/src/Advection/adapt_advection_order.jl +++ b/src/Advection/adapt_advection_order.jl @@ -115,8 +115,10 @@ function adapt_advection_order_x(advection::WENO{H, FT, XT, YT, ZT}, grid::Abstr end function adapt_advection_order_y(advection::WENO{H, FT, XT, YT, ZT}, grid::AbstractGrid, N::Int) where {H, FT, XT, YT, ZT} - - if N > H + + if N == 1 + return nothing + elseif N > H return advection else return new_weno_scheme(grid, N * 2 - 1, advection.bounds, YT) From 817d058699acba2316af14c97c421acc8d507e8d Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Mon, 26 Aug 2024 13:51:25 -0400 Subject: [PATCH 18/68] Update src/Advection/momentum_advection_operators.jl Co-authored-by: Gregory L. Wagner --- src/Advection/momentum_advection_operators.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/Advection/momentum_advection_operators.jl b/src/Advection/momentum_advection_operators.jl index 78b24dea8c..b475796146 100644 --- a/src/Advection/momentum_advection_operators.jl +++ b/src/Advection/momentum_advection_operators.jl @@ -94,8 +94,8 @@ end ##### Fallback advection fluxes! ##### -for flux_dir in (:Uu, :Vu, :Wu, :Uv, :Vv, :Wv, :Uw, :Vw, :Ww) - advective_momentum_flux = Symbol(:_advective_momentum_flux_, flux_dir) +for flux_type in (:Uu, :Vu, :Wu, :Uv, :Vv, :Wv, :Uw, :Vw, :Ww) + advective_momentum_flux = Symbol(:_advective_momentum_flux_, flux_type) @eval begin @inline $advective_momentum_flux(i, j, k, grid, ::Nothing, args...) = zero(grid) From 11355d6b37a7eb5fc37795223ebadc55ecec9a28 Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Mon, 26 Aug 2024 14:34:42 -0400 Subject: [PATCH 19/68] add topologies --- src/Advection/adapt_advection_order.jl | 33 ++++++++++++++++---------- 1 file changed, 20 insertions(+), 13 deletions(-) diff --git a/src/Advection/adapt_advection_order.jl b/src/Advection/adapt_advection_order.jl index 6767463a3a..77f735a9ec 100644 --- a/src/Advection/adapt_advection_order.jl +++ b/src/Advection/adapt_advection_order.jl @@ -1,3 +1,5 @@ +using Oceananigans.Grids: topology + """ adapt_advection_order(advection, grid::AbstractGrid) @@ -13,9 +15,9 @@ If the order of advection is changed in at least one direction, the adapted adve by this function is a `FluxFormAdvection`. """ function adapt_advection_order(advection, grid::AbstractGrid) - advection_x = adapt_advection_order_x(advection, grid, grid.Nx) - advection_y = adapt_advection_order_y(advection, grid, grid.Ny) - advection_z = adapt_advection_order_z(advection, grid, grid.Nz) + advection_x = adapt_advection_order_x(advection.x, topology(grid, 1), size(grid, 1), grid) + advection_y = adapt_advection_order_y(advection.y, topology(grid, 2), size(grid, 2), grid) + advection_z = adapt_advection_order_z(advection.z, topology(grid, 3), size(grid, 3), grid) # Check that we indeed changed the advection operator changed_x = advection_x != advection @@ -39,9 +41,9 @@ function adapt_advection_order(advection, grid::AbstractGrid) end function adapt_advection_order(advection::FluxFormAdvection, grid::AbstractGrid) - advection_x = adapt_advection_order_x(advection.x, grid, grid.Nx) - advection_y = adapt_advection_order_y(advection.y, grid, grid.Ny) - advection_z = adapt_advection_order_z(advection.z, grid, grid.Nz) + advection_x = adapt_advection_order_x(advection.x, topology(grid, 1), size(grid, 1), grid) + advection_y = adapt_advection_order_y(advection.y, topology(grid, 2), size(grid, 2), grid) + advection_z = adapt_advection_order_z(advection.z, topology(grid, 3), size(grid, 3), grid) # Check that we indeed changed the advection operator changed_x = advection_x != advection.x @@ -67,11 +69,16 @@ end # For the moment, we do not adapt the advection order for the VectorInvariant advection scheme adapt_advection_order(advection::VectorInvariant, grid::AbstractGrid) = advection +# We only need one halo in bounded directions! +adapt_advection_order_x(advection, topo, N, grid) = advection +adapt_advection_order_y(advection, topo, N, grid) = advection +adapt_advection_order_z(advection, topo, N, grid) = advection + ##### ##### Directional adapt advection order ##### -function adapt_advection_order_x(advection::Centered{H}, grid::AbstractGrid, N::Int) where H +function adapt_advection_order_x(advection::Centered{H}, ::Periodic, N::Int, grid::AbstractGrid) where H if N == 1 return nothing elseif N >= H @@ -81,7 +88,7 @@ function adapt_advection_order_x(advection::Centered{H}, grid::AbstractGrid, N:: end end -function adapt_advection_order_x(advection::UpwindBiased{H}, grid::AbstractGrid, N::Int) where H +function adapt_advection_order_x(advection::UpwindBiased{H}, ::Periodic, grid::AbstractGrid, N::Int) where H if N == 1 return nothing elseif N >= H @@ -91,8 +98,8 @@ function adapt_advection_order_x(advection::UpwindBiased{H}, grid::AbstractGrid, end end -adapt_advection_order_y(advection, grid::AbstractGrid, N::Int) = adapt_advection_order_x(advection, grid, N) -adapt_advection_order_z(advection, grid::AbstractGrid, N::Int) = adapt_advection_order_x(advection, grid, N) +adapt_advection_order_y(advection, topology::Periodic, grid::AbstractGrid, N::Int) = adapt_advection_order_x(advection, topology, grid, N) +adapt_advection_order_z(advection, topology::Periodic, grid::AbstractGrid, N::Int) = adapt_advection_order_x(advection, topology, grid, N) """ new_weno_scheme(grid, order, bounds, T) @@ -103,7 +110,7 @@ we rebuild the advection without passing the grid information, otherwise we use """ new_weno_scheme(grid, order, bounds, T) = ifelse(T == Nothing, WENO(; order, bounds), WENO(grid; order, bounds)) -function adapt_advection_order_x(advection::WENO{H, FT, XT, YT, ZT}, grid::AbstractGrid, N::Int) where {H, FT, XT, YT, ZT} +function adapt_advection_order_x(advection::WENO{H, FT, XT, YT, ZT}, ::Periodic, grid::AbstractGrid, N::Int) where {H, FT, XT, YT, ZT} if N == 1 return nothing @@ -114,7 +121,7 @@ function adapt_advection_order_x(advection::WENO{H, FT, XT, YT, ZT}, grid::Abstr end end -function adapt_advection_order_y(advection::WENO{H, FT, XT, YT, ZT}, grid::AbstractGrid, N::Int) where {H, FT, XT, YT, ZT} +function adapt_advection_order_y(advection::WENO{H, FT, XT, YT, ZT}, ::Periodic, grid::AbstractGrid, N::Int) where {H, FT, XT, YT, ZT} if N == 1 return nothing @@ -125,7 +132,7 @@ function adapt_advection_order_y(advection::WENO{H, FT, XT, YT, ZT}, grid::Abstr end end -function adapt_advection_order_z(advection::WENO{H, FT, XT, YT, ZT}, grid::AbstractGrid, N::Int) where {H, FT, XT, YT, ZT} +function adapt_advection_order_z(advection::WENO{H, FT, XT, YT, ZT}, ::Periodic, grid::AbstractGrid, N::Int) where {H, FT, XT, YT, ZT} if N == 1 return nothing From bc212e3703a59c5ee6a6d69e315b52d9fb595ed4 Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Mon, 26 Aug 2024 14:35:47 -0400 Subject: [PATCH 20/68] add topology argument --- src/Advection/adapt_advection_order.jl | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/src/Advection/adapt_advection_order.jl b/src/Advection/adapt_advection_order.jl index 77f735a9ec..f626e8fdeb 100644 --- a/src/Advection/adapt_advection_order.jl +++ b/src/Advection/adapt_advection_order.jl @@ -78,7 +78,7 @@ adapt_advection_order_z(advection, topo, N, grid) = advection ##### Directional adapt advection order ##### -function adapt_advection_order_x(advection::Centered{H}, ::Periodic, N::Int, grid::AbstractGrid) where H +function adapt_advection_order_x(advection::Centered{H}, topology, N::Int, grid::AbstractGrid) where H if N == 1 return nothing elseif N >= H @@ -88,7 +88,7 @@ function adapt_advection_order_x(advection::Centered{H}, ::Periodic, N::Int, gri end end -function adapt_advection_order_x(advection::UpwindBiased{H}, ::Periodic, grid::AbstractGrid, N::Int) where H +function adapt_advection_order_x(advection::UpwindBiased{H}, topology, grid::AbstractGrid, N::Int) where H if N == 1 return nothing elseif N >= H @@ -98,8 +98,8 @@ function adapt_advection_order_x(advection::UpwindBiased{H}, ::Periodic, grid::A end end -adapt_advection_order_y(advection, topology::Periodic, grid::AbstractGrid, N::Int) = adapt_advection_order_x(advection, topology, grid, N) -adapt_advection_order_z(advection, topology::Periodic, grid::AbstractGrid, N::Int) = adapt_advection_order_x(advection, topology, grid, N) +adapt_advection_order_y(advection, topology, grid::AbstractGrid, N::Int) = adapt_advection_order_x(advection, topology, grid, N) +adapt_advection_order_z(advection, topology, grid::AbstractGrid, N::Int) = adapt_advection_order_x(advection, topology, grid, N) """ new_weno_scheme(grid, order, bounds, T) @@ -110,7 +110,7 @@ we rebuild the advection without passing the grid information, otherwise we use """ new_weno_scheme(grid, order, bounds, T) = ifelse(T == Nothing, WENO(; order, bounds), WENO(grid; order, bounds)) -function adapt_advection_order_x(advection::WENO{H, FT, XT, YT, ZT}, ::Periodic, grid::AbstractGrid, N::Int) where {H, FT, XT, YT, ZT} +function adapt_advection_order_x(advection::WENO{H, FT, XT, YT, ZT}, topology, grid::AbstractGrid, N::Int) where {H, FT, XT, YT, ZT} if N == 1 return nothing @@ -121,7 +121,7 @@ function adapt_advection_order_x(advection::WENO{H, FT, XT, YT, ZT}, ::Periodic, end end -function adapt_advection_order_y(advection::WENO{H, FT, XT, YT, ZT}, ::Periodic, grid::AbstractGrid, N::Int) where {H, FT, XT, YT, ZT} +function adapt_advection_order_y(advection::WENO{H, FT, XT, YT, ZT}, topology, grid::AbstractGrid, N::Int) where {H, FT, XT, YT, ZT} if N == 1 return nothing @@ -132,7 +132,7 @@ function adapt_advection_order_y(advection::WENO{H, FT, XT, YT, ZT}, ::Periodic, end end -function adapt_advection_order_z(advection::WENO{H, FT, XT, YT, ZT}, ::Periodic, grid::AbstractGrid, N::Int) where {H, FT, XT, YT, ZT} +function adapt_advection_order_z(advection::WENO{H, FT, XT, YT, ZT}, topology, grid::AbstractGrid, N::Int) where {H, FT, XT, YT, ZT} if N == 1 return nothing From 261a31abf99403f95eb7c14f4527b4d526ad4e8c Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Mon, 26 Aug 2024 14:46:11 -0400 Subject: [PATCH 21/68] better adaptation --- src/Advection/adapt_advection_order.jl | 61 +++++++------------------- 1 file changed, 16 insertions(+), 45 deletions(-) diff --git a/src/Advection/adapt_advection_order.jl b/src/Advection/adapt_advection_order.jl index f626e8fdeb..f6f8ba6d5b 100644 --- a/src/Advection/adapt_advection_order.jl +++ b/src/Advection/adapt_advection_order.jl @@ -15,9 +15,9 @@ If the order of advection is changed in at least one direction, the adapted adve by this function is a `FluxFormAdvection`. """ function adapt_advection_order(advection, grid::AbstractGrid) - advection_x = adapt_advection_order_x(advection.x, topology(grid, 1), size(grid, 1), grid) - advection_y = adapt_advection_order_y(advection.y, topology(grid, 2), size(grid, 2), grid) - advection_z = adapt_advection_order_z(advection.z, topology(grid, 3), size(grid, 3), grid) + advection_x = adapt_advection_order(advection, topology(grid, 1), size(grid, 1), grid) + advection_y = adapt_advection_order(advection, topology(grid, 2), size(grid, 2), grid) + advection_z = adapt_advection_order(advection, topology(grid, 3), size(grid, 3), grid) # Check that we indeed changed the advection operator changed_x = advection_x != advection @@ -41,9 +41,9 @@ function adapt_advection_order(advection, grid::AbstractGrid) end function adapt_advection_order(advection::FluxFormAdvection, grid::AbstractGrid) - advection_x = adapt_advection_order_x(advection.x, topology(grid, 1), size(grid, 1), grid) - advection_y = adapt_advection_order_y(advection.y, topology(grid, 2), size(grid, 2), grid) - advection_z = adapt_advection_order_z(advection.z, topology(grid, 3), size(grid, 3), grid) + advection_x = adapt_advection_order(advection.x, topology(grid, 1), size(grid, 1), grid) + advection_y = adapt_advection_order(advection.y, topology(grid, 2), size(grid, 2), grid) + advection_z = adapt_advection_order(advection.z, topology(grid, 3), size(grid, 3), grid) # Check that we indeed changed the advection operator changed_x = advection_x != advection.x @@ -69,16 +69,11 @@ end # For the moment, we do not adapt the advection order for the VectorInvariant advection scheme adapt_advection_order(advection::VectorInvariant, grid::AbstractGrid) = advection -# We only need one halo in bounded directions! -adapt_advection_order_x(advection, topo, N, grid) = advection -adapt_advection_order_y(advection, topo, N, grid) = advection -adapt_advection_order_z(advection, topo, N, grid) = advection - ##### ##### Directional adapt advection order ##### -function adapt_advection_order_x(advection::Centered{H}, topology, N::Int, grid::AbstractGrid) where H +function adapt_advection_order(advection::Centered{H}, topology, N::Int, grid::AbstractGrid) where H if N == 1 return nothing elseif N >= H @@ -88,7 +83,7 @@ function adapt_advection_order_x(advection::Centered{H}, topology, N::Int, grid: end end -function adapt_advection_order_x(advection::UpwindBiased{H}, topology, grid::AbstractGrid, N::Int) where H +function adapt_advection_order(advection::UpwindBiased{H}, topology, N::Int, grid::AbstractGrid) where H if N == 1 return nothing elseif N >= H @@ -98,47 +93,23 @@ function adapt_advection_order_x(advection::UpwindBiased{H}, topology, grid::Abs end end -adapt_advection_order_y(advection, topology, grid::AbstractGrid, N::Int) = adapt_advection_order_x(advection, topology, grid, N) -adapt_advection_order_z(advection, topology, grid::AbstractGrid, N::Int) = adapt_advection_order_x(advection, topology, grid, N) - """ - new_weno_scheme(grid, order, bounds, T) + new_weno_scheme(grid, order, bounds, XT, YT, ZT) -Constructs a new WENO scheme based on the given parameters. `T` is the type of the weno coefficients. -A _non-stretched_ WENO scheme has `T` equal to `Nothing`. In case of a non-stretched WENO scheme, +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(grid, order, bounds, T) = ifelse(T == Nothing, WENO(; order, bounds), WENO(grid; order, bounds)) +new_weno_scheme(advection::WENO, grid, order, bounds, ::Type{Nothing}, ::Type{Nothing}, ::Type{Nothing},) = WENO(; order, bounds) +new_weno_scheme(advection::WENO, grid, order, bounds, XT, YT, ZT) = WENO(grid; order, bounds) -function adapt_advection_order_x(advection::WENO{H, FT, XT, YT, ZT}, topology, grid::AbstractGrid, N::Int) where {H, FT, XT, YT, ZT} +function adapt_advection_order(advection::WENO, topology, N::Int, grid::AbstractGrid) where {H, FT, XT, YT, ZT} if N == 1 return nothing elseif N >= H return advection else - return new_weno_scheme(grid, N * 2 - 1, advection.bounds, XT) - end -end - -function adapt_advection_order_y(advection::WENO{H, FT, XT, YT, ZT}, topology, grid::AbstractGrid, N::Int) where {H, FT, XT, YT, ZT} - - if N == 1 - return nothing - elseif N > H - return advection - else - return new_weno_scheme(grid, N * 2 - 1, advection.bounds, YT) + return new_weno_scheme(advection, grid, N * 2 - 1, advection.bounds, XT, YT, ZT) end -end - -function adapt_advection_order_z(advection::WENO{H, FT, XT, YT, ZT}, topology, grid::AbstractGrid, N::Int) where {H, FT, XT, YT, ZT} - - if N == 1 - return nothing - elseif N >= H - return advection - else - return new_weno_scheme(grid, N * 2 - 1, advection.bounds, XT) - end -end +end \ No newline at end of file From bcdb2c63cb348ed06d9f2d6ed6f4e18e8fe6f8cf Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Mon, 26 Aug 2024 14:46:37 -0400 Subject: [PATCH 22/68] formatting --- src/Advection/adapt_advection_order.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/Advection/adapt_advection_order.jl b/src/Advection/adapt_advection_order.jl index f6f8ba6d5b..163e4b9b46 100644 --- a/src/Advection/adapt_advection_order.jl +++ b/src/Advection/adapt_advection_order.jl @@ -100,8 +100,8 @@ Constructs a new WENO scheme based on the given parameters. `XT`, `YT`, and `ZT` 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(advection::WENO, grid, order, bounds, ::Type{Nothing}, ::Type{Nothing}, ::Type{Nothing},) = WENO(; order, bounds) -new_weno_scheme(advection::WENO, grid, order, bounds, XT, YT, ZT) = WENO(grid; order, bounds) +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, topology, N::Int, grid::AbstractGrid) where {H, FT, XT, YT, ZT} From d73b3d98ceb9a839534dd71b7bee298d46ae0d5b Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Mon, 26 Aug 2024 14:47:18 -0400 Subject: [PATCH 23/68] whoops --- src/Advection/adapt_advection_order.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/Advection/adapt_advection_order.jl b/src/Advection/adapt_advection_order.jl index 163e4b9b46..6ba848bbe4 100644 --- a/src/Advection/adapt_advection_order.jl +++ b/src/Advection/adapt_advection_order.jl @@ -103,7 +103,7 @@ we rebuild the advection without passing the grid information, otherwise we use 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, topology, N::Int, grid::AbstractGrid) where {H, FT, XT, YT, ZT} +function adapt_advection_order(advection::WENO{H, FT, XT, YT, ZT}, topology, N::Int, grid::AbstractGrid) where {H, FT, XT, YT, ZT} if N == 1 return nothing From d46eb4c68c74fd131b0c1f483ac7212fde546d97 Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Mon, 26 Aug 2024 15:49:08 -0400 Subject: [PATCH 24/68] remove topology and less metaprogramming --- src/Advection/adapt_advection_order.jl | 25 +++-- src/Advection/flux_form_advection.jl | 2 +- src/Advection/momentum_advection_operators.jl | 95 +++++++++++++++++-- src/Advection/tracer_advection_operators.jl | 39 ++++++-- 4 files changed, 130 insertions(+), 31 deletions(-) diff --git a/src/Advection/adapt_advection_order.jl b/src/Advection/adapt_advection_order.jl index 6ba848bbe4..df2057c0f3 100644 --- a/src/Advection/adapt_advection_order.jl +++ b/src/Advection/adapt_advection_order.jl @@ -15,9 +15,9 @@ If the order of advection is changed in at least one direction, the adapted adve by this function is a `FluxFormAdvection`. """ function adapt_advection_order(advection, grid::AbstractGrid) - advection_x = adapt_advection_order(advection, topology(grid, 1), size(grid, 1), grid) - advection_y = adapt_advection_order(advection, topology(grid, 2), size(grid, 2), grid) - advection_z = adapt_advection_order(advection, topology(grid, 3), size(grid, 3), grid) + advection_x = adapt_advection_order(advection, size(grid, 1), grid) + advection_y = adapt_advection_order(advection, size(grid, 2), grid) + advection_z = adapt_advection_order(advection, size(grid, 3), grid) # Check that we indeed changed the advection operator changed_x = advection_x != advection @@ -41,9 +41,9 @@ function adapt_advection_order(advection, grid::AbstractGrid) end function adapt_advection_order(advection::FluxFormAdvection, grid::AbstractGrid) - advection_x = adapt_advection_order(advection.x, topology(grid, 1), size(grid, 1), grid) - advection_y = adapt_advection_order(advection.y, topology(grid, 2), size(grid, 2), grid) - advection_z = adapt_advection_order(advection.z, topology(grid, 3), size(grid, 3), grid) + advection_x = adapt_advection_order(advection.x, size(grid, 1), grid) + advection_y = adapt_advection_order(advection.y, size(grid, 2), grid) + advection_z = adapt_advection_order(advection.z, size(grid, 3), grid) # Check that we indeed changed the advection operator changed_x = advection_x != advection.x @@ -73,8 +73,8 @@ adapt_advection_order(advection::VectorInvariant, grid::AbstractGrid) = advectio ##### Directional adapt advection order ##### -function adapt_advection_order(advection::Centered{H}, topology, N::Int, grid::AbstractGrid) where H - if N == 1 +function adapt_advection_order(advection::Centered{H}, N::Int, grid::AbstractGrid) where H + if N == 1 && H != 1 return nothing elseif N >= H return advection @@ -83,8 +83,8 @@ function adapt_advection_order(advection::Centered{H}, topology, N::Int, grid::A end end -function adapt_advection_order(advection::UpwindBiased{H}, topology, N::Int, grid::AbstractGrid) where H - if N == 1 +function adapt_advection_order(advection::UpwindBiased{H}, N::Int, grid::AbstractGrid) where H + if N == 1 && H != 1 return nothing elseif N >= H return advection @@ -103,9 +103,8 @@ we rebuild the advection without passing the grid information, otherwise we use 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{H, FT, XT, YT, ZT}, topology, N::Int, grid::AbstractGrid) where {H, FT, XT, YT, ZT} - - if N == 1 +function adapt_advection_order(advection::WENO{H, FT, XT, YT, ZT}, N::Int, grid::AbstractGrid) where {H, FT, XT, YT, ZT} + if N == 1 && H != 1 return nothing elseif N >= H return advection diff --git a/src/Advection/flux_form_advection.jl b/src/Advection/flux_form_advection.jl index aa1931a53d..faa35fb668 100644 --- a/src/Advection/flux_form_advection.jl +++ b/src/Advection/flux_form_advection.jl @@ -55,4 +55,4 @@ Adapt.adapt_structure(to, scheme::FluxFormAdvection{N, FT}) where {N, FT} = @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...) \ No newline at end of file +@inline _advective_momentum_flux_Ww(i, j, k, grid, advection::FluxFormAdvection, args...) = _advective_momentum_flux_Ww(i, j, k, grid, advection.z, args...) diff --git a/src/Advection/momentum_advection_operators.jl b/src/Advection/momentum_advection_operators.jl index b475796146..673d31b362 100644 --- a/src/Advection/momentum_advection_operators.jl +++ b/src/Advection/momentum_advection_operators.jl @@ -94,15 +94,92 @@ end ##### Fallback advection fluxes! ##### -for flux_type in (:Uu, :Vu, :Wu, :Uv, :Vv, :Wv, :Uw, :Vw, :Ww) - advective_momentum_flux = Symbol(:_advective_momentum_flux_, flux_type) - +# 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) + +# Fallback for `nothing` advection and `ZeroField` tracers and velocities +@inline _advective_momentum_flux_Uu(i, j, k, grid, ::Nothing, ::ZeroField, ::ZeroField) = zero(grid) +@inline _advective_momentum_flux_Uv(i, j, k, grid, ::Nothing, ::ZeroField, ::ZeroField) = zero(grid) +@inline _advective_momentum_flux_Uw(i, j, k, grid, ::Nothing, ::ZeroField, ::ZeroField) = zero(grid) + +@inline _advective_momentum_flux_Vu(i, j, k, grid, ::Nothing, ::ZeroField, ::ZeroField) = zero(grid) +@inline _advective_momentum_flux_Vv(i, j, k, grid, ::Nothing, ::ZeroField, ::ZeroField) = zero(grid) +@inline _advective_momentum_flux_Vw(i, j, k, grid, ::Nothing, ::ZeroField, ::ZeroField) = zero(grid) + +@inline _advective_momentum_flux_Wu(i, j, k, grid, ::Nothing, ::ZeroField, ::ZeroField) = zero(grid) +@inline _advective_momentum_flux_Wv(i, j, k, grid, ::Nothing, ::ZeroField, ::ZeroField) = zero(grid) +@inline _advective_momentum_flux_Ww(i, j, k, grid, ::Nothing, ::ZeroField, ::ZeroField) = zero(grid) + +@inline _advective_momentum_flux_Uu(i, j, k, grid, ::Nothing, U, ::ZeroField) = zero(grid) +@inline _advective_momentum_flux_Uv(i, j, k, grid, ::Nothing, U, ::ZeroField) = zero(grid) +@inline _advective_momentum_flux_Uw(i, j, k, grid, ::Nothing, U, ::ZeroField) = zero(grid) +@inline _advective_momentum_flux_Uu(i, j, k, grid, ::Nothing, ::ZeroField, u) = zero(grid) +@inline _advective_momentum_flux_Uv(i, j, k, grid, ::Nothing, ::ZeroField, v) = zero(grid) +@inline _advective_momentum_flux_Uw(i, j, k, grid, ::Nothing, ::ZeroField, w) = zero(grid) + +@inline _advective_momentum_flux_Vu(i, j, k, grid, ::Nothing, V, ::ZeroField) = zero(grid) +@inline _advective_momentum_flux_Vv(i, j, k, grid, ::Nothing, V, ::ZeroField) = zero(grid) +@inline _advective_momentum_flux_Vw(i, j, k, grid, ::Nothing, V, ::ZeroField) = zero(grid) +@inline _advective_momentum_flux_Vu(i, j, k, grid, ::Nothing, ::ZeroField, u) = zero(grid) +@inline _advective_momentum_flux_Vv(i, j, k, grid, ::Nothing, ::ZeroField, v) = zero(grid) +@inline _advective_momentum_flux_Vw(i, j, k, grid, ::Nothing, ::ZeroField, w) = zero(grid) + +@inline _advective_momentum_flux_Wu(i, j, k, grid, ::Nothing, W, ::ZeroField) = zero(grid) +@inline _advective_momentum_flux_Wv(i, j, k, grid, ::Nothing, W, ::ZeroField) = zero(grid) +@inline _advective_momentum_flux_Ww(i, j, k, grid, ::Nothing, W, ::ZeroField) = zero(grid) +@inline _advective_momentum_flux_Wu(i, j, k, grid, ::Nothing, ::ZeroField, u) = zero(grid) +@inline _advective_momentum_flux_Wv(i, j, k, grid, ::Nothing, ::ZeroField, v) = zero(grid) +@inline _advective_momentum_flux_Ww(i, j, k, grid, ::Nothing, ::ZeroField, w) = zero(grid) + +for scheme in (:UpwindBiased, :Centered, :WENO, :FluxFormAdvection) @eval begin - @inline $advective_momentum_flux(i, j, k, grid, ::Nothing, args...) = zero(grid) - @inline $advective_momentum_flux(i, j, k, grid, ::Nothing, ::ZeroField, ::ZeroField) = zero(grid) - @inline $advective_momentum_flux(i, j, k, grid, ::Nothing, U, ::ZeroField) = zero(grid) - @inline $advective_momentum_flux(i, j, k, grid, scheme, ::ZeroField, ::ZeroField) = zero(grid) - @inline $advective_momentum_flux(i, j, k, grid, scheme, U, ::ZeroField) = zero(grid) - @inline $advective_momentum_flux(i, j, k, grid, scheme, ::ZeroField, u) = zero(grid) + # Fallback for `ZeroField` tracers and velocities + @inline _advective_momentum_flux_Uu(i, j, k, grid, ::$Scheme, ::ZeroField, ::ZeroField) = zero(grid) + @inline _advective_momentum_flux_Uv(i, j, k, grid, ::$Scheme, ::ZeroField, ::ZeroField) = zero(grid) + @inline _advective_momentum_flux_Uw(i, j, k, grid, ::$Scheme, ::ZeroField, ::ZeroField) = zero(grid) + + @inline _advective_momentum_flux_Vu(i, j, k, grid, ::$Scheme, ::ZeroField, ::ZeroField) = zero(grid) + @inline _advective_momentum_flux_Vv(i, j, k, grid, ::$Scheme, ::ZeroField, ::ZeroField) = zero(grid) + @inline _advective_momentum_flux_Vw(i, j, k, grid, ::$Scheme, ::ZeroField, ::ZeroField) = zero(grid) + + @inline _advective_momentum_flux_Wu(i, j, k, grid, ::$Scheme, ::ZeroField, ::ZeroField) = zero(grid) + @inline _advective_momentum_flux_Wv(i, j, k, grid, ::$Scheme, ::ZeroField, ::ZeroField) = zero(grid) + @inline _advective_momentum_flux_Ww(i, j, k, grid, ::$Scheme, ::ZeroField, ::ZeroField) = zero(grid) + + # Fallback for `ZeroField` tracers + @inline _advective_momentum_flux_Uu(i, j, k, grid, ::$Scheme, U, ::ZeroField) = zero(grid) + @inline _advective_momentum_flux_Uv(i, j, k, grid, ::$Scheme, U, ::ZeroField) = zero(grid) + @inline _advective_momentum_flux_Uw(i, j, k, grid, ::$Scheme, U, ::ZeroField) = zero(grid) + + @inline _advective_momentum_flux_Vu(i, j, k, grid, ::$Scheme, V, ::ZeroField) = zero(grid) + @inline _advective_momentum_flux_Vv(i, j, k, grid, ::$Scheme, V, ::ZeroField) = zero(grid) + @inline _advective_momentum_flux_Vw(i, j, k, grid, ::$Scheme, V, ::ZeroField) = zero(grid) + + @inline _advective_momentum_flux_Wu(i, j, k, grid, ::$Scheme, W, ::ZeroField) = zero(grid) + @inline _advective_momentum_flux_Wv(i, j, k, grid, ::$Scheme, W, ::ZeroField) = zero(grid) + @inline _advective_momentum_flux_Ww(i, j, k, grid, ::$Scheme, W, ::ZeroField) = zero(grid) + + # Fallback for `ZeroField` velocities + @inline _advective_momentum_flux_Uu(i, j, k, grid, ::$Scheme, ::ZeroField, u) = zero(grid) + @inline _advective_momentum_flux_Uv(i, j, k, grid, ::$Scheme, ::ZeroField, v) = zero(grid) + @inline _advective_momentum_flux_Uw(i, j, k, grid, ::$Scheme, ::ZeroField, w) = zero(grid) + + @inline _advective_momentum_flux_Vu(i, j, k, grid, ::$Scheme, ::ZeroField, u) = zero(grid) + @inline _advective_momentum_flux_Vv(i, j, k, grid, ::$Scheme, ::ZeroField, v) = zero(grid) + @inline _advective_momentum_flux_Vw(i, j, k, grid, ::$Scheme, ::ZeroField, w) = zero(grid) + + @inline _advective_momentum_flux_Wu(i, j, k, grid, ::$Scheme, ::ZeroField, u) = zero(grid) + @inline _advective_momentum_flux_Wv(i, j, k, grid, ::$Scheme, ::ZeroField, v) = zero(grid) + @inline _advective_momentum_flux_Ww(i, j, k, grid, ::$Scheme, ::ZeroField, w) = zero(grid) end end \ No newline at end of file diff --git a/src/Advection/tracer_advection_operators.jl b/src/Advection/tracer_advection_operators.jl index 239628b481..bac920b72d 100644 --- a/src/Advection/tracer_advection_operators.jl +++ b/src/Advection/tracer_advection_operators.jl @@ -7,16 +7,39 @@ ##### Fallback tracer fluxes! ##### -for flux_dir in (:x, :y, :z) - advective_tracer_flux = Symbol(:_advective_tracer_flux_, flux_dir) +# 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) + +for scheme in (:UpwindBiased, :Centered, :WENO, :FluxFormAdvection) @eval begin - @inline $advective_tracer_flux(i, j, k, grid, ::Nothing, args...) = zero(grid) - @inline $advective_tracer_flux(i, j, k, grid, ::Nothing, ::ZeroField, ::ZeroField) = zero(grid) - @inline $advective_tracer_flux(i, j, k, grid, ::Nothing, U, ::ZeroField) = zero(grid) - @inline $advective_tracer_flux(i, j, k, grid, scheme, ::ZeroField, ::ZeroField) = zero(grid) - @inline $advective_tracer_flux(i, j, k, grid, scheme, U, ::ZeroField) = zero(grid) - @inline $advective_tracer_flux(i, j, k, grid, scheme, ::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) end end From 1b028f30241243bfb2e36511e813c3c36a634f47 Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Mon, 26 Aug 2024 17:26:48 -0400 Subject: [PATCH 25/68] add advection --- src/Advection/momentum_advection_operators.jl | 2 +- src/Advection/tracer_advection_operators.jl | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/src/Advection/momentum_advection_operators.jl b/src/Advection/momentum_advection_operators.jl index 673d31b362..e0548b4443 100644 --- a/src/Advection/momentum_advection_operators.jl +++ b/src/Advection/momentum_advection_operators.jl @@ -141,7 +141,7 @@ end @inline _advective_momentum_flux_Wv(i, j, k, grid, ::Nothing, ::ZeroField, v) = zero(grid) @inline _advective_momentum_flux_Ww(i, j, k, grid, ::Nothing, ::ZeroField, w) = zero(grid) -for scheme in (:UpwindBiased, :Centered, :WENO, :FluxFormAdvection) +for Scheme in (:UpwindBiased, :Centered, :WENO, :FluxFormAdvection) @eval begin # Fallback for `ZeroField` tracers and velocities @inline _advective_momentum_flux_Uu(i, j, k, grid, ::$Scheme, ::ZeroField, ::ZeroField) = zero(grid) diff --git a/src/Advection/tracer_advection_operators.jl b/src/Advection/tracer_advection_operators.jl index bac920b72d..c104d43ca9 100644 --- a/src/Advection/tracer_advection_operators.jl +++ b/src/Advection/tracer_advection_operators.jl @@ -24,7 +24,7 @@ @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) -for scheme in (:UpwindBiased, :Centered, :WENO, :FluxFormAdvection) +for Scheme in (:UpwindBiased, :Centered, :WENO, :FluxFormAdvection) @eval begin # Fallback for `ZeroField` tracers and velocities @inline _advective_tracer_flux_x(i, j, k, grid, ::$Scheme, ::ZeroField, ::ZeroField) = zero(grid) From 94568812fc053ee13ecdd6b978ecc48360144a91 Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Mon, 26 Aug 2024 17:59:27 -0400 Subject: [PATCH 26/68] dont adapt nothing --- src/Advection/adapt_advection_order.jl | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/Advection/adapt_advection_order.jl b/src/Advection/adapt_advection_order.jl index df2057c0f3..f1a34bf670 100644 --- a/src/Advection/adapt_advection_order.jl +++ b/src/Advection/adapt_advection_order.jl @@ -68,6 +68,8 @@ end # 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 From 3bc14be029349778138578da16b11a8d57dbcdbf Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Mon, 26 Aug 2024 20:07:03 -0400 Subject: [PATCH 27/68] retain same behavior as before --- src/Advection/adapt_advection_order.jl | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/Advection/adapt_advection_order.jl b/src/Advection/adapt_advection_order.jl index f1a34bf670..2f04b3135a 100644 --- a/src/Advection/adapt_advection_order.jl +++ b/src/Advection/adapt_advection_order.jl @@ -77,7 +77,7 @@ adapt_advection_order(advection::Nothing, N::Int, grid::AbstractGrid) = nothing function adapt_advection_order(advection::Centered{H}, N::Int, grid::AbstractGrid) where H if N == 1 && H != 1 - return nothing + return Centered(; order = 2) elseif N >= H return advection else @@ -87,7 +87,7 @@ end function adapt_advection_order(advection::UpwindBiased{H}, N::Int, grid::AbstractGrid) where H if N == 1 && H != 1 - return nothing + return UpwindBiased(; order = 1) elseif N >= H return advection else @@ -107,7 +107,7 @@ new_weno_scheme(::WENO, grid, order, bounds, XT, YT, ZT) function adapt_advection_order(advection::WENO{H, FT, XT, YT, ZT}, N::Int, grid::AbstractGrid) where {H, FT, XT, YT, ZT} if N == 1 && H != 1 - return nothing + return UpwindBiased(; order = 1) elseif N >= H return advection else From ce30f43df7ea11c7aa13109fdaee25333bced96d Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Mon, 26 Aug 2024 20:09:45 -0400 Subject: [PATCH 28/68] fix multi region advection --- src/MultiRegion/multi_region_models.jl | 7 ++++++- 1 file changed, 6 insertions(+), 1 deletion(-) diff --git a/src/MultiRegion/multi_region_models.jl b/src/MultiRegion/multi_region_models.jl index 3dabe07602..960a03e73e 100644 --- a/src/MultiRegion/multi_region_models.jl +++ b/src/MultiRegion/multi_region_models.jl @@ -9,12 +9,17 @@ using Oceananigans.Advection: OnlySelfUpwinding, CrossAndSelfUpwinding using Oceananigans.ImmersedBoundaries: GridFittedBottom, PartialCellBottom, GridFittedBoundary using Oceananigans.Solvers: PreconditionedConjugateGradientSolver -import Oceananigans.Advection: WENO, cell_advection_timescale +import Oceananigans.Advection: WENO, cell_advection_timescale, adapt_advection_order import Oceananigans.Models.HydrostaticFreeSurfaceModels: build_implicit_step_solver, validate_tracer_advection import Oceananigans.TurbulenceClosures: implicit_diffusion_solver const MultiRegionModel = HydrostaticFreeSurfaceModel{<:Any, <:Any, <:AbstractArchitecture, <:Any, <:MultiRegionGrids} +function adapt_advection_order(advection::MultiRegionObject, grid::MultiRegionGrids) + @apply_regionally new_advection = adapt_advection_order(advection, grid) + return new_advection +end + # Utility to generate the inputs to complex `getregion`s function getregionalproperties(T, inner=true) type = getglobal(@__MODULE__, T) From cd9691169133dcd2cf3804bf54fd561db25aec1e Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Tue, 27 Aug 2024 09:06:33 -0400 Subject: [PATCH 29/68] bugfix scalar diffusivities parameter types --- .../scalar_biharmonic_diffusivity.jl | 6 +++--- .../scalar_diffusivity.jl | 6 +++--- 2 files changed, 6 insertions(+), 6 deletions(-) diff --git a/src/TurbulenceClosures/turbulence_closure_implementations/scalar_biharmonic_diffusivity.jl b/src/TurbulenceClosures/turbulence_closure_implementations/scalar_biharmonic_diffusivity.jl index 677b4ebbc3..9aff3a615f 100644 --- a/src/TurbulenceClosures/turbulence_closure_implementations/scalar_biharmonic_diffusivity.jl +++ b/src/TurbulenceClosures/turbulence_closure_implementations/scalar_biharmonic_diffusivity.jl @@ -6,10 +6,10 @@ using Oceananigans.Utils: prettysummary Holds viscosity and diffusivities for models with prescribed isotropic diffusivities. """ -struct ScalarBiharmonicDiffusivity{F, V, K, N} <: AbstractScalarBiharmonicDiffusivity{F, N} +struct ScalarBiharmonicDiffusivity{F, N, V, K} <: AbstractScalarBiharmonicDiffusivity{F, N} ν :: V κ :: K - ScalarBiharmonicDiffusivity{F, N}(ν::V, κ::K) where {F, V, K, N} = new{F, V, K, N}(ν, κ) + ScalarBiharmonicDiffusivity{F, N}(ν::V, κ::K) where {F, V, K, N} = new{F, N, V, K}(ν, κ) end # Aliases that allow specify the floating type, assuming that the discretization is Explicit in time @@ -82,7 +82,7 @@ function ScalarBiharmonicDiffusivity(formulation = ThreeDimensionalFormulation() return ScalarBiharmonicDiffusivity{typeof(formulation), required_halo_size}(ν, κ) end -function with_tracers(tracers, closure::ScalarBiharmonicDiffusivity{F, N}) where {F, N} +function with_tracers(tracers, closure::ScalarBiharmonicDiffusivity{F, N, V, K}) where {F, N, V, K} κ = tracer_diffusivities(tracers, closure.κ) return ScalarBiharmonicDiffusivity{F, N}(closure.ν, κ) end diff --git a/src/TurbulenceClosures/turbulence_closure_implementations/scalar_diffusivity.jl b/src/TurbulenceClosures/turbulence_closure_implementations/scalar_diffusivity.jl index c5f4a71de2..8466b02528 100644 --- a/src/TurbulenceClosures/turbulence_closure_implementations/scalar_diffusivity.jl +++ b/src/TurbulenceClosures/turbulence_closure_implementations/scalar_diffusivity.jl @@ -3,10 +3,10 @@ using Oceananigans.Utils: prettysummary import Adapt import Oceananigans.Grids: required_halo_size_x, required_halo_size_y, required_halo_size_z -struct ScalarDiffusivity{TD, F, V, K, N} <: AbstractScalarDiffusivity{TD, F, N} +struct ScalarDiffusivity{TD, N, F, V, K} <: AbstractScalarDiffusivity{TD, F, N} ν :: V κ :: K - ScalarDiffusivity{TD, F, N}(ν::V, κ::K) where {TD, F, V, K, N} = new{TD, F, V, K, N}(ν, κ) + ScalarDiffusivity{TD, F, N}(ν::V, κ::K) where {TD, F, V, K, N} = new{TD, F, N, V, K}(ν, κ) end """ @@ -173,7 +173,7 @@ Shorthand for a `ScalarDiffusivity` with `HorizontalDivergenceFormulation()`. Se HorizontalScalarDiffusivity(FT::DataType; kwargs...) = ScalarDiffusivity(ExplicitTimeDiscretization(), HorizontalFormulation(), FT; kwargs...) HorizontalDivergenceScalarDiffusivity(FT::DataType; kwargs...) = ScalarDiffusivity(ExplicitTimeDiscretization(), HorizontalDivergenceFormulation(), FT; kwargs...) -@inline function with_tracers(tracers, closure::ScalarDiffusivity{TD, F, N}) where {TD, F, N} +@inline function with_tracers(tracers, closure::ScalarDiffusivity{TD, F, N, V, K}) where {TD, F, N, V, K} κ = tracer_diffusivities(tracers, closure.κ) return ScalarDiffusivity{TD, F, N}(closure.ν, κ) end From 7191393240ed6a77dcd97cb9f92c2444c57390e5 Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Tue, 27 Aug 2024 09:09:18 -0400 Subject: [PATCH 30/68] bugfix --- .../turbulence_closure_implementations/scalar_diffusivity.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/TurbulenceClosures/turbulence_closure_implementations/scalar_diffusivity.jl b/src/TurbulenceClosures/turbulence_closure_implementations/scalar_diffusivity.jl index 8466b02528..cd6b0b697d 100644 --- a/src/TurbulenceClosures/turbulence_closure_implementations/scalar_diffusivity.jl +++ b/src/TurbulenceClosures/turbulence_closure_implementations/scalar_diffusivity.jl @@ -3,7 +3,7 @@ using Oceananigans.Utils: prettysummary import Adapt import Oceananigans.Grids: required_halo_size_x, required_halo_size_y, required_halo_size_z -struct ScalarDiffusivity{TD, N, F, V, K} <: AbstractScalarDiffusivity{TD, F, N} +struct ScalarDiffusivity{TD, F, N, V, K} <: AbstractScalarDiffusivity{TD, F, N} ν :: V κ :: K ScalarDiffusivity{TD, F, N}(ν::V, κ::K) where {TD, F, V, K, N} = new{TD, F, N, V, K}(ν, κ) From e6a489ca6cdf2528ad6a7bb33c63765a8a1b087e Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Tue, 27 Aug 2024 10:21:09 -0400 Subject: [PATCH 31/68] enzyme error --- .../turbulence_closure_implementations/scalar_diffusivity.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/TurbulenceClosures/turbulence_closure_implementations/scalar_diffusivity.jl b/src/TurbulenceClosures/turbulence_closure_implementations/scalar_diffusivity.jl index cd6b0b697d..7227209cb7 100644 --- a/src/TurbulenceClosures/turbulence_closure_implementations/scalar_diffusivity.jl +++ b/src/TurbulenceClosures/turbulence_closure_implementations/scalar_diffusivity.jl @@ -6,7 +6,7 @@ import Oceananigans.Grids: required_halo_size_x, required_halo_size_y, required_ struct ScalarDiffusivity{TD, F, N, V, K} <: AbstractScalarDiffusivity{TD, F, N} ν :: V κ :: K - ScalarDiffusivity{TD, F, N}(ν::V, κ::K) where {TD, F, V, K, N} = new{TD, F, N, V, K}(ν, κ) + ScalarDiffusivity{TD, F, N}(ν::V, κ::K) where {TD, F, N, V, K} = new{TD, F, N, V, K}(ν, κ) end """ @@ -173,7 +173,7 @@ Shorthand for a `ScalarDiffusivity` with `HorizontalDivergenceFormulation()`. Se HorizontalScalarDiffusivity(FT::DataType; kwargs...) = ScalarDiffusivity(ExplicitTimeDiscretization(), HorizontalFormulation(), FT; kwargs...) HorizontalDivergenceScalarDiffusivity(FT::DataType; kwargs...) = ScalarDiffusivity(ExplicitTimeDiscretization(), HorizontalDivergenceFormulation(), FT; kwargs...) -@inline function with_tracers(tracers, closure::ScalarDiffusivity{TD, F, N, V, K}) where {TD, F, N, V, K} +@inline function with_tracers(tracers, closure::ScalarDiffusivity{TD, F, N}) where {TD, F, N} κ = tracer_diffusivities(tracers, closure.κ) return ScalarDiffusivity{TD, F, N}(closure.ν, κ) end From e041ea6e6e25b448a9ccead79b9c057e8599eaa1 Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Tue, 3 Sep 2024 21:54:32 -0400 Subject: [PATCH 32/68] make sure constructors are type-stable... --- .../isopycnal_skew_symmetric_diffusivity.jl | 10 +++++----- .../scalar_biharmonic_diffusivity.jl | 7 +++++-- .../scalar_diffusivity.jl | 7 +++++-- 3 files changed, 15 insertions(+), 9 deletions(-) diff --git a/src/TurbulenceClosures/turbulence_closure_implementations/isopycnal_skew_symmetric_diffusivity.jl b/src/TurbulenceClosures/turbulence_closure_implementations/isopycnal_skew_symmetric_diffusivity.jl index e79a1a76d1..3b36933f81 100644 --- a/src/TurbulenceClosures/turbulence_closure_implementations/isopycnal_skew_symmetric_diffusivity.jl +++ b/src/TurbulenceClosures/turbulence_closure_implementations/isopycnal_skew_symmetric_diffusivity.jl @@ -37,15 +37,15 @@ function IsopycnalSkewSymmetricDiffusivity(time_disc::TD = VerticallyImplicitTim κ_symmetric = 0, isopycnal_tensor = SmallSlopeIsopycnalTensor(), slope_limiter = FluxTapering(1e-2), - required_halo_size = 1) where TD + required_halo_size::Val{N} = Val(1)) where {TD, N} isopycnal_tensor isa SmallSlopeIsopycnalTensor || error("Only isopycnal_tensor=SmallSlopeIsopycnalTensor() is currently supported.") - return IsopycnalSkewSymmetricDiffusivity{TD, required_halo_size}(convert_diffusivity(FT, κ_skew), - convert_diffusivity(FT, κ_symmetric), - isopycnal_tensor, - slope_limiter) + return IsopycnalSkewSymmetricDiffusivity{TD, N}(convert_diffusivity(FT, κ_skew), + convert_diffusivity(FT, κ_symmetric), + isopycnal_tensor, + slope_limiter) end IsopycnalSkewSymmetricDiffusivity(FT::DataType; kw...) = diff --git a/src/TurbulenceClosures/turbulence_closure_implementations/scalar_biharmonic_diffusivity.jl b/src/TurbulenceClosures/turbulence_closure_implementations/scalar_biharmonic_diffusivity.jl index 9aff3a615f..ae2702916f 100644 --- a/src/TurbulenceClosures/turbulence_closure_implementations/scalar_biharmonic_diffusivity.jl +++ b/src/TurbulenceClosures/turbulence_closure_implementations/scalar_biharmonic_diffusivity.jl @@ -50,6 +50,9 @@ Keyword arguments * `discrete_form`: `Boolean`; default: `false`. +* `required_halo_size = Val(2)`: `Val(i::Int)` where `i` is the required halo size for the closure. + change only if using a function for `ν` or `κ` that requires a halo size larger than 1 to compute. + When prescribing the viscosities or diffusivities as functions, depending on the value of keyword argument `discrete_form`, the constructor expects: @@ -75,11 +78,11 @@ function ScalarBiharmonicDiffusivity(formulation = ThreeDimensionalFormulation() discrete_form = false, loc = (nothing, nothing, nothing), parameters = nothing, - required_halo_size = 2) + required_halo_size::Val{N} = Val(2)) where N ν = convert_diffusivity(FT, ν; discrete_form, loc, parameters) κ = convert_diffusivity(FT, κ; discrete_form, loc, parameters) - return ScalarBiharmonicDiffusivity{typeof(formulation), required_halo_size}(ν, κ) + return ScalarBiharmonicDiffusivity{typeof(formulation), N}(ν, κ) end function with_tracers(tracers, closure::ScalarBiharmonicDiffusivity{F, N, V, K}) where {F, N, V, K} diff --git a/src/TurbulenceClosures/turbulence_closure_implementations/scalar_diffusivity.jl b/src/TurbulenceClosures/turbulence_closure_implementations/scalar_diffusivity.jl index 7227209cb7..a036dbca61 100644 --- a/src/TurbulenceClosures/turbulence_closure_implementations/scalar_diffusivity.jl +++ b/src/TurbulenceClosures/turbulence_closure_implementations/scalar_diffusivity.jl @@ -63,6 +63,9 @@ value of keyword argument `discrete_form`, the constructor expects: - with `loc = (ℓx, ℓy, ℓz)` and specified `parameters`: functions of `(i, j, k, grid, clock, fields, parameters)`. +* `required_halo_size = Val(1)`: `Val(i::Int)` where `i` is the required halo size for the closure. + change only if using a function for `ν` or `κ` that requires a halo size larger than 1 to compute. + * `parameters`: `NamedTuple` with parameters used by the functions that compute viscosity and/or diffusivity; default: `nothing`. @@ -116,7 +119,7 @@ function ScalarDiffusivity(time_discretization=ExplicitTimeDiscretization(), discrete_form = false, loc = (nothing, nothing, nothing), parameters = nothing, - required_halo_size = 1) + required_halo_size::Val{N} = Val(1)) where N if formulation == HorizontalFormulation() && time_discretization == VerticallyImplicitTimeDiscretization() throw(ArgumentError("VerticallyImplicitTimeDiscretization is only supported for \ @@ -126,7 +129,7 @@ function ScalarDiffusivity(time_discretization=ExplicitTimeDiscretization(), κ = convert_diffusivity(FT, κ; discrete_form, loc, parameters) ν = convert_diffusivity(FT, ν; discrete_form, loc, parameters) - return ScalarDiffusivity{typeof(time_discretization), typeof(formulation), required_halo_size}(ν, κ) + return ScalarDiffusivity{typeof(time_discretization), typeof(formulation), N}(ν, κ) end # Explicit default From a9489556bdb7ddf1b8731112d5da6790c475a73d Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Tue, 3 Sep 2024 22:12:05 -0400 Subject: [PATCH 33/68] force required_halo_size to be an integer --- .../scalar_biharmonic_diffusivity.jl | 4 ++-- .../turbulence_closure_implementations/scalar_diffusivity.jl | 5 +++-- 2 files changed, 5 insertions(+), 4 deletions(-) diff --git a/src/TurbulenceClosures/turbulence_closure_implementations/scalar_biharmonic_diffusivity.jl b/src/TurbulenceClosures/turbulence_closure_implementations/scalar_biharmonic_diffusivity.jl index ae2702916f..341ebb220c 100644 --- a/src/TurbulenceClosures/turbulence_closure_implementations/scalar_biharmonic_diffusivity.jl +++ b/src/TurbulenceClosures/turbulence_closure_implementations/scalar_biharmonic_diffusivity.jl @@ -78,11 +78,11 @@ function ScalarBiharmonicDiffusivity(formulation = ThreeDimensionalFormulation() discrete_form = false, loc = (nothing, nothing, nothing), parameters = nothing, - required_halo_size::Val{N} = Val(2)) where N + required_halo_size::Int = 2) ν = convert_diffusivity(FT, ν; discrete_form, loc, parameters) κ = convert_diffusivity(FT, κ; discrete_form, loc, parameters) - return ScalarBiharmonicDiffusivity{typeof(formulation), N}(ν, κ) + return ScalarBiharmonicDiffusivity{typeof(formulation), required_halo_size}(ν, κ) end function with_tracers(tracers, closure::ScalarBiharmonicDiffusivity{F, N, V, K}) where {F, N, V, K} diff --git a/src/TurbulenceClosures/turbulence_closure_implementations/scalar_diffusivity.jl b/src/TurbulenceClosures/turbulence_closure_implementations/scalar_diffusivity.jl index a036dbca61..ebcddcd00c 100644 --- a/src/TurbulenceClosures/turbulence_closure_implementations/scalar_diffusivity.jl +++ b/src/TurbulenceClosures/turbulence_closure_implementations/scalar_diffusivity.jl @@ -65,6 +65,7 @@ value of keyword argument `discrete_form`, the constructor expects: * `required_halo_size = Val(1)`: `Val(i::Int)` where `i` is the required halo size for the closure. change only if using a function for `ν` or `κ` that requires a halo size larger than 1 to compute. + The required halo size has to be provided as a `Val` type, e.g., `Val(2)` if the closure requires a halo size of 2. * `parameters`: `NamedTuple` with parameters used by the functions that compute viscosity and/or diffusivity; default: `nothing`. @@ -119,7 +120,7 @@ function ScalarDiffusivity(time_discretization=ExplicitTimeDiscretization(), discrete_form = false, loc = (nothing, nothing, nothing), parameters = nothing, - required_halo_size::Val{N} = Val(1)) where N + required_halo_size::Int = 1) if formulation == HorizontalFormulation() && time_discretization == VerticallyImplicitTimeDiscretization() throw(ArgumentError("VerticallyImplicitTimeDiscretization is only supported for \ @@ -129,7 +130,7 @@ function ScalarDiffusivity(time_discretization=ExplicitTimeDiscretization(), κ = convert_diffusivity(FT, κ; discrete_form, loc, parameters) ν = convert_diffusivity(FT, ν; discrete_form, loc, parameters) - return ScalarDiffusivity{typeof(time_discretization), typeof(formulation), N}(ν, κ) + return ScalarDiffusivity{typeof(time_discretization), typeof(formulation), required_halo_size}(ν, κ) end # Explicit default From 1ffe23993863102f9134193aa35428665396a9ce Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Tue, 3 Sep 2024 22:14:17 -0400 Subject: [PATCH 34/68] add better comment --- .../scalar_biharmonic_diffusivity.jl | 2 +- .../turbulence_closure_implementations/scalar_diffusivity.jl | 3 +-- 2 files changed, 2 insertions(+), 3 deletions(-) diff --git a/src/TurbulenceClosures/turbulence_closure_implementations/scalar_biharmonic_diffusivity.jl b/src/TurbulenceClosures/turbulence_closure_implementations/scalar_biharmonic_diffusivity.jl index 341ebb220c..e35ecc649c 100644 --- a/src/TurbulenceClosures/turbulence_closure_implementations/scalar_biharmonic_diffusivity.jl +++ b/src/TurbulenceClosures/turbulence_closure_implementations/scalar_biharmonic_diffusivity.jl @@ -50,7 +50,7 @@ Keyword arguments * `discrete_form`: `Boolean`; default: `false`. -* `required_halo_size = Val(2)`: `Val(i::Int)` where `i` is the required halo size for the closure. +* `required_halo_size = 2`: the required halo size for the closure. This value should be an integer. change only if using a function for `ν` or `κ` that requires a halo size larger than 1 to compute. When prescribing the viscosities or diffusivities as functions, depending on the diff --git a/src/TurbulenceClosures/turbulence_closure_implementations/scalar_diffusivity.jl b/src/TurbulenceClosures/turbulence_closure_implementations/scalar_diffusivity.jl index ebcddcd00c..2e4acf9795 100644 --- a/src/TurbulenceClosures/turbulence_closure_implementations/scalar_diffusivity.jl +++ b/src/TurbulenceClosures/turbulence_closure_implementations/scalar_diffusivity.jl @@ -63,9 +63,8 @@ value of keyword argument `discrete_form`, the constructor expects: - with `loc = (ℓx, ℓy, ℓz)` and specified `parameters`: functions of `(i, j, k, grid, clock, fields, parameters)`. -* `required_halo_size = Val(1)`: `Val(i::Int)` where `i` is the required halo size for the closure. +* `required_halo_size = Val(1)`: the required halo size for the closure. This value should be an integer. change only if using a function for `ν` or `κ` that requires a halo size larger than 1 to compute. - The required halo size has to be provided as a `Val` type, e.g., `Val(2)` if the closure requires a halo size of 2. * `parameters`: `NamedTuple` with parameters used by the functions that compute viscosity and/or diffusivity; default: `nothing`. From 39407f6c9d77437aa6012b858523ac97d5123969 Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Tue, 3 Sep 2024 22:19:55 -0400 Subject: [PATCH 35/68] add a couple of unit tests --- .../scalar_diffusivity.jl | 2 +- test/test_turbulence_closures.jl | 20 +++++++++++++++++-- 2 files changed, 19 insertions(+), 3 deletions(-) diff --git a/src/TurbulenceClosures/turbulence_closure_implementations/scalar_diffusivity.jl b/src/TurbulenceClosures/turbulence_closure_implementations/scalar_diffusivity.jl index 2e4acf9795..2d54b16cd8 100644 --- a/src/TurbulenceClosures/turbulence_closure_implementations/scalar_diffusivity.jl +++ b/src/TurbulenceClosures/turbulence_closure_implementations/scalar_diffusivity.jl @@ -63,7 +63,7 @@ value of keyword argument `discrete_form`, the constructor expects: - with `loc = (ℓx, ℓy, ℓz)` and specified `parameters`: functions of `(i, j, k, grid, clock, fields, parameters)`. -* `required_halo_size = Val(1)`: the required halo size for the closure. This value should be an integer. +* `required_halo_size = 1`: the required halo size for the closure. This value should be an integer. change only if using a function for `ν` or `κ` that requires a halo size larger than 1 to compute. * `parameters`: `NamedTuple` with parameters used by the functions diff --git a/test/test_turbulence_closures.jl b/test/test_turbulence_closures.jl index 747f60171c..031d97af97 100644 --- a/test/test_turbulence_closures.jl +++ b/test/test_turbulence_closures.jl @@ -1,8 +1,9 @@ include("dependencies_for_runtests.jl") -using Oceananigans.TurbulenceClosures: CATKEVerticalDiffusivity, RiBasedVerticalDiffusivity +using Oceananigans.TurbulenceClosures: CATKEVerticalDiffusivity, RiBasedVerticalDiffusivity, DiscreteDiffusionFunction -using Oceananigans.TurbulenceClosures: viscosity_location, diffusivity_location +using Oceananigans.TurbulenceClosures: viscosity_location, diffusivity_location, + required_halo_size_x, required_halo_size_y, required_halo_size_z using Oceananigans.TurbulenceClosures: diffusive_flux_x, diffusive_flux_y, diffusive_flux_z, viscous_flux_ux, viscous_flux_uy, viscous_flux_uz, @@ -248,6 +249,21 @@ end @test closure.κ.T == T(κ) run_constant_isotropic_diffusivity_fluxdiv_tests(T) end + + @info " Testing ScalarDiffusivity with different halo requirements..." + closure = ScalarDiffusivity(ν=0.3) + @test required_halo_size_x(closure) == 1 + @test required_halo_size_y(closure) == 1 + @test required_halo_size_z(closure) == 1 + + @inline ν(i, j, k, grid, ℓx, ℓy, ℓz, clock, fields) = ℑxᶠᵃᵃ(i, j, k, grid, ℑxᶜᵃᵃ, fields.u) + closure = ScalarDiffusivity(; ν, discrete_form=true, required_halo_size=2) + + @test closure.ν isa DiscreteDiffusionFunction + @test required_halo_size_x(closure) == 2 + @test required_halo_size_y(closure) == 2 + @test required_halo_size_z(closure) == 2 + end @testset "HorizontalScalarDiffusivity" begin From b4f9e0fac9a61c9fdf3905ba6637f33965e3c77e Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Tue, 3 Sep 2024 22:20:45 -0400 Subject: [PATCH 36/68] test biharmonic diffusivity --- test/test_turbulence_closures.jl | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/test/test_turbulence_closures.jl b/test/test_turbulence_closures.jl index 031d97af97..c6336c49dc 100644 --- a/test/test_turbulence_closures.jl +++ b/test/test_turbulence_closures.jl @@ -256,6 +256,11 @@ end @test required_halo_size_y(closure) == 1 @test required_halo_size_z(closure) == 1 + closure = ScalarBiharmonicDiffusivity(ν=0.3) + @test required_halo_size_x(closure) == 2 + @test required_halo_size_y(closure) == 2 + @test required_halo_size_z(closure) == 2 + @inline ν(i, j, k, grid, ℓx, ℓy, ℓz, clock, fields) = ℑxᶠᵃᵃ(i, j, k, grid, ℑxᶜᵃᵃ, fields.u) closure = ScalarDiffusivity(; ν, discrete_form=true, required_halo_size=2) From 48112938fb2d8b923b940db384ca405d9f82073b Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Tue, 3 Sep 2024 22:48:05 -0400 Subject: [PATCH 37/68] fix enzyme tests --- test/test_enzyme.jl | 15 +++++++++------ 1 file changed, 9 insertions(+), 6 deletions(-) diff --git a/test/test_enzyme.jl b/test/test_enzyme.jl index ac457cae6e..ce56c0a9d3 100644 --- a/test/test_enzyme.jl +++ b/test/test_enzyme.jl @@ -28,10 +28,8 @@ const maximum_diffusivity = 100 Change diffusivity of model to `diffusivity`. """ function set_diffusivity!(model, diffusivity) - closure = VerticalScalarDiffusivity(; κ=diffusivity) - names = tuple(:c) # tracernames(model.tracers) - closure = with_tracers(names, closure) - model.closure = closure + closure = model.closure + fill!(closure.κ, diffusivity) return nothing end @@ -202,7 +200,10 @@ end topology = (Periodic, Periodic, Bounded) grid = RectilinearGrid(size=(Nx, Ny, Nz); x, y, z, topology) - diffusion = VerticalScalarDiffusivity(κ=0.1) + κ = ZFaceField(grid) + fill!(κ, 0.1) + + diffusion = VerticalScalarDiffusivity(; κ) u = XFaceField(grid) v = YFaceField(grid) @@ -248,7 +249,9 @@ end # Now for real amplitude = 1.0 - κ = 1.0 + κ = ZFaceField(grid) + fill!(κ, 1.0) + dmodel = Enzyme.make_zero(model) set_diffusivity!(dmodel, 0) From 2e02dfe5ce414ba3925599b0754c2d5c7f667d6e Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Tue, 3 Sep 2024 22:55:48 -0400 Subject: [PATCH 38/68] correct adaptbetter explanation --- src/Advection/adapt_advection_order.jl | 17 ++++++----------- 1 file changed, 6 insertions(+), 11 deletions(-) diff --git a/src/Advection/adapt_advection_order.jl b/src/Advection/adapt_advection_order.jl index 2f04b3135a..ecf11898a3 100644 --- a/src/Advection/adapt_advection_order.jl +++ b/src/Advection/adapt_advection_order.jl @@ -4,8 +4,9 @@ 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 `nothing`. -A high order advection sheme is reduced to a lower order advection scheme if the grid has fewer points in that 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. @@ -76,9 +77,7 @@ adapt_advection_order(advection::Nothing, N::Int, grid::AbstractGrid) = nothing ##### function adapt_advection_order(advection::Centered{H}, N::Int, grid::AbstractGrid) where H - if N == 1 && H != 1 - return Centered(; order = 2) - elseif N >= H + if N >= H return advection else return Centered(; order = N * 2) @@ -86,9 +85,7 @@ function adapt_advection_order(advection::Centered{H}, N::Int, grid::AbstractGri end function adapt_advection_order(advection::UpwindBiased{H}, N::Int, grid::AbstractGrid) where H - if N == 1 && H != 1 - return UpwindBiased(; order = 1) - elseif N >= H + if N >= H return advection else return UpwindBiased(; order = N * 2 - 1) @@ -106,9 +103,7 @@ new_weno_scheme(::WENO, grid, order, bounds, ::Type{Nothing}, ::Type{Nothing}, : new_weno_scheme(::WENO, grid, order, bounds, XT, YT, ZT) = WENO(grid; order, bounds) function adapt_advection_order(advection::WENO{H, FT, XT, YT, ZT}, N::Int, grid::AbstractGrid) where {H, FT, XT, YT, ZT} - if N == 1 && H != 1 - return UpwindBiased(; order = 1) - elseif N >= H + if N >= H return advection else return new_weno_scheme(advection, grid, N * 2 - 1, advection.bounds, XT, YT, ZT) From 4a45a958af09e0232613a1de16e5c845ebdf9385 Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Tue, 3 Sep 2024 23:00:06 -0400 Subject: [PATCH 39/68] add a test for adapt_advection --- test/test_nonhydrostatic_models.jl | 11 +++++++++++ 1 file changed, 11 insertions(+) diff --git a/test/test_nonhydrostatic_models.jl b/test/test_nonhydrostatic_models.jl index dedc95c9ce..78a754bd18 100644 --- a/test/test_nonhydrostatic_models.jl +++ b/test/test_nonhydrostatic_models.jl @@ -66,6 +66,17 @@ include("dependencies_for_runtests.jl") model = NonhydrostaticModel(closure=ScalarBiharmonicDiffusivity(), grid=funny_grid) @test model.grid.Hx == 2 && model.grid.Hy == 3 && model.grid.Hz == 4 + + + @info " Testing adjustment of advection schemes in NonhydrostaticModel constructor..." + small_grid = RectilinearGrid(size=(4, 2, 4), extent=(1, 2, 3), halo=(1, 1, 1)) + + # Model ensures that halos are at least of size 1 + model = NonhydrostaticModel(grid=small_grid, advection=WENO()) + @test model.advection isa FluxFormAdvection + @test required_halo_size(model.advection.x) == 2 + @test required_halo_size(model.advection.y) == 1 + @test required_halo_size(model.advection.z) == 2 end @testset "Model construction with single tracer and nothing tracer" begin From 2e33b5697e3de60eefa3c291ca0591bb57244570 Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Tue, 3 Sep 2024 23:00:28 -0400 Subject: [PATCH 40/68] whoops wrong test --- test/test_nonhydrostatic_models.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/test/test_nonhydrostatic_models.jl b/test/test_nonhydrostatic_models.jl index 78a754bd18..f72b23a337 100644 --- a/test/test_nonhydrostatic_models.jl +++ b/test/test_nonhydrostatic_models.jl @@ -74,9 +74,9 @@ include("dependencies_for_runtests.jl") # Model ensures that halos are at least of size 1 model = NonhydrostaticModel(grid=small_grid, advection=WENO()) @test model.advection isa FluxFormAdvection - @test required_halo_size(model.advection.x) == 2 + @test required_halo_size(model.advection.x) == 3 @test required_halo_size(model.advection.y) == 1 - @test required_halo_size(model.advection.z) == 2 + @test required_halo_size(model.advection.z) == 3 end @testset "Model construction with single tracer and nothing tracer" begin From ca876fab64a6e5d57f7d77cce516711ee8dbaf32 Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Tue, 3 Sep 2024 23:03:06 -0400 Subject: [PATCH 41/68] revert changes --- .../isopycnal_skew_symmetric_diffusivity.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/TurbulenceClosures/turbulence_closure_implementations/isopycnal_skew_symmetric_diffusivity.jl b/src/TurbulenceClosures/turbulence_closure_implementations/isopycnal_skew_symmetric_diffusivity.jl index 3b36933f81..78967ae51d 100644 --- a/src/TurbulenceClosures/turbulence_closure_implementations/isopycnal_skew_symmetric_diffusivity.jl +++ b/src/TurbulenceClosures/turbulence_closure_implementations/isopycnal_skew_symmetric_diffusivity.jl @@ -37,12 +37,12 @@ function IsopycnalSkewSymmetricDiffusivity(time_disc::TD = VerticallyImplicitTim κ_symmetric = 0, isopycnal_tensor = SmallSlopeIsopycnalTensor(), slope_limiter = FluxTapering(1e-2), - required_halo_size::Val{N} = Val(1)) where {TD, N} + required_halo_size::Int = 1) where {TD, N} isopycnal_tensor isa SmallSlopeIsopycnalTensor || error("Only isopycnal_tensor=SmallSlopeIsopycnalTensor() is currently supported.") - return IsopycnalSkewSymmetricDiffusivity{TD, N}(convert_diffusivity(FT, κ_skew), + return IsopycnalSkewSymmetricDiffusivity{TD, required_halo_size}(convert_diffusivity(FT, κ_skew), convert_diffusivity(FT, κ_symmetric), isopycnal_tensor, slope_limiter) From b3a94d8f7be3e82ea8ed3dbb343c786cc89ebfcc Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Tue, 3 Sep 2024 23:08:51 -0400 Subject: [PATCH 42/68] add couple more tests --- test/test_nonhydrostatic_models.jl | 14 ++++++++++++++ 1 file changed, 14 insertions(+) diff --git a/test/test_nonhydrostatic_models.jl b/test/test_nonhydrostatic_models.jl index f72b23a337..523f11327e 100644 --- a/test/test_nonhydrostatic_models.jl +++ b/test/test_nonhydrostatic_models.jl @@ -77,6 +77,20 @@ include("dependencies_for_runtests.jl") @test required_halo_size(model.advection.x) == 3 @test required_halo_size(model.advection.y) == 1 @test required_halo_size(model.advection.z) == 3 + + + model = NonhydrostaticModel(grid=small_grid, advection=UpwindBiased(; order = 9)) + @test model.advection isa FluxFormAdvection + @test required_halo_size(model.advection.x) == 4 + @test required_halo_size(model.advection.y) == 1 + @test required_halo_size(model.advection.z) == 4 + + + model = NonhydrostaticModel(grid=small_grid, advection=Centered(; order = 11)) + @test model.advection isa FluxFormAdvection + @test required_halo_size(model.advection.x) == 4 + @test required_halo_size(model.advection.y) == 1 + @test required_halo_size(model.advection.z) == 4 end @testset "Model construction with single tracer and nothing tracer" begin From e1085a48a70daa8583ea28449041f2f476496ce4 Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Tue, 3 Sep 2024 23:10:20 -0400 Subject: [PATCH 43/68] fix tests --- test/test_nonhydrostatic_models.jl | 20 +++++++++++--------- 1 file changed, 11 insertions(+), 9 deletions(-) diff --git a/test/test_nonhydrostatic_models.jl b/test/test_nonhydrostatic_models.jl index 523f11327e..a30d74f2ba 100644 --- a/test/test_nonhydrostatic_models.jl +++ b/test/test_nonhydrostatic_models.jl @@ -1,5 +1,7 @@ include("dependencies_for_runtests.jl") +using Oceananigans.Grids: required_halo_size_x, required_halo_size_y, required_halo_size_z + @testset "Models" begin @info "Testing models..." @@ -74,23 +76,23 @@ include("dependencies_for_runtests.jl") # Model ensures that halos are at least of size 1 model = NonhydrostaticModel(grid=small_grid, advection=WENO()) @test model.advection isa FluxFormAdvection - @test required_halo_size(model.advection.x) == 3 - @test required_halo_size(model.advection.y) == 1 - @test required_halo_size(model.advection.z) == 3 + @test required_halo_size_x(model.advection) == 3 + @test required_halo_size_y(model.advection) == 1 + @test required_halo_size_z(model.advection) == 3 model = NonhydrostaticModel(grid=small_grid, advection=UpwindBiased(; order = 9)) @test model.advection isa FluxFormAdvection - @test required_halo_size(model.advection.x) == 4 - @test required_halo_size(model.advection.y) == 1 - @test required_halo_size(model.advection.z) == 4 + @test required_halo_size_x(model.advection) == 4 + @test required_halo_size_y(model.advection) == 1 + @test required_halo_size_z(model.advection) == 4 model = NonhydrostaticModel(grid=small_grid, advection=Centered(; order = 11)) @test model.advection isa FluxFormAdvection - @test required_halo_size(model.advection.x) == 4 - @test required_halo_size(model.advection.y) == 1 - @test required_halo_size(model.advection.z) == 4 + @test required_halo_size_x(model.advection) == 4 + @test required_halo_size_y(model.advection) == 1 + @test required_halo_size_z(model.advection) == 4 end @testset "Model construction with single tracer and nothing tracer" begin From 4a781689b626bcff8512cb1c51e4fd2bb9c61d68 Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Tue, 3 Sep 2024 23:10:51 -0400 Subject: [PATCH 44/68] without the extra space --- test/test_nonhydrostatic_models.jl | 2 -- 1 file changed, 2 deletions(-) diff --git a/test/test_nonhydrostatic_models.jl b/test/test_nonhydrostatic_models.jl index a30d74f2ba..a24de43480 100644 --- a/test/test_nonhydrostatic_models.jl +++ b/test/test_nonhydrostatic_models.jl @@ -80,14 +80,12 @@ using Oceananigans.Grids: required_halo_size_x, required_halo_size_y, required_h @test required_halo_size_y(model.advection) == 1 @test required_halo_size_z(model.advection) == 3 - model = NonhydrostaticModel(grid=small_grid, advection=UpwindBiased(; order = 9)) @test model.advection isa FluxFormAdvection @test required_halo_size_x(model.advection) == 4 @test required_halo_size_y(model.advection) == 1 @test required_halo_size_z(model.advection) == 4 - model = NonhydrostaticModel(grid=small_grid, advection=Centered(; order = 11)) @test model.advection isa FluxFormAdvection @test required_halo_size_x(model.advection) == 4 From f5b778928f8e5c6593400a46b6ceb6ef1e91adb6 Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Tue, 3 Sep 2024 23:43:01 -0400 Subject: [PATCH 45/68] fix enzyme test --- .../isopycnal_skew_symmetric_diffusivity.jl | 2 +- test/test_enzyme.jl | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/src/TurbulenceClosures/turbulence_closure_implementations/isopycnal_skew_symmetric_diffusivity.jl b/src/TurbulenceClosures/turbulence_closure_implementations/isopycnal_skew_symmetric_diffusivity.jl index 78967ae51d..87f6ceff5f 100644 --- a/src/TurbulenceClosures/turbulence_closure_implementations/isopycnal_skew_symmetric_diffusivity.jl +++ b/src/TurbulenceClosures/turbulence_closure_implementations/isopycnal_skew_symmetric_diffusivity.jl @@ -37,7 +37,7 @@ function IsopycnalSkewSymmetricDiffusivity(time_disc::TD = VerticallyImplicitTim κ_symmetric = 0, isopycnal_tensor = SmallSlopeIsopycnalTensor(), slope_limiter = FluxTapering(1e-2), - required_halo_size::Int = 1) where {TD, N} + required_halo_size::Int = 1) where TD isopycnal_tensor isa SmallSlopeIsopycnalTensor || error("Only isopycnal_tensor=SmallSlopeIsopycnalTensor() is currently supported.") diff --git a/test/test_enzyme.jl b/test/test_enzyme.jl index ce56c0a9d3..be4aa983d0 100644 --- a/test/test_enzyme.jl +++ b/test/test_enzyme.jl @@ -29,7 +29,7 @@ Change diffusivity of model to `diffusivity`. """ function set_diffusivity!(model, diffusivity) closure = model.closure - fill!(closure.κ, diffusivity) + fill!(closure.κ.c, diffusivity) return nothing end From 98abf58ad30a9d93032d6f5f9d87ce219b589653 Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Wed, 4 Sep 2024 06:36:37 -0400 Subject: [PATCH 46/68] fix enzyme tests --- test/test_enzyme.jl | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/test/test_enzyme.jl b/test/test_enzyme.jl index be4aa983d0..404960e6d0 100644 --- a/test/test_enzyme.jl +++ b/test/test_enzyme.jl @@ -249,8 +249,9 @@ end # Now for real amplitude = 1.0 + diffusivity = 1.0 κ = ZFaceField(grid) - fill!(κ, 1.0) + fill!(κ, diffusivity) dmodel = Enzyme.make_zero(model) set_diffusivity!(dmodel, 0) @@ -259,7 +260,7 @@ end stable_diffusion!, Duplicated(model, dmodel), Const(amplitude), - Active(κ)) + Active(diffusivity)) @info """ \n Enzyme computed $dc²_dκ From 9f35ebec11c8e4bd0391d8508a9602ce4d39b9dd Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Wed, 4 Sep 2024 06:38:03 -0400 Subject: [PATCH 47/68] fix nonhydrostatic tests --- test/test_nonhydrostatic_models.jl | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/test/test_nonhydrostatic_models.jl b/test/test_nonhydrostatic_models.jl index a24de43480..51c9f23502 100644 --- a/test/test_nonhydrostatic_models.jl +++ b/test/test_nonhydrostatic_models.jl @@ -77,19 +77,19 @@ using Oceananigans.Grids: required_halo_size_x, required_halo_size_y, required_h model = NonhydrostaticModel(grid=small_grid, advection=WENO()) @test model.advection isa FluxFormAdvection @test required_halo_size_x(model.advection) == 3 - @test required_halo_size_y(model.advection) == 1 + @test required_halo_size_y(model.advection) == 2 @test required_halo_size_z(model.advection) == 3 model = NonhydrostaticModel(grid=small_grid, advection=UpwindBiased(; order = 9)) @test model.advection isa FluxFormAdvection @test required_halo_size_x(model.advection) == 4 - @test required_halo_size_y(model.advection) == 1 + @test required_halo_size_y(model.advection) == 2 @test required_halo_size_z(model.advection) == 4 - model = NonhydrostaticModel(grid=small_grid, advection=Centered(; order = 11)) + model = NonhydrostaticModel(grid=small_grid, advection=Centered(; order = 10)) @test model.advection isa FluxFormAdvection @test required_halo_size_x(model.advection) == 4 - @test required_halo_size_y(model.advection) == 1 + @test required_halo_size_y(model.advection) == 2 @test required_halo_size_z(model.advection) == 4 end From cc98d5e3fe34f85f838b0b1215861594c1c22d4f Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Wed, 4 Sep 2024 06:38:18 -0400 Subject: [PATCH 48/68] rmove one newline --- test/test_nonhydrostatic_models.jl | 1 - 1 file changed, 1 deletion(-) diff --git a/test/test_nonhydrostatic_models.jl b/test/test_nonhydrostatic_models.jl index 51c9f23502..bd10653bfd 100644 --- a/test/test_nonhydrostatic_models.jl +++ b/test/test_nonhydrostatic_models.jl @@ -69,7 +69,6 @@ using Oceananigans.Grids: required_halo_size_x, required_halo_size_y, required_h model = NonhydrostaticModel(closure=ScalarBiharmonicDiffusivity(), grid=funny_grid) @test model.grid.Hx == 2 && model.grid.Hy == 3 && model.grid.Hz == 4 - @info " Testing adjustment of advection schemes in NonhydrostaticModel constructor..." small_grid = RectilinearGrid(size=(4, 2, 4), extent=(1, 2, 3), halo=(1, 1, 1)) From ab0827ea690f554f3b3bf4836acbb540146f8f20 Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Wed, 4 Sep 2024 07:42:13 -0400 Subject: [PATCH 49/68] simplify test --- test/test_enzyme.jl | 2 -- 1 file changed, 2 deletions(-) diff --git a/test/test_enzyme.jl b/test/test_enzyme.jl index 404960e6d0..b6ad1cad9d 100644 --- a/test/test_enzyme.jl +++ b/test/test_enzyme.jl @@ -250,8 +250,6 @@ end # Now for real amplitude = 1.0 diffusivity = 1.0 - κ = ZFaceField(grid) - fill!(κ, diffusivity) dmodel = Enzyme.make_zero(model) set_diffusivity!(dmodel, 0) From 9c2103a287ffb0ed0ebaca048be8e963075ef4b7 Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Fri, 6 Sep 2024 07:48:23 -0400 Subject: [PATCH 50/68] Update src/Models/HydrostaticFreeSurfaceModels/hydrostatic_free_surface_model.jl Co-authored-by: Gregory L. Wagner --- .../hydrostatic_free_surface_model.jl | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/src/Models/HydrostaticFreeSurfaceModels/hydrostatic_free_surface_model.jl b/src/Models/HydrostaticFreeSurfaceModels/hydrostatic_free_surface_model.jl index 9094f912d1..6b5a03fd00 100644 --- a/src/Models/HydrostaticFreeSurfaceModels/hydrostatic_free_surface_model.jl +++ b/src/Models/HydrostaticFreeSurfaceModels/hydrostatic_free_surface_model.jl @@ -130,9 +130,7 @@ function HydrostaticFreeSurfaceModel(; grid, tracers = tupleit(tracers) # supports tracers=:c keyword argument (for example) - # Adjust advection scheme to be valid on a particular grid size. i.e. if the grid size - # is smaller than the advection order, reduce the order of the advection in that particular - # direction + # Reduce the advection order in directions that do not have enough grid points momentum_advection = adapt_advection_order(momentum_advection, grid) tracer_advection = adapt_advection_order(tracer_advection, grid) From 2da6007ec425715f2841534ff42c2dc2c0677849 Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Fri, 6 Sep 2024 08:17:37 -0400 Subject: [PATCH 51/68] remove extra dispatches --- src/Advection/momentum_advection_operators.jl | 93 +------------------ src/Advection/tracer_advection_operators.jl | 33 +------ src/ImmersedBoundaries/conditional_fluxes.jl | 40 ++++---- 3 files changed, 24 insertions(+), 142 deletions(-) diff --git a/src/Advection/momentum_advection_operators.jl b/src/Advection/momentum_advection_operators.jl index e0548b4443..eda8f78423 100644 --- a/src/Advection/momentum_advection_operators.jl +++ b/src/Advection/momentum_advection_operators.jl @@ -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) @@ -105,81 +97,4 @@ end @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) - -# Fallback for `nothing` advection and `ZeroField` tracers and velocities -@inline _advective_momentum_flux_Uu(i, j, k, grid, ::Nothing, ::ZeroField, ::ZeroField) = zero(grid) -@inline _advective_momentum_flux_Uv(i, j, k, grid, ::Nothing, ::ZeroField, ::ZeroField) = zero(grid) -@inline _advective_momentum_flux_Uw(i, j, k, grid, ::Nothing, ::ZeroField, ::ZeroField) = zero(grid) - -@inline _advective_momentum_flux_Vu(i, j, k, grid, ::Nothing, ::ZeroField, ::ZeroField) = zero(grid) -@inline _advective_momentum_flux_Vv(i, j, k, grid, ::Nothing, ::ZeroField, ::ZeroField) = zero(grid) -@inline _advective_momentum_flux_Vw(i, j, k, grid, ::Nothing, ::ZeroField, ::ZeroField) = zero(grid) - -@inline _advective_momentum_flux_Wu(i, j, k, grid, ::Nothing, ::ZeroField, ::ZeroField) = zero(grid) -@inline _advective_momentum_flux_Wv(i, j, k, grid, ::Nothing, ::ZeroField, ::ZeroField) = zero(grid) -@inline _advective_momentum_flux_Ww(i, j, k, grid, ::Nothing, ::ZeroField, ::ZeroField) = zero(grid) - -@inline _advective_momentum_flux_Uu(i, j, k, grid, ::Nothing, U, ::ZeroField) = zero(grid) -@inline _advective_momentum_flux_Uv(i, j, k, grid, ::Nothing, U, ::ZeroField) = zero(grid) -@inline _advective_momentum_flux_Uw(i, j, k, grid, ::Nothing, U, ::ZeroField) = zero(grid) -@inline _advective_momentum_flux_Uu(i, j, k, grid, ::Nothing, ::ZeroField, u) = zero(grid) -@inline _advective_momentum_flux_Uv(i, j, k, grid, ::Nothing, ::ZeroField, v) = zero(grid) -@inline _advective_momentum_flux_Uw(i, j, k, grid, ::Nothing, ::ZeroField, w) = zero(grid) - -@inline _advective_momentum_flux_Vu(i, j, k, grid, ::Nothing, V, ::ZeroField) = zero(grid) -@inline _advective_momentum_flux_Vv(i, j, k, grid, ::Nothing, V, ::ZeroField) = zero(grid) -@inline _advective_momentum_flux_Vw(i, j, k, grid, ::Nothing, V, ::ZeroField) = zero(grid) -@inline _advective_momentum_flux_Vu(i, j, k, grid, ::Nothing, ::ZeroField, u) = zero(grid) -@inline _advective_momentum_flux_Vv(i, j, k, grid, ::Nothing, ::ZeroField, v) = zero(grid) -@inline _advective_momentum_flux_Vw(i, j, k, grid, ::Nothing, ::ZeroField, w) = zero(grid) - -@inline _advective_momentum_flux_Wu(i, j, k, grid, ::Nothing, W, ::ZeroField) = zero(grid) -@inline _advective_momentum_flux_Wv(i, j, k, grid, ::Nothing, W, ::ZeroField) = zero(grid) -@inline _advective_momentum_flux_Ww(i, j, k, grid, ::Nothing, W, ::ZeroField) = zero(grid) -@inline _advective_momentum_flux_Wu(i, j, k, grid, ::Nothing, ::ZeroField, u) = zero(grid) -@inline _advective_momentum_flux_Wv(i, j, k, grid, ::Nothing, ::ZeroField, v) = zero(grid) -@inline _advective_momentum_flux_Ww(i, j, k, grid, ::Nothing, ::ZeroField, w) = zero(grid) - -for Scheme in (:UpwindBiased, :Centered, :WENO, :FluxFormAdvection) - @eval begin - # Fallback for `ZeroField` tracers and velocities - @inline _advective_momentum_flux_Uu(i, j, k, grid, ::$Scheme, ::ZeroField, ::ZeroField) = zero(grid) - @inline _advective_momentum_flux_Uv(i, j, k, grid, ::$Scheme, ::ZeroField, ::ZeroField) = zero(grid) - @inline _advective_momentum_flux_Uw(i, j, k, grid, ::$Scheme, ::ZeroField, ::ZeroField) = zero(grid) - - @inline _advective_momentum_flux_Vu(i, j, k, grid, ::$Scheme, ::ZeroField, ::ZeroField) = zero(grid) - @inline _advective_momentum_flux_Vv(i, j, k, grid, ::$Scheme, ::ZeroField, ::ZeroField) = zero(grid) - @inline _advective_momentum_flux_Vw(i, j, k, grid, ::$Scheme, ::ZeroField, ::ZeroField) = zero(grid) - - @inline _advective_momentum_flux_Wu(i, j, k, grid, ::$Scheme, ::ZeroField, ::ZeroField) = zero(grid) - @inline _advective_momentum_flux_Wv(i, j, k, grid, ::$Scheme, ::ZeroField, ::ZeroField) = zero(grid) - @inline _advective_momentum_flux_Ww(i, j, k, grid, ::$Scheme, ::ZeroField, ::ZeroField) = zero(grid) - - # Fallback for `ZeroField` tracers - @inline _advective_momentum_flux_Uu(i, j, k, grid, ::$Scheme, U, ::ZeroField) = zero(grid) - @inline _advective_momentum_flux_Uv(i, j, k, grid, ::$Scheme, U, ::ZeroField) = zero(grid) - @inline _advective_momentum_flux_Uw(i, j, k, grid, ::$Scheme, U, ::ZeroField) = zero(grid) - - @inline _advective_momentum_flux_Vu(i, j, k, grid, ::$Scheme, V, ::ZeroField) = zero(grid) - @inline _advective_momentum_flux_Vv(i, j, k, grid, ::$Scheme, V, ::ZeroField) = zero(grid) - @inline _advective_momentum_flux_Vw(i, j, k, grid, ::$Scheme, V, ::ZeroField) = zero(grid) - - @inline _advective_momentum_flux_Wu(i, j, k, grid, ::$Scheme, W, ::ZeroField) = zero(grid) - @inline _advective_momentum_flux_Wv(i, j, k, grid, ::$Scheme, W, ::ZeroField) = zero(grid) - @inline _advective_momentum_flux_Ww(i, j, k, grid, ::$Scheme, W, ::ZeroField) = zero(grid) - - # Fallback for `ZeroField` velocities - @inline _advective_momentum_flux_Uu(i, j, k, grid, ::$Scheme, ::ZeroField, u) = zero(grid) - @inline _advective_momentum_flux_Uv(i, j, k, grid, ::$Scheme, ::ZeroField, v) = zero(grid) - @inline _advective_momentum_flux_Uw(i, j, k, grid, ::$Scheme, ::ZeroField, w) = zero(grid) - - @inline _advective_momentum_flux_Vu(i, j, k, grid, ::$Scheme, ::ZeroField, u) = zero(grid) - @inline _advective_momentum_flux_Vv(i, j, k, grid, ::$Scheme, ::ZeroField, v) = zero(grid) - @inline _advective_momentum_flux_Vw(i, j, k, grid, ::$Scheme, ::ZeroField, w) = zero(grid) - - @inline _advective_momentum_flux_Wu(i, j, k, grid, ::$Scheme, ::ZeroField, u) = zero(grid) - @inline _advective_momentum_flux_Wv(i, j, k, grid, ::$Scheme, ::ZeroField, v) = zero(grid) - @inline _advective_momentum_flux_Ww(i, j, k, grid, ::$Scheme, ::ZeroField, w) = zero(grid) - end -end \ No newline at end of file +@inline _advective_momentum_flux_Ww(i, j, k, grid, ::Nothing, args...) = zero(grid) \ No newline at end of file diff --git a/src/Advection/tracer_advection_operators.jl b/src/Advection/tracer_advection_operators.jl index c104d43ca9..c22dee4cd9 100644 --- a/src/Advection/tracer_advection_operators.jl +++ b/src/Advection/tracer_advection_operators.jl @@ -12,37 +12,6 @@ @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) - -for Scheme in (:UpwindBiased, :Centered, :WENO, :FluxFormAdvection) - @eval begin - # 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) - end -end - ##### ##### Tracer advection operator ##### @@ -67,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) diff --git a/src/ImmersedBoundaries/conditional_fluxes.jl b/src/ImmersedBoundaries/conditional_fluxes.jl index 03be5fb0f9..fd2d45fb8b 100644 --- a/src/ImmersedBoundaries/conditional_fluxes.jl +++ b/src/ImmersedBoundaries/conditional_fluxes.jl @@ -124,31 +124,27 @@ end # Fallback for `nothing` Advection -for flux_dir in (:x, :y, :z) - advective_tracer_flux = Symbol(:_advective_tracer_flux_, flux_dir) +# dx(uu), dy(vu), dz(wu) +# ccc, ffc, fcf +@inline _advective_momentum_flux_Uu(i, j, k, ibg::IBG, ::Nothing, args...) = zero(ibg) +@inline _advective_momentum_flux_Vu(i, j, k, ibg::IBG, ::Nothing, args...) = zero(ibg) +@inline _advective_momentum_flux_Wu(i, j, k, ibg::IBG, ::Nothing, args...) = zero(ibg) - @eval begin - @inline $advective_tracer_flux(i, j, k, ibg::IBG, ::Nothing, args...) = zero(ibg) - @inline $advective_tracer_flux(i, j, k, ibg::IBG, ::Nothing, ::ZeroField, ::ZeroField) = zero(ibg) - @inline $advective_tracer_flux(i, j, k, ibg::IBG, ::Nothing, U, ::ZeroField) = zero(ibg) - @inline $advective_tracer_flux(i, j, k, ibg::IBG, scheme, ::ZeroField, ::ZeroField) = zero(ibg) - @inline $advective_tracer_flux(i, j, k, ibg::IBG, scheme, U, ::ZeroField) = zero(ibg) - @inline $advective_tracer_flux(i, j, k, ibg::IBG, scheme, ::ZeroField, c) = zero(ibg) - end -end +# dx(uv), dy(vv), dz(wv) +# ffc, ccc, cff +@inline _advective_momentum_flux_Uv(i, j, k, ibg::IBG, ::Nothing, args...) = zero(ibg) +@inline _advective_momentum_flux_Vv(i, j, k, ibg::IBG, ::Nothing, args...) = zero(ibg) +@inline _advective_momentum_flux_Wv(i, j, k, ibg::IBG, ::Nothing, args...) = zero(ibg) -for flux_dir in (:Uu, :Vu, :Wu, :Uv, :Vv, :Wv, :Uw, :Vw, :Ww) - advective_momentum_flux = Symbol(:_advective_momentum_flux_, flux_dir) +# dx(uw), dy(vw), dz(ww) +# fcf, cff, ccc +@inline _advective_momentum_flux_Uw(i, j, k, ibg::IBG, ::Nothing, args...) = zero(ibg) +@inline _advective_momentum_flux_Vw(i, j, k, ibg::IBG, ::Nothing, args...) = zero(ibg) +@inline _advective_momentum_flux_Ww(i, j, k, ibg::IBG, ::Nothing, args...) = zero(ibg) - @eval begin - @inline $advective_momentum_flux(i, j, k, ibg::IBG, ::Nothing, args...) = zero(ibg) - @inline $advective_momentum_flux(i, j, k, ibg::IBG, ::Nothing, ::ZeroField, ::ZeroField) = zero(ibg) - @inline $advective_momentum_flux(i, j, k, ibg::IBG, ::Nothing, U, ::ZeroField) = zero(ibg) - @inline $advective_momentum_flux(i, j, k, ibg::IBG, scheme, ::ZeroField, ::ZeroField) = zero(ibg) - @inline $advective_momentum_flux(i, j, k, ibg::IBG, scheme, U, ::ZeroField) = zero(ibg) - @inline $advective_momentum_flux(i, j, k, ibg::IBG, scheme, ::ZeroField, u) = zero(ibg) - end -end +@inline _advective_tracer_flux_x(i, j, k, ibg::IBG, ::Nothing, args...) = zero(ibg) +@inline _advective_tracer_flux_y(i, j, k, ibg::IBG, ::Nothing, args...) = zero(ibg) +@inline _advective_tracer_flux_z(i, j, k, ibg::IBG, ::Nothing, args...) = zero(ibg) ##### ##### "Boundary-aware" reconstruct From da88f62df3edfffef5b4c89e4e0ec419ef6b23f5 Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Mon, 16 Sep 2024 19:11:21 -0600 Subject: [PATCH 52/68] relaunch the tests --- .../compute_hydrostatic_free_surface_boundary_tendencies.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/Models/HydrostaticFreeSurfaceModels/compute_hydrostatic_free_surface_boundary_tendencies.jl b/src/Models/HydrostaticFreeSurfaceModels/compute_hydrostatic_free_surface_boundary_tendencies.jl index 4096cf0972..db8ba9839f 100644 --- a/src/Models/HydrostaticFreeSurfaceModels/compute_hydrostatic_free_surface_boundary_tendencies.jl +++ b/src/Models/HydrostaticFreeSurfaceModels/compute_hydrostatic_free_surface_boundary_tendencies.jl @@ -17,7 +17,7 @@ function compute_boundary_tendencies!(model::HydrostaticFreeSurfaceModel) p_parameters = boundary_p_kernel_parameters(grid, arch) κ_parameters = boundary_κ_kernel_parameters(grid, model.closure, arch) - # We need new values for `w`, `p` and `κ` + # Compute new values for `w`, `p` and `κ` on the perifery compute_auxiliaries!(model; w_parameters, p_parameters, κ_parameters) # parameters for communicating North / South / East / West side From 0d2e7bb06d5f9dc7c6e5b7fe5bd4c726f9817280 Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Tue, 8 Oct 2024 09:45:21 +0200 Subject: [PATCH 53/68] change comment --- src/Advection/adapt_advection_order.jl | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/src/Advection/adapt_advection_order.jl b/src/Advection/adapt_advection_order.jl index ecf11898a3..f27d428ec7 100644 --- a/src/Advection/adapt_advection_order.jl +++ b/src/Advection/adapt_advection_order.jl @@ -29,13 +29,13 @@ function adapt_advection_order(advection, grid::AbstractGrid) changed_advection = any((changed_x, changed_y, changed_z)) if changed_x - @info "User-defined advection scheme $(summary(advection)) reduced to $(summary(new_advection.x)) in the x-direction to comply with grid-size limitations." + @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 "User-defined advection scheme $(summary(advection)) reduced to $(summary(new_advection.y)) in the y-direction to comply with grid-size limitations." + @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 "User-defined advection scheme $(summary(advection)) reduced to $(summary(new_advection.z)) in the z-direction to comply with grid-size limitations." + @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) @@ -55,13 +55,13 @@ function adapt_advection_order(advection::FluxFormAdvection, grid::AbstractGrid) changed_advection = any((changed_x, changed_y, changed_z)) if changed_x - @info "User-defined advection scheme $(summary(advection)) reduced to $(summary(new_advection.x)) in the x-direction to comply with grid-size limitations." + @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 "User-defined advection scheme $(summary(advection)) reduced to $(summary(new_advection.y)) in the y-direction to comply with grid-size limitations." + @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 "User-defined advection scheme $(summary(advection)) reduced to $(summary(new_advection.z)) in the z-direction to comply with grid-size limitations." + @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) From b6d65e273b5f289bb1ded84693f2699f8babd49a Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Tue, 8 Oct 2024 09:46:51 +0200 Subject: [PATCH 54/68] switch to B --- src/Advection/adapt_advection_order.jl | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/src/Advection/adapt_advection_order.jl b/src/Advection/adapt_advection_order.jl index f27d428ec7..288723ffa4 100644 --- a/src/Advection/adapt_advection_order.jl +++ b/src/Advection/adapt_advection_order.jl @@ -76,16 +76,16 @@ adapt_advection_order(advection::Nothing, N::Int, grid::AbstractGrid) = nothing ##### Directional adapt advection order ##### -function adapt_advection_order(advection::Centered{H}, N::Int, grid::AbstractGrid) where H - if N >= H +function adapt_advection_order(advection::Centered{B}, N::Int, grid::AbstractGrid) where B + if N >= B return advection else return Centered(; order = N * 2) end end -function adapt_advection_order(advection::UpwindBiased{H}, N::Int, grid::AbstractGrid) where H - if N >= H +function adapt_advection_order(advection::UpwindBiased{B}, N::Int, grid::AbstractGrid) where B + if N >= B return advection else return UpwindBiased(; order = N * 2 - 1) @@ -102,8 +102,8 @@ we rebuild the advection without passing the grid information, otherwise we use 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{H, FT, XT, YT, ZT}, N::Int, grid::AbstractGrid) where {H, FT, XT, YT, ZT} - if N >= H +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, N * 2 - 1, advection.bounds, XT, YT, ZT) From 147988464ff644785aac85b761a6dbdb594d1935 Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Tue, 8 Oct 2024 10:00:04 +0200 Subject: [PATCH 55/68] new syntax for KernelParameters --- ...static_free_surface_boundary_tendencies.jl | 18 +++--- ...pute_nonhydrostatic_boundary_tendencies.jl | 61 ++++++++----------- 2 files changed, 32 insertions(+), 47 deletions(-) diff --git a/src/Models/HydrostaticFreeSurfaceModels/compute_hydrostatic_free_surface_boundary_tendencies.jl b/src/Models/HydrostaticFreeSurfaceModels/compute_hydrostatic_free_surface_boundary_tendencies.jl index 2ea6abc0a4..271e5562b6 100644 --- a/src/Models/HydrostaticFreeSurfaceModels/compute_hydrostatic_free_surface_boundary_tendencies.jl +++ b/src/Models/HydrostaticFreeSurfaceModels/compute_hydrostatic_free_surface_boundary_tendencies.jl @@ -56,19 +56,15 @@ function boundary_w_kernel_parameters(grid, arch) Nx, Ny, _ = size(grid) Hx, Hy, _ = halo_size(grid) - Sx = (Hx, Ny+2) - Sy = (Nx+2, Hy) - # Offsets in tangential direction are == -1 to # cover the required corners - Oxᴸ = (-Hx+1, -1) - Oyᴸ = (-1, -Hy+1) - Oxᴿ = (Nx-1, -1) - Oyᴿ = (-1, Ny-1) + param_west = (-Hx+1:0, 0:Ny+1, 1:Nz) + param_east = (Nx+1:Nx+Hx-1, 0:Ny+1, 1:Nz) + param_north = (0:Nx+1, -Hy+1:0, 1:Nz) + param_south = (0:Nx+1, Ny+1:Ny+Hy-1, 1:Nz) - sizes = (Sx, Sy, Sx, Sy) - offs = (Oxᴸ, Oyᴸ, Oxᴿ, Oyᴿ) - - return boundary_parameters(sizes, offs, grid, arch) + params = (param_west, param_east, param_north, param_south) + + return boundary_parameters(params, grid, arch) end diff --git a/src/Models/NonhydrostaticModels/compute_nonhydrostatic_boundary_tendencies.jl b/src/Models/NonhydrostaticModels/compute_nonhydrostatic_boundary_tendencies.jl index f6a5cbf1a9..80294d828c 100644 --- a/src/Models/NonhydrostaticModels/compute_nonhydrostatic_boundary_tendencies.jl +++ b/src/Models/NonhydrostaticModels/compute_nonhydrostatic_boundary_tendencies.jl @@ -28,36 +28,29 @@ end function boundary_tendency_kernel_parameters(grid, arch) Nx, Ny, Nz = size(grid) Hx, Hy, _ = halo_size(grid) - - Sx = (Hx, Ny, Nz) - Sy = (Nx, Hy, Nz) - Oᴸ = (0, 0, 0) - Oxᴿ = (Nx-Hx, 0, 0) - Oyᴿ = (0, Ny-Hy, 0) + param_west = (1:Hx, 1:Ny, 1:Nz) + param_east = (Nx-Hx+1:Nx, 1:Ny, 1:Nz) + param_north = (1:Nx, 1:Hy, 1:Nz) + param_south = (1:Nx, Ny-Hy+1:Ny, 1:Nz) - sizes = (Sx, Sy, Sx, Sy) - offs = (Oᴸ, Oᴸ, Oxᴿ, Oyᴿ) + params = (param_west, param_east, param_north, param_south) - return boundary_parameters(sizes, offs, grid, arch) + return boundary_parameters(params, grid, arch) end # p needs computing in the range 0 : 0 and N + 1 : N + 1 function boundary_p_kernel_parameters(grid, arch) Nx, Ny, _ = size(grid) - Sx = (1, Ny) - Sy = (Nx, 1) + param_west = (0:0, 1:Ny, 1:Nz) + param_east = (Nx+1:Nx+1, 1:Ny, 1:Nz) + param_north = (1:Nx, 0:0, 1:Nz) + param_south = (1:Nx, Ny+1:Ny+1, 1:Nz) - Oxᴸ = (-1, 0) - Oyᴸ = (0, -1) - Oxᴿ = (Nx, 0) - Oyᴿ = (0, Ny) + params = (param_west, param_east, param_north, param_south) - sizes = (Sx, Sy, Sx, Sy) - offs = (Oxᴸ, Oyᴸ, Oxᴿ, Oyᴿ) - - return boundary_parameters(sizes, offs, grid, arch) + return boundary_parameters(params, grid, arch) end # diffusivities need recomputing in the range 0 : B and N - B + 1 : N + 1 @@ -67,32 +60,28 @@ function boundary_κ_kernel_parameters(grid, closure, arch) Bx = required_halo_size_x(closure) By = required_halo_size_y(closure) - Sx = (Bx+1, Ny, Nz) - Sy = (Nx, By+1, Nz) - - Oxᴸ = (-1, 0, 0) - Oyᴸ = (0, -1, 0) - Oxᴿ = (Nx-Bx, 0, 0) - Oyᴿ = (0, Ny-By, 0) + param_west = (0:Bx, 1:Ny, 1:Nz) + param_east = (Nx-Bx+1:Nx+1, 1:Ny, 1:Nz) + param_north = (1:Nx, 0:By, 1:Nz) + param_south = (1:Nx, Ny-By+1:Ny+1, 1:Nz) - sizes = (Sx, Sy, Sx, Sy) - offs = (Oxᴸ, Oyᴸ, Oxᴿ, Oyᴿ) + params = (param_west, param_east, param_north, param_south) - return boundary_parameters(sizes, offs, grid, arch) + return boundary_parameters(params, grid, arch) end # Recompute only on communicating sides -function boundary_parameters(S, O, grid, arch) +function boundary_parameters(parameters, grid, arch) Rx, Ry, _ = arch.ranks Tx, Ty, _ = topology(grid) - include_xᴸ = !isa(grid, XFlatGrid) && (Rx != 1) && !(Tx == RightConnected) - include_yᴸ = !isa(grid, YFlatGrid) && (Ry != 1) && !(Ty == RightConnected) - include_xᴿ = !isa(grid, XFlatGrid) && (Rx != 1) && !(Tx == LeftConnected) - include_yᴿ = !isa(grid, YFlatGrid) && (Ry != 1) && !(Ty == LeftConnected) + include_west = !isa(grid, XFlatGrid) && (Rx != 1) && !(Tx == RightConnected) + include_east = !isa(grid, XFlatGrid) && (Rx != 1) && !(Tx == LeftConnected) + include_south = !isa(grid, YFlatGrid) && (Ry != 1) && !(Ty == RightConnected) + include_north = !isa(grid, YFlatGrid) && (Ry != 1) && !(Ty == LeftConnected) - include_side = (include_xᴸ, include_yᴸ, include_xᴿ, include_yᴿ) + include_side = (include_west, include_east, include_south, include_north) - return Tuple(KernelParameters(S[i], O[i]) for i in findall(include_side)) + return Tuple(KernelParameters(parameters[i]) for i in findall(include_side)) end From e741ce89ed24a01acbe96c98190010dc778d38ec Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Tue, 8 Oct 2024 10:56:24 +0200 Subject: [PATCH 56/68] bugfix --- ...static_free_surface_boundary_tendencies.jl | 10 ++++----- ...pute_nonhydrostatic_boundary_tendencies.jl | 22 +++++++++---------- 2 files changed, 16 insertions(+), 16 deletions(-) diff --git a/src/Models/HydrostaticFreeSurfaceModels/compute_hydrostatic_free_surface_boundary_tendencies.jl b/src/Models/HydrostaticFreeSurfaceModels/compute_hydrostatic_free_surface_boundary_tendencies.jl index 271e5562b6..791340c3a0 100644 --- a/src/Models/HydrostaticFreeSurfaceModels/compute_hydrostatic_free_surface_boundary_tendencies.jl +++ b/src/Models/HydrostaticFreeSurfaceModels/compute_hydrostatic_free_surface_boundary_tendencies.jl @@ -58,12 +58,12 @@ function boundary_w_kernel_parameters(grid, arch) # Offsets in tangential direction are == -1 to # cover the required corners - param_west = (-Hx+1:0, 0:Ny+1, 1:Nz) - param_east = (Nx+1:Nx+Hx-1, 0:Ny+1, 1:Nz) - param_north = (0:Nx+1, -Hy+1:0, 1:Nz) - param_south = (0:Nx+1, Ny+1:Ny+Hy-1, 1:Nz) + param_west = (-Hx+1:0, 0:Ny+1) + param_east = (Nx+1:Nx+Hx-1, 0:Ny+1) + param_south = (0:Nx+1, -Hy+1:0) + param_north = (0:Nx+1, Ny+1:Ny+Hy-1) - params = (param_west, param_east, param_north, param_south) + params = (param_west, param_east, param_south, param_north) return boundary_parameters(params, grid, arch) end diff --git a/src/Models/NonhydrostaticModels/compute_nonhydrostatic_boundary_tendencies.jl b/src/Models/NonhydrostaticModels/compute_nonhydrostatic_boundary_tendencies.jl index 80294d828c..b481b7cdb8 100644 --- a/src/Models/NonhydrostaticModels/compute_nonhydrostatic_boundary_tendencies.jl +++ b/src/Models/NonhydrostaticModels/compute_nonhydrostatic_boundary_tendencies.jl @@ -31,10 +31,10 @@ function boundary_tendency_kernel_parameters(grid, arch) param_west = (1:Hx, 1:Ny, 1:Nz) param_east = (Nx-Hx+1:Nx, 1:Ny, 1:Nz) - param_north = (1:Nx, 1:Hy, 1:Nz) - param_south = (1:Nx, Ny-Hy+1:Ny, 1:Nz) + param_south = (1:Nx, 1:Hy, 1:Nz) + param_north = (1:Nx, Ny-Hy+1:Ny, 1:Nz) - params = (param_west, param_east, param_north, param_south) + params = (param_west, param_east, param_south, param_north) return boundary_parameters(params, grid, arch) end @@ -43,12 +43,12 @@ end function boundary_p_kernel_parameters(grid, arch) Nx, Ny, _ = size(grid) - param_west = (0:0, 1:Ny, 1:Nz) - param_east = (Nx+1:Nx+1, 1:Ny, 1:Nz) - param_north = (1:Nx, 0:0, 1:Nz) - param_south = (1:Nx, Ny+1:Ny+1, 1:Nz) + param_west = (0:0, 1:Ny) + param_east = (Nx+1:Nx+1, 1:Ny) + param_south = (1:Nx, 0:0) + param_north = (1:Nx, Ny+1:Ny+1) - params = (param_west, param_east, param_north, param_south) + params = (param_west, param_east, param_south, param_north) return boundary_parameters(params, grid, arch) end @@ -62,10 +62,10 @@ function boundary_κ_kernel_parameters(grid, closure, arch) param_west = (0:Bx, 1:Ny, 1:Nz) param_east = (Nx-Bx+1:Nx+1, 1:Ny, 1:Nz) - param_north = (1:Nx, 0:By, 1:Nz) - param_south = (1:Nx, Ny-By+1:Ny+1, 1:Nz) + param_south = (1:Nx, 0:By, 1:Nz) + param_north = (1:Nx, Ny-By+1:Ny+1, 1:Nz) - params = (param_west, param_east, param_north, param_south) + params = (param_west, param_east, param_south, param_north) return boundary_parameters(params, grid, arch) end From b03dc9d237d22e2765ece1662ef93825a974f48c Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Tue, 8 Oct 2024 13:04:24 +0200 Subject: [PATCH 57/68] another bugfix --- .../compute_nonhydrostatic_boundary_tendencies.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/Models/NonhydrostaticModels/compute_nonhydrostatic_boundary_tendencies.jl b/src/Models/NonhydrostaticModels/compute_nonhydrostatic_boundary_tendencies.jl index b481b7cdb8..cfe3691075 100644 --- a/src/Models/NonhydrostaticModels/compute_nonhydrostatic_boundary_tendencies.jl +++ b/src/Models/NonhydrostaticModels/compute_nonhydrostatic_boundary_tendencies.jl @@ -82,6 +82,6 @@ function boundary_parameters(parameters, grid, arch) include_side = (include_west, include_east, include_south, include_north) - return Tuple(KernelParameters(parameters[i]) for i in findall(include_side)) + return Tuple(KernelParameters(parameters[i]...) for i in findall(include_side)) end From 7cc91496c32d31573aacf1d7f63cf87812a5c3b3 Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Tue, 8 Oct 2024 15:01:35 +0200 Subject: [PATCH 58/68] make sure the range is correct --- ...ompute_hydrostatic_free_surface_boundary_tendencies.jl | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/Models/HydrostaticFreeSurfaceModels/compute_hydrostatic_free_surface_boundary_tendencies.jl b/src/Models/HydrostaticFreeSurfaceModels/compute_hydrostatic_free_surface_boundary_tendencies.jl index 791340c3a0..bbc13a43d1 100644 --- a/src/Models/HydrostaticFreeSurfaceModels/compute_hydrostatic_free_surface_boundary_tendencies.jl +++ b/src/Models/HydrostaticFreeSurfaceModels/compute_hydrostatic_free_surface_boundary_tendencies.jl @@ -58,10 +58,10 @@ function boundary_w_kernel_parameters(grid, arch) # Offsets in tangential direction are == -1 to # cover the required corners - param_west = (-Hx+1:0, 0:Ny+1) - param_east = (Nx+1:Nx+Hx-1, 0:Ny+1) - param_south = (0:Nx+1, -Hy+1:0) - param_north = (0:Nx+1, Ny+1:Ny+Hy-1) + param_west = (-Hx+2:1, 0:Ny+1) + param_east = (Nx:Nx+Hx-1, 0:Ny+1) + param_south = (0:Nx+1, -Hy+2:1) + param_north = (0:Nx+1, Ny:Ny+Hy-1) params = (param_west, param_east, param_south, param_north) From c90934449e50d1c0494fc4e544518b737a64e27e Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Wed, 9 Oct 2024 17:04:11 +0200 Subject: [PATCH 59/68] Update src/Advection/adapt_advection_order.jl Co-authored-by: Gregory L. Wagner --- src/Advection/adapt_advection_order.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/Advection/adapt_advection_order.jl b/src/Advection/adapt_advection_order.jl index 288723ffa4..70479af3f9 100644 --- a/src/Advection/adapt_advection_order.jl +++ b/src/Advection/adapt_advection_order.jl @@ -80,7 +80,7 @@ function adapt_advection_order(advection::Centered{B}, N::Int, grid::AbstractGri if N >= B return advection else - return Centered(; order = N * 2) + return Centered(; order = 2N) end end From 58699525370f34861621a7ae2e8291ae0ad0c2c8 Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Wed, 9 Oct 2024 17:04:44 +0200 Subject: [PATCH 60/68] Update src/Advection/adapt_advection_order.jl Co-authored-by: Gregory L. Wagner --- src/Advection/adapt_advection_order.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/Advection/adapt_advection_order.jl b/src/Advection/adapt_advection_order.jl index 70479af3f9..90765e429e 100644 --- a/src/Advection/adapt_advection_order.jl +++ b/src/Advection/adapt_advection_order.jl @@ -106,6 +106,6 @@ function adapt_advection_order(advection::WENO{B, FT, XT, YT, ZT}, N::Int, grid: if N >= B return advection else - return new_weno_scheme(advection, grid, N * 2 - 1, advection.bounds, XT, YT, ZT) + return new_weno_scheme(advection, grid, 2N - 1, advection.bounds, XT, YT, ZT) end end \ No newline at end of file From 13c21aab36d247cc0bf75288cc6aa8074fd1b1a5 Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Wed, 9 Oct 2024 17:04:49 +0200 Subject: [PATCH 61/68] Update src/Advection/adapt_advection_order.jl Co-authored-by: Gregory L. Wagner --- src/Advection/adapt_advection_order.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/Advection/adapt_advection_order.jl b/src/Advection/adapt_advection_order.jl index 90765e429e..08c57e4aa7 100644 --- a/src/Advection/adapt_advection_order.jl +++ b/src/Advection/adapt_advection_order.jl @@ -88,7 +88,7 @@ function adapt_advection_order(advection::UpwindBiased{B}, N::Int, grid::Abstrac if N >= B return advection else - return UpwindBiased(; order = N * 2 - 1) + return UpwindBiased(; order = 2N - 1) end end From 96c5edab812e965a09dad953e297fbd4e49a3d5e Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Wed, 9 Oct 2024 17:34:22 +0200 Subject: [PATCH 62/68] remove comments --- src/Advection/adapt_advection_order.jl | 41 ++++++------------- .../scalar_biharmonic_diffusivity.jl | 24 +---------- .../scalar_diffusivity.jl | 24 +---------- 3 files changed, 14 insertions(+), 75 deletions(-) diff --git a/src/Advection/adapt_advection_order.jl b/src/Advection/adapt_advection_order.jl index 288723ffa4..311e4f6242 100644 --- a/src/Advection/adapt_advection_order.jl +++ b/src/Advection/adapt_advection_order.jl @@ -16,35 +16,9 @@ If the order of advection is changed in at least one direction, the adapted adve by this function is a `FluxFormAdvection`. """ function adapt_advection_order(advection, grid::AbstractGrid) - advection_x = adapt_advection_order(advection, size(grid, 1), grid) - advection_y = adapt_advection_order(advection, size(grid, 2), grid) - advection_z = adapt_advection_order(advection, size(grid, 3), grid) - - # Check that we indeed changed the advection operator - changed_x = advection_x != advection - changed_y = advection_y != advection - changed_z = advection_z != advection - - new_advection = FluxFormAdvection(advection_x, advection_y, 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 - -function adapt_advection_order(advection::FluxFormAdvection, grid::AbstractGrid) - advection_x = adapt_advection_order(advection.x, size(grid, 1), grid) - advection_y = adapt_advection_order(advection.y, size(grid, 2), grid) - advection_z = adapt_advection_order(advection.z, size(grid, 3), grid) + advection_x = adapt_advection_order(x_advection(advection), size(grid, 1), grid) + advection_y = adapt_advection_order(y_advection(advection), size(grid, 2), grid) + advection_z = adapt_advection_order(z_advection(advection), size(grid, 3), grid) # Check that we indeed changed the advection operator changed_x = advection_x != advection.x @@ -67,6 +41,15 @@ function adapt_advection_order(advection::FluxFormAdvection, grid::AbstractGrid) 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 diff --git a/src/TurbulenceClosures/turbulence_closure_implementations/scalar_biharmonic_diffusivity.jl b/src/TurbulenceClosures/turbulence_closure_implementations/scalar_biharmonic_diffusivity.jl index e35ecc649c..6a446659f0 100644 --- a/src/TurbulenceClosures/turbulence_closure_implementations/scalar_biharmonic_diffusivity.jl +++ b/src/TurbulenceClosures/turbulence_closure_implementations/scalar_biharmonic_diffusivity.jl @@ -119,26 +119,4 @@ function on_architecture(to, closure::ScalarBiharmonicDiffusivity{F, <:Any, <:An ν = on_architecture(to, closure.ν) κ = on_architecture(to, closure.κ) return ScalarBiharmonicDiffusivity{F, N}(ν, κ) -end - -# We do not have at the moment the ability to distinguish between the halos required in the x, y and z -# direction for a general closure since the closure could contain any user-defined function - -# In the end we probably want something like this: -#= -required_halo_size_x(::ScalarBiharmonicDiffusivity{<:Any, ThreeDimensionalFormulation, N}) where N = N -required_halo_size_y(::ScalarBiharmonicDiffusivity{<:Any, ThreeDimensionalFormulation, N}) where N = N -required_halo_size_z(::ScalarBiharmonicDiffusivity{<:Any, ThreeDimensionalFormulation, N}) where N = N - -required_halo_size_x(::ScalarBiharmonicDiffusivity{<:Any, VerticalFormulation, N}) where N = 0 -required_halo_size_y(::ScalarBiharmonicDiffusivity{<:Any, VerticalFormulation, N}) where N = 0 -required_halo_size_z(::ScalarBiharmonicDiffusivity{<:Any, VerticalFormulation, N}) where N = N - -required_halo_size_x(::ScalarBiharmonicDiffusivity{<:Any, HorizontalFormulation, N}) where N = N -required_halo_size_y(::ScalarBiharmonicDiffusivity{<:Any, HorizontalFormulation, N}) where N = N -required_halo_size_z(::ScalarBiharmonicDiffusivity{<:Any, HorizontalFormulation, N}) where N = 0 - -required_halo_size_x(::ScalarBiharmonicDiffusivity{<:Any, HorizontalDivergenceFormulation, N}) where N = N -required_halo_size_y(::ScalarBiharmonicDiffusivity{<:Any, HorizontalDivergenceFormulation, N}) where N = N -required_halo_size_z(::ScalarBiharmonicDiffusivity{<:Any, HorizontalDivergenceFormulation, N}) where N = 0 -=# \ No newline at end of file +end \ No newline at end of file diff --git a/src/TurbulenceClosures/turbulence_closure_implementations/scalar_diffusivity.jl b/src/TurbulenceClosures/turbulence_closure_implementations/scalar_diffusivity.jl index 2d54b16cd8..3b8b683f77 100644 --- a/src/TurbulenceClosures/turbulence_closure_implementations/scalar_diffusivity.jl +++ b/src/TurbulenceClosures/turbulence_closure_implementations/scalar_diffusivity.jl @@ -219,26 +219,4 @@ function on_architecture(to, closure::ScalarDiffusivity{TD, F, <:Any, <:Any, N}) ν = on_architecture(to, closure.ν) κ = on_architecture(to, closure.κ) return ScalarDiffusivity{TD, F, N}(ν, κ) -end - -# We do not have at the moment the ability to distinguish between the halos required in the x, y and z -# direction for a general closure since the closure could contain any user-defined function - -# In the end we probably want something like this: -#= -required_halo_size_x(::ScalarDiffusivity{<:Any, ThreeDimensionalFormulation, N}) where N = N -required_halo_size_y(::ScalarDiffusivity{<:Any, ThreeDimensionalFormulation, N}) where N = N -required_halo_size_z(::ScalarDiffusivity{<:Any, ThreeDimensionalFormulation, N}) where N = N - -required_halo_size_x(::ScalarDiffusivity{<:Any, VerticalFormulation, N}) where N = 0 -required_halo_size_y(::ScalarDiffusivity{<:Any, VerticalFormulation, N}) where N = 0 -required_halo_size_z(::ScalarDiffusivity{<:Any, VerticalFormulation, N}) where N = N - -required_halo_size_x(::ScalarDiffusivity{<:Any, HorizontalFormulation, N}) where N = N -required_halo_size_y(::ScalarDiffusivity{<:Any, HorizontalFormulation, N}) where N = N -required_halo_size_z(::ScalarDiffusivity{<:Any, HorizontalFormulation, N}) where N = 0 - -required_halo_size_x(::ScalarDiffusivity{<:Any, HorizontalDivergenceFormulation, N}) where N = N -required_halo_size_y(::ScalarDiffusivity{<:Any, HorizontalDivergenceFormulation, N}) where N = N -required_halo_size_z(::ScalarDiffusivity{<:Any, HorizontalDivergenceFormulation, N}) where N = 0 -=# \ No newline at end of file +end \ No newline at end of file From da68014b116fd7e1d4d37e5ec6c82cbdd7f10a21 Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Wed, 9 Oct 2024 19:06:28 +0200 Subject: [PATCH 63/68] this should be type stable for the moment --- .../scalar_biharmonic_diffusivity.jl | 5 +++++ .../scalar_diffusivity.jl | 12 +++++++++--- test/test_enzyme.jl | 16 +++++++--------- 3 files changed, 21 insertions(+), 12 deletions(-) diff --git a/src/TurbulenceClosures/turbulence_closure_implementations/scalar_biharmonic_diffusivity.jl b/src/TurbulenceClosures/turbulence_closure_implementations/scalar_biharmonic_diffusivity.jl index 6a446659f0..17df25ccc5 100644 --- a/src/TurbulenceClosures/turbulence_closure_implementations/scalar_biharmonic_diffusivity.jl +++ b/src/TurbulenceClosures/turbulence_closure_implementations/scalar_biharmonic_diffusivity.jl @@ -82,6 +82,11 @@ function ScalarBiharmonicDiffusivity(formulation = ThreeDimensionalFormulation() ν = convert_diffusivity(FT, ν; discrete_form, loc, parameters) κ = convert_diffusivity(FT, κ; discrete_form, loc, parameters) + + if ν isa Number && κ isa Number + return ScalarBiharmonicDiffusivity{typeof(formulation), 2}(ν, κ) + end + return ScalarBiharmonicDiffusivity{typeof(formulation), required_halo_size}(ν, κ) end diff --git a/src/TurbulenceClosures/turbulence_closure_implementations/scalar_diffusivity.jl b/src/TurbulenceClosures/turbulence_closure_implementations/scalar_diffusivity.jl index 3b8b683f77..cbe1b4cbb8 100644 --- a/src/TurbulenceClosures/turbulence_closure_implementations/scalar_diffusivity.jl +++ b/src/TurbulenceClosures/turbulence_closure_implementations/scalar_diffusivity.jl @@ -114,12 +114,13 @@ ScalarDiffusivity{ExplicitTimeDiscretization}(ν=0.0, κ=Oceananigans.Turbulence ``` """ function ScalarDiffusivity(time_discretization=ExplicitTimeDiscretization(), - formulation=ThreeDimensionalFormulation(), FT=Float64; + formulation=ThreeDimensionalFormulation(), + FT=Float64, + required_halo_size::Int = 1; ν=0, κ=0, discrete_form = false, loc = (nothing, nothing, nothing), - parameters = nothing, - required_halo_size::Int = 1) + parameters = nothing) if formulation == HorizontalFormulation() && time_discretization == VerticallyImplicitTimeDiscretization() throw(ArgumentError("VerticallyImplicitTimeDiscretization is only supported for \ @@ -129,6 +130,11 @@ function ScalarDiffusivity(time_discretization=ExplicitTimeDiscretization(), κ = convert_diffusivity(FT, κ; discrete_form, loc, parameters) ν = convert_diffusivity(FT, ν; discrete_form, loc, parameters) + # Force a type-stable constructor if ν and κ are numbers + if ν isa Number && κ isa Number + return ScalarDiffusivity{typeof(time_discretization), typeof(formulation), 1}(ν, κ) + end + return ScalarDiffusivity{typeof(time_discretization), typeof(formulation), required_halo_size}(ν, κ) end diff --git a/test/test_enzyme.jl b/test/test_enzyme.jl index b5f027a42a..00ff0b90ba 100644 --- a/test/test_enzyme.jl +++ b/test/test_enzyme.jl @@ -16,8 +16,10 @@ const maximum_diffusivity = 100 Change diffusivity of model to `diffusivity`. """ function set_diffusivity!(model, diffusivity) - closure = model.closure - fill!(closure.κ.c, diffusivity) + closure = VerticalScalarDiffusivity(; κ=diffusivity) + names = tuple(:c) # tracernames(model.tracers) + closure = with_tracers(names, closure) + model.closure = closure return nothing end @@ -167,10 +169,7 @@ end topology = (Periodic, Periodic, Bounded) grid = RectilinearGrid(size=(Nx, Ny, Nz); x, y, z, topology) - κ = ZFaceField(grid) - fill!(κ, 0.1) - - diffusion = VerticalScalarDiffusivity(; κ) + diffusion = VerticalScalarDiffusivity(κ=0.1) u = XFaceField(grid) v = YFaceField(grid) @@ -215,8 +214,7 @@ end # Now for real amplitude = 1.0 - diffusivity = 1.0 - + κ = 1.0 dmodel = Enzyme.make_zero(model) set_diffusivity!(dmodel, 0) @@ -224,7 +222,7 @@ end stable_diffusion!, Duplicated(model, dmodel), Const(amplitude), - Active(diffusivity)) + Active(κ)) @info """ \n Enzyme computed $dc²_dκ From 9f3fffcc6981c74a2cff86775ad5fcb3cd065062 Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Thu, 10 Oct 2024 09:45:45 +0200 Subject: [PATCH 64/68] bugfix --- src/Advection/adapt_advection_order.jl | 18 +++++++++++------- 1 file changed, 11 insertions(+), 7 deletions(-) diff --git a/src/Advection/adapt_advection_order.jl b/src/Advection/adapt_advection_order.jl index 241a5b18bb..41f3f60c9a 100644 --- a/src/Advection/adapt_advection_order.jl +++ b/src/Advection/adapt_advection_order.jl @@ -16,16 +16,20 @@ If the order of advection is changed in at least one direction, the adapted adve by this function is a `FluxFormAdvection`. """ function adapt_advection_order(advection, grid::AbstractGrid) - advection_x = adapt_advection_order(x_advection(advection), size(grid, 1), grid) - advection_y = adapt_advection_order(y_advection(advection), size(grid, 2), grid) - advection_z = adapt_advection_order(z_advection(advection), size(grid, 3), grid) + 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 = advection_x != advection.x - changed_y = advection_y != advection.y - changed_z = advection_z != advection.z + changed_x = new_advection_x != advection_x + changed_y = new_advection_y != advection_y + changed_z = new_advection_z != advection_z - new_advection = FluxFormAdvection(advection_x, advection_y, 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 From 8bbaf03d41adfe4beeebbe1cc0eb5fa5a1ca1613 Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Thu, 10 Oct 2024 14:36:13 +0200 Subject: [PATCH 65/68] bugfix --- .../scalar_biharmonic_diffusivity.jl | 3 +++ .../scalar_diffusivity.jl | 10 ++++++---- 2 files changed, 9 insertions(+), 4 deletions(-) diff --git a/src/TurbulenceClosures/turbulence_closure_implementations/scalar_biharmonic_diffusivity.jl b/src/TurbulenceClosures/turbulence_closure_implementations/scalar_biharmonic_diffusivity.jl index 17df25ccc5..8b6bae5c7e 100644 --- a/src/TurbulenceClosures/turbulence_closure_implementations/scalar_biharmonic_diffusivity.jl +++ b/src/TurbulenceClosures/turbulence_closure_implementations/scalar_biharmonic_diffusivity.jl @@ -83,6 +83,9 @@ function ScalarBiharmonicDiffusivity(formulation = ThreeDimensionalFormulation() ν = convert_diffusivity(FT, ν; discrete_form, loc, parameters) κ = convert_diffusivity(FT, κ; discrete_form, loc, parameters) + # Force a type-stable constructor if ν and κ are numbers + # This particular short-circuiting of the constructor is necessary to perform parameter + # estimation of the diffusivity coefficients using autodiff. if ν isa Number && κ isa Number return ScalarBiharmonicDiffusivity{typeof(formulation), 2}(ν, κ) end diff --git a/src/TurbulenceClosures/turbulence_closure_implementations/scalar_diffusivity.jl b/src/TurbulenceClosures/turbulence_closure_implementations/scalar_diffusivity.jl index cbe1b4cbb8..851c1a9253 100644 --- a/src/TurbulenceClosures/turbulence_closure_implementations/scalar_diffusivity.jl +++ b/src/TurbulenceClosures/turbulence_closure_implementations/scalar_diffusivity.jl @@ -115,15 +115,15 @@ ScalarDiffusivity{ExplicitTimeDiscretization}(ν=0.0, κ=Oceananigans.Turbulence """ function ScalarDiffusivity(time_discretization=ExplicitTimeDiscretization(), formulation=ThreeDimensionalFormulation(), - FT=Float64, - required_halo_size::Int = 1; + FT=Float64; ν=0, κ=0, discrete_form = false, loc = (nothing, nothing, nothing), - parameters = nothing) + parameters = nothing, + required_halo_size::Int = 1) if formulation == HorizontalFormulation() && time_discretization == VerticallyImplicitTimeDiscretization() - throw(ArgumentError("VerticallyImplicitTimeDiscretization is only supported for \ + throw(ArgumentError("VerticallyImplicitTimeDiscretization is only supported for \ `VerticalFormulation` or `ThreeDimensionalFormulation`")) end @@ -131,6 +131,8 @@ function ScalarDiffusivity(time_discretization=ExplicitTimeDiscretization(), ν = convert_diffusivity(FT, ν; discrete_form, loc, parameters) # Force a type-stable constructor if ν and κ are numbers + # This particular short-circuiting of the constructor is necessary to perform parameter + # estimation of the diffusivity coefficients using autodiff. if ν isa Number && κ isa Number return ScalarDiffusivity{typeof(time_discretization), typeof(formulation), 1}(ν, κ) end From 402d0314db3e1607ee372497aa4f1a86641ffa74 Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Thu, 10 Oct 2024 14:36:39 +0200 Subject: [PATCH 66/68] better comment --- .../scalar_biharmonic_diffusivity.jl | 2 +- .../turbulence_closure_implementations/scalar_diffusivity.jl | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/src/TurbulenceClosures/turbulence_closure_implementations/scalar_biharmonic_diffusivity.jl b/src/TurbulenceClosures/turbulence_closure_implementations/scalar_biharmonic_diffusivity.jl index 8b6bae5c7e..003fabbeea 100644 --- a/src/TurbulenceClosures/turbulence_closure_implementations/scalar_biharmonic_diffusivity.jl +++ b/src/TurbulenceClosures/turbulence_closure_implementations/scalar_biharmonic_diffusivity.jl @@ -84,7 +84,7 @@ function ScalarBiharmonicDiffusivity(formulation = ThreeDimensionalFormulation() κ = convert_diffusivity(FT, κ; discrete_form, loc, parameters) # Force a type-stable constructor if ν and κ are numbers - # This particular short-circuiting of the constructor is necessary to perform parameter + # This particular short-circuiting of the required_halo_size kwargs is necessary to perform parameter # estimation of the diffusivity coefficients using autodiff. if ν isa Number && κ isa Number return ScalarBiharmonicDiffusivity{typeof(formulation), 2}(ν, κ) diff --git a/src/TurbulenceClosures/turbulence_closure_implementations/scalar_diffusivity.jl b/src/TurbulenceClosures/turbulence_closure_implementations/scalar_diffusivity.jl index 851c1a9253..a4f85814b8 100644 --- a/src/TurbulenceClosures/turbulence_closure_implementations/scalar_diffusivity.jl +++ b/src/TurbulenceClosures/turbulence_closure_implementations/scalar_diffusivity.jl @@ -131,7 +131,7 @@ function ScalarDiffusivity(time_discretization=ExplicitTimeDiscretization(), ν = convert_diffusivity(FT, ν; discrete_form, loc, parameters) # Force a type-stable constructor if ν and κ are numbers - # This particular short-circuiting of the constructor is necessary to perform parameter + # This particular short-circuiting of the required_halo_size kwargs is necessary to perform parameter # estimation of the diffusivity coefficients using autodiff. if ν isa Number && κ isa Number return ScalarDiffusivity{typeof(time_discretization), typeof(formulation), 1}(ν, κ) From 40e5dfae865473e4185b2ffdff19d9d82a51311e Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Fri, 11 Oct 2024 09:18:30 +0200 Subject: [PATCH 67/68] fix the docs --- docs/src/model_setup/buoyancy_and_equation_of_state.md | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/docs/src/model_setup/buoyancy_and_equation_of_state.md b/docs/src/model_setup/buoyancy_and_equation_of_state.md index c547374072..42b670078e 100644 --- a/docs/src/model_setup/buoyancy_and_equation_of_state.md +++ b/docs/src/model_setup/buoyancy_and_equation_of_state.md @@ -64,7 +64,7 @@ HydrostaticFreeSurfaceModel{CPU, RectilinearGrid}(time = 0 seconds, iteration = ├── free surface: ImplicitFreeSurface with gravitational acceleration 9.80665 m s⁻² │ └── solver: FFTImplicitFreeSurfaceSolver ├── advection scheme: -│ └── momentum: Centered reconstruction order 2 +│ └── momentum: Vector Invariant, Dimension-by-dimension reconstruction └── coriolis: Nothing ``` @@ -98,7 +98,7 @@ HydrostaticFreeSurfaceModel{CPU, RectilinearGrid}(time = 0 seconds, iteration = ├── free surface: ImplicitFreeSurface with gravitational acceleration 9.80665 m s⁻² │ └── solver: FFTImplicitFreeSurfaceSolver ├── advection scheme: -│ ├── momentum: Centered reconstruction order 2 +│ ├── momentum: Vector Invariant, Dimension-by-dimension reconstruction │ └── b: Centered reconstruction order 2 └── coriolis: Nothing ``` @@ -139,7 +139,7 @@ HydrostaticFreeSurfaceModel{CPU, RectilinearGrid}(time = 0 seconds, iteration = ├── free surface: ImplicitFreeSurface with gravitational acceleration 9.80665 m s⁻² │ └── solver: FFTImplicitFreeSurfaceSolver ├── advection scheme: -│ ├── momentum: Centered reconstruction order 2 +│ ├── momentum: Vector Invariant, Dimension-by-dimension reconstruction │ ├── T: Centered reconstruction order 2 │ └── S: Centered reconstruction order 2 └── coriolis: Nothing From 92141102111b572936260e4ab995c1c3a26eb67f Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Fri, 11 Oct 2024 13:48:40 +0200 Subject: [PATCH 68/68] whoops --- docs/src/model_setup/buoyancy_and_equation_of_state.md | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/docs/src/model_setup/buoyancy_and_equation_of_state.md b/docs/src/model_setup/buoyancy_and_equation_of_state.md index 42b670078e..c547374072 100644 --- a/docs/src/model_setup/buoyancy_and_equation_of_state.md +++ b/docs/src/model_setup/buoyancy_and_equation_of_state.md @@ -64,7 +64,7 @@ HydrostaticFreeSurfaceModel{CPU, RectilinearGrid}(time = 0 seconds, iteration = ├── free surface: ImplicitFreeSurface with gravitational acceleration 9.80665 m s⁻² │ └── solver: FFTImplicitFreeSurfaceSolver ├── advection scheme: -│ └── momentum: Vector Invariant, Dimension-by-dimension reconstruction +│ └── momentum: Centered reconstruction order 2 └── coriolis: Nothing ``` @@ -98,7 +98,7 @@ HydrostaticFreeSurfaceModel{CPU, RectilinearGrid}(time = 0 seconds, iteration = ├── free surface: ImplicitFreeSurface with gravitational acceleration 9.80665 m s⁻² │ └── solver: FFTImplicitFreeSurfaceSolver ├── advection scheme: -│ ├── momentum: Vector Invariant, Dimension-by-dimension reconstruction +│ ├── momentum: Centered reconstruction order 2 │ └── b: Centered reconstruction order 2 └── coriolis: Nothing ``` @@ -139,7 +139,7 @@ HydrostaticFreeSurfaceModel{CPU, RectilinearGrid}(time = 0 seconds, iteration = ├── free surface: ImplicitFreeSurface with gravitational acceleration 9.80665 m s⁻² │ └── solver: FFTImplicitFreeSurfaceSolver ├── advection scheme: -│ ├── momentum: Vector Invariant, Dimension-by-dimension reconstruction +│ ├── momentum: Centered reconstruction order 2 │ ├── T: Centered reconstruction order 2 │ └── S: Centered reconstruction order 2 └── coriolis: Nothing