Rayleigh-Taylor instability

This example simulates a three-dimensional Rayleigh-Taylor instability: a heavy fluid initially resting above a light one, which is unstable under gravity. All of the kinetic energy is released from the potential energy stored in the unstable stratification, with no fluxes through the boundaries. We run it as a large-eddy simulation (LES) and use Oceanostics to close the volume-integrated sub-filter-scale kinetic-energy budget, that is, the budget of the kinetic energy carried by the scales that a low-pass filter removes from the flow.

Before starting, make sure you have the required packages installed for this example, which can be done with

using Pkg
pkg"add Oceananigans, Oceanostics, CairoMakie, NCDatasets"

Parameters and grid

using Oceananigans

We work with nondimensional quantities. We take the domain height H as the length scale and the buoyancy jump Δb across the initial interface as the buoyancy scale, so the free-fall velocity U = √(Δb H) and the free-fall time τ = √(H / Δb) follow. Setting H = Δb = 1 makes U = τ = 1, and time is measured in free-fall units:

H  = 1     # domain height (length scale)
Δb = 1     # buoyancy jump across the interface (buoyancy scale)
U  = sqrt(Δb * H)   # free-fall velocity scale
τ  = sqrt(H / Δb)   # free-fall time scale

1.0

We use a cuboid that is periodic in the two horizontal directions and bounded in the vertical, with the unstable interface at mid-depth z = 0:

N = 48
grid = RectilinearGrid(size=(N, N÷2, N), x=(-H/2, H/2), y=(-H/4, H/4), z=(-H/2, H/2),
                       topology=(Periodic, Periodic, Bounded))

48×24×48 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo
├── Periodic x ∈ [-0.5, 0.5)   regularly spaced with Δx=0.0208333
├── Periodic y ∈ [-0.25, 0.25) regularly spaced with Δy=0.0208333
└── Bounded  z ∈ [-0.5, 0.5]   regularly spaced with Δz=0.0208333

Closure and model

We model the stresses with a Large Eddy Simulation (LES) closure:

closure = SmagorinskyLilly(C=0.3) # A large C increases eddy viscosity, keeping this very coarse example well behaved

Smagorinsky closure with
├── coefficient = LillyCoefficient(smagorinsky = 0.3, reduction_factor = 1.0)
└── Pr = 1.0

We build a NonhydrostaticModel with a fourth-order centered advection scheme, a third-order Runge-Kutta timestepper, and a buoyancy b as the active tracer. A centered scheme is non-dissipative and adds no numerical dissipation of its own, so essentially all of the dissipation comes from the LES closure.

model = NonhydrostaticModel(grid; timestepper = :RungeKutta3,
                            advection = Centered(order=4),
                            closure,
                            buoyancy = BuoyancyTracer(), tracers = :b)

NonhydrostaticModel{CPU, RectilinearGrid}(time = 0 seconds, iteration = 0)
├── grid: 48×24×48 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo
├── timestepper: RungeKutta3TimeStepper
├── advection scheme: Centered(order=4)
├── tracers: b
├── closure: Smagorinsky with coefficient = LillyCoefficient(smagorinsky = 0.3, reduction_factor = 1.0), Pr=(b = 1.0,)
├── buoyancy: BuoyancyTracer with ĝ = NegativeZDirection()
└── coriolis: Nothing

Initial condition

The initial buoyancy is a hyperbolic-tangent profile that is heavy on top: it decreases from +Δb/2 at the bottom to −Δb/2 at the top, so ∂b/∂z < 0 and the stratification is unstable. The interface is thin (its half-thickness δ is about one grid spacing) and we perturb it with small-amplitude random noise localized to the interface, which seeds a broad band of horizontal wavelengths and produces a multi-mode, turbulent instability rather than a single growing bubble. We seed the random number generator so the run is reproducible:

δ = 0.02 * H                    # interface half-thickness
b₀(z) = -(Δb / 2) * tanh(z / δ) # +Δb/2 at the bottom, −Δb/2 at the top

using Random
Random.seed!(43)
noise_amplitude = 1e-2
bᵢ(x, y, z) = b₀(z) + noise_amplitude * Δb * randn() * exp(-(z / δ)^2)

set!(model, b=bᵢ)

Simulation

We create an adaptive-time-step simulation. The initial step is set conservatively from the horizontal grid spacing and the free-fall velocity, and a TimeStepWizard adapts it as the mixing layer accelerates:

Δx = minimum_xspacing(grid)
simulation = Simulation(model, Δt = 0.1 * Δx / U, stop_time = 7τ)
conjure_time_step_wizard!(simulation, IterationInterval(2), cfl=0.8, max_change=1.1)

Model diagnostics

We report progress with the TimedMessenger

using Oceanostics
progress = ProgressMessengers.TimedMessenger()
add_callback!(simulation, progress, IterationInterval(100))

Sub-filter-scale kinetic-energy budget

Rayleigh-Taylor turbulence converts potential energy into kinetic energy (KE) across a wide range of scales, so we follow it with a coarse-graining analysis in the spirit of Aluie et al. (2018). A low-pass filter (overbar) splits each field into a filtered part and a sub-filter remainder. Here we budget the kinetic energy carried by the scales the filter removes,

\[K^s = \tfrac{1}{2}\,\tau^{ii} , \qquad \tau^{ij} = \overline{u^i u^j} - \overline{u}^i\,\overline{u}^j ,\]

where $\tau^{ij}$ is the sub-filter stress (subfilter_stress_tensor), so $K^s$ itself is computed by SubFilterKineticEnergy. We apply a Gaussian filter of width ℓ in the two horizontal directions, which are statistically homogeneous; the vertical direction is left unfiltered.

Volume integrated (the transport terms integrate to zero, because the flow is doubly periodic and w = 0 with free slip at the top and bottom), the budget reads

\[\frac{d}{dt} \int K^s\, dV = \int \Pi_K\, dV + \int \tau(w, b)\, dV - \int \varepsilon_K^s\, dV ,\]

with two sources and one sink:

• $\Pi_K = -\tau^{ij}\overline{S}^{ij}$ is the cross-scale kinetic-energy flux (KineticEnergyCrossScaleFlux), the rate at which the filtered scales hand kinetic energy down to the sub-filter scales.
• $\tau(w, b) = \overline{wb} - \overline{w}\,\overline{b}$ is the sub-filter buoyancy flux (a subfilter_covariance), which converts sub-filter potential energy into sub-filter kinetic energy.
• $\varepsilon_K^s = \overline{\varepsilon} - \varepsilon^{\ell}$ is the sub-filter dissipation (SubFilterKineticEnergyDissipationRate): the filtered total dissipation $\varepsilon$ (KineticEnergyDissipationRate) minus the dissipation $\varepsilon^{\ell}$ of the filtered flow (FilteredKineticEnergyDissipationRate). For a constant viscosity it reduces to $2\nu[\overline{S^{ij}S^{ij}} - \overline{S}^{ij}\overline{S}^{ij}] \ge 0$, a strictly positive sink; with an LES closure it is the dissipation that the modeled stress carries out on the sub-filter scales.

In the code below, $K^s$, $\Pi_K$, $\tau(w,b)$ and $\varepsilon_K^s$ are written Kˢ, Πₖ, wbˢ and εˢ.

using Oceananigans.AbstractOperations: @at

ℓ  = 8 * Δx                          # filter scale (full width at half maximum of the Gaussian kernel)
σℓ = ℓ / (2 * sqrt(2 * log(2)))      # corresponding Gaussian standard deviation
gfilter = GaussianFilter(dims=(1, 2), σ=σℓ)

u, v, w = model.velocities
b = model.tracers.b

Kˢ  = SubFilterKineticEnergy(model, gfilter)     # sub-filter kinetic energy ½τⁱⁱ
Πₖ  = KineticEnergyCrossScaleFlux(model, gfilter)  # cross-scale flux from filtered scales
wbˢ = subfilter_covariance(w, b, gfilter)          # sub-filter buoyancy flux τ(w, b)
εˢ  = SubFilterKineticEnergyDissipationRate(model, gfilter)  # sub-filter dissipation ε̄ − εˡ

SubFilterKineticEnergyDissipationRate KernelFunctionOperation at (Center, Center, Center)
├── grid: 48×24×48 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo
├── kernel_function: subfilter_ke_dissipation_rate_ccc (generic function with 1 method)
└── arguments: ("Oceananigans.AbstractOperations.BinaryOperation",)
└── computes: sub-filter kinetic energy dissipation rate  εˢ = filter(ε) - εˡ

The budget needs only the (cheap) volume integrals of these terms:

∫Kˢ  = Integral(Kˢ)
∫Πₖ  = Integral(Πₖ)
∫wbˢ = Integral(wbˢ)
∫εˢ  = Integral(εˢ)

Integral of BinaryOperation at (Center, Center, Center) over dims (1, 2, 3)
└── operand: BinaryOperation at (Center, Center, Center)
    └── grid: 48×24×48 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo

For the movie we also keep the coarse-grained kinetic energy $\overline{K} = \tfrac{1}{2}\,\overline{u}_i\overline{u}_i$ (FilteredKineticEnergy), the filtered counterpart of $K^s$. Together the two show how the filter splits the flow's kinetic energy between the scales it keeps and the scales it removes:

# `FilteredKineticEnergy` materializes the filtered velocities internally, so the multi-direction filter
# runs on its fast staged path (see the filter performance notes and `check_filter_staging`).
Kˡ = FilteredKineticEnergy(model, gfilter)  # kinetic energy of the coarse-grained (filtered) flow

FilteredKineticEnergy KernelFunctionOperation at (Center, Center, Center)
├── grid: 48×24×48 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo
├── kernel_function: filtered_kinetic_energy_ccc (generic function with 1 method)
└── arguments: ("Field", "Field", "Field")
└── computes: kinetic energy of the filtered flow  Kˡ = ½ūᵢūᵢ

Output

We use two NetCDF writers. A snapshot writer stores vertical (x–z) slices of the buoyancy b, the cross-scale flux Πₖ and the two kinetic energies Kˡ and Kˢ, at a fixed y index (the flow is periodic and statistically homogeneous in y, so the particular plane makes no difference), while a budget writer stores only the integrated scalars on ConsecutiveIterations(TimeInterval(τ/5)), which takes a second sample one model step after each output time. That lets us finite-difference ∫Kˢ across that single step to estimate d/dt, exactly as in the Kelvin-Helmholtz example.

using NCDatasets
filename = joinpath(@__DIR__, "rayleigh_taylor_instability")

simulation.output_writers[:fields] = NetCDFWriter(model, (; b, Πₖ, Kˡ, Kˢ),
                                                  filename = filename,
                                                  schedule = TimeInterval(τ / 5),
                                                  indices = (:, 1, :),
                                                  overwrite_existing = true)

simulation.output_writers[:budget] = NetCDFWriter(model, (; ∫Kˢ, ∫Πₖ, ∫wbˢ, ∫εˢ),
                                                  filename = filename * "_budget",
                                                  schedule = ConsecutiveIterations(TimeInterval(τ / 5)),
                                                  overwrite_existing = true)

NetCDFWriter scheduled on ConsecutiveIterations(TimeInterval(200 ms), 1):
├── filepath: rayleigh_taylor_instability_budget.nc
├── dimensions: time(0), y_afa(24), x_faa(48), x_caa(48), y_aca(24), z_aaf(49), z_aac(48)
├── 4 outputs: (∫wbˢ, ∫Πₖ, ∫Kˢ, ∫εˢ)
├── array_type: Array{Float32}
├── file_splitting: NoFileSplitting
└── file size: 33.7 KiB

Run the simulation and process results

To run the simulation:

run!(simulation)

[ Info: Initializing simulation...
┌ Info: iter =      0,  [000.00%] time = 0 seconds,  Δt = 2.292 ms,  walltime = 1.809 minutes,  walltime / timestep = 0 seconds
└       |u⃗|ₘₐₓ = [0.00e+00,  0.00e+00,  0.00e+00] m/s,  advective CFL = 0,  diffusive CFL = 0,  νₘₐₓ = 0 m²/s
[ Info:     ... simulation initialization complete (25.107 seconds)
[ Info: Executing initial time step...
[ Info:     ... initial time step complete (3.258 seconds).
┌ Info: iter =    100,  [040.00%] time = 2.800 seconds,  Δt = 45.487 ms,  walltime = 2.671 minutes,  walltime / timestep = 516.956 ms
└       |u⃗|ₘₐₓ = [1.01e-01,  9.81e-02,  2.86e-01] m/s,  advective CFL = 0.8,  diffusive CFL = 0.051,  νₘₐₓ = 0.00049 m²/s
┌ Info: iter =    200,  [076.16%] time = 5.331 seconds,  Δt = 18.767 ms,  walltime = 3.081 minutes,  walltime / timestep = 246.093 ms
└       |u⃗|ₘₐₓ = [3.39e-01,  3.42e-01,  7.03e-01] m/s,  advective CFL = 0.8,  diffusive CFL = 0.037,  νₘₐₓ = 0.00086 m²/s
┌ Info: iter =    300,  [099.53%] time = 6.967 seconds,  Δt = 14.831 ms,  walltime = 3.360 minutes,  walltime / timestep = 167.657 ms
└       |u⃗|ₘₐₓ = [4.90e-01,  4.62e-01,  7.47e-01] m/s,  advective CFL = 0.8,  diffusive CFL = 0.03,  νₘₐₓ = 0.00089 m²/s
[ Info: Simulation is stopping after running for 1.584 minutes.
[ Info: Simulation time 7 seconds equals or exceeds stop time 7 seconds.

We read the four field slices, and their x–z coordinates, back with NCDataset (all of them live at cell centers, so their coordinates are x_caa and z_aac); the singleton y dimension of the slice is dropped:

using CairoMakie

ds = NCDataset(simulation.output_writers[:fields].filepath)
times = ds["time"][:]
x_caa = ds["x_caa"][:]
z_aac = ds["z_aac"][:]
b_arr = ds["b"][:, 1, :, :]
Π_arr = ds["Πₖ"][:, 1, :, :]
Kˡ_arr = ds["Kˡ"][:, 1, :, :]
Kˢ_arr = ds["Kˢ"][:, 1, :, :]
close(ds)

closed Dataset

The integrated budget scalars come in consecutive-iteration pairs (2k-1, 2k); a one-step finite difference inside each pair gives d(∫Kˢ)/dt, and each budget term is evaluated at the pair midpoint. The residual measures how well the sub-filter-scale budget closes.

bud_filepath = simulation.output_writers[:budget].filepath
ds_bud = NCDataset(bud_filepath)
times_bud = ds_bud["time"][:]
∫Kˢ_t = ds_bud["∫Kˢ"][:]
∫Πₖ_t = ds_bud["∫Πₖ"][:]
∫wbˢ_t = ds_bud["∫wbˢ"][:]
∫εˢ_t = ds_bud["∫εˢ"][:]
close(ds_bud)

i1 = 1:2:length(times_bud)-1   # primary snapshots
i2 = 2:2:length(times_bud)     # consecutive-iteration snapshots
Δt_pair = times_bud[i2] .- times_bud[i1]
t_pair = @. 0.5 * (times_bud[i1] + times_bud[i2])

dKˢdt    = (∫Kˢ_t[i2] .- ∫Kˢ_t[i1]) ./ Δt_pair
Πₖ_pair  = @. 0.5 * (∫Πₖ_t[i1] + ∫Πₖ_t[i2])
wbˢ_pair = @. 0.5 * (∫wbˢ_t[i1] + ∫wbˢ_t[i2])
εˢ_pair  = @. 0.5 * (∫εˢ_t[i1] + ∫εˢ_t[i2])

35-element Vector{Float64}:
 5.95359235713432e-16
 1.087107737163251e-9
 1.6115127010607466e-8
 1.3989108538225992e-7
 1.0053045116364956e-6
 5.913985660299659e-6
 2.537825275794603e-5
 7.691980863455683e-5
 0.0001773710537236184
 0.0003089088131673634
 0.0004670735215768218
 0.0006454169051721692
 0.0008495457004755735
 0.001091537531465292
 0.0013851976254954934
 0.0017412814777344465
 0.002166382735595107
 0.002652437426149845
 0.0031752674840390682
 0.0037090927362442017
 0.004239959642291069
 0.004766242112964392
 0.005304264836013317
 0.005912557244300842
 0.006579696200788021
 0.007150683552026749
 0.007710915058851242
 0.0084513109177351
 0.009435515850782394
 0.0103079192340374
 0.01102607510983944
 0.011534986086189747
 0.011679408140480518
 0.011499028652906418
 0.011203847825527191

Residual in sum-to-zero form: the negative tendency plus the sources, so the plotted curves add to it

resid = @. -dKˢdt + Πₖ_pair + wbˢ_pair - εˢ_pair

Plotting

We build the figure in three rows: the vertical slices of the buoyancy b (the spikes and bubbles) and of the cross-scale flux Πₖ on top; the two kinetic energies that the filter separates, Kˡ and Kˢ, in the middle; and the volume-integrated sub-filter-scale kinetic-energy budget at the bottom.

set_theme!(Theme(fontsize=18))
fig = Figure(size=(1000, 1200))

n = Observable(1)

axb = Axis(fig[2, 1]; title="buoyancy, b", xlabel="x", ylabel="z", aspect=1)
axΠ = Axis(fig[2, 3]; title="cross-scale KE flux, Πₖ", xlabel="x", ylabel="z", aspect=1)

blim = 0.8 * maximum(abs, b_arr[:, :, 1])
Πlim = 0.5 * maximum(abs, Π_arr)

bₙ = @lift b_arr[:, :, $n]
Πₙ = @lift Π_arr[:, :, $n]

hmb = heatmap!(axb, x_caa, z_aac, bₙ; colormap=:balance, colorrange=(-blim, blim))
Colorbar(fig[2, 2], hmb)

hmΠ = heatmap!(axΠ, x_caa, z_aac, Πₙ; colormap=:balance, colorrange=(-Πlim, Πlim))
Colorbar(fig[2, 4], hmΠ)

Makie.Colorbar()

The middle row splits the kinetic energy across the filter scale: on the left the coarse-grained (filtered) energy Kˡ, on the right the sub-filter energy Kˢ whose budget the bottom panel closes. Both are non-negative, so they get a sequential colormap, and each gets its own colour scale because the two differ by orders of magnitude.

axKˡ = Axis(fig[3, 1]; title="coarse-grained KE, Kˡ", xlabel="x", ylabel="z", aspect=1)
axKˢ = Axis(fig[3, 3]; title="sub-filter-scale KE, Kˢ", xlabel="x", ylabel="z", aspect=1)

Kˡlim = maximum(Kˡ_arr)
Kˢlim = maximum(Kˢ_arr)

Kˡₙ = @lift Kˡ_arr[:, :, $n]
Kˢₙ = @lift Kˢ_arr[:, :, $n]

hmKˡ = heatmap!(axKˡ, x_caa, z_aac, Kˡₙ; colormap=:magma, colorrange=(0, Kˡlim))
Colorbar(fig[3, 2], hmKˡ)

hmKˢ = heatmap!(axKˢ, x_caa, z_aac, Kˢₙ; colormap=:magma, colorrange=(0, Kˢlim))
Colorbar(fig[3, 4], hmKˢ)

Makie.Colorbar()

The bottom panel shows the volume-integrated sub-filter-scale kinetic-energy budget. We plot the negative tendency −d(∫Kˢ)/dt together with the two sources that feed it, the cross-scale flux ∫Πₖ dV handed down from the filtered scales and the sub-filter buoyancy flux ∫τ(w,b) dV, and the single sink that drains it, the sub-filter dissipation −∫εˢ dV. With the tendency negated, the four curves sum to the residual.

ax_bud = Axis(fig[4, 1:4]; xlabel="time [free-fall units]", title="Sub-filter-scale KE budget")
lines!(ax_bud, t_pair ./ τ, -dKˢdt,   label="−d(∫Kˢ)/dt")
lines!(ax_bud, t_pair ./ τ, Πₖ_pair,  label="∫Πₖ dV  (flux from filtered scales)")
lines!(ax_bud, t_pair ./ τ, wbˢ_pair, label="∫τ(w,b) dV  (sub-filter buoyancy flux)")
lines!(ax_bud, t_pair ./ τ, -εˢ_pair, label="−∫εˢ dV  (sub-filter dissipation)")
lines!(ax_bud, t_pair ./ τ, resid,    label="residual", color=:black, linestyle=:dash)
axislegend(ax_bud; position=:lt, labelsize=10)

vlines!(ax_bud, @lift(times[$n] / τ), color=:black, linestyle=:dash)

title = @lift "Rayleigh-Taylor instability, t = " * string(round(times[$n] / τ, digits=2)) * " τ"
fig[1, 1:4] = Label(fig, title, fontsize=22, tellwidth=false)

@info "Animating..."
record(fig, "rayleigh_taylor_instability.mp4", 1:length(times), framerate=12) do i
    n[] = i
end

"rayleigh_taylor_instability.mp4"

As the heavy fluid falls in spikes and the light fluid rises in bubbles, the flow rolls up and breaks into a turbulent mixing layer. The bottom panel shows the volume-integrated SFS KE budget. The sub-filter buoyancy flux ∫τ(w,b) dV and dissipation ∫εˢ dV are the dominant terms, with the tendency −d(∫Kˢ)/dt in third place. The residual is small, which shows that the budget closes well and that the coarse-graining analysis is consistent with the simulation.

This page was generated using Literate.jl.