Lock Exchange Example
This is a simple example of a lock exchange problem on a slope. This example demonstrates:
- How to set up a 2D grid with a sloping bottom with an immersed boundary
- Initializing a hydrostatic free surface model
- Including variable density initial conditions
- Saving outputs of the simulation
- Creating an animation with CairoMakie
The Lock Exchange Problem
This use case is a basic example where there are fluids of two different densities (due to temperature, salinity, etc.) that are separated by a ‘lock’ at time $t=0$, and we calculate the evolution of how these fluids interact as time progresses. This lock exchange implementation can be a representation of scenarios where water of different salinities or temperatures meet and form sharp density gradients. For example, in estuaries or in the Denmark Strait overflow. Solutions of this problem describe how the fluids interact with each other as time evolves and can be described by hydrostatic Boussinesq equations. In this example, we use buoyancy as a tracer; see the Boussinesq approximation section for more details on the Boussinesq approximation.
Install dependencies
First let's make sure we have all required packages installed.
using Pkg
pkg"add Oceananigans, CairoMakie"Import Required Packages
using Oceananigans
using Oceananigans.Units
using Oceananigans.Grids
using Oceananigans.Diagnostics
using Oceananigans.OutputWriters
using Oceananigans.Operators
using Printf
using CairoMakieSet up 2D Rectilinear Grid
Set resolution of the simulation grid
Nx, Nz = 128, 64(128, 64)Set grid size
L = 8kilometers # horizontal length
H = 50meters # depth50.0Set domain: We wrap z into a MutableVerticalDiscretization. This allows the HydrostaticFreeSurface model (which we construct further down) to use a time-evolving ZStarCoordinate free-surface–following vertical coordinate.
x = (0, L)
z = MutableVerticalDiscretization((-50, 0))MutableVerticalDiscretization with reference interfaces r:
(-50, 0)Initialize the grid: Additional details on the grid set up can be found in Grids section of the documentation
underlying_grid = RectilinearGrid(; size=(Nx, Nz), x, z, halo=(5, 5),
topology=(Bounded, Flat, Bounded))128×1×64 RectilinearGrid{Float64, Bounded, Flat, Bounded} on CPU with 5×0×5 halo
├── Bounded x ∈ [0.0, 8000.0] regularly spaced with Δx=62.5
├── Flat y
└── Bounded z ∈ [-50.0, 0.0] regularly spaced with mutable Δr=0.78125Add a slope at the bottom of the grid
h_left = -H
h_right = -H/2
slope = (h_right - h_left) / L
bottom(x) = h_left + slope * xbottom (generic function with 1 method)Use an immersed boundary with grid fitted bottom to describe the sloped bottom of the domain
grid = ImmersedBoundaryGrid(underlying_grid, PartialCellBottom(bottom))128×1×64 ImmersedBoundaryGrid{Float64, Bounded, Flat, Bounded} on CPU with 5×0×5 halo:
├── immersed_boundary: PartialCellBottom(mean(zb)=-35.3401, min(zb)=-49.9023, max(zb)=0.0, ϵ=0.2)
├── underlying_grid: 128×1×64 RectilinearGrid{Float64, Bounded, Flat, Bounded} on CPU with 5×0×5 halo
├── Bounded x ∈ [0.0, 8000.0] regularly spaced with Δx=62.5
├── Flat y
└── Bounded z ∈ [-50.0, 0.0] regularly spaced with mutable Δr=0.78125Initialize the model
- Want to use a hydrostatic model since horizontal motion may be more significant than vertical motion
- Tracers act as markers within the fluid to track movement and dispersion
- Vertical closure
CATKEVerticalDiffusivityhandles small-scale vertical turbulence - Weighted Essentially Non-Oscillatory (WENO) methods are useful for capturing sharp changes in density
model = HydrostaticFreeSurfaceModel(; grid,
tracers = :b,
buoyancy = BuoyancyTracer(),
closure = CATKEVerticalDiffusivity(),
momentum_advection = WENO(order=5),
tracer_advection = WENO(order=7),
free_surface = SplitExplicitFreeSurface(grid; substeps=10))HydrostaticFreeSurfaceModel{CPU, ImmersedBoundaryGrid}(time = 0 seconds, iteration = 0)
├── grid: 128×1×64 ImmersedBoundaryGrid{Float64, Bounded, Flat, Bounded} on CPU with 5×0×5 halo
├── timestepper: QuasiAdamsBashforth2TimeStepper
├── tracers: (b, e)
├── closure: CATKEVerticalDiffusivity{VerticallyImplicitTimeDiscretization}
├── buoyancy: BuoyancyTracer with ĝ = NegativeZDirection()
├── free surface: SplitExplicitFreeSurface with gravitational acceleration 9.80665 m s⁻²
│ └── substepping: FixedSubstepNumber(7)
├── advection scheme:
│ ├── momentum: WENO{3, Float64, Float32}(order=5)
│ ├── b: WENO{4, Float64, Float32}(order=7)
│ └── e: WENO{4, Float64, Float32}(order=7)
├── vertical_coordinate: ZStarCoordinate
└── coriolis: NothingSet variable density initial conditions
Set initial conditions for lock exchange with different buoyancies.
bᵢ(x, z) = x > L/2 ? 0.01 : 0.06
set!(model, b=bᵢ)Construct a Simulation
Fast wave speeds make the equations stiff, so the CFL condition restricts the timestep to adequately small values to maintain numerical stability.
Set the timesteps
Δt = 1second
stop_time = 5hours
simulation = Simulation(model; Δt, stop_time)Simulation of HydrostaticFreeSurfaceModel{CPU, ImmersedBoundaryGrid}(time = 0 seconds, iteration = 0)
├── Next time step: 1 second
├── run_wall_time: 0 seconds
├── run_wall_time / iteration: NaN days
├── stop_time: 5 hours
├── stop_iteration: Inf
├── wall_time_limit: Inf
├── minimum_relative_step: 0.0
├── callbacks: OrderedDict with 4 entries:
│ ├── stop_time_exceeded => Callback of stop_time_exceeded on IterationInterval(1)
│ ├── stop_iteration_exceeded => Callback of stop_iteration_exceeded on IterationInterval(1)
│ ├── wall_time_limit_exceeded => Callback of wall_time_limit_exceeded on IterationInterval(1)
│ └── nan_checker => Callback of NaNChecker for u on IterationInterval(100)
└── output_writers: OrderedDict with no entriesThe TimeStepWizard is incorporated in the simulation via the conjure_time_step_wizard! helper function and it ensures stable time-stepping with a Courant-Freidrichs-Lewy (CFL) number of 0.35.
conjure_time_step_wizard!(simulation, cfl=0.35)Track Simulation Progress
Wall clock represents the real world time as opposed to simulation time
wall_clock = Ref(time_ns())Base.RefValue{UInt64}(0x003ebe3aedc7407b)Define callback function to show how the simulation is progressing alongside with some flow statistics.
function progress(sim)
elapsed = 1e-9 * (time_ns() - wall_clock[])
msg = @sprintf("Iter: %6d, time: %s, Δt: %s, wall: %s, max|w| = %6.3e m s⁻¹",
iteration(sim),
prettytime(sim),
prettytime(sim.Δt),
prettytime(elapsed),
maximum(abs, sim.model.velocities.w))
wall_clock[] = time_ns()
@info msg
return nothing
end
add_callback!(simulation, progress, name = :progress, TimeInterval(30minutes))Add output writer
Here, we construct a JLD2Writer to save some output every save_interval.
b = model.tracers.b
e = model.tracers.e
u, v, w = model.velocities
N² = ∂z(b)
filename = "lock_exchange.jld2"
save_interval = 2minutes
simulation.output_writers[:fields] = JLD2Writer(model, (; b, e, u, N²);
filename = filename,
schedule = TimeInterval(save_interval),
overwrite_existing = true)JLD2Writer scheduled on TimeInterval(2 minutes):
├── filepath: lock_exchange.jld2
├── 4 outputs: (b, e, u, N²)
├── array_type: Array{Float32}
├── including: [:grid, :coriolis, :buoyancy, :closure]
├── file_splitting: NoFileSplitting
└── file size: 93.8 KiBRun Simulation
run!(simulation)
@info "Simulation finished. Output saved to $(filename)"[ Info: Initializing simulation...
[ Info: Iter: 0, time: 0 seconds, Δt: 1.100 seconds, wall: 43.931 seconds, max|w| = 0.000e+00 m s⁻¹
[ Info: ... simulation initialization complete (32.903 seconds)
[ Info: Executing initial time step...
[ Info: ... initial time step complete (48.394 seconds).
[ Info: Iter: 823, time: 30 minutes, Δt: 2.060 seconds, wall: 1.003 minutes, max|w| = 1.273e-01 m s⁻¹
[ Info: Iter: 1747, time: 1 hour, Δt: 1.914 seconds, wall: 11.476 seconds, max|w| = 1.553e-01 m s⁻¹
[ Info: Iter: 2791, time: 1.500 hours, Δt: 1.504 seconds, wall: 12.796 seconds, max|w| = 1.808e-01 m s⁻¹
[ Info: Iter: 3901, time: 2 hours, Δt: 4.168 seconds, wall: 13.843 seconds, max|w| = 6.461e-02 m s⁻¹
[ Info: Iter: 4360, time: 2.500 hours, Δt: 3.609 seconds, wall: 5.893 seconds, max|w| = 7.482e-02 m s⁻¹
[ Info: Iter: 4907, time: 3 hours, Δt: 3.972 seconds, wall: 7.020 seconds, max|w| = 7.278e-02 m s⁻¹
[ Info: Iter: 5309, time: 3.500 hours, Δt: 6.631 seconds, wall: 5.049 seconds, max|w| = 3.137e-02 m s⁻¹
[ Info: Iter: 5595, time: 4 hours, Δt: 8.418 seconds, wall: 3.638 seconds, max|w| = 3.366e-02 m s⁻¹
[ Info: Iter: 5799, time: 4.500 hours, Δt: 7.208 seconds, wall: 2.662 seconds, max|w| = 4.000e-02 m s⁻¹
[ Info: Simulation is stopping after running for 2.641 minutes.
[ Info: Simulation time 5 hours equals or exceeds stop time 5 hours.
[ Info: Iter: 6150, time: 5 hours, Δt: 4.429 seconds, wall: 4.491 seconds, max|w| = 6.054e-02 m s⁻¹
[ Info: Simulation finished. Output saved to lock_exchange.jld2
Load Saved TimeSeries Values
ut = FieldTimeSeries(filename, "u")
N²t = FieldTimeSeries(filename, "N²")
bt = FieldTimeSeries(filename, "b")
et = FieldTimeSeries(filename, "e")
times = bt.times0.0:120.0:18000.0Visualize Simulation Outputs
We use Makie's Observable to animate the data. To dive into how Observables work we refer to Makie.jl's Documentation.
n = Observable(1)
title = @lift @sprintf("t = %5.2f hours", times[$n] / hour)
un = @lift ut[$n]
N²n = @lift N²t[$n]
bn = @lift bt[$n]
en = @lift et[$n]For visualization color ranges (use last snapshot)
umax = maximum(abs, ut[end])
N²max = maximum(abs, N²t[end])
bmax = maximum(abs, bt[end])
emax = maximum(abs, et[end])Use snapshots to create Makie visualization for $b$, $e$, $u$, and $N^2$.
axis_kwargs = (xlabel = "x [m]", ylabel = "z [m]",
limits = ((0, L), (-H, 0)), titlesize = 18)
fig = Figure(size = (800, 900))
fig[1, :] = Label(fig, title, fontsize = 24, tellwidth = false)
nan_color = :grey
ax_b = Axis(fig[2, 1]; title = "b (buoyancy)", axis_kwargs...)
hm_b = heatmap!(ax_b, bn; nan_color, colorrange = (0, bmax), colormap = :thermal)
Colorbar(fig[2, 2], hm_b, label = "m s⁻²")
ax_e = Axis(fig[3, 1]; title = "e (turbulent kinetic energy)", axis_kwargs...)
hm_e = heatmap!(ax_e, en; nan_color, colorrange = (0, emax), colormap = :magma)
Colorbar(fig[3, 2], hm_e, label = "m² s⁻²")
ax_u = Axis(fig[4, 1]; title = "u (horizontal velocity)", axis_kwargs...)
hm_u = heatmap!(ax_u, un; nan_color, colorrange = (-umax, umax), colormap = :balance)
Colorbar(fig[4, 2], hm_u, label = "m s⁻¹")
ax_N² = Axis(fig[5, 1]; title = "N² (stratification)", axis_kwargs...)
hm_N² = heatmap!(ax_N², N²n; nan_color, colorrange = (-N²max/4, N²max), colormap = :haline)
Colorbar(fig[5, 2], hm_N², label = "s⁻²")
fig
Create Animation
frames = 1:length(times)
record(fig, "lock_exchange.mp4", frames; framerate = 8) do i
@info "Animating frame $i out of $(frames[end])"
n[] = i
end"lock_exchange.mp4"The visualization shows the time evolution of buoyancy $b$, turbulent kinetic energy $e$, stratification $N²$, and zonal velocity $u$. Initially, the two water masses are separated horizontally by a sharp density interface. As the flow evolves, gravity currents form and the dense fluid moves beneath the lighter fluid. This shows the characteristic transition from horizontal to vertical density separation for the lock exchange problem.
Julia version and environment information
This example was executed with the following version of Julia:
using InteractiveUtils: versioninfo
versioninfo()Julia Version 1.12.2
Commit ca9b6662be4 (2025-11-20 16:25 UTC)
Build Info:
Official https://julialang.org release
Platform Info:
OS: Linux (x86_64-linux-gnu)
CPU: 128 × AMD EPYC 9374F 32-Core Processor
WORD_SIZE: 64
LLVM: libLLVM-18.1.7 (ORCJIT, znver4)
GC: Built with stock GC
Threads: 1 default, 1 interactive, 1 GC (on 128 virtual cores)
Environment:
LD_LIBRARY_PATH =
JULIA_PKG_SERVER_REGISTRY_PREFERENCE = eager
JULIA_DEPOT_PATH = /var/lib/buildkite-agent/.julia-oceananigans
JULIA_PROJECT = /var/lib/buildkite-agent/Oceananigans.jl-27790/docs/
JULIA_VERSION = 1.12.2
JULIA_LOAD_PATH = @:@v#.#:@stdlib
JULIA_VERSION_ENZYME = 1.10.10
JULIA_PYTHONCALL_EXE = /var/lib/buildkite-agent/Oceananigans.jl-27790/docs/.CondaPkg/.pixi/envs/default/bin/python
JULIA_DEBUG = Literate
These were the top-level packages installed in the environment:
import Pkg
Pkg.status()Status `~/Oceananigans.jl-27790/docs/Project.toml`
[79e6a3ab] Adapt v4.4.0
[052768ef] CUDA v5.9.5
[13f3f980] CairoMakie v0.15.8
[e30172f5] Documenter v1.16.1
[daee34ce] DocumenterCitations v1.4.1
[033835bb] JLD2 v0.6.3
[98b081ad] Literate v2.21.0
[da04e1cc] MPI v0.20.23
[85f8d34a] NCDatasets v0.14.10
[9e8cae18] Oceananigans v0.103.0 `~/Oceananigans.jl-27790`
[f27b6e38] Polynomials v4.1.0
[6038ab10] Rotations v1.7.1
[d496a93d] SeawaterPolynomials v0.3.10
[09ab397b] StructArrays v0.7.2
[bdfc003b] TimesDates v0.3.3
[2e0b0046] XESMF v0.1.6
[b77e0a4c] InteractiveUtils v1.11.0
[37e2e46d] LinearAlgebra v1.12.0
[44cfe95a] Pkg v1.12.0
This page was generated using Literate.jl.