IncompressibleNavierStokes.jl

Appearance

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

Taylor-Green vortex - 3D ​

In this example we consider the Taylor-Green vortex.

julia
using CairoMakie
using IncompressibleNavierStokes

Floating point precision

julia
T = Float64

Backend ​

Running in 3D is heavier than in 2D. If you are running this on a CPU, consider using multiple threads by starting Julia with julia -t auto, or add "julia.NumThreads": "auto" to the settings in VSCode.

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

Setup ​

julia
n = 32
r = range(T(0), T(1), n + 1)
setup = Setup(; x = (r, r, r), Re = T(1e3), backend);
psolver = psolver_spectral(setup);
nothing #hide

Initial conditions

julia
U(dim, x, y, z) =
    if dim == 1
        sinpi(2x) * cospi(2y) * sinpi(2z) / 2
    elseif dim == 2
        -cospi(2x) * sinpi(2y) * sinpi(2z) / 2
    else
        zero(x)
    end
ustart = velocityfield(setup, U, psolver);
nothing #hide

Solve unsteady problem ​

julia
state, outputs = solve_unsteady(;
    setup,
    ustart,
    tlims = (T(0), T(1.0)),
    Δt = T(1e-3),
    processors = (
        # rtp = realtimeplotter(; setup, plot = fieldplot, nupdate = 10),
        ehist = realtimeplotter(;
            setup,
            plot = energy_history_plot,
            nupdate = 1,
            displayfig = false,
        ),
        espec = realtimeplotter(; setup, plot = energy_spectrum_plot, nupdate = 10),
        # anim = animator(; setup, path = "$outdir/solution.mkv", nupdate = 20),
        # vtk = vtk_writer(; setup, nupdate = 10, dir = outdir, filename = "solution"),
        log = timelogger(; nupdate = 100),
    ),
    psolver,
);
nothing #hide

Post-process ​

Energy history

julia
outputs.ehist

Energy spectrum

julia
outputs.espec

Copy-pasteable code ​

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

julia
using GLMakie
using IncompressibleNavierStokes

T = Float64

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

n = 32
r = range(T(0), T(1), n + 1)
setup = Setup(; x = (r, r, r), Re = T(1e3), backend);
psolver = psolver_spectral(setup);

U(dim, x, y, z) =
    if dim == 1
        sinpi(2x) * cospi(2y) * sinpi(2z) / 2
    elseif dim == 2
        -cospi(2x) * sinpi(2y) * sinpi(2z) / 2
    else
        zero(x)
    end
ustart = velocityfield(setup, U, psolver);

state, outputs = solve_unsteady(;
    setup,
    ustart,
    tlims = (T(0), T(1.0)),
    Δt = T(1e-3),
    processors = (
        # rtp = realtimeplotter(; setup, plot = fieldplot, nupdate = 10),
        ehist = realtimeplotter(;
            setup,
            plot = energy_history_plot,
            nupdate = 1,
            displayfig = false,
        ),
        espec = realtimeplotter(; setup, plot = energy_spectrum_plot, nupdate = 10),
        # anim = animator(; setup, path = "$outdir/solution.mkv", nupdate = 20),
        # vtk = vtk_writer(; setup, nupdate = 10, dir = outdir, filename = "solution"),
        log = timelogger(; nupdate = 100),
    ),
    psolver,
);

outputs.ehist

outputs.espec

This page was generated using Literate.jl.