Linear modes



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