-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathconvergence_plot_from_fem_order1.jl
66 lines (57 loc) · 3.27 KB
/
convergence_plot_from_fem_order1.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
# 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 first-order Raviart-Thomas FE discretizations
"""
# Load the package to use the generate_mesh function
include("meshing.jl")
# Load the computing.jl file to compute the relative errors
include("../computing.jl")
# Define the range of N values: logspace between 1e1 and 1e4
N_values = Integer.(round.(10.0.^range(log10(25), stop=log10(600), length=15)))
k=50.0 # wave number (2, 20, 100)
orderFE = 1 # order of the finite element
orderFE_postprocessing = 1 # order of the finite element for postprocessing
# Define the exact solution# Compute the analytical solution and the displacement error
p0 = 1.0 # Pressure amplitude
Zs = 1.0 # Specific surface impedance
pex(x) = p0 * exp(1im*k*x[2]) # =-divergence(uex(x))
uex(x) = VectorValue(0.0, 1im * p0 / k * exp(1im*k*x[2]))
grad_pex(x) = VectorValue(0.0, 1im * p0 * k * exp(1im*k*x[2])) # =k^2*uex(x)
exact_solution = Dict("ExactPressure"=>pex, "ExactDisplacement"=>uex, "ExactGradientPressure"=>grad_pex)
# Dictionaonary with the boundary data associated to each boundary tag
uD(x) = VectorValue(0.0, 0.0)
pD(x) = p0
boundary_data = Dict("RigidBoundary"=>uD, "WallPressure"=>pD, "AbsorbingBoundary"=>Zs)
# Define the arrays to store the complex-valued relative errors
error_L2_u = zeros(size(N_values, 1))
error_H1_u = zeros(size(N_values, 1))
error_L2_p = zeros(size(N_values, 1))
error_H1_p = zeros(size(N_values, 1))
error_L2_proj_DG = zeros(size(N_values, 1))
error_H1_proj_DG = zeros(size(N_values, 1))
error_L2_proj_CG = zeros(size(N_values, 1))
error_H1_proj_CG = zeros(size(N_values, 1))
error_L2_proj_H1 = zeros(size(N_values, 1))
error_H1_proj_H1 = zeros(size(N_values, 1))
error_L2_proj_H2_DG = zeros(size(N_values, 1))
error_H1_proj_H2_DG = zeros(size(N_values, 1))
# Perform sweeping in h_values with a fixed k value
for j in 1:size(N_values, 1)
# Compute the relative errors
error_L2_u[j], error_H1_u[j], error_L2_p[j], error_H1_p[j], error_L2_proj_DG[j], error_H1_proj_DG[j], error_L2_proj_CG[j], error_H1_proj_CG[j], error_L2_proj_H1[j], error_H1_proj_H1[j], error_L2_proj_H2_DG[j], error_H1_proj_H2_DG[j] = compute_from_fem(k, N_values[j], orderFE, orderFE_postprocessing, exact_solution, boundary_data)
end
h_values = 1.0./N_values
# Save errors to a JDL2 file
using JLD2
save("results/convergence_data_from_fem_k=$(Integer(k))_orderFE=$(orderFE)_orderFEpostprocessing=$(orderFE_postprocessing).jld2",
"h_values", h_values, "k", k, "orderFE", orderFE, "orderFE_postprocessing", orderFE_postprocessing,
"error_L2_u", error_L2_u, "error_H1_u", error_H1_u, "error_L2_p", error_L2_p, "error_H1_p", error_H1_p,
"error_L2_proj_DG", error_L2_proj_DG, "error_H1_proj_DG", error_H1_proj_DG, "error_L2_proj_CG", error_L2_proj_CG, "error_H1_proj_CG", error_H1_proj_CG,
"error_L2_proj_H1", error_L2_proj_H1, "error_H1_proj_H1", error_H1_proj_H1, "error_L2_proj_H2_DG", error_L2_proj_H2_DG, "error_H1_proj_H2_DG", error_H1_proj_H2_DG)
# Plot the relative errors L2 and H1 with latex labels and log scale in Makie
# include("plot_fem_order1.jl")