%config InlineBackend.figure_format='retina'
import numpy as np
import matplotlib.pyplot as plt
from IPython import display
2D MHD simulation¶
I show several problems used to test a MHD code in this post. It solves the equations for 2 dimension below.
$$\frac{{\partial {\bf{Q}}}}{{\partial t}} + \frac{{\partial {\bf{E}}}}{{\partial x}} + \frac{{\partial {\bf{F}}}}{{\partial y}} = \frac{{\partial {{\bf{E}}_v}}}{{\partial x}} + \frac{{\partial {{\bf{F}}_v}}}{{\partial y}},$$
$${\bf{U}} = \left( {\begin{array}{*{20}{c}} \rho \\ {\rho u}\\ {\rho v}\\ {\rho w}\\ e\\ {{B_x}}\\ {{B_y}}\\ {{B_z}} \end{array}} \right),\quad {\bf{E}} = \left( {\begin{array}{*{20}{c}} {\rho u}\\ {\rho {u^2} + P + {B^2}/2 - B_x^2}\\ {\rho uv - {B_x}{B_y}}\\ {\rho uw - {B_x}{B_z}}\\ {\left( {e + P + {B^2}/2} \right)u - {B_x}\left( {{\bf{u}} \cdot {\bf{B}}} \right)}\\ 0\\ {{B_y}u - {B_x}v}\\ {{B_z}u - {B_x}w} \end{array}} \right),\quad {\bf{F}} = \left( {\begin{array}{*{20}{c}} {\rho v}\\ {\rho uv - {B_x}{B_y}}\\ {\rho {v^2} + P + {B^2}/2 - B_y^2}\\ {\rho vw - {B_y}{B_z}}\\ {\left( {e + P + {B^2}/2} \right)v - {B_y}\left( {{\bf{u}} \cdot {\bf{B}}} \right)}\\ {{B_x}v - {B_y}u}\\ 0\\ {{B_z}v - {B_y}w} \end{array}} \right),\quad {{\bf{E}}_v} = \frac{{{M_\infty }}}{{{\mathop{\rm Re}\nolimits} }}\left( {\begin{array}{*{20}{c}} 0\\ {{\tau _{xx}}}\\ {{\tau _{xy}}}\\ {{\tau _{xz}}}\\ {{\beta _x}}\\ 0\\ 0\\ 0 \end{array}} \right) + \frac{{{M_\infty }}}{S}\left( {\begin{array}{*{20}{c}} 0\\ 0\\ 0\\ 0\\ {{\alpha _x}}\\ 0\\ {\eta {J_z}}\\ { - \eta {J_y}} \end{array}} \right),\quad {{\bf{F}}_v} = \frac{{{M_\infty }}}{{{\mathop{\rm Re}\nolimits} }}\left( {\begin{array}{*{20}{c}} 0\\ {{\tau _{yx}}}\\ {{\tau _{yy}}}\\ {{\tau _{yz}}}\\ {{\beta _y}}\\ 0\\ 0\\ 0 \end{array}} \right) + \frac{{{M_\infty }}}{S}\left( {\begin{array}{*{20}{c}} 0\\ 0\\ 0\\ 0\\ {{\alpha _y}}\\ { - \eta {J_z}}\\ 0\\ {\eta {J_x}} \end{array}} \right).$$
Here, $$\begin{array}{l} e = \frac{P}{{\gamma - 1}} + \rho \frac{{{v^2}}}{2} + \frac{{{B^2}}}{{2{\mu _0}}},\\ {\beta _x} = u{\tau _{xx}} + v{\tau _{xy}} + w{\tau _{xz}} - {q_x},\\ {\beta _y} = u{\tau _{yx}} + v{\tau _{yy}} + w{\tau _{yz}} - {q_y},\\ {\alpha _x} = \eta \left( {{B_y}{J_z} - {B_z}{J_y}} \right),\\ {\alpha _y} = \eta \left( {{B_z}{J_x} - {B_x}{J_z}} \right). \end{array}$$
The viscous stress tensor $\tau_{ij}$ (for a Newtonian fluid) is given by $${\tau _{ij}} = \mu \left( {\frac{{\partial {u_j}}}{{\partial {x_i}}} + \frac{{\partial {u_i}}}{{\partial {x_j}}}} \right) + \left( {\beta - \frac{2}{3}\mu } \right)\frac{{\partial {u_k}}}{{\partial {x_k}}}{\delta _{ij}},$$ where $\mu$ is the dynamic (shear) viscosity, $\beta$ is the bulk viscosity.
The heat flux vector $q_i$ is given by $${q_i} = - \kappa \frac{{\partial T}}{{\partial {x_i}}} = - \frac{1}{{\gamma - 1}}\frac{\mu }{{\Pr }}\frac{{\partial {a^2}}}{{\partial {x_i}}},$$ where $\kappa$ is the thermal conductivity, and $a$ is the speed of sound.
The current density vector $\bf{J}$ is defined as $${\bf{J}} = \nabla \times {\bf{B}} = \left[ {\begin{array}{*{20}{c}} {\frac{{\partial {B_z}}}{{\partial y}} - \frac{{\partial {B_y}}}{{\partial z}}}\\ {\frac{{\partial {B_x}}}{{\partial z}} - \frac{{\partial {B_z}}}{{\partial x}}}\\ {\frac{{\partial {B_y}}}{{\partial x}} - \frac{{\partial {B_x}}}{{\partial y}}} \end{array}} \right].$$
In this whole post, $\mu$ and $\kappa$ is non-active and $\eta$ and $\beta$ is active as a diffusion term to capture discontinuities.
In addition to these terms, I put the artificial mass diffusivity term to capture contact discontinuities.
These equations are solved by discritizing the space with the sixth-order compact scheme and the time with the third-order TVD RK scheme.
You can see this paper (Kawai2013JCP) as a reference of this post.
In the following part, 5 problems (the Orszag_Tang problem, The rotor problem, the blast wave problem, the KH instability and the two current sheet) are shown.
The first 3 problems can be compared with the paper above.
#This function is to plot all values.
def overview(it,nof,nx):
f1 =open('%s/data/rho%05d.d' %(nof,it),'rb')
f2 =open('%s/data/vx%05d.d' %(nof,it),'rb')
f3 =open('%s/data/vy%05d.d' %(nof,it),'rb')
f4 =open('%s/data/vz%05d.d' %(nof,it),'rb')
f5 =open('%s/data/pr%05d.d' %(nof,it),'rb')
f6 =open('%s/data/bx%05d.d' %(nof,it),'rb')
f7 =open('%s/data/by%05d.d' %(nof,it),'rb')
f8 =open('%s/data/bz%05d.d' %(nof,it),'rb')
f9 =open('%s/data/beta%05d.d' %(nof,it),'rb')
f10=open('%s/data/eta%05d.d' %(nof,it),'rb')
f11=open('%s/data/chi%05d.d' %(nof,it),'rb')
f12=open('%s/data/jx%05d.d' %(nof,it),'rb')
f13=open('%s/data/jy%05d.d' %(nof,it),'rb')
f14=open('%s/data/jz%05d.d' %(nof,it),'rb')
f15=open('%s/data/divu%05d.d' %(nof,it),'rb')
f16=open('%s/data/divb%05d.d' %(nof,it),'rb')
rho =np.fromfile(f1, dtype='float64').reshape(-1,nx)
vx =np.fromfile(f2, dtype='float64').reshape(-1,nx)
vy =np.fromfile(f3, dtype='float64').reshape(-1,nx)
vz =np.fromfile(f4, dtype='float64').reshape(-1,nx)
pr =np.fromfile(f5, dtype='float64').reshape(-1,nx)
bx =np.fromfile(f6, dtype='float64').reshape(-1,nx)
by =np.fromfile(f7, dtype='float64').reshape(-1,nx)
bz =np.fromfile(f8, dtype='float64').reshape(-1,nx)
beta=np.fromfile(f9, dtype='float64').reshape(-1,nx)
eta =np.fromfile(f10,dtype='float64').reshape(-1,nx)
chi =np.fromfile(f11,dtype='float64').reshape(-1,nx)
jx =np.fromfile(f12,dtype='float64').reshape(-1,nx)
jy =np.fromfile(f13,dtype='float64').reshape(-1,nx)
jz =np.fromfile(f14,dtype='float64').reshape(-1,nx)
divu=np.fromfile(f15,dtype='float64').reshape(-1,nx)
divb=np.fromfile(f16,dtype='float64').reshape(-1,nx)
display.clear_output(True)
plt.figure(figsize=(20,12))
plt.subplots_adjust(wspace=0.2, hspace=0.2)
plt.subplot(3,5,1) ;plt.imshow(rho, aspect='equal',origin='lower',cmap='terrain_r');plt.colorbar();plt.title(r'$\rho$'); plt.tick_params(labelbottom='off',labelleft='off')
plt.subplot(3,5,2) ;plt.imshow(vx, aspect='equal',origin='lower',cmap='terrain_r');plt.colorbar();plt.title(r'$v_x$'); plt.tick_params(labelbottom='off',labelleft='off')
plt.subplot(3,5,3) ;plt.imshow(vy, aspect='equal',origin='lower',cmap='terrain_r');plt.colorbar();plt.title(r'$v_y$'); plt.tick_params(labelbottom='off',labelleft='off')
plt.subplot(3,5,4) ;plt.imshow(vz, aspect='equal',origin='lower',cmap='terrain_r');plt.colorbar();plt.title(r'$v_z$'); plt.tick_params(labelbottom='off',labelleft='off')
plt.subplot(3,5,5) ;plt.imshow(pr, aspect='equal',origin='lower',cmap='terrain_r');plt.colorbar();plt.title(r'$Pr$'); plt.tick_params(labelbottom='off',labelleft='off')
plt.subplot(3,5,6) ;plt.imshow(bx, aspect='equal',origin='lower',cmap='terrain_r');plt.colorbar();plt.title(r'$B_x$'); plt.tick_params(labelbottom='off',labelleft='off')
plt.subplot(3,5,7) ;plt.imshow(by, aspect='equal',origin='lower',cmap='terrain_r');plt.colorbar();plt.title(r'$B_y$'); plt.tick_params(labelbottom='off',labelleft='off')
plt.subplot(3,5,8) ;plt.imshow(bz, aspect='equal',origin='lower',cmap='terrain_r');plt.colorbar();plt.title(r'$B_z$'); plt.tick_params(labelbottom='off',labelleft='off')
plt.subplot(3,5,9) ;plt.imshow(beta,aspect='equal',origin='lower',cmap='terrain_r');plt.colorbar();plt.title(r'$\beta$'); plt.tick_params(labelbottom='off',labelleft='off')
plt.subplot(3,5,10);plt.imshow(eta, aspect='equal',origin='lower',cmap='terrain_r');plt.colorbar();plt.title(r'$\eta$'); plt.tick_params(labelbottom='off',labelleft='off')
plt.subplot(3,5,11);plt.imshow(chi, aspect='equal',origin='lower',cmap='terrain_r');plt.colorbar();plt.title(r'$\chi$'); plt.tick_params(labelbottom='off',labelleft='off')
plt.subplot(3,5,12);plt.imshow(divb,aspect='equal',origin='lower',cmap='terrain_r');plt.colorbar();plt.title(r'$\nabla\cdot\bf{B}$'); plt.tick_params(labelbottom='off',labelleft='off')
plt.subplot(3,5,13);plt.imshow(jx, aspect='equal',origin='lower',cmap='terrain_r');plt.colorbar();plt.title(r'$j_x$'); plt.tick_params(labelbottom='off',labelleft='off')
plt.subplot(3,5,14);plt.imshow(jy, aspect='equal',origin='lower',cmap='terrain_r');plt.colorbar();plt.title(r'$j_y$'); plt.tick_params(labelbottom='off',labelleft='off')
plt.subplot(3,5,15);plt.imshow(jz, aspect='equal',origin='lower',cmap='terrain_r');plt.colorbar();plt.title(r'$j_z$'); plt.tick_params(labelbottom='off',labelleft='off')
#plt.suptitle('it=%06d' %it)
plt.show()
The Orszag_Tang Vortex Problem¶
The initial conditions have X-points in both the velocity and magnetic field: $\rho=2.778, u_x=-\sin y, u_y=\sin x, u_z=0, p=5/3, B_x=-\sin y, B_y=\sin 2x,$ and $B_z=0$ with $\gamma=5/3$.
The number of grid points is $512$ in the x- and y-direnctions.
The CFL number is fixed with $0.4$.
The computational domain extends $0\le x,y\le 2\pi$.
%%time
nof='Orszag_Tang'; nx=512
overview(0,nof,nx)
Here is plots at $t=3.14$.
%%time
nof='Orszag_Tang'; nx=512
overview(31,nof,nx)
The Rotor Problem¶
There is a large pressure jumpy with a very low-pressure region around the center of the rotor. Because of this large pressure jump, many conventional Godunov-type algorithms fail owing to the negative pressure.
The initial dense rotating disk with $\rho=10, u_x=-u_0(y-0.5)/r_0, u_y=-u_0(x-0.5)/r_0,$ and $u_z=0$ is placed at $r<r_0$.
Here, $r_0=1$ and $u_0=2$.
The ambient fluid is at rest with $\rho=1, u_x=u_y=0$ for $r>r_1=0.115$.
The linear taper is applied to the density and the toroidal velocity between the disk and the ambient fluid in the regions $r_0\le r\le r_1$, thus $\rho=1+9f, u_x=-fu_0(y-0.5)/r_0$, and $u_y=-fu_0(x-0.5)/r_0$ where $f=(r_1-r)/(r_1-r_0)$.
The uniform pressure, magnetic field, and specific heat ration are imposed, $p=1, B_x=5/\sqrt{4\pi}, B_y=B_z=0$, and $\gamma=1.4$.
The number of grid points is $512$ in the x- and y-direnctions.
The CFL number is fixed with $0.4$.
The computational domain extends $0\le x,y\le 1$.
%%time
nof='Rotor/'; nx=512
overview(0,nof,nx)
Here is the plots at $t=0.15$
%%time
nof='Rotor/'; nx=512
overview(15,nof,nx)
it=15
nx=512
gam=1.4
f0 =open('Rotor/data/rho%05d.d' %(it),'rb')
f1 =open('Rotor/data/vx%05d.d' %(it),'rb')
f2 =open('Rotor/data/vy%05d.d' %(it),'rb')
f3 =open('Rotor/data/vz%05d.d' %(it),'rb')
f4 =open('Rotor/data/pr%05d.d' %(it),'rb')
f5 =open('Rotor/data/bx%05d.d' %(it),'rb')
f6 =open('Rotor/data/by%05d.d' %(it),'rb')
f7 =open('Rotor/data/bz%05d.d' %(it),'rb')
f8 =open('Rotor/data/eta%05d.d' %(it),'rb')
display.clear_output(True)
rho =np.fromfile(f0, dtype='float64').reshape(-1,nx)
vx =np.fromfile(f1, dtype='float64').reshape(-1,nx)
vy =np.fromfile(f2, dtype='float64').reshape(-1,nx)
vz =np.fromfile(f3, dtype='float64').reshape(-1,nx)
pr =np.fromfile(f4, dtype='float64').reshape(-1,nx)
bx =np.fromfile(f5, dtype='float64').reshape(-1,nx)
by =np.fromfile(f6, dtype='float64').reshape(-1,nx)
bz =np.fromfile(f7, dtype='float64').reshape(-1,nx)
eta =np.fromfile(f8, dtype='float64').reshape(-1,nx)
mach=np.sqrt(vx**2+vy**2+vz**2)/np.sqrt(gam*pr/rho)
pm =0.5*(bx**2+by**2+bz**2)
plt.figure(figsize=(15,5))
plt.subplot(131);plt.contour(mach,20);plt.title('Mach number')
plt.subplot(132);plt.contour(pm,20);plt.title('magnetic pressure')
plt.subplot(133);plt.imshow(eta,aspect='auto',origin='lower');plt.title('Artificial magnetic resistivity')
plt.show()
The Blast Wave Problem¶
The uniform density, velocity, magnetic field, and specific heat ration are imposed, $\rho=1, u_x=u_y=u_z=0, B_x=B_y=10/\sqrt{2}, B_z=0$, and $\gamma=5/3$.
The initial centrally over-pressurized region with $p=100$ and $\beta=2$ is placed at $r<r_0$ and $r_0=0.110$.
The linear taper is applied to the pressure between the over-pressurized region and the outside low-pressure region $p=1$ at $r_0\le r \le r_1$ where $r_1=0.125$.
The number of grid points is $512$ in the x- and y-direnctions.
The CFL number is fixed with $0.4$.
The computational domain extends $-0.5\le x,y\le 0.5$.
%%time
nof='Blast_Wave'; nx=512
overview(0,nof,nx)
%%time
nof='Blast_Wave'; nx=512
overview(20,nof,nx)
it=20
nx=512
f0 =open('Blast_Wave/data/rho%05d.d' %(it),'rb')
f5 =open('Blast_Wave/data/bx%05d.d' %(it),'rb')
f6 =open('Blast_Wave/data/by%05d.d' %(it),'rb')
f7 =open('Blast_Wave/data/bz%05d.d' %(it),'rb')
f8 =open('Blast_Wave/data/eta%05d.d' %(it),'rb')
display.clear_output(True)
rho =np.fromfile(f0, dtype='float64').reshape(-1,nx)
bx =np.fromfile(f5, dtype='float64').reshape(-1,nx)
by =np.fromfile(f6, dtype='float64').reshape(-1,nx)
bz =np.fromfile(f7, dtype='float64').reshape(-1,nx)
eta =np.fromfile(f8, dtype='float64').reshape(-1,nx)
pm =0.5*(bx**2+by**2+bz**2)
intervalrho = np.linspace(0.17,3.34,20)
intervalpm = np.linspace(25,74,20)
plt.figure(figsize=(15,5))
plt.subplot(131);plt.contour(rho,intervalrho);plt.title('Density')
plt.subplot(132);plt.contour(pm, intervalpm) ;plt.title('Magnetic Pressure')
plt.subplot(133);plt.imshow(eta,aspect='auto',origin='lower');plt.title('Artificial magnetic resistivity')
plt.show()
Kelvin-Helmholtz¶
This is not a test problem. I put this case to show a vortex due to this instability.
The parameters can be referred in this paper (Matsumoto&Hoshino2004JGRL).
The grid number is $160\times400$.
I put pertubations for $v_y$ to facilitate the instability.
%%time
nof='KH_Instability'; nx=160
overview(0,nof,nx)
%%time
nof='KH_Instability'; nx=160
overview(40,nof,nx)
Two Current Sheet (Magnetic Reconnection)¶
This problem is shown to prove that this code can capture magnetic reconnections due to artificial magnetic resistivity automaticaly ditected.
The parameters and settings are written in this paper (Keppens2013PoP).
The grid size is $512\times512$.
%%time
nof='Current_Sheet'; nx=512
overview(0,nof,nx)
Here is the plots at $t=100$.
%%time
nof='Current_Sheet'; nx=512
overview(10,nof,nx)
