Hydrostatic lock exchange with CATKEVerticalDiffusivity
This example simulates a lock exchange problem on a slope. It 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
- Applying bottom drag boundary conditions on immersed boundaries using
BulkDrag - 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 Printf
using CairoMakieSet up 2D Rectilinear Grid
Set resolution of the simulation grid
Nx = 128
Nz = 64
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 = (-L/8, 7L/8)
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 ∈ [-1000.0, 7000.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)=-38.0791, min(zb)=-50.0, 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 ∈ [-1000.0, 7000.0] regularly spaced with Δx=62.5
├── Flat y
└── Bounded z ∈ [-50.0, 0.0] regularly spaced with mutable Δr=0.78125Set up a bottom drag boundary condition on the immersed boundary
To apply drag to the flow, we use BulkDrag, which implements quadratic drag proportional to Cᴰ |U| u, where Cᴰ is the drag coefficient and |U| = sqrt(u² + v²) is the horizontal speed. We use a drag coefficient Cᴰ = 0.002, a reasonable value for seafloor drag.
Cᴰ = 0.002
drag = BulkDrag(coefficient=Cᴰ)FluxBoundaryCondition: BulkDragFunction(QuadraticFormulation(), Nothing, Cᴰ=0.002)BulkDrag can be applied both to domain boundaries (like the bottom of a Bounded grid) and to immersed boundaries. Here we apply it to the immersed sloping bottom boundary:
u_bcs = FieldBoundaryConditions(bottom=drag, immersed=drag)Oceananigans.FieldBoundaryConditions, with boundary conditions
├── west: DefaultBoundaryCondition (FluxBoundaryCondition: Nothing)
├── east: DefaultBoundaryCondition (FluxBoundaryCondition: Nothing)
├── south: DefaultBoundaryCondition (FluxBoundaryCondition: Nothing)
├── north: DefaultBoundaryCondition (FluxBoundaryCondition: Nothing)
├── bottom: FluxBoundaryCondition: BulkDragFunction(QuadraticFormulation(), Nothing, Cᴰ=0.002)
├── top: DefaultBoundaryCondition (FluxBoundaryCondition: Nothing)
└── immersed: FluxBoundaryCondition: BulkDragFunction(QuadraticFormulation(), Nothing, Cᴰ=0.002)Note: in a 2D simulation with Flat in the y-direction, we don't need v boundary conditions.
Initialize 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),
boundary_conditions = (; u=u_bcs),
free_surface = SplitExplicitFreeSurface(grid; substeps=20))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(14)
├── 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 = 6hours
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: 6 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.3)Track Simulation Progress
# Wall clock represents the real world time as opposed to simulation time
wall_clock = Ref(time_ns())
# 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))
@info msg
wall_clock[] = time_ns()
return nothing
end
add_callback!(simulation, progress, IterationInterval(1000))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 = "hydrostatic_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: hydrostatic_lock_exchange.jld2
├── 4 outputs: (b, e, u, N²)
├── array_type: Array{Float32}
├── including: [:grid, :coriolis, :buoyancy, :closure]
├── file_splitting: NoFileSplitting
└── file size: 0 bytes (file not yet created)Run Simulation
run!(simulation)
@info "Simulation finished. Output saved to $(filename)"[ Info: Initializing simulation...
[ Info: Iter: 0, time: 0 seconds, Δt: 1.100 seconds, wall: 53.223 seconds, max|w| = 0.000e+00 m s⁻¹
[ Info: ... simulation initialization complete (48.426 seconds)
[ Info: Executing initial time step...
[ Info: ... initial time step complete (50.707 seconds).
[ Info: Iter: 1000, time: 31.860 minutes, Δt: 2.036 seconds, wall: 1.097 minutes, max|w| = 1.152e-01 m s⁻¹
[ Info: Iter: 2000, time: 1.004 hours, Δt: 1.522 seconds, wall: 13.604 seconds, max|w| = 1.525e-01 m s⁻¹
[ Info: Iter: 3000, time: 1.436 hours, Δt: 1.503 seconds, wall: 13.448 seconds, max|w| = 1.557e-01 m s⁻¹
[ Info: Iter: 4000, time: 1.833 hours, Δt: 1.658 seconds, wall: 13.563 seconds, max|w| = 1.414e-01 m s⁻¹
[ Info: Iter: 5000, time: 2.213 hours, Δt: 1.392 seconds, wall: 13.792 seconds, max|w| = 1.657e-01 m s⁻¹
[ Info: Iter: 6000, time: 2.674 hours, Δt: 2.601 seconds, wall: 13.640 seconds, max|w| = 8.337e-02 m s⁻¹
[ Info: Iter: 7000, time: 3.678 hours, Δt: 5.195 seconds, wall: 14.083 seconds, max|w| = 4.492e-02 m s⁻¹
[ Info: Iter: 8000, time: 4.343 hours, Δt: 2.048 seconds, wall: 13.967 seconds, max|w| = 1.070e-01 m s⁻¹
[ Info: Iter: 9000, time: 4.839 hours, Δt: 1.659 seconds, wall: 13.741 seconds, max|w| = 1.402e-01 m s⁻¹
[ Info: Iter: 10000, time: 5.442 hours, Δt: 2.708 seconds, wall: 13.648 seconds, max|w| = 8.528e-02 m s⁻¹
[ Info: Simulation is stopping after running for 4.081 minutes.
[ Info: Simulation time 6 hours equals or exceeds stop time 6 hours.
[ Info: Simulation finished. Output saved to hydrostatic_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:21600.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$.
nan_color = :grey
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)
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⁻²")
Nt = length(times)
record(fig, "hydrostatic_lock_exchange.mp4", 1:Nt; framerate = 8) do nn
@info "Animating frame $nn out of $Nt"
n[] = nn
end"hydrostatic_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-28008/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-28008/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-28008/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.1 `~/Oceananigans.jl-28008`
[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.