-
Notifications
You must be signed in to change notification settings - Fork 0
/
plot_exact_order2.jl
82 lines (77 loc) · 5.18 KB
/
plot_exact_order2.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
# Copyright (c) 2024 Ashwin. S. Nayak, Andrés Prieto, Daniel Fernández Comesaña
#
# This file is part of pressure-projections
#
# SPDX-License-Identifier: MIT
"""
Computation and plot of the error curves obtained from the post-processing projection techniques
when the displacement field is computed from the interpolation of the exact solution in second-order FE discretizations
"""
# Load variables from JDL2 file
using JLD2
file = "results/convergence_data_from_exact_k=50_orderFEpostprocessing=2.jld2"
vars = load(file)
k = vars["k"]
orderFE_postprocessing = vars["orderFE_postprocessing"]
N_values = 1.0./vars["h_values"]
error_L2_u = vars["error_L2_u"]
error_H1_u = vars["error_H1_u"]
error_L2_p = vars["error_L2_p"]
error_H1_p = vars["error_H1_p"]
error_L2_proj_DG = vars["error_L2_proj_DG"]
error_H1_proj_DG = vars["error_H1_proj_DG"]
error_L2_proj_CG = vars["error_L2_proj_CG"]
error_H1_proj_CG = vars["error_H1_proj_CG"]
error_L2_proj_H1 = vars["error_L2_proj_H1"]
error_H1_proj_H1 = vars["error_H1_proj_H1"]
error_L2_proj_H2_DG = vars["error_L2_proj_H2_DG"]
error_H1_proj_H2_DG = vars["error_H1_proj_H2_DG"]
# Plot the relative errors L2 and H1 with latex labels and log scale in Makie
using CairoMakie, LaTeXStrings
fig = Figure(size = (600, 400), fonts = (; regular="CMU Serif")) ## probably you need to install this font in your system
ax = Axis(fig[1, 1], xlabel = L"$\lambda/h$", ylabel = L"Relative error (%)$$", xtickalign = 1,
xticksize = 10, ytickalign = 1, yticksize = 10, xscale=log10, yscale=log10, #xticks = vcat([2*pi/k.*N_values[end]], 10 .^range(-3, stop=-1, length=3)), yticks = 10 .^range(-5, stop=2, length=8),
yminorticksvisible = true, yminorgridvisible = true, yminorticks = IntervalsBetween(9),
xminorticksvisible = true, xminorgridvisible = true, xminorticks = IntervalsBetween(9))
# Error u
scatterlines!(2*pi/k.*N_values, error_L2_u, color = :black, linestyle = :dash , label = L"$\mathrm{L}^2$-error $I_{h,\mathrm{RT}_{1}}$", marker = '△', markersize = 10)
scatterlines!(2*pi/k.*N_values, error_H1_u, color = :black, label = L"$\mathrm{H(div)}$-error $I_{h,\mathrm{RT}_{1}}$", marker = '△', markersize = 10)
# Error p
scatterlines!(2*pi/k.*N_values, error_L2_p, color = :green, linestyle = :dash, label = L"$\mathrm{L}^2$-error $I_{h,\mathrm{CG}_{2}}$", marker = '◇', markersize = 10)
scatterlines!(2*pi/k.*N_values, error_H1_p, color = :green, label = L"$\mathrm{H}^1$-error $I_{h,\mathrm{CG}_{2}}$", marker = '◇', markersize = 10)
# Error L2 DG projection
scatterlines!(2*pi/k.*N_values, error_L2_proj_DG, color = :orange, linestyle = :dash, label = L"$\mathrm{L}^2$-error $\Pi_{h,\mathrm{DG}_{1}}^{\mathrm{L}^2}$", marker = '▷', markersize = 10)
scatterlines!(2*pi/k.*N_values, error_H1_proj_DG, color = :orange, label = L"$\mathrm{H}^1$-error $\Pi_{h,\mathrm{DG}_{1}}^{\mathrm{L}^2}$", marker = '▷', markersize = 10)
# Error L2 CG projection
scatterlines!(2*pi/k.*N_values, error_L2_proj_CG, color = :blue, linestyle = :dash, label = L"$\mathrm{L}^2$-error $\Pi_{h,\mathrm{CG}_{2}}^{\mathrm{L}^2}$", marker = '□', markersize = 10)
scatterlines!(2*pi/k.*N_values, error_H1_proj_CG, color = :blue, label = L"$\mathrm{H}^1$-error $\Pi_{h,\mathrm{CG}_{2}}^{\mathrm{L}^2}$", marker = '□', markersize = 10)
# Error H1 CG projection
scatterlines!(2*pi/k.*N_values, error_L2_proj_H1, color = :red, linestyle = :dash, label = L"$\mathrm{L}^2$-error $\Pi_{h,\mathrm{CG}_{2}}^{\mathrm{H}^1}$", marker = '○', markersize = 10)
scatterlines!(2*pi/k.*N_values, error_H1_proj_H1, color = :red, label = L"$\mathrm{H}^1$-error $\Pi_{h,\mathrm{CG}_{2}}^{\mathrm{H}^1}$", marker = '○', markersize = 10)
# Error H2 DG projection
scatterlines!(2*pi/k.*N_values, error_L2_proj_H2_DG, color = :purple, linestyle = :dash, label = L"$\mathrm{L}^2$-error $\Pi_{h,\mathrm{DG}_{2}}^{\mathrm{H}^1}$", marker = '☆', markersize = 10)
scatterlines!(2*pi/k.*N_values, error_H1_proj_H2_DG, color = :purple, label = L"$\mathrm{H}^1$-error $\Pi_{h,\mathrm{DG}_{2}}^{\mathrm{H}^1}$", marker = '☆', markersize = 10)
# axislegend(; position = :rb, nbanks = 1, framecolor = (:grey, 0.5))
Legend(fig[1, 2], ax, halign = :right, valign = :top, labelsize = 17, markersize = 10, markerstrokewidth = 0.5)
xlims!(ax, 2*pi/k.*N_values[end], 2*pi/k.*N_values[4])
ylims!(ax, 1e-3, 1e2)
# Plots triangles for order of convergence
order=3
x0=10.0^1.4; y0=3e-3
x1=10.0; y1=y0*(x0/x1)^order
lines!([(x0, y0), (x1, y0), (x1, y1), (x0, y0)], color=:black, linewidth=0.5)
text!([(10.0^((log10(x0)+log10(x1))/2), 10.0^(0.8*log10(y0)+0.2*log10(y1)))], text=L"O(h^3)", color=:black, fontsize=10)
# Triangle inverted
order=2
x0=10.0^1.4; y0=1.5
x1=10.0; y1=y0*(x1/x0)^order
lines!([(x0, y0), (x1, y0), (x0, y1), (x0, y0)], color=:black, linewidth=0.5)
text!([(10.0^(0.95*log10(x0)+0.05*log10(x1)), 10.0^(0.5*log10(y1)+0.5*log10(y0)))], text=L"O(h^2)", color=:black, fontsize=10)
# Triangle inverted
order=1
x0=10.0^1.4; y0=15.0
x1=10.0; y1=y0*(x1/x0)^order
lines!([(x0, y0), (x1, y0), (x0, y1), (x0, y0)], color=:black, linewidth=0.5)
text!([(10.0^(0.95*log10(x0)+0.05*log10(x1)), 10.0^(0.6*log10(y1)+0.4*log10(y0)))], text=L"O(h)", color=:black, fontsize=10)
fig
save("results/convergence_plot_from_exact_k=$(Integer(k))_orderFEpostprocessing=$(orderFE_postprocessing).pdf", fig)