<jupyter_start><jupyter_text>MHD Equation with CentPy in 1D  Import packages<jupyter_code># Install the centpy package
!pip install centpy
# Import numpy and centpy for the solution 
import numpy as np
import centpy
# Imports functions from matplotlib and setup for the animation
import matplotlib.pyplot as plt
from matplotlib import animation
from IPython.display import HTML<jupyter_output><empty_output><jupyter_text>Equation We solve the equations of ideal magnetohydrodynamics in 1D \begin{equation} \partial_t \begin{bmatrix} \rho \\ \rho v_x \\ \rho v_y \\ \rho v_z \\ B_y \\ B_z \\ E \end{bmatrix} + \partial_x \begin{bmatrix} \rho v_x \\ \rho v_x^2 + p^* - B_x^2 \\ \rho v_x v_y - B_x B_y \\\rho v_x v_z - B_x B_z \\ B_y v_x - B_x v_y \\ B_z v_x - B_x v_z \\(E+p^*) v_x - B_x (B_x v_x + B_y v_y + B_z v_Z) \end{bmatrix} = 0 \end{equation}where the total pressure is given by \begin{equation}p^* = p + \frac{1}{2} (B_x^2 + B_y^2 + B_z^2)\end{equation}with the equation of state\begin{equation}p = (\gamma-1) \left(E-\frac{1}{2} \rho (v_x^2+v_y^2+v_z^2) - \frac{1}{2}(B_x^2 + B_y^2 + B_z^2)\right), \qquad \gamma=2.0\end{equation}The solution is computed on the domain $(x,t)\in([-1,1]\times[0,0.2])$ with initial data for a *Brio-Wu shock tube*:\begin{equation}(\rho, v_x, v_y, v_z, B_y, B_z, p)_{t=0} = \begin{cases}(1,0,0,0,1,0,1) & \text{if} & -1<x\leq 0 \\(0.125, 0, 0, 0, -1, 0, 0.1) & \text{if} & \ \ 0<x<1\end{cases}\end{equation}and Dirichlet boundary data set by initial data on each boundary. The solution is computed using 400 cells and CFL number 0.475.<jupyter_code>pars = centpy.Pars1d(x_init=-1., x_final=1., t_final=0.2, dt_out=0.002, J=400, cfl=0.475, scheme="fd2")
pars.B1 = 0.75
# MHD equation
class MHD1d(centpy.Equation1d):
    def pressure(self, u):
        return u[:, 6] - 0.5*((u[:, 1] ** 2 + u[:, 2] ** 2 + u[:, 3] ** 2)/u[:, 0]) 
        - 0.5 * (self.B1 ** 2 + u[:, 4] ** 2 + u[:, 5] ** 2)

    def initial_data(self):
        u = np.zeros((self.J + 4, 7))
        midpoint = int(self.J / 2) + 2

        # Left side
        u[:midpoint, 0] = 1.0
        u[:midpoint, 1] = 0.0
        u[:midpoint, 2] = 0.0
        u[:midpoint, 3] = 0.0
        u[:midpoint, 4] = 1.0
        u[:midpoint, 5] = 0.0
        u[:midpoint, 6] = 1.0 + 25.0 / 32.0

        # Right side
        u[midpoint:, 0] = 0.125
        u[midpoint:, 1] = 0.0
        u[midpoint:, 2] = 0.0
        u[midpoint:, 3] = 0.0
        u[midpoint:, 4] = -1.0
        u[midpoint:, 5] = 0.0
        u[midpoint:, 6] = 0.1 + 25.0 / 32.0

        return u

    def boundary_conditions(self, u):

        left_v = [1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 1.0 + 25.0 / 32.0]
        right_v = [0.125, 0.0, 0.0, 0.0, -1.0, 0.0, 0.1 + 25.0 / 32]

        if self.odd:
            u[0] = left_v
            u[-1] = right_v
            u[-2] = right_v
        else:
            u[0] = left_v
            u[1] = left_v
            u[-1] = right_v

    def flux_x(self, u):
        f = np.zeros_like(u)

        B1 = self.B1

        p_star = self.pressure(u) + 0.5 * (B1 ** 2 + u[:, 4] ** 2 + u[:, 5] ** 2)

        f[:, 0] = u[:, 1]
        f[:, 1] = u[:, 1] ** 2 / u[:, 0] + p_star
        f[:, 2] = u[:, 1] * u[:, 2] / u[:, 0] - B1 * u[:, 4]
        f[:, 3] = u[:, 1] * u[:, 3] / u[:, 0] - B1 * u[:, 5]
        f[:, 4] = u[:, 1] * u[:, 4] / u[:, 0] - B1 * u[:, 2] / u[:, 0]
        f[:, 5] = u[:, 1] * u[:, 5] / u[:, 0] - B1 * u[:, 3] / u[:, 0]
        f[:, 6] = (u[:, 6] + p_star) * (u[:, 1] / u[:, 0]) - B1 * (
            B1 * u[:, 1] + u[:, 2] * u[:, 4] + u[:, 3] * u[:, 5]
        ) / u[:, 0]

        return f

    def spectral_radius_x(self, u):
        rho = u[:, 0]
        u_x = u[:, 1] / rho
        p = self.pressure(u)
        A = 2.0 * p / rho
        B = (self.B1 ** 2 + u[:, 4] ** 2 + u[:, 5] ** 2) / rho
        cf = np.sqrt(
            0.5 * (A + B + np.sqrt((A + B) ** 2 - 4.0 * A * self.B1 ** 2 / rho))
        )
        return np.abs(u_x) + cf<jupyter_output><empty_output><jupyter_text>Solution<jupyter_code>eqn = MHD1d(pars)
soln = centpy.Solver1d(eqn)
soln.solve()<jupyter_output><empty_output><jupyter_text>Animation<jupyter_code># Animation 
j0 = slice(2,-2)

fig = plt.figure(figsize=(12,6))
ax1=fig.add_subplot(1,3,1)
ax2=fig.add_subplot(2,3,2)
ax3=fig.add_subplot(2,3,3)
ax4=fig.add_subplot(2,3,5)
ax5=fig.add_subplot(2,3,6)

line_u=[]
for ax in [ax1,ax2,ax3,ax4,ax5]:
  ax.set_xlim(-1, 1)
  line_u.append(ax.plot([], [], linewidth=1, marker='o', markersize=2)[0])

ax1.set_xlabel('x')
ax4.set_xlabel('x')
ax5.set_xlabel('x')

ax2.set_xticks([])
ax3.set_xticks([])
ax2.set_yticks([])
ax3.set_yticks([])
ax4.set_yticks([])
ax5.set_yticks([])

ax1.set_ylabel(r'$\rho$', fontsize=12)
ax2.set_ylabel(r'$v_x$', fontsize=12)
ax3.set_ylabel(r'$v_y$', fontsize=12)
ax4.set_ylabel(r'$B_y$', fontsize=12)
ax5.set_ylabel(r'$p$', fontsize=12)

# Primitive variables
rho=soln.u_n[:,j0,0]
v_x = soln.u_n[:,j0,1]/soln.u_n[:,j0,0]
v_y = soln.u_n[:,j0,2]/soln.u_n[:,j0,0]
B_y = soln.u_n[:,j0,4]

ax1.set_ylim(np.min(rho), np.max(rho))
ax2.set_ylim(np.min(v_x), np.max(v_x))
ax3.set_ylim(np.min(v_y), np.max(v_y))
ax4.set_ylim(np.min(B_y), np.max(B_y))
ax5.set_ylim(0.5, 2.)

def animate(i):
    p = eqn.pressure(soln.u_n[i,j0,:])
    line_u[0].set_data(soln.x[j0], rho[i])
    line_u[1].set_data(soln.x[j0], v_x[i])
    line_u[2].set_data(soln.x[j0], v_y[i])
    line_u[3].set_data(soln.x[j0], B_y[i])
    line_u[4].set_data(soln.x[j0], p)
    
plt.close()
anim = animation.FuncAnimation(fig, animate, frames=soln.Nt, interval=100, blit=False);
HTML(anim.to_html5_video())<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.21574239713774598, 'fraction_numerical': 0.06690518783542039, 'mean_word_length': 2.666229508196721, 'pattern_counts': {'":': 0, '<': 20, '<?xml version=': 0, '>': 17, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 5, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '4', 'dir': 'jupyter-structured-clean-dedup', 'id': '454966', 'license': 'MIT', 'n_tokens_mistral': 2782, 'n_tokens_neox': 2635, 'n_words': 558, 'path': 'notebooks/MHD_1d.ipynb', 'repo_name': 'olekravchenko/centpy', 'size': '5590'}