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

output = "output/LidDrivenCavity2D"
"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
Float32

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
Array

Here we choose a moderate Reynolds number. Note how we pass the floating point type.

Re = T(1_000)
1000.0f0

Non-zero Dirichlet boundary conditions are specified as plain Julia functions. Note that time derivatives are required.

boundary_conditions = (
    # x left, x right
    (DirichletBC(), DirichletBC()),

    # y bottom, y top
    (
        DirichletBC(),
        DirichletBC(
            (dim, x, y, t) -> dim() == 1 ? one(x) : zero(x),
            (dim, x, y, t) -> zero(x),
        ),
    ),
)
((DirichletBC{Nothing, Nothing}(nothing, nothing), DirichletBC{Nothing, Nothing}(nothing, nothing)), (DirichletBC{Nothing, Nothing}(nothing, nothing), DirichletBC{Main.var"#1#3", Main.var"#2#4"}(Main.var"#1#3"(), Main.var"#2#4"())))

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)
y = cosine_grid(lims..., n)
plotgrid(x, y)
Example block output

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, y; boundary_conditions, Re, ArrayType);

The initial conditions are provided in function. The value dim() determines the velocity component.

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

Solve problems

Problems can be solved.

The solve_steady_state function is for computing a state where the right hand side of the momentum equation is zero.

# u, p = solve_steady_state(setup, u₀, p₀)
nothing

For this test case, the same steady state may be obtained by solving an unsteady problem for a sufficiently long time.

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 = "$output/solution.mkv", nupdate = 20),
    # vtk = vtk_writer(; setup, nupdate = 100, dir = output, filename = "solution"),
    # field = fieldsaver(; setup, nupdate = 10),
    log = timelogger(; nupdate = 1000),
);

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);
Iteration 1000	t = 0.999991	Δt = 0.000999987	umax = 1
Iteration 2000	t = 2.00004	Δt = 0.00100005	umax = 1
Iteration 3000	t = 2.99996	Δt = 0.000999928	umax = 1
Iteration 4000	t = 3.99989	Δt = 0.000999928	umax = 1
Iteration 5000	t = 4.99982	Δt = 0.000999928	umax = 1
Iteration 6000	t = 5.99975	Δt = 0.000999928	umax = 1
Iteration 7000	t = 6.99968	Δt = 0.000999928	umax = 1
Iteration 8000	t = 7.9996	Δt = 0.000999928	umax = 1
Iteration 9000	t = 9.00001	Δt = 0.0010004	umax = 1
Iteration 10000	t = 10.0004	Δt = 0.0010004	umax = 1

Post-process

We may visualize or export the computed fields

Export fields to VTK. The file output/solution.vti may be opened for visualization in ParaView. This is particularly useful for inspecting results from 3D simulations.

save_vtk(setup, state.u, state.t, "$output/solution")
1-element Vector{String}:
 "output/LidDrivenCavity2D/solution.vtr"

Plot pressure

fieldplot(state; setup, fieldname = :pressure)
Example block output

Plot velocity. Note the time stamp used for computing boundary conditions, if any.

fieldplot(state; setup, fieldname = :velocitynorm)
Example block output

Plot vorticity

fieldplot(state; setup, fieldname = :vorticity)
Example block output

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.