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 backend allows for moving arrays to a different device such as a GPU.

Note: For GPUs, single precision is preferred.

backend = CPU()
# using CUDA; backend = CUDABackend()

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, backend);
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

Copy-pasteable code

Below is the full code for this example stripped of comments and output.

using GLMakie
using IncompressibleNavierStokes

outdir = joinpath(@__DIR__, "output", "LidDrivenCavity2D")

T = Float32
# T = Float16

backend = CPU()
# using CUDA; backend = CUDABackend()

Re = T(1_000)

U = (T(1), T(0))
boundary_conditions = (
    # x left, x right
    (DirichletBC(), DirichletBC()),

    # y bottom, y top
    (DirichletBC(), DirichletBC(U)),
)

n = 32
lims = T(0), T(1)
x = cosine_grid(lims..., n), cosine_grid(lims..., n)
plotgrid(x...)

setup = Setup(; x, boundary_conditions, Re, backend);

ustart = velocityfield(setup, (dim, x, y) -> zero(x));

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),
);

tlims = (T(0), T(10))
state, outputs = solve_unsteady(; setup, ustart, tlims, Δt = T(1e-3), processors);

save_vtk(state; setup, filename = joinpath(outdir, "solution"))

fieldplot(state; setup, fieldname = :pressure)

fieldplot(state; setup, fieldname = :velocitynorm)

fieldplot(state; setup, fieldname = :vorticity)

# outputs.rtp
# outputs.ehist
# outputs.espec
# outputs.anim
# outputs.vtk
# outputs.field
outputs.log

---

This page was generated using Literate.jl.