Enforce total buoyancy flux BC in tilted geometry example (#3581)

* Enforce total buoyancy flux BC in tilted geometry example

This example includes no explicit modification to the BCs on buoyancy,
meaning that it defaults to a no-flux BC (or `FluxBoundaryCondition()`)
on buoyancy.

Following equation (6) of Wenegrat and Thomas (2020), we instead enforce
a no-normal diffusive flux boundary condition on the *total* buoyancy, i.e.
the perturbation plus the constantly-stratified `BackgroundField`. Because
Oceananigans.jl does not allow diffusion to act on `BackgroundField`s, the
background part of the no-flux BC is missing. However, we here can enforce it
by specifying a perturbation flux BC that matches the implied background flux.

Other minor changes:
- Added equally-spaced buoyancy surfaces to movie panels
- Guess timestep by minimum of advective and diffusive timescales,
instead of just the advective timescale

* Update examples/tilted_bottom_boundary_layer.jl

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

* Imeplement Greg's suggestions

* Add equations to explain the choice of bottom BC on buoyancy

* More intentional phrasing to explain perturbation BC

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

* More transparent assignment of bottom BC

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

* Rename background field to match mathematical notation

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

* Rename background field to match mathematical notation, follow-up

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

---------

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

examples/tilted_bottom_boundary_layer.jl

@@ -95,9 +95,22 @@ coriolis = ConstantCartesianCoriolis(f = 1e-4, rotation_axis = ĝ)
 # _perturbations_ away from the constant density stratification by imposing
 # a constant stratification as a `BackgroundField`,
 
-B_field = BackgroundField(constant_stratification, parameters=(; ĝ, N² = 1e-5))
+N² = 1e-5 # s⁻² # background vertical buoyancy gradient
+B∞_field = BackgroundField(constant_stratification, parameters=(; ĝ, N² = N²))
 
-# where ``N^2 = 10^{-5} \rm{s}^{-2}`` is the background buoyancy gradient.
+# We choose to impose a bottom boundary condition of zero *total* diffusive buoyancy
+# flux across the seafloor,
+# ```math
+# ∂_z B = ∂_z b + N^{2} \cos{\theta} = 0.
+# ```
+# This shows that to impose a no-flux boundary condition on the total buoyancy field ``B``, we must apply a boundary condition to the perturbation buoyancy ``b``,
+# ```math
+# ∂_z b = - N^{2} \cos{\theta}.
+#```
+
+∂z_b_bottom = - N² * cosd(θ)
+negative_background_diffusive_flux = GradientBoundaryCondition(∂z_b_bottom)
+b_bcs = FieldBoundaryConditions(bottom = negative_background_diffusive_flux)
 
 # ## Bottom drag
 #
@@ -128,14 +141,16 @@ v_bcs = FieldBoundaryConditions(bottom = drag_bc_v)
 # and a constant viscosity and diffusivity. Here we use a smallish value
 # of ``10^{-4} \, \rm{m}^2\, \rm{s}^{-1}``.
 
-closure = ScalarDiffusivity(ν=1e-4, κ=1e-4)
+ν = 1e-4
+κ = 1e-4
+closure = ScalarDiffusivity(ν=ν, κ=κ)
 
 model = NonhydrostaticModel(; grid, buoyancy, coriolis, closure,
                             timestepper = :RungeKutta3,
                             advection = UpwindBiasedFifthOrder(),
                             tracers = :b,
-                            boundary_conditions = (u=u_bcs, v=v_bcs),
-                            background_fields = (; b=B_field))
+                            boundary_conditions = (u=u_bcs, v=v_bcs, b=b_bcs),
+                            background_fields = (; b=B∞_field))
 
 # Let's introduce a bit of random noise at the bottom of the domain to speed up the onset of
 # turbulence:
@@ -146,9 +161,11 @@ set!(model, u=noise, w=noise)
 # ## Create and run a simulation
 #
 # We are now ready to create the simulation. We begin by setting the initial time step
-# conservatively, based on the smallest grid size of our domain.
+# conservatively, based on the smallest grid size of our domain and either an advective
+# or diffusive time scaling, depending on which is shorter.
 
-simulation = Simulation(model, Δt = 0.5 * minimum_zspacing(grid) / V∞, stop_time = 1day)
+Δt₀ = 0.5 * minimum([minimum_zspacing(grid) / V∞, minimum_zspacing(grid)^2/κ])
+simulation = Simulation(model, Δt = Δt₀, stop_time = 1day)
 
 # We use a `TimeStepWizard` to adapt our time-step,
 
@@ -199,6 +216,7 @@ run!(simulation)
 
 using NCDatasets, CairoMakie
 
+xb, yb, zb = nodes(B)
 xω, yω, zω = nodes(ωy)
 xv, yv, zv = nodes(V)
 
@@ -220,14 +238,17 @@ ax_v = Axis(fig[3, 1]; title = "Along-slope velocity (v)", axis_kwargs...)
 n = Observable(1)
 
 ωy = @lift ds["ωy"][:, 1, :, $n]
+B = @lift ds["B"][:, 1, :, $n]
 hm_ω = heatmap!(ax_ω, xω, zω, ωy, colorrange = (-0.015, +0.015), colormap = :balance)
 Colorbar(fig[2, 2], hm_ω; label = "s⁻¹")
+ct_b = contour!(ax_ω, xb, zb, B, levels=-1e-3:0.5e-4:1e-3, color=:black)
 
 V = @lift ds["V"][:, 1, :, $n]
 V_max = @lift maximum(abs, ds["V"][:, 1, :, $n])
 
 hm_v = heatmap!(ax_v, xv, zv, V, colorrange = (-V∞, +V∞), colormap = :balance)
 Colorbar(fig[3, 2], hm_v; label = "m s⁻¹")
+ct_b = contour!(ax_v, xb, zb, B, levels=-1e-3:0.5e-4:1e-3, color=:black)
 
 times = collect(ds["time"])
 title = @lift "t = " * string(prettytime(times[$n]))