Skip to content

Convergence study: Taylor-Green vortex (2D)

In this example we consider the Taylor-Green vortex. In 2D, it has an analytical solution, given by

u1(x,y,t)=sin(x)cos(y)e2t/Reu2(x,y,t)=+cos(x)sin(y)e2t/Re

This allows us to test the convergence of our solver.

julia
using CairoMakie
using IncompressibleNavierStokes
using LinearAlgebra

Output directory

julia
outdir = joinpath(@__DIR__, "output", "TaylorGreenVortex2D")
ispath(outdir) || mkpath(outdir)
"/home/runner/work/IncompressibleNavierStokes.jl/IncompressibleNavierStokes.jl/docs/build/examples/generated/output/TaylorGreenVortex2D"

Convergence

julia
"""
Compare numerical solution with analytical solution at final time.
"""
function compute_convergence(; D, nlist, lims, Re, tlims, Δt, uref, ArrayType = Array)
    T = typeof(lims[1])
    e = zeros(T, length(nlist))
    for (i, n) in enumerate(nlist)
        @info "Computing error for n = $n"
        x = ntuple-> LinRange(lims..., n + 1), D)
        setup = Setup(; x, Re, ArrayType)
        psolver = psolver_spectral(setup)
        ustart = velocityfield(
            setup,
            (dim, x...) -> uref(dim, x..., tlims[1]),
            tlims[1];
            psolver,
        )
        ut = velocityfield(
            setup,
            (dim, x...) -> uref(dim, x..., tlims[2]),
            tlims[2];
            psolver,
            doproject = false,
        )
        (; u, t), outputs = solve_unsteady(; setup, ustart, tlims, Δt, psolver)
        (; Ip) = setup.grid
        a, b = T(0), T(0)
        for α = 1:D
            a += sum(abs2, u[α][Ip] - ut[α][Ip])
            b += sum(abs2, ut[α][Ip])
        end
        e[i] = sqrt(a) / sqrt(b)
    end
    e
end
Main.compute_convergence

Analytical solution for 2D Taylor-Green vortex

julia
solution(Re) =
    (dim, x, y, t) -> (dim == 1 ? -sin(x) * cos(y) : cos(x) * sin(y)) * exp(-2t / Re)
solution (generic function with 1 method)

Compute error for different resolutions

julia
Re = 2.0e3
nlist = [2, 4, 8, 16, 32, 64, 128, 256]
e = compute_convergence(;
    D = 2,
    nlist,
    lims = (0.0, ),
    Re,
    tlims = (0.0, 2.0),
    Δt = 0.01,
    uref = solution(Re),
)
8-element Vector{Float64}:
 0.4242656052086717
 0.000378932838635171
 0.00010072266508321756
 2.557066539066516e-5
 6.417292500660421e-6
 1.6058662108541681e-6
 4.015630574104089e-7
 1.0039679673157669e-7

Plot convergence

julia
fig = Figure()
ax = Axis(
    fig[1, 1];
    xscale = log10,
    yscale = log10,
    xticks = nlist,
    xlabel = "n",
    title = "Relative error",
)
scatterlines!(ax, nlist, e; label = "Data")
lines!(ax, collect(extrema(nlist)), n -> n^-2.0; linestyle = :dash, label = "n^-2")
axislegend(ax)
fig

Save figure

julia
save(joinpath(outdir, "convergence.png"), fig)
CairoMakie.Screen{IMAGE}

This page was generated using Literate.jl.