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 subfilter-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 Pkgpkg"add Oceananigans, Oceanostics, CairoMakie, NCDatasets"Parameters and grid
using OceananigansWe 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 scale1.0We 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 = 48grid = 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.0208333Closure 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 behavedSmagorinsky closure with
├── coefficient = LillyCoefficient(smagorinsky = 0.3, reduction_factor = 1.0)
└── Pr = 1.0We 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:
│ ├── momentum: Centered(order=4)
│ └── b: 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: NothingInitial 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-thicknessb₀(z) = -(Δb / 2) * tanh(z / δ) # +Δb/2 at the bottom, −Δb/2 at the topusing RandomRandom.seed!(43)noise_amplitude = 1e-2bᵢ(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 Oceanosticsprogress = ProgressMessengers.TimedMessenger()add_callback!(simulation, progress, IterationInterval(100))Subfilter-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 filtering analysis in the spirit of Aluie et al. (2018). A low-pass filter (overbar) splits each field into a filtered part and a subfilter remainder. Here we budget the kinetic energy carried by the scales the filter removes,
\[e_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 subfilter stress (subfilter_stress_tensor), so $e_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 e_k^s\, \mathrm{d}V = \int \Pi_k\, \mathrm{d}V + \int \tau(w, b)\, \mathrm{d}V - \int \varepsilon_k^s\, \mathrm{d}V ,\]
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 subfilter scales. - $\tau(w, b) = \overline{wb} - \overline{w}\,\overline{b}$ is the subfilter buoyancy flux (a
subfilter_covariance), which converts subfilter potential energy into subfilter kinetic energy. The filter here is purely horizontal, so this equals the $\tau^l(w, b_r)$ of the subfilter kinetic energy budget: the reference profile $b^\star(z)$ is constant along the filtered directions and drops out of the covariance. - $\varepsilon_k^s = \overline{\varepsilon_k} - \varepsilon_k^l$ is the subfilter dissipation (
SubFilterKineticEnergyDissipationRate): the filtered total dissipation $\varepsilon_k$ (KineticEnergyDissipationRate) minus the dissipation $\varepsilon_k^l$ 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 subfilter scales.
In the code below, $e_k^s$, $\Pi_k$, $\tau(w,b)$ and $\varepsilon_k^s$ are written eₖˢ, Πₖ, 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 deviationgfilter = GaussianFilter(dims=(1, 2), σ=σℓ)u, v, w = model.velocitiesb = model.tracers.beₖˢ = SubFilterKineticEnergy(model, gfilter) # subfilter kinetic energy ½τᵢᵢΠₖ = KineticEnergyCrossScaleFlux(model, gfilter) # cross-scale flux from filtered scaleswbˢ = subfilter_covariance(w, b, gfilter) # subfilter buoyancy flux τ(w, b)εₖˢ = SubFilterKineticEnergyDissipationRate(model, gfilter) # subfilter 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: subfilter kinetic energy dissipation rate filter(εₖ) - εₖˡThe budget needs only the (cheap) volume integrals of these terms:
∫eₖˢ = Integral(eₖˢ)∫Πₖ = 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 haloFor the movie we also keep the filtered kinetic energy $e_k^l = \tfrac{1}{2}\,\overline{u}_i\overline{u}_i$ (FilteredKineticEnergy), the filtered counterpart of $e_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`).eₖˡ = FilteredKineticEnergy(model, gfilter) # kinetic energy of the filtered flowFilteredKineticEnergy 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 ½ūᵢūᵢ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 eₖˡ and eₖˢ, 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 ∫eₖˢ across that single step to estimate d/dt, exactly as in the Kelvin-Helmholtz example.
using NCDatasetsfilename = joinpath(@__DIR__, "rayleigh_taylor_instability")simulation.output_writers[:fields] = NetCDFWriter(model, (; b, Πₖ, eₖˡ, eₖˢ), filename = filename, schedule = TimeInterval(τ / 5), indices = (:, 1, :), overwrite_existing = true)simulation.output_writers[:budget] = NetCDFWriter(model, (; ∫eₖˢ, ∫Πₖ, ∫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ˢ, ∫Πₖ, ∫eₖˢ, ∫εₖˢ)
├── array_type: Array{Float32}
├── file_splitting: NoFileSplitting
└── file size: (file not yet created)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 = 2.887 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 (1.058 minutes)
[ Info: Executing initial time step...
[ Info: ... initial time step complete (3.399 seconds).
┌ Info: iter = 100, [040.00%] time = 2.800 seconds, Δt = 45.487 ms, walltime = 3.838 minutes, walltime / timestep = 570.482 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 = 4.363 minutes, walltime / timestep = 314.800 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 = 4.714 minutes, walltime / timestep = 210.620 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 2.525 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 CairoMakieds = 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, :, :]eₖˡ_arr = ds["eₖˡ"][:, 1, :, :]eₖˢ_arr = ds["eₖˢ"][:, 1, :, :]close(ds)closed DatasetThe integrated budget scalars come in consecutive-iteration pairs (2k-1, 2k); a one-step finite difference inside each pair gives d(∫eₖˢ)/dt, and each budget term is evaluated at the pair midpoint. The residual measures how well the subfilter-scale budget closes.
bud_filepath = simulation.output_writers[:budget].filepathds_bud = NCDataset(bud_filepath)times_bud = ds_bud["time"][:]∫eₖˢ_t = ds_bud["∫eₖˢ"][:]∫Πₖ_t = ds_bud["∫Πₖ"][:]∫wbˢ_t = ds_bud["∫wbˢ"][:]∫εₖˢ_t = ds_bud["∫εₖˢ"][:]close(ds_bud)i1 = 1:2:length(times_bud)-1 # primary snapshotsi2 = 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])deₖˢdt = (∫eₖˢ_t[i2] .- ∫eₖˢ_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]);Residual in sum-to-zero form: the negative tendency plus the sources, so the plotted curves add to it
resid = @. -deₖˢdt + Πₖ_pair + wbˢ_pair - εₖˢ_pairPlotting
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, eₖˡ and eₖˢ, in the middle; and the volume-integrated subfilter-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 kinetic energy 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 filtered energy eₖˡ, on the right the subfilter energy eₖˢ 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.
axeₖˡ = Axis(fig[3, 1]; title="filtered kinetic energy, eₖˡ", xlabel="x", ylabel="z", aspect=1)axeₖˢ = Axis(fig[3, 3]; title="subfilter kinetic energy, eₖˢ", xlabel="x", ylabel="z", aspect=1)eₖˡlim = maximum(eₖˡ_arr)eₖˢlim = maximum(eₖˢ_arr)eₖˡₙ = @lift eₖˡ_arr[:, :, $n]eₖˢₙ = @lift eₖˢ_arr[:, :, $n]hmeₖˡ = heatmap!(axeₖˡ, x_caa, z_aac, eₖˡₙ; colormap=:magma, colorrange=(0, eₖˡlim))Colorbar(fig[3, 2], hmeₖˡ)hmeₖˢ = heatmap!(axeₖˢ, x_caa, z_aac, eₖˢₙ; colormap=:magma, colorrange=(0, eₖˢlim))Colorbar(fig[3, 4], hmeₖˢ)Makie.Colorbar()The bottom panel shows the volume-integrated subfilter-scale kinetic-energy budget. We plot the negative tendency −d(∫eₖˢ)/dt together with the two sources that feed it, the cross-scale flux ∫Πₖ dV handed down from the filtered scales and the subfilter buoyancy flux ∫τ(w,b) dV, and the single sink that drains it, the subfilter 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="Subfilter kinetic energy budget")lines!(ax_bud, t_pair ./ τ, -deₖˢdt, label="−d(∫eₖˢ)/dt")lines!(ax_bud, t_pair ./ τ, Πₖ_pair, label="∫Πₖ dV (flux from filtered scales)")lines!(ax_bud, t_pair ./ τ, wbˢ_pair, label="∫τ(w,b) dV (subfilter buoyancy flux)")lines!(ax_bud, t_pair ./ τ, -εₖˢ_pair, label="−∫εₖˢ dV (subfilter 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[] = iend"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 subfilter buoyancy flux ∫τ(w,b) dV and dissipation ∫εₖˢ dV are the dominant terms, with the tendency −d(∫eₖˢ)/dt in third place. The residual is small, which shows that the budget closes well and that the filtering analysis is consistent with the simulation.
This page was generated using Literate.jl.