Dont adapt Flat direction advection scheme plus better summary for FluxFormAdvection (#3908)

* Dont adapt Flat directions plus better summary for FluxFormAdvection

* Update docs

---------

Co-authored-by: Navid C. Constantinou <[email protected]>

Author: glwagner

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

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)

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

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

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,23 +76,26 @@ 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
 
 function adapt_advection_order(advection::UpwindBiased{B}, N::Int, grid::AbstractGrid) where B
     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
+end
+

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",

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))

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",

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

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",

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)

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()