forked from agdestein/IncompressibleNavierStokes.jl
-
Notifications
You must be signed in to change notification settings - Fork 0
/
TaylorGreenVortex2D.jl
100 lines (91 loc) · 2.5 KB
/
TaylorGreenVortex2D.jl
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
# # Convergence study: Taylor-Green vortex (2D)
#
# In this example we consider the Taylor-Green vortex.
# In 2D, it has an analytical solution, given by
#
# ```math
# \begin{split}
# u^1(x, y, t) & = - \sin(x) \cos(y) \mathrm{e}^{-2 t / Re} \\
# u^2(x, y, t) & = + \cos(x) \sin(y) \mathrm{e}^{-2 t / Re}
# \end{split}
# ```
#
# This allows us to test the convergence of our solver.
#md using CairoMakie
using GLMakie #!md
using IncompressibleNavierStokes
using LinearAlgebra
# Output directory
outdir = joinpath(@__DIR__, "output", "TaylorGreenVortex2D")
ispath(outdir) || mkpath(outdir)
# Convergence
"""
Compare numerical solution with analytical solution at final time.
"""
function compute_convergence(; D, nlist, lims, Re, tlims, Δt, uref, backend = CPU())
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, backend)
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 = sum(abs2, u[Ip, :] - ut[Ip, :])
b = sum(abs2, ut[Ip, :])
e[i] = sqrt(a) / sqrt(b)
end
e
end
# Analytical solution for 2D Taylor-Green vortex
solution(Re) =
(dim, x, y, t) -> (dim == 1 ? -sin(x) * cos(y) : cos(x) * sin(y)) * exp(-2t / Re)
# Compute error for different resolutions
Re = 2.0e3
nlist = [2, 4, 8, 16, 32, 64, 128, 256]
e = compute_convergence(;
D = 2,
nlist,
lims = (0.0, 2π),
Re,
tlims = (0.0, 2.0),
Δt = 0.01,
uref = solution(Re),
)
# Plot convergence
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
save(joinpath(outdir, "convergence.png"), fig)
#md # ## Copy-pasteable code
#md #
#md # Below is the full code for this example stripped of comments and output.
#md #
#md # ```julia
#md # CODE_CONTENT
#md # ```