Page MenuHomec4science
Files F93999428

convtest.py
No OneTemporary
Actions

File Metadata

Created
Tue, Dec 3, 03:36

convtest.py
View Options

#!/usr/bin/python
import matplotlib
import numpy as np
import matplotlib.pyplot as plt
#Full explicit values
dt_exp = [ 5e-2,
1e-2,
5e-3,
1e-3,
5e-4 ]
values_exp = [ -5.126686497443573E-003,
-5.126672399036673E-003,
-5.126672374797059E-003,
-5.126672373157391E-003,
-5.126672373154985E-003 ]
# IMEX values
dt_imex = [ 1e-2,
8e-3,
5e-3,
2e-3,
1e-3 ]
values_imex = [ -5.126672537958238E-003,
-5.126672440535181E-003,
-5.126672383407881E-003,
-5.126672373416452E-003,
-5.126672373171177E-003 ]
fig, axis = plt.subplots()
axis.set_title("Convergence test")
axis.set_xlabel("dt")
axis.set_ylabel("|value - ref|")
dt_exp = np.array(dt_exp)
values_exp = np.array(values_exp)
dt_imex = np.array(dt_imex)
values_imex = np.array(values_imex)
# full explicit precise reference
ref = -5.126672373154763E-003
# imex reference
#ref_imex = -5.126672373154302E-003
err_exp = np.abs(values_exp - ref)
err_imex = np.abs(values_imex - ref)
# curve fits
coeff_exp = np.polyfit(np.log(dt_exp), np.log(err_exp), 1)
coeff_imex = np.polyfit(np.log(dt_imex), np.log(err_imex), 1)
#print("Errors full explicit: ", list(err_exp))
#print("Errors IMEX: ", list(err_imex))
print("Convergence order (full explicit):", coeff_exp[0])
print("Stability coefficient (full explicit):", np.exp(coeff_exp[1]))
print("Convergence order (IMEX):", coeff_imex[0])
print("Stability coefficient (IMEX):", np.exp(coeff_imex[1]))
def fit_exp(dt):
return np.exp(coeff_exp[0] * np.log(dt) + coeff_exp[1])
def fit_imex(dt):
return np.exp(coeff_imex[0] * np.log(dt) + coeff_imex[1])
# plots
axis.loglog(dt_exp, err_exp, 'b*', label='Full explicit')
axis.loglog(dt_exp, fit_exp(dt_exp), 'k-')
axis.loglog(dt_imex, err_imex, 'r*', label='IMEX')
axis.loglog(dt_imex, fit_imex(dt_imex), 'k-')
axis.legend()
# Easy test
dt_easy = [1e-2,
8e-3,
6.25e-3,
5e-3,
2.5e-3,
2e-3,
1e-3]
values_easy = [ -0.623588259377137,
-0.623588255920246,
-0.623588254409333,
-0.623588253877998,
-0.623588253531279,
-0.623588253517578,
-0.623588253508655]
fig_, axis_ = plt.subplots()
axis_.set_title("Easy convergence test")
axis_.set_xlabel("dt")
axis_.set_ylabel("|value - ref|")
dt_easy = np.array(dt_easy)
values_easy = np.array(values_easy)
# easy reference, -0.6235882520493409
#ref_easy = -0.623588253508096
ref_easy = -0.623588253508061
u = 5.0
#ref_easy = - np.sin(2 * np.pi * np.sqrt(2) * u / 18.0)
print(ref_easy)
err_easy = np.abs(values_easy - ref_easy)
for err in err_easy:
print(err)
axis_.loglog(dt_easy, err_easy, 'b*')
coeff_easy = np.polyfit(np.log(dt_easy), np.log(err_easy), 1)
def fit_easy(dt):
return np.exp(coeff_easy[0] * np.log(dt) + coeff_easy[1])
print("Convergence order (easy):", coeff_easy[0])
print("Stability coefficient (easy):", np.exp(coeff_easy[1]))
axis_.loglog(dt_easy, fit_easy(dt_easy), 'k-')
# show all plots
plt.show()

Event Timeline

c4science · Help