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, and the budget of the available potential energy carried by the same scales.

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 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 = 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.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-resolved

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:
│   ├── 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 ĝ = 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-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 = 8τ)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(u_i, u_i) , \qquad \tau(u_i, u_j) = \overline{u_i u_j} - \overline{u}_i\,\overline{u}_j ,\]

where $\tau(u_i, u_j)$ 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 \tau(w, b_r)\, \mathrm{d}V + \int \Pi_k\, \mathrm{d}V - \int \varepsilon_k^s\, \mathrm{d}V ,\]

with two sources and one sink:

$\tau(w, b_r) = \overline{w b_r} - \overline{w}\,\overline{b_r}$ is the subfilter buoyancy flux, which converts subfilter available potential energy into subfilter kinetic energy. Here $b_r = b - b^\star(z)$ is the buoyancy measured from the sorted reference profile $b^\star(z)$ of the APE budget below, which shares this term and computes it (SubFilterAvailablePotentialToKineticEnergyConversion). The filter is purely horizontal, so $b^\star(z)$ is constant along the filtered directions and drops out of the covariance: $\tau(w, b_r) = \tau(w, b)$.
$\Pi_k = -\tau(u_i, u_j)\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.
$\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$ and $\varepsilon_k^s$ are written eₖˢ, Πₖ and εₖˢ.

ℓ  = 0.2                         # filter scale (full width at half maximum of the Gaussian kernel)σℓ = ℓ / (2 * sqrt(2 * log(2)))  # corresponding Gaussian standard deviationgfilter = GaussianFilter(dims=(1, 2), σ=σℓ)b = model.tracers.beₖˢ = SubFilterKineticEnergy(model, gfilter)       # subfilter kinetic energy ½τᵢᵢΠₖ  = KineticEnergyCrossScaleFlux(model, gfilter)  # cross-scale flux from filtered scalesεₖˢ = 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: ("BinaryOperation",)
└── computes: subfilter kinetic energy dissipation rate  filter(εₖ) - εₖˡ

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

∫eₖˢ = Integral(eₖˢ)∫Πₖ  = Integral(Πₖ)∫εₖˢ = Integral(εₖˢ)∂ₜ∫eₖˢ = TimeDerivative(∫eₖˢ)

TimeDerivative of 1×1×1 Field{Nothing, Nothing, Nothing} reduced over dims = (1, 2, 3) on RectilinearGrid on CPU

Subfilter available potential energy budget

The same filter splits the local available potential energy (APE) $e_a$ of the Available potential energy equation, measured against the sorted reference profile $b^\star(z, t)$ of the full buoyancy, into the APE of the filtered buoyancy, $e_a^l$, and the subfilter remainder $e_a^s = \overline{e_a} - e_a^l$ (SubFilterAvailablePotentialEnergy), which a horizontal filter keeps non-negative. Volume integrated, its transport terms vanish as well, since the top and bottom are also insulating, and the budget reads

\[\frac{d}{dt} \int e_a^s\, \mathrm{d}V = - \int \tau(w, b_r)\, \mathrm{d}V + \int \Pi_a\, \mathrm{d}V - \int \varepsilon_a^s\, \mathrm{d}V + \int R^s\, \mathrm{d}V ,\]

where

$\tau(w, b_r) = \overline{w b_r} - \overline{w}\, \overline{b_r}$ converts subfilter APE into subfilter kinetic energy (SubFilterAvailablePotentialToKineticEnergyConversion); it is the conversion term of the KE budget above, with the opposite sign.
$\Pi_a$ is the cross-scale APE flux (AvailablePotentialEnergyCrossScaleFlux), the rate at which the filtered scales hand APE down to the subfilter scales.
$\varepsilon_a^s = \overline{\varepsilon_a} - \varepsilon_a^l$ is the subfilter APE dissipation (SubFilterAvailablePotentialEnergyDissipationRate), the APE that the closure's buoyancy diffusion removes from the subfilter scales.
$R^s = \overline{R} - R^l$ is the subfilter reference tendency (SubFilterAvailablePotentialEnergyReferenceTendency). Mixing reshapes the sorted profile, which changes the APE of every parcel at a rate $R$. That rate integrates to zero over the domain while its filtered counterpart $R^l$ does not, so $R^s$ moves APE between the filtered and the subfilter scales.

In the code below these are written eₐˢ, wbᵣˢ, Πₐ, εₐˢ and Rˢ. All of them are measured against one sorted column of the buoyancy, which ProfileLookup shares between them, and $R^s$ integrates the time derivative of its profile. Reading a TimeDerivative does not update it, so a callback does, on the iteration before each output and at the output, the same step over which the budget writer differences $\int e_a^s\, \mathrm{d}V$:

z✶_column = reference_height(model, method=VerticalSort())lookup    = ProfileLookup(z✶_column)∂ₜb✶ = TimeDerivative(reference_buoyancy(z✶_column))simulation.callbacks[:∂ₜb✶] = Callback(∂ₜb✶, PrecedingIterations(TimeInterval(τ / 5)))eₐˢ  = SubFilterAvailablePotentialEnergy(model, gfilter; method=lookup)wbᵣˢ = SubFilterAvailablePotentialToKineticEnergyConversion(model, gfilter; method=lookup)Πₐ   = AvailablePotentialEnergyCrossScaleFlux(model, gfilter; method=lookup)εₐˢ  = SubFilterAvailablePotentialEnergyDissipationRate(model, gfilter; method=lookup)Rˢ   = SubFilterAvailablePotentialEnergyReferenceTendency(model, gfilter, ∂ₜb✶; method=lookup)∫wbᵣˢ = Integral(wbᵣˢ)∫Πₐ   = Integral(Πₐ)∫εₐˢ  = Integral(εₐˢ)∫Rˢ   = Integral(Rˢ)∂ₜ∫eₐˢ = TimeDerivative(Integral(eₐˢ))

TimeDerivative of 1×1×1 Field{Nothing, Nothing, Nothing} reduced over dims = (1, 2, 3) on RectilinearGrid on CPU

Output

A single NetCDF writer stores, every τ/5, vertical (x–z) slices at a fixed y index (the flow is periodic and statistically homogeneous in y, so the particular plane makes no difference) of the buoyancy b and of five budget terms: the two cross-scale fluxes Πₖ and Πₐ, the two subfilter dissipations εₖˢ and εₐˢ, and the reference tendency Rˢ. It also stores the tendencies and volume integrals of the two budgets. The indices that cut the slices leave the volume integrals whole, since those have no spatial dimension left to slice.

using NCDatasetsfilename = joinpath(@__DIR__, "rayleigh_taylor_instability")outputs = (; b, Πₖ, εₖˢ, Rˢ, Πₐ, εₐˢ, ∂ₜ∫eₖˢ, ∫wbᵣˢ, ∫Πₖ, ∫εₖˢ, ∂ₜ∫eₐˢ, ∫Πₐ, ∫εₐˢ, ∫Rˢ)simulation.output_writers[:nc] = NetCDFWriter(model, outputs,                                              filename = filename,                                              schedule = TimeInterval(τ / 5),                                              indices = (:, 1, :),                                              overwrite_files = true)

NetCDFWriter scheduled on TimeInterval(200 ms):
├── filepath: rayleigh_taylor_instability.nc
├── dimensions: time(0), x_faa(48), x_caa(48), y_afa(1), y_aca(1), z_aaf(49), z_aac(48)
├── 14 outputs: (b, Πₖ, εₖˢ, Rˢ, Πₐ, εₐˢ, ∂ₜ∫eₖˢ, ∫wbᵣˢ, ∫Πₖ, ∫εₖˢ, ∂ₜ∫eₐˢ, ∫Πₐ, ∫εₐˢ, ∫Rˢ)
├── 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.282 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 (40.312 seconds)
[ Info: Executing initial time step...
[ Info:     ... initial time step complete (778.650 ms).
┌ Info: iter =    100,  [034.85%] time = 2.788 seconds,  Δt = 43.053 ms,  walltime = 2.913 minutes,  walltime / timestep = 378.016 ms
└       |u⃗|ₘₐₓ = [8.94e-02,  1.03e-01,  3.41e-01] m/s,  advective CFL = 0.8,  diffusive CFL = 0.05,  νₘₐₓ = 0.0005 m²/s
┌ Info: iter =    200,  [066.41%] time = 5.313 seconds,  Δt = 22.522 ms,  walltime = 3.232 minutes,  walltime / timestep = 191.820 ms
└       |u⃗|ₘₐₓ = [3.53e-01,  3.30e-01,  5.92e-01] m/s,  advective CFL = 0.8,  diffusive CFL = 0.048,  νₘₐₓ = 0.00093 m²/s
┌ Info: iter =    300,  [088.52%] time = 7.081 seconds,  Δt = 16.260 ms,  walltime = 3.478 minutes,  walltime / timestep = 147.544 ms
└       |u⃗|ₘₐₓ = [5.27e-01,  6.01e-01,  6.69e-01] m/s,  advective CFL = 0.8,  diffusive CFL = 0.044,  νₘₐₓ = 0.0012 m²/s
[ Info: Simulation is stopping after running for 1.677 minutes.
[ Info: Simulation time 8 seconds equals or exceeds stop time 8 seconds.

We load the output back as FieldTimeSeries. Each snapshot of a slice, such as b_ts[n], is a Field that carries its grid and location, so it plots directly as an x–z heatmap:

using CairoMakiefilepath = simulation.output_writers[:nc].filepathb_ts   = FieldTimeSeries(filepath, "b")Πₖ_ts  = FieldTimeSeries(filepath, "Πₖ")εₖˢ_ts = FieldTimeSeries(filepath, "εₖˢ")Rˢ_ts  = FieldTimeSeries(filepath, "Rˢ")Πₐ_ts  = FieldTimeSeries(filepath, "Πₐ")εₐˢ_ts = FieldTimeSeries(filepath, "εₐˢ")times  = b_ts.times

Every record carries both tendencies and their budget terms at the same time, and each volume integral comes back as a FieldTimeSeries of a single point. The four that the budget plots show as stored, ∫τ(w,bᵣ) dV, ∫Πₖ dV, ∫Πₐ dV and ∫Rˢ dV, go to lines! as they are. The negated curves and the residuals need arithmetic, which a FieldTimeSeries does not support, so interior turns each term into a vector over the records from the second on: the first record of a tendency has no earlier state to difference against and is written as zero. The residuals measure how well the two subfilter budgets close.

∫wbᵣˢ_ts = FieldTimeSeries(filepath, "∫wbᵣˢ")∫Πₖ_ts   = FieldTimeSeries(filepath, "∫Πₖ")∫Πₐ_ts   = FieldTimeSeries(filepath, "∫Πₐ")∫Rˢ_ts   = FieldTimeSeries(filepath, "∫Rˢ")nb    = 2:length(times)t_bud = times[nb]budget_term(series) = interior(series, 1, 1, 1, nb)deₖˢdt   = budget_term(FieldTimeSeries(filepath, "∂ₜ∫eₖˢ"))wbᵣˢ_bud = budget_term(∫wbᵣˢ_ts)Πₖ_bud   = budget_term(∫Πₖ_ts)εₖˢ_bud  = budget_term(FieldTimeSeries(filepath, "∫εₖˢ"))deₐˢdt   = budget_term(FieldTimeSeries(filepath, "∂ₜ∫eₐˢ"))Πₐ_bud   = budget_term(∫Πₐ_ts)εₐˢ_bud  = budget_term(FieldTimeSeries(filepath, "∫εₐˢ"))Rˢ_bud   = budget_term(∫Rˢ_ts)

Residuals in sum-to-zero form: the negative tendency plus the right-hand side, so the plotted curves add to it

eₖˢ_resid = @. -deₖˢdt + wbᵣˢ_bud + Πₖ_bud - εₖˢ_budeₐˢ_resid = @. -deₐˢdt - wbᵣˢ_bud + Πₐ_bud - εₐˢ_bud + Rˢ_bud

Plotting

We build the figure in four rows. The top two hold vertical slices, three panels each: the buoyancy b (the spikes and bubbles) and the KE terms Πₖ and εₖˢ in the first row, and the APE terms Rˢ, Πₐ and εₐˢ in the second. The bottom two rows hold the two volume-integrated subfilter budgets.

set_theme!(Theme(fontsize=18))fig = Figure(size=(1500, 1600))n = Observable(1)

Observable(1)

Every slice gets a diverging colormap centred on zero. The two cross-scale fluxes share one colour range, and so do the two dissipations, so the KE and APE terms stacked in each column compare directly; the buoyancy and Rˢ get ranges of their own.

blim = 0.8 * maximum(abs, b_ts[1])Rlim = 0.5 * maximum(abs, Rˢ_ts)Πlim = 0.5 * max(maximum(abs, Πₖ_ts), maximum(abs, Πₐ_ts))εlim = 0.5 * max(maximum(abs, εₖˢ_ts), maximum(abs, εₐˢ_ts))panels = [(1, 1, b_ts,   "buoyancy, b",                      blim),          (1, 2, Πₖ_ts,  "cross-scale KE flux, Πₖ",          Πlim),          (1, 3, εₖˢ_ts, "subfilter KE dissipation, εₖˢ",    εlim),          (2, 1, Rˢ_ts,  "subfilter reference tendency, Rˢ", Rlim),          (2, 2, Πₐ_ts,  "cross-scale APE flux, Πₐ",         Πlim),          (2, 3, εₐˢ_ts, "subfilter APE dissipation, εₐˢ",   εlim)]for (row, column, series, label, lim) in panels    ax = Axis(fig[row + 1, 2column - 1]; title=label, xlabel="x", ylabel="z", aspect=1)    hm = heatmap!(ax, @lift(series[$n]); colormap=:balance, colorrange=(-lim, lim))    Colorbar(fig[row + 1, 2column], hm)end

The fourth row 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 subfilter buoyancy flux ∫τ(w,bᵣ) dV and the cross-scale flux ∫Πₖ dV handed down from the filtered scales, 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:6]; xlabel="time [free-fall units]", title="Subfilter kinetic energy budget")lines!(ax_bud, t_bud, -deₖˢdt,   label="−d(∫eₖˢ)/dt")lines!(ax_bud, ∫wbᵣˢ_ts,         label="∫τ(w,bᵣ) dV  (conversion from subfilter APE)")lines!(ax_bud, ∫Πₖ_ts,           label="∫Πₖ dV  (flux from filtered scales)")lines!(ax_bud, t_bud, -εₖˢ_bud,  label="−∫εₖˢ dV  (subfilter dissipation)")lines!(ax_bud, t_bud, eₖˢ_resid, label="residual", color=:black, linestyle=:dash)axislegend(ax_bud; position=:lt, labelsize=12)

Makie.Legend()

The bottom row shows the subfilter APE budget in the same form: the negative tendency −d(∫eₐˢ)/dt, the conversion −∫τ(w,bᵣ) dV into subfilter kinetic energy, the cross-scale flux ∫Πₐ dV, the subfilter dissipation −∫εₐˢ dV and the reference tendency ∫Rˢ dV.

ax_ape = Axis(fig[5, 1:6]; xlabel="time [free-fall units]", title="Subfilter available potential energy budget")lines!(ax_ape, t_bud, -deₐˢdt,   label="−d(∫eₐˢ)/dt")lines!(ax_ape, t_bud, -wbᵣˢ_bud, label="−∫τ(w,bᵣ) dV  (conversion to subfilter KE)")lines!(ax_ape, ∫Πₐ_ts,           label="∫Πₐ dV  (flux from filtered scales)")lines!(ax_ape, t_bud, -εₐˢ_bud,  label="−∫εₐˢ dV  (subfilter dissipation)")lines!(ax_ape, ∫Rˢ_ts,           label="∫Rˢ dV  (reference tendency)")lines!(ax_ape, t_bud, eₐˢ_resid, label="residual", color=:black, linestyle=:dash)axislegend(ax_ape; position=:lt, labelsize=12)for ax in (ax_bud, ax_ape)    vlines!(ax, @lift(times[$n]), color=:black, linestyle=:dash)endtitle = @lift "Rayleigh-Taylor instability, t = " * string(round(times[$n] / τ, digits=2)) * " τ"fig[1, 1:6] = Label(fig, title, fontsize=22, tellwidth=false)@info "Animating..."record(fig, "rayleigh_taylor_instability.mp4", 1:length(times), framerate=11) 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. In the subfilter 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. In the subfilter APE budget, the cross-scale flux ∫Πₐ dV is the main source. About half of what it supplies is converted into subfilter kinetic energy, through the same ∫τ(w,bᵣ) dV that feeds the KE budget, and most of the rest is dissipated. The reference tendency ∫Rˢ dV is the smallest of the four terms on the right-hand side, but still many times the residual. Both residuals are small, which shows that the budgets close well and that the filtering analysis is consistent with the simulation.

This page was generated using Literate.jl.