DUBeat 1.0.1
High-order discontinuous Galerkin methods and applications to cardiac electrophysiology
Loading...
Searching...
No Matches
generate_convergence_plots.py
Go to the documentation of this file.
1#!/usr/bin/env python3
2
16
17
20
21# Author: Federica Botta <federica.botta@mail.polimi.it>.
22# Author: Matteo Calafà <matteo.calafa@mail.polimi.it>.
23
24
25#
26# Generate plots for the convergence analysis on DG schemes.
27# N.B. requires pandas and matplotlib.
28#
29
30
31import numpy as np
32from matplotlib import pyplot as plt
33import pandas as pd
34
35
36
37errors = pd.read_csv("../errors_Laplace_3D_u.data", sep='\t')
38errors = errors[["l_inf","l_2","h_1","DG"]]
39errors = np.array(errors)
40errors = errors[errors[:,0]!='x',:].astype(float)
41n = errors.shape[0]
42
43
44if (n > 0):
45
46 x = np.logspace(-2, -n, num=n, base=2.0)
47
48
49 k = errors[0,0]
50 plt.figure()
51 ax = plt.axes()
52 plt.loglog(x, errors[:,0], 'tab:blue', marker = 'o', linestyle = '-', label='$L^\infty$ error')
53 plt.loglog(x, k*np.power(x,1)/np.power(x[0],1), 'k', linestyle = '-.', label='$h^{-1}$')
54 plt.loglog(x, k*np.power(x,2)/np.power(x[0],2), 'k', linestyle = '--', label='$h^{-2}$')
55 plt.loglog(x, k*np.power(x,3)/np.power(x[0],3), 'k', linestyle = 'dotted', label = '$h^{-3}$')
56 plt.xlabel('h')
57 plt.ylabel('$L^\infty$ error')
58 plt.legend()
59 ax.set_title('Convergence test ($L^\infty$ error)')
60 plt.savefig("Linf")
61
62
63
64 k = errors[0,1]
65 plt.figure()
66 ax = plt.axes()
67 plt.loglog(x, errors[:,1], 'tab:blue', marker = 'o', linestyle = '-', label='$L^2$ error')
68 plt.loglog(x, k*np.power(x,1)/np.power(x[0],1), 'k', linestyle = '-.', label='$h^{-1}$')
69 plt.loglog(x, k*np.power(x,2)/np.power(x[0],2), 'k', linestyle = '--', label='$h^{-2}$')
70 plt.loglog(x, k*np.power(x,3)/np.power(x[0],3), 'k', linestyle = 'dotted', label = '$h^{-3}$')
71 plt.xlabel('h')
72 plt.ylabel('$L^2$ error')
73 plt.legend()
74 ax.set_title('Convergence test ($L^2$ error)')
75 plt.savefig("L2")
76
77
78 k = errors[0,2]
79 plt.figure()
80 ax = plt.axes()
81 plt.loglog(x, errors[:,2], 'tab:blue', marker = 'o', linestyle = '-', label='$H^1$ error')
82 plt.loglog(x, k*np.power(x,1)/np.power(x[0],1), 'k', linestyle = '-.', label='$h^{-1}$')
83 plt.loglog(x, k*np.power(x,2)/np.power(x[0],2), 'k', linestyle = '--', label='$h^{-2}$')
84 plt.loglog(x, k*np.power(x,3)/np.power(x[0],3), 'k', linestyle = 'dotted', label = '$h^{-3}$')
85 plt.xlabel('h')
86 plt.ylabel('$H^1$ error')
87 plt.legend()
88 ax.set_title('Convergence test ($H^1$ error)')
89 plt.savefig("H1")
90
91
92 k = errors[0,3]
93 plt.figure()
94 ax = plt.axes()
95 plt.loglog(x, errors[:,3], 'tab:blue', marker = 'o', linestyle = '-', label='$DG$ error')
96 plt.loglog(x, k*np.power(x,1)/np.power(x[0],1), 'k', linestyle = '-.', label='$h^{-1}$')
97 plt.loglog(x, k*np.power(x,2)/np.power(x[0],2), 'k', linestyle = '--', label='$h^{-2}$')
98 plt.loglog(x, k*np.power(x,3)/np.power(x[0],3), 'k', linestyle = 'dotted', label = '$h^{-3}$')
99 plt.xlabel('h')
100 plt.ylabel('$DG$ error')
101 plt.legend()
102 ax.set_title('Convergence test ($DG$ error)')
103 plt.savefig("DG")