Decaying Homogeneous Isotropic Turbulence - 2D 
In this example we consider decaying homogeneous isotropic turbulence, similar to the cases considered in [1] and [2]. The initial velocity field is created randomly, but with a specific energy spectrum. Due to viscous dissipation, the turbulent features eventually group to form larger visible eddies.
Packages 
We just need IncompressibleNavierStokes and a Makie plotting backend.
julia
using CairoMakie
using IncompressibleNavierStokesSetup
julia
n = 256
ax = LinRange(0.0, 1.0, n + 1)
setup = Setup(; x = (ax, ax), Re = 4e3);
ustart = random_field(setup, 0.0);Solve unsteady problem
julia
state, outputs = solve_unsteady(;
    setup,
    ustart,
    tlims = (0.0, 1.0),
    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),
    ),
);[ Info: t = 0.103785	Δt = 0.001	umax = 3.4	itertime = 0.027
[ Info: t = 0.205464	Δt = 0.0012	umax = 3	itertime = 0.016
[ Info: t = 0.331177	Δt = 0.0013	umax = 2.7	itertime = 0.016
[ Info: t = 0.460566	Δt = 0.0012	umax = 2.9	itertime = 0.016
[ Info: t = 0.5867	Δt = 0.0013	umax = 2.8	itertime = 0.016
[ Info: t = 0.7122	Δt = 0.0015	umax = 2.4	itertime = 0.016
[ Info: t = 0.865263	Δt = 0.0014	umax = 2.4	itertime = 0.016Post-process 
We may visualize or export the computed fields
Energy history
julia
outputs.ehist
Energy spectrum
julia
outputs.espec
Plot field
julia
fieldplot(state; setup)
Copy-pasteable code 
Below is the full code for this example stripped of comments and output.
julia
using GLMakie
using IncompressibleNavierStokes
n = 256
ax = LinRange(0.0, 1.0, n + 1)
setup = Setup(; x = (ax, ax), Re = 4e3);
ustart = random_field(setup, 0.0);
state, outputs = solve_unsteady(;
    setup,
    ustart,
    tlims = (0.0, 1.0),
    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.ehist
outputs.espec
fieldplot(state; setup)This page was generated using Literate.jl.