From fb2c670626e38e3bffe298c485f95625fb1d83be Mon Sep 17 00:00:00 2001 From: Henri Drake Date: Wed, 8 May 2024 14:23:03 -0700 Subject: [PATCH] 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 * 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 * More transparent assignment of bottom BC Co-authored-by: Gregory L. Wagner * Rename background field to match mathematical notation Co-authored-by: Gregory L. Wagner * Rename background field to match mathematical notation, follow-up Co-authored-by: Gregory L. Wagner --------- Co-authored-by: Henri Drake Co-authored-by: Gregory L. Wagner Co-authored-by: Gregory L. Wagner --- examples/tilted_bottom_boundary_layer.jl | 35 +++++++++++++++++++----- 1 file changed, 28 insertions(+), 7 deletions(-) diff --git a/examples/tilted_bottom_boundary_layer.jl b/examples/tilted_bottom_boundary_layer.jl index 62a24ababb..8211a8a6c2 100644 --- a/examples/tilted_bottom_boundary_layer.jl +++ b/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]))