In [4]:
import h5py
import glob
import numpy as np
import pylab as pl
%matplotlib inline
# Set plot parameters to make beautiful plots
pl.rcParams['figure.figsize'] = 12, 7.5
pl.rcParams['lines.linewidth'] = 1.5
pl.rcParams['font.family'] = 'serif'
pl.rcParams['font.weight'] = 'bold'
pl.rcParams['font.size'] = 30
pl.rcParams['font.sans-serif'] = 'serif'
pl.rcParams['text.usetex'] = True
pl.rcParams['axes.linewidth'] = 1.5
pl.rcParams['axes.titlesize'] = 'medium'
pl.rcParams['axes.labelsize'] = 'medium'
pl.rcParams['xtick.major.size'] = 8
pl.rcParams['xtick.minor.size'] = 4
pl.rcParams['xtick.major.pad'] = 8
pl.rcParams['xtick.minor.pad'] = 8
pl.rcParams['xtick.color'] = 'k'
pl.rcParams['xtick.labelsize'] = 'medium'
pl.rcParams['xtick.direction'] = 'in'
pl.rcParams['ytick.major.size'] = 8
pl.rcParams['ytick.minor.size'] = 4
pl.rcParams['ytick.major.pad'] = 8
pl.rcParams['ytick.minor.pad'] = 8
pl.rcParams['ytick.color'] = 'k'
pl.rcParams['ytick.labelsize'] = 'medium'
pl.rcParams['ytick.direction'] = 'in'
In [5]:
N = np.array([32, 64, 128, 256, 512])
errors_rho = [3.988654e-12, 4.969112e-13, 5.100655e-14, 9.139936e-15, 1.210917e-15]
errors_u = [1.169794e-11, 1.228705e-12, 1.224386e-13, 1.150570e-14, 1.428690e-15]
errors_u1 = [2.550613e-13, 4.841199e-14, 6.637352e-15, 5.871516e-16, 6.168086e-17]
errors_u2 = [1.426849e-12, 1.326205e-13, 9.958794e-15, 1.233547e-15, 2.263401e-16]
errors_q = [6.370257e-12, 5.447082e-13, 4.694195e-14, 4.554252e-15, 7.538579e-16]
errors_dP = [4.509557e-12, 4.659176e-13, 4.400123e-14, 4.295960e-15, 5.914910e-16]
errors_B1 = [4.857199e-13, 7.734114e-14, 8.072464e-15, 5.271674e-16, 9.312974e-17]
errors_B2 = [3.446792e-13, 4.281455e-14, 4.445065e-15, 3.543071e-16, 5.441856e-17]
In [6]:
pl.figure(figsize=(10, 24))
pl.subplots_adjust(hspace=0.1)
axes1 = pl.subplot(8, 1, 1)
pl.title('$L_1\;\mathrm{norm}\;\mathrm{of}\;\mathrm{error}$')
pl.loglog(N, errors_rho, 'o-', label='$\mathtt{grim}$', color='blue', linewidth=1.5, basex=2)
pl.loglog(N, 1e-8*np.array(N)**(-2.), color='black', linewidth=3.5, linestyle='--', label='$O(\\Delta X_1^2, \\Delta X_2^2)$', basex=2)
pl.legend(fontsize='small', loc='lower left', frameon=False)
pl.ylabel('$\\rho$')
pl.yticks([1e-15, 1e-13, 1e-11])
pl.xticks([])
pl.xlim(xmin=24)
pl.xlim(xmax=600)
pl.ylim(ymax=2e-10)
axes1 = pl.subplot(8, 1, 2)
pl.loglog(N, errors_u, 'o-', label='$\mathtt{grim}$', color='blue', linewidth=1.5, basex=2)
pl.loglog(N, 3.5e-8*np.array(N)**(-2.), color='black', linewidth=3.5, linestyle='--', label='$O(\\Delta X_1^2, \\Delta X_2^2)$', basex=2)
pl.ylabel('$u$')
#pl.yticks([])
pl.yticks([1e-15, 1e-13, 1e-11])
pl.xticks([])
pl.xlim(xmin=24)
pl.xlim(xmax=600)
pl.ylim(ymax=4e-10)
axes1 = pl.subplot(8, 1, 3)
pl.loglog(N, errors_u1, 'o-', label='$\mathtt{grim}$', color='blue', linewidth=1.5, basex=2)
pl.loglog(N, 1.5e-9*np.array(N)**(-2.), color='black', linewidth=3.5, linestyle='--', label='$O(\\Delta X_1^2, \\Delta X_2^2)$', basex=2)
pl.ylabel('$u^1$')
#pl.yticks([])
pl.yticks([1e-15, 1e-13, 1e-11])
pl.xticks([])
pl.xlim(xmin=24)
pl.xlim(xmax=600)
pl.ylim(ymax=2.5e-11)
axes1 = pl.subplot(8, 1, 4)
pl.loglog(N, errors_u2, 'o-', label='$\mathtt{grim}$', color='blue', linewidth=1.5, basex=2)
pl.loglog(N, 3.5e-9*np.array(N)**(-2.), color='black', linewidth=3.5, linestyle='--', label='$O(\\Delta X_1^2, \\Delta X_2^2)$', basex=2)
pl.ylabel('$u^2$')
#pl.yticks([])
pl.yticks([1e-15, 1e-13, 1e-11])
pl.xticks([])
pl.xlim(xmin=24)
pl.xlim(xmax=600)
pl.ylim(ymax=1e-10)
axes1 = pl.subplot(8, 1, 5)
pl.loglog(N, errors_B1, 'o-', label='$\mathtt{grim}$', color='blue', linewidth=1.5, basex=2)
pl.loglog(N, 2e-9*np.array(N)**(-2.), color='black', linewidth=3.5, linestyle='--', label='$O(\\Delta X_1^2, \\Delta X_2^2)$', basex=2)
pl.ylabel('$B^1$')
#pl.yticks([])
pl.yticks([1e-15, 1e-13, 1e-11])
pl.xticks([])
pl.xlim(xmin=24)
pl.xlim(xmax=600)
pl.ylim(ymax=1e-10)
axes1 = pl.subplot(8, 1, 6)
pl.loglog(N, errors_B2, 'o-', label='$\mathtt{grim}$', color='blue', linewidth=1.5, basex=2)
pl.loglog(N, 1e-9*np.array(N)**(-2.), color='black', linewidth=3.5, linestyle='--', label='$O(\\Delta X_1^2, \\Delta X_2^2)$', basex=2)
pl.ylabel('$B^2$')
#pl.yticks([])
pl.yticks([1e-15, 1e-13, 1e-11])
pl.xticks([])
pl.xlim(xmin=24)
pl.xlim(xmax=600)
pl.ylim(ymax=1e-10)
axes1 = pl.subplot(8, 1, 7)
pl.loglog(N, errors_q, 'o-', label='$\mathtt{grim}$', color='blue', linewidth=1.5, basex=2)
pl.loglog(N, 7e-9*np.array(N)**(-2.), color='black', linewidth=3.5, linestyle='--', label='$O(\\Delta X_1^2, \\Delta X_2^2)$', basex=2)
pl.ylabel('$q$')
#pl.yticks([])
pl.yticks([1e-15, 1e-13, 1e-11])
pl.xticks([32, 64, 128, 256, 512])
pl.xticks([])
pl.xlim(xmin=24)
pl.xlim(xmax=600)
pl.ylim(ymax=3e-10)
axes1 = pl.subplot(8, 1, 8)
pl.loglog(N, errors_dP, 'o-', label='$\mathtt{grim}$', color='blue', linewidth=1.5, basex=2)
pl.loglog(N, 7e-9*np.array(N)**(-2.), color='black', linewidth=3.5, linestyle='--', label='$O(\\Delta X_1^2, \\Delta X_2^2)$', basex=2)
pl.ylabel('$\\Delta P$')
#pl.yticks([])
pl.yticks([1e-15, 1e-13, 1e-11])
#pl.xticks([32, 64, 128, 256, 512])
#pl.xticks([1, 10, 100])
pl.xlim(xmin=24)
pl.xlim(xmax=600)
pl.ylim(ymax=3e-10)
pl.xlabel('$\mathrm{Resolution}$')
pl.savefig('linear_modes_convergence_full_emhd_new.png')
In [16]:
amplitude = 1e-8
k1 = 2.*np.pi
k2 = 4.*np.pi
rho0 = 1.
u0 = 2.
u10 = 0.
u20 = 0.
u30 = 0.
B10 = 0.1
B20 = 0.3
B30 = 0.
q0 = 0.
dp0 = 0.
omega = -0.5533585207638141 - 3.6262571286888425*1j
delta_rho = -0.518522524082246 - 0.1792647678001878*1j
delta_u = 0.5516170736393801
delta_u1 = 0.008463122479547733 + 0.011862022608466213*1j
delta_u2 = -0.16175466371870767 - 0.03482808082360316*1j
delta_u3 = 0.
delta_B1 = -0.05973794979640759 - 0.03351707506150907*1j
delta_B2 = 0.029868974898203816 + 0.01675853753075452*1j
delta_B3 = 0.
delta_q = 0.5233486841539436 + 0.04767672501939603*1j
delta_dp = 0.2909106062057661 + 0.02159452055336522*1j
"""
omega = -0.5236801536105606 - 2.2074259426146234*1j
delta_rho = 0.8245277155993238
delta_u = -0.2655263062409864 + 0.11295573948161661*1j
delta_u1 = 0.01870577997997143 - 0.07463791562430622*1j
delta_u2 = 0.13548478237175182 + 0.0029582968057185223*1j
delta_u3 = 0.
delta_B1 = 0.013767510382234251 + 0.13241920582661107*1j
delta_B2 = -0.006883755191117101 - 0.06620960291330555*1j
delta_B3 = 0.
delta_q = -0.38685433587246587 + 0.11071197829033866*1j
delta_dp = -0.1599729569151908 + 0.054267588035646866*1j
"""
mode = np.exp(1j*k1*x + 1j*k2*y + omega*final_time)
theory_rho = amplitude * (delta_rho * mode).real
theory_u = amplitude * (delta_u * mode).real
theory_u1 = amplitude * (delta_u1 * mode).real
theory_u2 = amplitude * (delta_u2 * mode).real
theory_u3 = amplitude * (delta_u3 * mode).real
theory_B1 = amplitude * (delta_B1 * mode).real
theory_B2 = amplitude * (delta_B2 * mode).real
theory_B3 = amplitude * (delta_B3 * mode).real
theory_q = amplitude * (delta_q * mode).real
theory_dp = amplitude * (delta_dp * mode).real
numerical_rho = (final_data[:, :, RHO] - rho0)
numerical_u = (final_data[:, :, UU] - u0)
numerical_u1 = (final_data[:, :, U1] - u10)
numerical_u2 = (final_data[:, :, U2] - u20)
numerical_u3 = (final_data[:, :, U3] - u30)
numerical_B1 = (final_data[:, :, B1] - B10)
numerical_B2 = (final_data[:, :, B2] - B20)
numerical_B3 = (final_data[:, :, B3] - B30)
numerical_q = (final_data[:, :, PHI] - q0)
numerical_dp = (final_data[:, :, PSI] - dp0)
res = N1 * N2
err_rho = simps(simps(abs(theory_rho - numerical_rho), axis=0)) / res
err_u = simps(simps(abs(theory_u - numerical_u), axis=0)) / res
err_u1 = simps(simps(abs(theory_u1 - numerical_u1), axis=0)) / res
err_u2 = simps(simps(abs(theory_u2 - numerical_u2), axis=0)) / res
err_u3 = simps(simps(abs(theory_u3 - numerical_u3), axis=0)) / res
err_B1 = simps(simps(abs(theory_B1 - numerical_B1), axis=0)) / res
err_B2 = simps(simps(abs(theory_B2 - numerical_B2), axis=0)) / res
err_B3 = simps(simps(abs(theory_B3 - numerical_B3), axis=0)) / res
err_q = simps(simps(abs(theory_q - numerical_q), axis=0)) / res
err_dp = simps(simps(abs(theory_dp - numerical_dp), axis=0)) / res