Note: Output is not generated for this example (to save resources on GitHub).
Tutorial: Lid-Driven Cavity - 2D
In this example we consider a box with a moving lid. The velocity is initially at rest. The solution should reach at steady state equilibrium after a certain time. The same steady state should be obtained when solving a steady state problem.
We start by loading packages. A Makie plotting backend is needed
for plotting. GLMakie
creates an interactive window (useful for real-time plotting), but does not work when building this example on GitHub. CairoMakie
makes high-quality static vector-graphics plots.
using CairoMakie
using IncompressibleNavierStokes
Case name for saving results
outdir = joinpath(@__DIR__, "output", "LidDrivenCavity2D")
The code allows for using different floating point number types, including single precision (Float32
) and double precision (Float64
). On the CPU, the speed is not really different, but double precision uses twice as much memory as single precision. When running on the GPU, single precision is preferred. Half precision (Float16
) is also an option, but then the values should be scaled judiciously to avoid vanishing digits when applying differential operators of the form "right minus left divided by small distance".
Note how floating point type hygiene is enforced in the following using T
to avoid mixing different precisions.
T = Float64
T = Float32
# T = Float16
We can also choose to do the computations on a different device. By default, the computations are performed on the host (CPU). An optional ArrayType
allows for moving arrays to a different device such as a GPU.
Note: For GPUs, single precision is preferred.
ArrayType = Array
# using CUDA; ArrayType = CuArray
# using AMDGPU; ArrayType = ROCArray
# using oneAPI; ArrayType = oneArray
# using Metal; ArrayType = MtlArray
Here we choose a moderate Reynolds number. Note how we pass the floating point type.
Re = T(1_000)
Non-zero Dirichlet boundary conditions are specified as plain Julia functions.
U = (T(1), T(0))
boundary_conditions = (
# x left, x right
(DirichletBC(), DirichletBC()),
# y bottom, y top
(DirichletBC(), DirichletBC(U)),
)
We create a two-dimensional domain with a box of size [1, 1]
. The grid is created as a Cartesian product between two vectors. We add a refinement near the walls.
n = 32
lims = T(0), T(1)
x = cosine_grid(lims..., n), cosine_grid(lims..., n)
plotgrid(x...)
We can now build the setup and assemble operators. A 3D setup is built if we also provide a vector of z-coordinates.
setup = Setup(; x, boundary_conditions, Re, ArrayType);
nothing #hide
The initial conditions are provided in function. The value dim()
determines the velocity component.
ustart = velocityfield(setup, (dim, x, y) -> zero(x));
nothing #hide
Iteration processors are called after every nupdate
time steps. This can be useful for logging, plotting, or saving results. Their respective outputs are later returned by solve_unsteady
.
processors = (
# rtp = realtimeplotter(; setup, plot = fieldplot, nupdate = 50),
# ehist = realtimeplotter(; setup, plot = energy_history_plot, nupdate = 10),
# espec = realtimeplotter(; setup, plot = energy_spectrum_plot, nupdate = 10),
# anim = animator(; setup, path = "$outdir/solution.mkv", nupdate = 20),
# vtk = vtk_writer(; setup, nupdate = 100, dir = outdir, filename = "solution"),
# field = fieldsaver(; setup, nupdate = 10),
log = timelogger(; nupdate = 1000),
);
nothing #hide
By default, a standard fourth order Runge-Kutta method is used. If we don't provide the time step explicitly, an adaptive time step is used.
tlims = (T(0), T(10))
state, outputs = solve_unsteady(; setup, ustart, tlims, Δt = T(1e-3), processors);
nothing #hide
Post-process
We may visualize or export the computed fields
Export fields to VTK. The file outdir/solution.vti
may be opened for visualization in ParaView. This is particularly useful for inspecting results from 3D simulations.
save_vtk(state; setup, filename = joinpath(outdir, "solution"))
Plot pressure
fieldplot(state; setup, fieldname = :pressure)
Plot velocity. Note the time stamp used for computing boundary conditions, if any.
fieldplot(state; setup, fieldname = :velocitynorm)
Plot vorticity
fieldplot(state; setup, fieldname = :vorticity)
In addition, the named tuple outputs
contains quantities from our processors. The logger returns nothing.
# outputs.rtp
# outputs.ehist
# outputs.espec
# outputs.anim
# outputs.vtk
# outputs.field
outputs.log
This page was generated using Literate.jl.