Note: Output is not generated for this example (to save resources on GitHub).

Kolmogorov flow (2D)

The Kolmogorov flow in a periodic box $\mathrm{\Omega }=[0,1{]}^2$ is initiated via the force field

$$f(x,y)=\left(\begin{array}{c}\sin (\pi ky)\\ 0\end{array}\right)$$

where k is the wavenumber where energy is injected.

Packages

We just need the IncompressibleNavierStokes and a Makie plotting package.

using CairoMakie
using IncompressibleNavierStokes
# using CUDA

Problem setup

Define a uniform grid with a steady body force field.

n = 256
axis = range(0.0, 1.0, n + 1)
setup = Setup(;
    x = (axis, axis),
    boundary_conditions = (;
        u = ((PeriodicBC(), PeriodicBC()), (PeriodicBC(), PeriodicBC()))
    ),
    # backend = CUDABackend(),
);
u = random_field(setup, 0.0; A = 1e-2);
nothing #hide

This is the right-hand side force in the momentum equation By default, it is just navierstokes!. Here we add a pre-computed body force.

function force!(f, state, t; setup, cache, viscosity)
    navierstokes!(f, state, t; setup, cache, viscosity)
    f.u .+= cache.bodyforce
end

Tell IncompressibleNavierStokes how to prepare the cache for force!. The cache is created before time stepping begins.

function IncompressibleNavierStokes.get_cache(::typeof(force!), setup)
    f(dim, x, y) = (dim == 1) * 5 * sinpi(8 * y)
    bodyforce = velocityfield(setup, f; doproject = false)
    (; bodyforce)
end

Plot body force

Since the force is steady, it is just stored as a field.

let
    (; bodyforce) = IncompressibleNavierStokes.get_cache(force!, setup)
    bodyforce[:, :, 1] |> Array |> heatmap
end

Solve unsteady problem

state, outputs = solve_unsteady(;
    setup,
    force!,
    start = (; u),
    tlims = (0.0, 2.0),
    params = (; viscosity = 5e-4),
    processors = (
        rtp = realtimeplotter(; setup, nupdate = 10),
        ehist = realtimeplotter(;
            setup,
            plot = energy_history_plot,
            nupdate = 10,
            displayfig = false,
        ),
        espec = realtimeplotter(;
            setup,
            plot = energy_spectrum_plot,
            nupdate = 10,
            displayfig = false,
        ),
        log = timelogger(; nupdate = 100),
    ),
);
nothing #hide

Field plot

outputs.rtp

Energy history

outputs.ehist

Energy spectrum

outputs.espec

Copy-pasteable code

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

using WGLMakie
using IncompressibleNavierStokes
# using CUDA

n = 256
axis = range(0.0, 1.0, n + 1)
setup = Setup(;
    x = (axis, axis),
    boundary_conditions = (;
        u = ((PeriodicBC(), PeriodicBC()), (PeriodicBC(), PeriodicBC()))
    ),
    # backend = CUDABackend(),
);
u = random_field(setup, 0.0; A = 1e-2);

function force!(f, state, t; setup, cache, viscosity)
    navierstokes!(f, state, t; setup, cache, viscosity)
    f.u .+= cache.bodyforce
end

function IncompressibleNavierStokes.get_cache(::typeof(force!), setup)
    f(dim, x, y) = (dim == 1) * 5 * sinpi(8 * y)
    bodyforce = velocityfield(setup, f; doproject = false)
    (; bodyforce)
end

let
    (; bodyforce) = IncompressibleNavierStokes.get_cache(force!, setup)
    bodyforce[:, :, 1] |> Array |> heatmap
end

state, outputs = solve_unsteady(;
    setup,
    force!,
    start = (; u),
    tlims = (0.0, 2.0),
    params = (; viscosity = 5e-4),
    processors = (
        rtp = realtimeplotter(; setup, nupdate = 10),
        ehist = realtimeplotter(;
            setup,
            plot = energy_history_plot,
            nupdate = 10,
            displayfig = false,
        ),
        espec = realtimeplotter(;
            setup,
            plot = energy_spectrum_plot,
            nupdate = 10,
            displayfig = false,
        ),
        log = timelogger(; nupdate = 100),
    ),
);

outputs.rtp

outputs.ehist

outputs.espec

---

This page was generated using Literate.jl.