diff --git a/docs/make.jl b/docs/make.jl
index e8adfd0958..2b1b627f51 100644
--- a/docs/make.jl
+++ b/docs/make.jl
@@ -99,9 +99,9 @@ numerical_pages = [
]
verification_pages = [
- "Taylor-Green vortex" => "verification/taylor_green_vortex.md",
- "Stratified Couette flow" => "verification/stratified_couette_flow.md",
- "Convergence tests" => "verification/convergence_tests.md"
+ "Convergence tests" => "verification/convergence_tests.md",
+ "Lid-driven cavity" => "verification/lid_driven_cavity.md",
+ "Stratified Couette flow" => "verification/stratified_couette_flow.md"
]
appendix_pages = [
diff --git a/docs/src/verification/lid_driven_cavity.md b/docs/src/verification/lid_driven_cavity.md
index 338f1efab0..d0c3d51e00 100644
--- a/docs/src/verification/lid_driven_cavity.md
+++ b/docs/src/verification/lid_driven_cavity.md
@@ -2,8 +2,96 @@
The lid-driven cavity test problem has been used for a long time as a simple verification test for computational fluid
dynamics codes. First described by [Burggraf66](@cite), the fluid is contained in a square cavity with Dirchlet boundary
-conditions on all four sides. The top wall moves with some velocity ``U`` while the other three walls are stationary.
+conditions on all four sides. The top wall moves with velocity ``U = 1`` while the other three walls are stationary.
The solution reaches a laminar steady-state whose properties can be compared with a huge amount of existing data. The
canonical database is given by [Ghia82](@cite) who report detailed information on the velocity fields as well as the
streamline and vorticity contours at various Reynolds numbers. More accurate data is reported by [Botella98](@cite),
[Erturk05](@cite), and [Bruneau06](@cite).
+
+Below are lid-driven cavity simulation results from Oceananigans.jl compared with [Ghia82](@cite) for Re = 100, 400,
+1000, 3200, 5000, 7500, and 10000.
+
+## Re = 100
+
+```@raw html
+
+
+```
+
+## Re = 400
+
+```@raw html
+
+```
+
+## Re = 1,000
+
+```@raw html
+
+```
+
+## Re = 3,200
+
+```@raw html
+
+```
+
+## Re = 5,000
+
+```@raw html
+
+```
+
+## Re = 7,500
+
+```@raw html
+
+```
+
+## Re = 10,000
+
+```@raw html
+
+```
diff --git a/docs/src/verification/taylor_green_vortex.md b/docs/src/verification/taylor_green_vortex.md
deleted file mode 100644
index 79e32ce2f3..0000000000
--- a/docs/src/verification/taylor_green_vortex.md
+++ /dev/null
@@ -1,13 +0,0 @@
-# Taylor-Green vortex
-An exact solution to the two-dimensional incompressible Navier-Stokes equations is given by [Taylor37](@cite) describing
-the unsteady flow of a vortex decaying under viscous dissipation. The viscous terms balance the time derivatives while
-the nonlinear advection terms balance the pressure gradient term. We use the doubly-periodic solution described by
-(p. 310) [Hesthaven07](@cite)
-
-```math
-\begin{aligned}
- u(x, y, t) &= -\sin(2\pi y) e^{-4\pi^2\nu t} \\
- v(x, y, t) &= \sin(2\pi x) e^{-4\pi^2\nu t} \\
- p(x, y, t) &= -\cos(2\pi x) \cos(2\pi y) e^{-8\pi^2\nu t}
-\end{aligned}
-```
diff --git a/verification/lid_driven_cavity/lid_driven_cavity.jl b/verification/lid_driven_cavity/lid_driven_cavity.jl
new file mode 100644
index 0000000000..2a1654a900
--- /dev/null
+++ b/verification/lid_driven_cavity/lid_driven_cavity.jl
@@ -0,0 +1,97 @@
+using Printf
+using Logging
+
+using Oceananigans
+using Oceananigans: Face, Cell
+using Oceananigans.Diagnostics
+using Oceananigans.OutputWriters
+using Oceananigans.AbstractOperations
+
+Logging.global_logger(OceananigansLogger())
+
+function simulate_lid_driven_cavity(; Re, N, end_time)
+ topology = (Flat, Bounded, Bounded)
+ domain = (x=(0, 1), y=(0, 1), z=(0, 1))
+ grid = RegularCartesianGrid(topology=topology, size=(1, N, N); domain...)
+
+ v_bcs = VVelocityBoundaryConditions(grid,
+ top = ValueBoundaryCondition(1.0),
+ bottom = ValueBoundaryCondition(0.0)
+ )
+
+ w_bcs = WVelocityBoundaryConditions(grid,
+ north = ValueBoundaryCondition(0.0),
+ south = ValueBoundaryCondition(0.0)
+ )
+
+ model = IncompressibleModel(
+ grid = grid,
+ buoyancy = nothing,
+ tracers = nothing,
+ coriolis = nothing,
+ boundary_conditions = (v=v_bcs, w=w_bcs),
+ closure = IsotropicDiffusivity(ν=1/Re)
+ )
+
+ u, v, w = model.velocities
+ ζ_op = ∂y(w) - ∂z(v)
+ ζ = Field(Cell, Face, Face, model.architecture, model.grid, TracerBoundaryConditions(grid))
+ ζ_computation = Computation(ζ_op, ζ)
+
+ fields = Dict(
+ "v" => model.velocities.v,
+ "w" => model.velocities.w,
+ "ζ" => model -> ζ_computation(model)
+ )
+
+ dims = Dict("ζ" => ("xC", "yF", "zF"))
+ global_attributes = Dict("Re" => Re)
+ output_attributes = Dict("ζ" => Dict("longname" => "vorticity", "units" => "1/s"))
+
+ field_output_writer =
+ NetCDFOutputWriter(model, fields, filename="lid_driven_cavity_Re$Re.nc", time_interval=0.1,
+ global_attributes=global_attributes, output_attributes=output_attributes,
+ dimensions=dims)
+
+ max_Δt = 0.25 * model.grid.Δy^2 * Re / 2 # Make sure not to violate diffusive CFL.
+ wizard = TimeStepWizard(cfl=0.1, Δt=1e-6, max_change=1.1, max_Δt=max_Δt)
+
+ cfl = AdvectiveCFL(wizard)
+ dcfl = DiffusiveCFL(wizard)
+
+ simulation = Simulation(model, Δt=wizard, stop_time=end_time, progress=print_progress,
+ iteration_interval=20, parameters=(cfl=cfl, dcfl=dcfl))
+
+ simulation.output_writers[:fields] = field_output_writer
+
+ run!(simulation)
+
+ return simulation
+end
+
+function print_progress(simulation)
+ model = simulation.model
+ cfl, dcfl = simulation.parameters
+
+ # Calculate simulation progress in %.
+ progress = 100 * (model.clock.time / simulation.stop_time)
+
+ # Find maximum velocities.
+ vmax = maximum(abs, interior(model.velocities.v))
+ wmax = maximum(abs, interior(model.velocities.w))
+
+ i, t = model.clock.iteration, model.clock.time
+ @info @sprintf("[%06.2f%%] i: %d, t: %.3f, U_max: (%.2e, %.2e), CFL: %.2e, dCFL: %.2e, next Δt: %.2e",
+ progress, i, t, vmax, wmax, cfl(model), dcfl(model), simulation.Δt.Δt)
+
+ return nothing
+end
+
+ simulate_lid_driven_cavity(Re=100, N=128, end_time=15)
+ simulate_lid_driven_cavity(Re=400, N=128, end_time=20)
+ simulate_lid_driven_cavity(Re=1000, N=128, end_time=25)
+ simulate_lid_driven_cavity(Re=3200, N=128, end_time=50)
+ simulate_lid_driven_cavity(Re=5000, N=256, end_time=50)
+ simulate_lid_driven_cavity(Re=7500, N=256, end_time=75)
+ simulate_lid_driven_cavity(Re=10000, N=256, end_time=100)
+
diff --git a/verification/lid_driven_cavity/plot_lid_driven_cavity.py b/verification/lid_driven_cavity/plot_lid_driven_cavity.py
new file mode 100644
index 0000000000..4256b25d69
--- /dev/null
+++ b/verification/lid_driven_cavity/plot_lid_driven_cavity.py
@@ -0,0 +1,138 @@
+import joblib
+import numpy as np
+import xarray as xr
+import matplotlib.pyplot as plt
+import matplotlib.colors as colors
+import cmocean
+import ffmpeg
+
+from numpy import nan
+from xarray.ufuncs import fabs
+
+#####
+##### Data from tables 1, 2, and 4 of Ghia et al. (1982).
+#####
+
+j_Ghia = [1, 8, 9, 10, 14, 23, 37, 59, 65, 80, 95, 110, 123, 124, 125, 126, 129]
+y_Ghia = [0.0000, 0.0547, 0.0625, 0.0703, 0.1016, 0.1719, 0.2813, 0.4531, 0.5000, 0.6172, 0.7344, 0.8516, 0.9531, 0.9609, 0.9688, 0.9766, 1.0000]
+
+v_Ghia = {
+ 100: [0.00000, -0.03717, -0.04192, -0.04775, -0.06434, -0.10150, -0.15662, -0.21090, -0.20581, -0.13641, 0.00332, 0.23151, 0.68717, 0.73722, 0.78871, 0.84123, 1.00000],
+ 400: [0.00000, -0.08186, -0.09266, -0.10338, -0.14612, -0.24299, -0.32726, -0.17119, -0.11477, 0.02135, 0.16256, 0.29093, 0.55892, 0.61756, 0.68439, 0.75837, 1.00000],
+ 1000: [0.00000, -0.18109, -0.20196, -0.22220, -0.29730, -0.38289, -0.27805, -0.10648, -0.06080, 0.05702, 0.18719, 0.33304, 0.46604, 0.51117, 0.57492, 0.65928, 1.00000],
+ 3200: [0.00000, -0.32407, -0.35344, -0.37827, -0.41933, -0.34323, -0.24427, -0.86636, -0.04272, 0.07156, 0.19791, 0.34682, 0.46101, 0.46547, 0.48296, 0.53236, 1.00000],
+ 5000: [0.00000, -0.41165, -0.42901, -0.43643, -0.40435, -0.33050, -0.22855, -0.07404, -0.03039, 0.08183, 0.20087, 0.33556, 0.46036, 0.45992, 0.46120, 0.48223, 1.00000],
+ 7500: [0.00000, -0.43154, -0.43590, -0.43025, -0.38324, -0.32393, -0.23176, -0.07503, -0.03800, 0.08342, 0.20591, 0.34228, 0.47167, 0.47323, 0.47048, 0.47244, 1.00000],
+ 10000: [0.00000, -0.42735, -0.42537, -0.41657, -0.38000, -0.32709, -0.23186, -0.07540, 0.03111, 0.08344, 0.20673, 0.34635, 0.47804, 0.48070, 0.47783, 0.47221, 1.00000]
+}
+
+k_Ghia = [1, 9, 10, 11, 13, 21, 30, 31, 65, 104, 111, 117, 122, 123, 124, 125, 129]
+z_Ghia = [0.0000, 0.0625, 0.0703, 0.0781, 0.0938, 0.1563, 0.2266, 0.2344, 0.5000, 0.8047, 0.8594, 0.9063, 0.9453, 0.9531, 0.9609, 0.9688, 1.0000]
+
+w_Ghia = {
+ 100: [0.00000, 0.09233, 0.10091, 0.10890, 0.12317, 0.16077, 0.17507, 0.17527, 0.05454, -0.24533, -0.22445, -0.16914, -0.10313, -0.08864, -0.07391, -0.05906, 0.00000],
+ 400: [0.00000, 0.18360, 0.19713, 0.20920, 0.22965, 0.28124, 0.30203, 0.30174, 0.05186, -0.38598, -0.44993, -0.23827, -0.22847, -0.19254, -0.15663, -0.12146, 0.00000],
+ 1000: [0.00000, 0.27485, 0.29012, 0.30353, 0.32627, 0.37095, 0.33075, 0.32235, 0.02526, -0.31966, -0.42665, -0.51550, -0.39188, -0.33714, -0.27669, -0.21388, 0.00000],
+ 3200: [0.00000, 0.39560, 0.40917, 0.41906, 0.42768, 0.37119, 0.29030, 0.28188, 0.00999, -0.31184, -0.37401, -0.44307, -0.54053, -0.52357, -0.47425, -0.39017, 0.00000],
+ 5000: [0.00000, 0.42447, 0.43329, 0.43648, 0.42951, 0.35368, 0.28066, 0.27280, 0.00945, -0.30018, -0.36214, -0.41442, -0.52876, -0.55408, -0.55069, -0.49774, 0.00000],
+ 7500: [0.00000, 0.43979, 0.44030, 0.43564, 0.41824, 0.35060, 0.28117, 0.27348, 0.00824, -0.30448, -0.36213, -0.41050, -0.48590, -0.52347, -0.55216, -0.53858, 0.00000],
+ 10000: [0.00000, 0.43983, 0.43733, 0.43124, 0.41487, 0.35070, 0.28003, 0.27224, 0.00831, -0.30719, -0.36737, -0.41496, -0.45863, -0.49099, -0.52987, -0.54302, 0.00000]
+}
+
+y_ζ_Ghia = [0.0000, 0.0625, 0.1250, 0.1875, 0.2500, 0.3125, 0.3750, 0.4375, 0.5000, 0.5625, 0.6250, 0.6875, 0.7500, 0.8125, 0.8750, 0.9375, 1.0000]
+
+ζ_Ghia = {
+ 100: [nan, 40.0110, 22.5378, 16.2862, 12.7844, 10.4199, 8.69628, 7.43218, 6.57451, 6.13973, 6.18946, 6.82674, 8.22110, 10.7414, 15.6591, 30.7923, nan],
+ 400: [nan, 53.6863, 34.6351, 26.5825, 21.0985, 16.8900, 13.7040, 11.4537, 10.0545, 9.38889, 9.34599, 9.88979, 11.2018, 13.9068, 19.6859, 35.0773, nan],
+ 1000: [nan, 42.1124, 23.8707, 18.3120, 16.0458, 14.8061, 14.1374, 14.0928, 14.8901, 16.8350, 20.2666, 25.4341, 32.2953, 40.5437, 51.0557, 75.5980, nan][::-1],
+ 3200: [nan, 49.9664, 34.2327, 30.4779, 27.9514, 25.8572, 24.4639, 24.1457, 25.3889, 28.9413, 35.8795, 47.1443, 61.7864, 75.6401, 89.3391, 126.670, nan][::-1],
+ 5000: [nan, 56.7091, 41.3394, 38.0436, 35.3504, 33.0486, 31.5791, 31.3793, 33.0115, 37.3609, 45.8622, 60.0065, 77.9509, 91.5682, 103.436, 146.702, nan][::-1],
+ 7500: [nan, 61.4046, 50.0769, 46.8901, 43.5641, 40.6123, 38.6951, 38.3834, 40.3982, 45.9128, 56.9345, 75.6334, 98.2364, 111.115, 125.131, 180.927, nan][::-1],
+ 10000: [nan, 66.0352, 57.7756, 54.3725, 50.3792, 46.8672, 44.6303, 44.3287, 46.8271, 53.5905, 67.1400, 90.0231, 116.275, 127.928, 145.073, 209.452, nan][::-1]
+}
+
+def plot_lid_driven_cavity_frame(Re, n):
+ ds = xr.open_dataset(f"lid_driven_cavity_Re{Re}.nc")
+ Ny = ds.yC.size
+ Nz = ds.zC.size
+
+ fig, axes = plt.subplots(nrows=3, ncols=2, figsize=(16, 16), dpi=200)
+ plt.subplots_adjust(hspace=0.25)
+ fig.suptitle(f"Lid-driven cavity, Re = {Re}, t = {ds.time[n].values:.2f}", fontsize=16)
+
+ ax_v_line, ax_v_mesh = axes[0, 0], axes[0, 1]
+ ax_w_line, ax_w_mesh = axes[1, 0], axes[1, 1]
+ ax_ζ_line, ax_ζ_mesh = axes[2, 0], axes[2, 1]
+
+ v_line = ds.v.isel(time=n, yF=Ny//2)
+ ax_v_line.plot(y_Ghia, v_Ghia[Re], label="Ghia et al. (1982)", color="tab:blue", linestyle="", marker="o", fillstyle="none")
+ ax_v_line.plot(ds.zC, v_line.values.flatten(), label="Oceananigans.jl", color="tab:blue")
+ ax_v_line.legend(loc="lower left", bbox_to_anchor=(0, 1.01, 1, 0.2), ncol=2, frameon=False)
+ ax_v_line.set_xlabel("z")
+ ax_v_line.set_ylabel("v")
+ ax_v_line.set_xlim([0, 1])
+
+ w_line = ds.w.isel(time=n, zF=Nz//2)
+ ax_w_line.plot(z_Ghia, w_Ghia[Re], label="Ghia et al. (1982)", color="tab:orange", linestyle="", marker="o", fillstyle="none")
+ ax_w_line.plot(ds.yC, w_line.values.flatten(), label="Oceananigans.jl", color="tab:orange")
+ ax_w_line.legend(loc="lower left", bbox_to_anchor=(0, 1.01, 1, 0.2), ncol=2, frameon=False)
+ ax_w_line.set_xlabel("y")
+ ax_w_line.set_ylabel("w")
+ ax_w_line.set_xlim([0, 1])
+
+ v = ds.v.isel(time=n).squeeze()
+ img_v = v.plot.pcolormesh(ax=ax_v_mesh, vmin=-1, vmax=1, cmap=cmocean.cm.balance, add_colorbar=False)
+ fig.colorbar(img_v, ax=ax_v_mesh, extend="both")
+ ax_v_mesh.axvline(x=0.5, color="tab:blue", alpha=0.5)
+ ax_v_mesh.set_title("v-velocity")
+ ax_v_mesh.set_xlabel("y")
+ ax_v_mesh.set_ylabel("z")
+ ax_v_mesh.set_aspect("equal")
+
+ w = ds.w.isel(time=n).squeeze()
+ img_w = w.plot.pcolormesh(ax=ax_w_mesh, vmin=-1, vmax=1, cmap=cmocean.cm.balance, add_colorbar=False)
+ fig.colorbar(img_w, ax=ax_w_mesh, extend="both")
+ ax_w_mesh.axhline(y=0.5, color="tab:orange", alpha=0.5)
+ ax_w_mesh.set_title("w-velocity")
+ ax_w_mesh.set_xlabel("y")
+ ax_w_mesh.set_ylabel("z")
+ ax_w_mesh.set_aspect("equal")
+
+ ζ_line = fabs(ds.ζ.isel(time=n, zF=Nz))
+ ax_ζ_line.plot(y_ζ_Ghia, ζ_Ghia[Re], label="Ghia et al. (1982)", color="tab:purple", linestyle="", marker="o", fillstyle="none")
+ ax_ζ_line.plot(ds.yF, ζ_line.values.flatten(), label="Oceananigans.jl", color="tab:purple")
+ ax_ζ_line.legend(loc="lower left", bbox_to_anchor=(0, 1.01, 1, 0.2), ncol=2, frameon=False)
+ ax_ζ_line.set_xlabel("y")
+ ax_ζ_line.set_ylabel("vorticity $|\zeta$|")
+ ax_ζ_line.set_xlim([0, 1])
+ ax_ζ_line.set_ylim(bottom=0)
+
+ ζ = ds.ζ.isel(time=n).squeeze()
+ img_ζ = ζ.plot.pcolormesh(ax=ax_ζ_mesh, cmap=cmocean.cm.curl, extend="both", add_colorbar=False,
+ norm=colors.SymLogNorm(base=10, linthresh=1, vmin=-1e2, vmax=1e2))
+ fig.colorbar(img_ζ, ax=ax_ζ_mesh, extend="both")
+ ax_ζ_mesh.axhline(y=1, color="tab:purple", alpha=1.0)
+ ax_ζ_mesh.set_title("vorticity")
+ ax_ζ_mesh.set_xlabel("y")
+ ax_ζ_mesh.set_ylabel("z")
+ ax_ζ_mesh.set_aspect("equal")
+
+ print(f"Saving lid-driven cavity Re={Re} frame {n}/{ds.time.size-1}...")
+ plt.savefig(f"lid_driven_cavity_Re{Re}_{n:05d}.png")
+ plt.close("all")
+
+for Re in [100, 400, 1000, 3200, 5000, 7500, 10000]:
+ ds = xr.open_dataset(f"lid_driven_cavity_Re{Re}.nc")
+
+ joblib.Parallel(n_jobs=-1)(
+ joblib.delayed(plot_lid_driven_cavity_frame)(Re, n)
+ for n in range(ds.time.size)
+ )
+
+ (
+ ffmpeg
+ .input(f"lid_driven_cavity_Re{Re}_%05d.png", framerate=30)
+ .output(f"lid_driven_cavity_Re{Re}.mp4", crf=15, pix_fmt='yuv420p')
+ .overwrite_output()
+ .run()
+ )