diff --git a/docs/src/model_setup/boundary_conditions.md b/docs/src/model_setup/boundary_conditions.md index a8c0656012..2a3ebc60d5 100644 --- a/docs/src/model_setup/boundary_conditions.md +++ b/docs/src/model_setup/boundary_conditions.md @@ -43,7 +43,7 @@ julia> model = NonhydrostaticModel(; grid, boundary_conditions=(u=no_slip_field_ NonhydrostaticModel{CPU, RectilinearGrid}(time = 0 seconds, iteration = 0) ├── grid: 16×16×16 RectilinearGrid{Float64, Periodic, Bounded, Bounded} on CPU with 3×3×3 halo ├── timestepper: RungeKutta3TimeStepper -├── advection scheme: Centered reconstruction order 2 +├── advection scheme: Centered(order=2) ├── tracers: () ├── closure: Nothing ├── buoyancy: Nothing @@ -405,7 +405,7 @@ julia> model = NonhydrostaticModel(grid=grid, boundary_conditions=(u=u_bcs, c=c_ NonhydrostaticModel{CPU, RectilinearGrid}(time = 0 seconds, iteration = 0) ├── grid: 16×16×16 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo ├── timestepper: RungeKutta3TimeStepper -├── advection scheme: Centered reconstruction order 2 +├── advection scheme: Centered(order=2) ├── tracers: c ├── closure: Nothing ├── buoyancy: Nothing 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..d21df10c0f 100644 --- a/docs/src/model_setup/buoyancy_and_equation_of_state.md +++ b/docs/src/model_setup/buoyancy_and_equation_of_state.md @@ -29,7 +29,7 @@ julia> model = NonhydrostaticModel(; grid, buoyancy=nothing) NonhydrostaticModel{CPU, RectilinearGrid}(time = 0 seconds, iteration = 0) ├── grid: 8×8×8 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo ├── timestepper: RungeKutta3TimeStepper -├── advection scheme: Centered reconstruction order 2 +├── advection scheme: Centered(order=2) ├── tracers: () ├── closure: Nothing ├── buoyancy: Nothing @@ -44,7 +44,7 @@ julia> model = NonhydrostaticModel(; grid) NonhydrostaticModel{CPU, RectilinearGrid}(time = 0 seconds, iteration = 0) ├── grid: 8×8×8 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo ├── timestepper: RungeKutta3TimeStepper -├── advection scheme: Centered reconstruction order 2 +├── advection scheme: Centered(order=2) ├── tracers: () ├── closure: Nothing ├── buoyancy: Nothing @@ -78,7 +78,7 @@ julia> model = NonhydrostaticModel(; grid, buoyancy=BuoyancyTracer(), tracers=:b NonhydrostaticModel{CPU, RectilinearGrid}(time = 0 seconds, iteration = 0) ├── grid: 8×8×8 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo ├── timestepper: RungeKutta3TimeStepper -├── advection scheme: Centered reconstruction order 2 +├── advection scheme: Centered(order=2) ├── tracers: b ├── closure: Nothing ├── buoyancy: BuoyancyTracer with ĝ = NegativeZDirection() @@ -99,7 +99,7 @@ HydrostaticFreeSurfaceModel{CPU, RectilinearGrid}(time = 0 seconds, iteration = │ └── solver: FFTImplicitFreeSurfaceSolver ├── advection scheme: │ ├── momentum: Vector Invariant, Dimension-by-dimension reconstruction -│ └── b: Centered reconstruction order 2 +│ └── b: Centered(order=2) └── coriolis: Nothing ``` @@ -119,7 +119,7 @@ julia> model = NonhydrostaticModel(; grid, buoyancy=SeawaterBuoyancy(), tracers= NonhydrostaticModel{CPU, RectilinearGrid}(time = 0 seconds, iteration = 0) ├── grid: 8×8×8 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo ├── timestepper: RungeKutta3TimeStepper -├── advection scheme: Centered reconstruction order 2 +├── advection scheme: Centered(order=2) ├── tracers: (T, S) ├── closure: Nothing ├── buoyancy: SeawaterBuoyancy with g=9.80665 and LinearEquationOfState(thermal_expansion=0.000167, haline_contraction=0.00078) with ĝ = NegativeZDirection() @@ -140,8 +140,8 @@ HydrostaticFreeSurfaceModel{CPU, RectilinearGrid}(time = 0 seconds, iteration = │ └── solver: FFTImplicitFreeSurfaceSolver ├── advection scheme: │ ├── momentum: Vector Invariant, Dimension-by-dimension reconstruction -│ ├── T: Centered reconstruction order 2 -│ └── S: Centered reconstruction order 2 +│ ├── T: Centered(order=2) +│ └── S: Centered(order=2) └── coriolis: Nothing ``` @@ -158,7 +158,7 @@ julia> model = NonhydrostaticModel(; grid, buoyancy, tracers=(:T, :S)) NonhydrostaticModel{CPU, RectilinearGrid}(time = 0 seconds, iteration = 0) ├── grid: 8×8×8 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo ├── timestepper: RungeKutta3TimeStepper -├── advection scheme: Centered reconstruction order 2 +├── advection scheme: Centered(order=2) ├── tracers: (T, S) ├── closure: Nothing ├── buoyancy: SeawaterBuoyancy with g=1.3 and LinearEquationOfState(thermal_expansion=0.000167, haline_contraction=0.00078) with ĝ = NegativeZDirection() @@ -239,7 +239,7 @@ julia> model = NonhydrostaticModel(; grid, buoyancy, tracers=:b) NonhydrostaticModel{CPU, RectilinearGrid}(time = 0 seconds, iteration = 0) ├── grid: 8×8×8 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo ├── timestepper: RungeKutta3TimeStepper -├── advection scheme: Centered reconstruction order 2 +├── advection scheme: Centered(order=2) ├── tracers: b ├── closure: Nothing ├── buoyancy: BuoyancyTracer with ĝ = (0.0, 0.707107, 0.707107) diff --git a/docs/src/model_setup/lagrangian_particles.md b/docs/src/model_setup/lagrangian_particles.md index 01e8f8fbe4..be10c534b5 100644 --- a/docs/src/model_setup/lagrangian_particles.md +++ b/docs/src/model_setup/lagrangian_particles.md @@ -44,7 +44,7 @@ model = NonhydrostaticModel(grid=grid, particles=lagrangian_particles) NonhydrostaticModel{CPU, RectilinearGrid}(time = 0 seconds, iteration = 0) ├── grid: 10×10×10 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo ├── timestepper: RungeKutta3TimeStepper -├── advection scheme: Centered reconstruction order 2 +├── advection scheme: Centered(order=2) ├── tracers: () ├── closure: Nothing ├── buoyancy: Nothing diff --git a/docs/src/model_setup/tracers.md b/docs/src/model_setup/tracers.md index 947cb6e9ff..2aa1668352 100644 --- a/docs/src/model_setup/tracers.md +++ b/docs/src/model_setup/tracers.md @@ -16,7 +16,7 @@ julia> model = NonhydrostaticModel(; grid) NonhydrostaticModel{CPU, RectilinearGrid}(time = 0 seconds, iteration = 0) ├── grid: 16×16×16 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo ├── timestepper: RungeKutta3TimeStepper -├── advection scheme: Centered reconstruction order 2 +├── advection scheme: Centered(order=2) ├── tracers: () ├── closure: Nothing ├── buoyancy: Nothing @@ -31,7 +31,7 @@ julia> model = NonhydrostaticModel(; grid, tracers=(:T, :S)) NonhydrostaticModel{CPU, RectilinearGrid}(time = 0 seconds, iteration = 0) ├── grid: 16×16×16 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo ├── timestepper: RungeKutta3TimeStepper -├── advection scheme: Centered reconstruction order 2 +├── advection scheme: Centered(order=2) ├── tracers: (T, S) ├── closure: Nothing ├── buoyancy: Nothing @@ -66,7 +66,7 @@ julia> model = NonhydrostaticModel(; grid, tracers=(:T, :S, :C₁, :CO₂, :nitr NonhydrostaticModel{CPU, RectilinearGrid}(time = 0 seconds, iteration = 0) ├── grid: 16×16×16 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo ├── timestepper: RungeKutta3TimeStepper -├── advection scheme: Centered reconstruction order 2 +├── advection scheme: Centered(order=2) ├── tracers: (T, S, C₁, CO₂, nitrogen) ├── closure: Nothing ├── buoyancy: Nothing diff --git a/src/Advection/adapt_advection_order.jl b/src/Advection/adapt_advection_order.jl index 41f3f60c9a..21c993e0a3 100644 --- a/src/Advection/adapt_advection_order.jl +++ b/src/Advection/adapt_advection_order.jl @@ -20,9 +20,14 @@ function adapt_advection_order(advection, grid::AbstractGrid) 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) + tx = topology(grid, 1)() + ty = topology(grid, 2)() + tz = topology(grid, 3)() + + # Note: no adaptation in Flat directions. + new_advection_x = adapt_advection_order(tx, advection_x, size(grid, 1), grid) + new_advection_y = adapt_advection_order(ty, advection_y, size(grid, 2), grid) + new_advection_z = adapt_advection_order(tz, advection_z, size(grid, 3), grid) # Check that we indeed changed the advection operator changed_x = new_advection_x != advection_x @@ -45,7 +50,6 @@ function adapt_advection_order(advection, 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 @@ -54,10 +58,15 @@ x_advection(advection) = advection y_advection(advection) = advection z_advection(advection) = advection +# Don't adapt advection in Flat directions +adapt_advection_order(topo, advection, N, grid) = adapt_advection_order(advection, N, grid) +adapt_advection_order(::Flat, advection, N, grid) = advection + # For the moment, we do not adapt the advection order for the VectorInvariant advection scheme adapt_advection_order(advection::VectorInvariant, grid::AbstractGrid) = advection adapt_advection_order(advection::Nothing, grid::AbstractGrid) = nothing adapt_advection_order(advection::Nothing, N::Int, grid::AbstractGrid) = nothing +adapt_advection_order(advection, N, grid) = advection ##### ##### Directional adapt advection order @@ -67,7 +76,7 @@ function adapt_advection_order(advection::Centered{B}, N::Int, grid::AbstractGri if N >= B return advection else - return Centered(; order = 2N) + return Centered(; order=2N) end end @@ -75,15 +84,18 @@ function adapt_advection_order(advection::UpwindBiased{B}, N::Int, grid::Abstrac if N >= B return advection else - return UpwindBiased(; order = 2N - 1) + return UpwindBiased(; order=2N-1) end end """ new_weno_scheme(grid, order, bounds, XT, YT, ZT) -Constructs a new WENO scheme based on the given parameters. `XT`, `YT`, and `ZT` is the type of the precomputed weno coefficients in the -x-direction, y-direction and z-direction. A _non-stretched_ WENO scheme has `T` equal to `Nothing` everywhere. In case of a non-stretched WENO scheme, +Constructs a new WENO scheme based on the given parameters. +`XT`, `YT`, and `ZT` is the type of the precomputed weno coefficients in the +x-direction, y-direction and z-direction. +A _non-stretched_ WENO scheme has `T` equal to `Nothing` everywhere. +In case of a non-stretched WENO scheme, we rebuild the advection without passing the grid information, otherwise we use the grid to account for stretched directions. """ new_weno_scheme(::WENO, grid, order, bounds, ::Type{Nothing}, ::Type{Nothing}, ::Type{Nothing},) = WENO(; order, bounds) @@ -93,6 +105,7 @@ 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, 2N - 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 +end + diff --git a/src/Advection/centered_reconstruction.jl b/src/Advection/centered_reconstruction.jl index 5884d37b5a..403343564b 100644 --- a/src/Advection/centered_reconstruction.jl +++ b/src/Advection/centered_reconstruction.jl @@ -58,7 +58,7 @@ function Centered(FT::DataType = Float64; grid = nothing, order = 2) return Centered{N, FT}(coefficients..., buffer_scheme) end -Base.summary(a::Centered{N}) where N = string("Centered reconstruction order ", N*2) +Base.summary(a::Centered{N}) where N = string("Centered(order=", 2N, ")") Base.show(io::IO, a::Centered{N, FT, XT, YT, ZT}) where {N, FT, XT, YT, ZT} = print(io, summary(a), " \n", diff --git a/src/Advection/flux_form_advection.jl b/src/Advection/flux_form_advection.jl index 962f462ecb..6e4de07af8 100644 --- a/src/Advection/flux_form_advection.jl +++ b/src/Advection/flux_form_advection.jl @@ -26,8 +26,13 @@ function FluxFormAdvection(x_advection, y_advection, z_advection) return FluxFormAdvection{H, FT}(x_advection, y_advection, z_advection) end +Base.summary(scheme::FluxFormAdvection) = string("FluxFormAdvection(x=", + summary(scheme.x), ", y=", + summary(scheme.y), ", z=", + summary(scheme.z), ")") + Base.show(io::IO, scheme::FluxFormAdvection) = - print(io, "FluxFormAdvection with reconstructions: ", " \n", + print(io, "FluxFormAdvection with direction-based reconstructions:", " \n", " ├── x: ", summary(scheme.x), "\n", " ├── y: ", summary(scheme.y), "\n", " └── z: ", summary(scheme.z)) diff --git a/src/Advection/upwind_biased_reconstruction.jl b/src/Advection/upwind_biased_reconstruction.jl index eb3c6805ec..8269112353 100644 --- a/src/Advection/upwind_biased_reconstruction.jl +++ b/src/Advection/upwind_biased_reconstruction.jl @@ -65,7 +65,7 @@ function UpwindBiased(FT::DataType = Float64; grid = nothing, order = 3) return UpwindBiased{N, FT}(coefficients..., buffer_scheme, advecting_velocity_scheme) end -Base.summary(a::UpwindBiased{N}) where N = string("Upwind Biased reconstruction order ", N*2-1) +Base.summary(a::UpwindBiased{N}) where N = string("UpwindBiased(order=", 2N-1, ")") Base.show(io::IO, a::UpwindBiased{N, FT, XT, YT, ZT}) where {N, FT, XT, YT, ZT} = print(io, summary(a), " \n", diff --git a/src/Advection/vector_invariant_advection.jl b/src/Advection/vector_invariant_advection.jl index fd8e4ff92f..5d5becb2f5 100644 --- a/src/Advection/vector_invariant_advection.jl +++ b/src/Advection/vector_invariant_advection.jl @@ -88,16 +88,16 @@ Vector Invariant, Dimension-by-dimension reconstruction julia> using Oceananigans julia> VectorInvariant(vorticity_scheme = WENO(), vertical_scheme = WENO(order = 3)) -Vector Invariant, Dimension-by-dimension reconstruction - Vorticity flux scheme: - ├── WENO reconstruction order 5 +Vector Invariant, Dimension-by-dimension reconstruction + Vorticity flux scheme: + ├── WENO(order=5) └── smoothness ζ: Oceananigans.Advection.VelocityStencil() - Vertical advection / Divergence flux scheme: - ├── WENO reconstruction order 3 - └── upwinding treatment: OnlySelfUpwinding - KE gradient and Divergence flux cross terms reconstruction: - └── Centered reconstruction order 2 - Smoothness measures: + Vertical advection / Divergence flux scheme: + ├── WENO(order=3) + └── upwinding treatment: OnlySelfUpwinding + KE gradient and Divergence flux cross terms reconstruction: + └── Centered(order=2) + Smoothness measures: ├── smoothness δU: FunctionStencil f = divergence_smoothness ├── smoothness δV: FunctionStencil f = divergence_smoothness ├── smoothness δu²: FunctionStencil f = u_smoothness diff --git a/src/Advection/weno_reconstruction.jl b/src/Advection/weno_reconstruction.jl index f91c633557..d95a1a2c4f 100644 --- a/src/Advection/weno_reconstruction.jl +++ b/src/Advection/weno_reconstruction.jl @@ -58,14 +58,14 @@ Examples julia> using Oceananigans julia> WENO() -WENO reconstruction order 5 - Boundary scheme: - └── WENO reconstruction order 3 - Symmetric scheme: - └── Centered reconstruction order 4 +WENO(order=5) + Boundary scheme: + └── WENO(order=3) + Symmetric scheme: + └── Centered(order=4) Directions: - ├── X regular - ├── Y regular + ├── X regular + ├── Y regular └── Z regular ``` @@ -82,14 +82,14 @@ julia> grid = RectilinearGrid(size = (Nx, Nz), halo = (4, 4), topology=(Periodic x = (0, Lx), z = chebychev_spaced_z_faces); julia> WENO(grid; order=7) -WENO reconstruction order 7 - Boundary scheme: - └── WENO reconstruction order 5 - Symmetric scheme: - └── Centered reconstruction order 6 +WENO(order=7) + Boundary scheme: + └── WENO(order=5) + Symmetric scheme: + └── Centered(order=6) Directions: - ├── X regular - ├── Y regular + ├── X regular + ├── Y regular └── Z stretched ``` """ @@ -105,14 +105,13 @@ function WENO(FT::DataType=Float64; mod(order, 2) == 0 && throw(ArgumentError("WENO reconstruction scheme is defined only for odd orders")) if order < 3 - # WENO(order = 1) is equivalent to UpwindBiased(order = 1) - return UpwindBiased(FT; order = 1) + # WENO(order=1) is equivalent to UpwindBiased(order=1) + return UpwindBiased(FT; order=1) else - N = Int((order + 1) ÷ 2) - + N = Int((order + 1) ÷ 2) weno_coefficients = compute_reconstruction_coefficients(grid, FT, :WENO; order = N) advecting_velocity_scheme = Centered(FT; grid, order = order - 1) - buffer_scheme = WENO(FT; grid, order = order - 2, bounds) + buffer_scheme = WENO(FT; grid, order=order-2, bounds) end return WENO{N, FT}(weno_coefficients..., bounds, buffer_scheme, advecting_velocity_scheme) @@ -127,7 +126,7 @@ WENOFifthOrder(grid=nothing, FT::DataType=Float64; kwargs...) = WENO(grid, FT; # Flavours of WENO const PositiveWENO = WENO{<:Any, <:Any, <:Any, <:Any, <:Any, <:Tuple} -Base.summary(a::WENO{N}) where N = string("WENO reconstruction order ", N*2-1) +Base.summary(a::WENO{N}) where N = string("WENO(order=", N*2-1, ")") Base.show(io::IO, a::WENO{N, FT, RX, RY, RZ, PP}) where {N, FT, RX, RY, RZ, PP} = print(io, summary(a), " \n", diff --git a/src/BuoyancyModels/buoyancy.jl b/src/BuoyancyModels/buoyancy.jl index 3754cb4a4d..188a7876a4 100644 --- a/src/BuoyancyModels/buoyancy.jl +++ b/src/BuoyancyModels/buoyancy.jl @@ -32,7 +32,7 @@ model = NonhydrostaticModel(; grid, buoyancy, tracers=:b) NonhydrostaticModel{CPU, RectilinearGrid}(time = 0 seconds, iteration = 0) ├── grid: 1×8×8 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 1×3×3 halo ├── timestepper: RungeKutta3TimeStepper -├── advection scheme: Centered reconstruction order 2 +├── advection scheme: Centered(order=2) ├── tracers: b ├── closure: Nothing ├── buoyancy: BuoyancyTracer with ĝ = (0.0, -0.707107, -0.707107) diff --git a/src/Models/seawater_density.jl b/src/Models/seawater_density.jl index 6dbfdefb69..c5e723f101 100644 --- a/src/Models/seawater_density.jl +++ b/src/Models/seawater_density.jl @@ -82,7 +82,7 @@ julia> model = NonhydrostaticModel(; grid, buoyancy, tracers) NonhydrostaticModel{CPU, RectilinearGrid}(time = 0 seconds, iteration = 0) ├── grid: 1×1×100 RectilinearGrid{Float64, Flat, Flat, Bounded} on CPU with 0×0×3 halo ├── timestepper: RungeKutta3TimeStepper -├── advection scheme: Centered reconstruction order 2 +├── advection scheme: Centered(order=2) ├── tracers: (T, S) ├── closure: Nothing ├── buoyancy: SeawaterBuoyancy with g=9.80665 and BoussinesqEquationOfState{Float64} with ĝ = NegativeZDirection()