34.2 Magnetohydrodynamics Test Problems

The magnetohydrodynamics (MHD) test problems provided in this release can be found in source/Simulation/SimulationMain/magnetoHD/. In order to set up an MHD problem, users need to specify the magnetoHD path in a setup script. For instance, the BrioWu problem can be configured by typing ./setup magnetoHD/BrioWu -auto -1d.

34.2.1 Brio-Wu MHD Shock Tube

The Brio-Wu MHD shock tube problem (Brio and Wu, 1988), magnetoHD/BrioWu, is a coplanar magnetohydrodynamic counterpart of the hydrodynamic Sod problem (Sec:SimulationSod). The initial left and right states are given by $\rho_l=1$, $u_l=v_l=0$, $p_l=1$, $(B_y)_l=1$; and $\rho_r=0.125$, $u_r=v_r=0$, $p_r=0.1$, $(B_y)_r=-1$. In addition, $B_x=0.75$ and $\gamma=2$. This is a good problem to test wave properties of a particular MHD solver, because it involves two fast rarefaction waves, a slow compound wave, a contact discontinuity and a slow shock wave.

The conventional 800 point solution to this problem computed with FLASH is presented in Figure 34.35, Figure 34.36, Figure 34.37, Figure 34.38, and Figure 34.39 . The figures show the distribution of density, normal and tangential velocity components, tangential magnetic field component and pressure at $t=0.1$ (in non-dimensional units). As can be seen, the code accurately and sharply resolves all waves present in the solution. There is a small undershoot in the solution at $x\approx0.44$, which results from a discontinuity-enhancing monotonized centered gradient limiting function (LeVeque 1997). This undershoot can be easily removed if a less aggressive limiter, e.g. a minmod or a van Leer limiter, is used instead. This, however, will degrade the sharp resolution of other discontinuities.

The directionally splitting 8Wave MHD solver with a second-order MUSCL-Hancock scheme (setup with +8wave) was used for the results shown in this simulation. The StaggeredMesh MHD solver (setup with +usm) can also be used for this Brio-Wu problem in one- and two-dimensions. However, in the latter case, the StaggeredMesh solver only supports non-rotated setups for which a shock normal is parallel to the $x$-axis that initially intersects that axis at $x=0.5$ (halfway across a box with unit dimensions). This limitation occurs in the StaggeredMesh scheme because the currently released version of the FLASH code does not truly support physically correct boundary conditions for this rotated shock geometry.

Figure 34.35: Density profile for the Brio-Wu shock tube problem.

Figure 34.36: Pressure profile for the Brio-Wu shock tube problem.

Figure 34.37: Tangential magnetic field profile for the Brio-Wu shock tube problem.

Figure 34.38: Normal velocity profile for the Brio-Wu shock tube problem.

Figure 34.39: Tangential velocity profile for the Brio-Wu shock tube problem.

34.2.1.1 Slowly moving shocks (SMS) issues in the Brio-Wu problem

Figure 34.40 clearly demonstrates that the conventional PPM reconstruction method fails to preserve monotonicities, shedding oscillations especially in the plateau near strong discontinuities such as the contact and right going slow MHD shock. In Fig. 34.41, Mach numbers are plotted with varying strengths of the transverse field $B_y$. The oscillations increase with an increase of $B_y$; the reason being that stronger $B_y$ introduces more transverse effect that resists shock propagation in the $x$ direction causing the shock to move slowly. This effect is clearly seen in the locations of the shock fronts (right going slow MHD shocks in this case), which remain closer to the initial location $x=0.5$ when $B_y$ is stronger.

Figure: Density (thick red) and Mach number (thick blue) at $t=0.1$ of the Brio-Wu test with the conventional PPM's MC limiter. Thin black curves represent reference solutions using the PLM of MUSCL-Hancock scheme. Severe numerical oscillations are evident in the solution using the conventional PPM reconstruction on 400 grid points.

Figure: Mach numbers at $t=0.1$ with varying $B_y$ from 0 to 1 with the conventional PPM's MC limiter. Curves in red, blue, green, purple, black, and cyan represent $B_y=0, 0.2, 0.4, 0.6, 0.8,$ and $1$, respectively. Severe numerical oscillations are evident in the solution using the conventional PPM reconstruction on 400 grid points.

Results from using the upwind biased slope limiter for PPM are illustrated in Figures 34.42 - 34.45 . The oscillation shedding found in the conventional PPM (e.g., Fig. 34.40) are significantly reduced in both density, and Mach number profiles. The overall qualitative solution behavior of using the upwind PPM approach, shown in Fig. 34.42, is very much similar to those from the 5th order WENO method as illustrated in Fig. 34.43. As shown in the closeup views in Fig. 34.44 and 34.45, the solutions with the upwind PPM slope limiter (blue curves) outperforms the conventional PPM method (red curves), and compares well with the WENO scheme (green curves). In fact, the upwind approach shows the most flat density plateau of all methods. Note that the SMS issue can be observed regardless of the choice of Riemann solvers and the dissipation mechanism in PPM (e.g., even with flattening on). 400 grid points were used.

Figure: Density (red) and Mach number (blue) at $t=0.1$ of the Brio-Wu test using the upwind biased PPM limiter.

Figure: Density (red) and Mach number (blue) at $t=0.1$ of the Brio-Wu test using the 5th order WENO scheme.

Figure: Density closeup at $t=0.1$ of the Brio-Wu test. The conventional PPM in red curve, the upwind PPM in blue curve, and the WENO scheme in green curve.

Figure: Mach number closeup at $t=0.1$ of the Brio-Wu test. The conventional PPM in red curve, the upwind PPM in blue curve, and the WENO scheme in green curve.

34.2.2 Orszag-Tang MHD Vortex

The Orszag-Tang MHD vortex problem (Orszag and Tang, 1979), magnetoHD/OrszagTang, is a simple two-dimensional problem that has become a classic test for MHD codes. In this problem a simple, non-random initial condition is imposed at time $t=0$

$${\bf V} = V_0\left(-\sin(2\pi y),\sin(2\pi x), 0\right), {\bf B} = B_0\left(-\sin(2\pi y),\sin(4\pi x), 0\right), (x,y)=\in[0,1]^2,$$ (34.17)

where $B_0$ is chosen so that the ratio of the gas pressure to the RMS magnetic pressure is equal to $2\gamma$. In this setup the initial density, the speed of sound and $V_0$ are set to unity; therefore, the initial pressure $p_0=1/\gamma$ and $B_0=1/\gamma$.

As the evolution time increases, the vortex flow pattern becomes increasingly complicated due to the nonlinear interactions of waves. A highly resolved simulation of this problem should produce two-dimensional MHD turbulence. Figure 34.46 and Figure 34.47 shows density and magnetic field contours at $t=0.5$. As one can observe, the flow pattern at this time is already quite complicated. A number of strong waves have formed and passed through each other, creating turbulent flow features at all spatial scales.

The results were obtained using the directionally splitting 8Wave MHD solver for this Orszag-Tang problem.

Figure: Density contours in the Oszag-Tang MHD vortex problem at $t=0.5$.

Figure: Magnetic field contours in the Oszag-Tang MHD vortex problem at $t=0.5$.

The 3D version are also shown in the below solved using the unsplit staggered mesh MHD solver.

Figure: Density plots of a 3D version of the Orszag-Tang problem on a $128^3$ uniform grid resolution. (a) Density at $t=0.2$ (b) Density at $t=0.5$ (c) Density at $t=0.7$ (d) Density at $t=1.0$.

[a] [b] [c] [d]

34.2.3 Magnetized Accretion Torus

The magnetized accretion torus problem is based on the global magneto-rotational instability (MRI) simulations of Hawley (2000). It can be found under magnetoHD/Torus. We consider a magnetized torus of constant angular momentum ( $\Omega\propto r^{-2}$), inside a normalized Paczynsky-Wiita pseudo-Newtonean gravitational potential (Paczynsky&Wiita, 1980) of the form $\Phi=-1/(R-1)$, where $R=\sqrt{r^2+z^2}$.

The cylindrical computational domain, with lengths normalized at $r_0$, extends from 1.5 $r_0$ to 15.5 $r_0$ in the radial direction and from -7 $r_0$ to 7 $r_0$ in the $z$ direction. For this specific simulation we use seven levels of refinement, linear reconstruction with vanLeer limiter and the HLLC Riemann solver. The boundary conditions are set to outflow, except for the leftmost side where a diode-like condition is applied.

Assuming an adiabatic equation of state, the initial density profile of the torus is given by

$$\displaystyle \frac{\Gamma P}{(\Gamma-1)\rho} = C - \Phi -\frac{l_K^2}{2r^2},$$ (34.18)

where the specific heats ratio is $\Gamma=5/3$, $l_K$ is the Keplerian angular momentum at the pressure maximum ( $r_{Pmax}=4.7 r_0$) and C is an integration constant that specifies the outer surface of the torus, given its inner radius ( $r_{in}=3 r_0$). The initial poloidal magnetic field is computed using the $\phi$ component of the vector potential, $A_{\phi} \propto \textrm{max}(\rho-5, 0)$, and normalized so as the initial minimum value of the plasma $\beta=2P/{\bf B}^2$ is equal to $10^2$. The resulting field follows the contours of density, i.e. torus-like nested loops, and is embedded well within the torus.

We allow the system to evolve for $t=150 t_0$. The torus is MRI unstable and after approximately one revolution accretion sets in. The strong shear generates an azimuthal field component and the angular momentum is redistributed. Due to the instability, fillamentary structures form at the torus surface that account for its rich morphology (Figure 34.49). These results can be promtly compared to those in Hawley (2000) and Mignone et al. (2007).

Figure: Left: 3D rendering of the axisymmetric torus evolution. Right: Density logarithm of the magnetized accretion torus after 150 $t_0$. Superposed are the AMR levels and the mesh.

34.2.4 Magnetized Noh Z-pinch

In this test we consider the magnetized version of the classical Noh problem (Noh, 1987). It consists of a cylindrically symmetric implosion of a pressure-less gas: the gas stagnates at the symmetry axis and creates an outward moving accretion shock which propagates at a constant velocity. The self-similar analytic solution of this problem has been widely used as benchmark test for hydrodynamic codes, especially those targeting implosion.

Recently, Velikovich et al. (2012) have extended the original test problem to include an embedded azimuthal magnetic field, in accordance with the Z-pinch physics. In their study, they present a family of exact analytic solutions for which the values of primitive variables are finite everywhere, providing an excellent benchmark for Z-pinch and general MHD codes.

We perform this test in both cartesian and cylindrical geometries. The tests can be found respectively in magnetoHD/Noh and magnetoHD/NohCylindrical. For brevity we describe only the cylindrical initialization. The cartesian can be easily recovered by projecting the solution onto the X-Y plane. A 3T MHD version of this test can also be found under magnetoHD/unitTest/NohCylindricalRagelike.

The simulation is initialized in a computational box that spans $[0, 3]$ cm in the $r$ and $z$ directions, in cylindrical ($r-z$) geometry. The leftmost boundary condition is set to axisymmetry, whereas the remaining boundaries are set to outflow (zero gradient). The initial condition in primitive variables is defined as

$$\displaystyle \rho = 3.1831\times10^{-5} r^2 g/cm^3,$$

$$\displaystyle {\bf v} = (-3.24101\times10^7, 0, 0) \textrm{cm/s},$$

$$\displaystyle {\bf B} = (0, 6.35584\times10^5 r, 0) \textrm{gauss},$$

$$\displaystyle P = C {\bf B}^2.$$ (34.19)

Since Godunov-type codes cannot run with zero pressure, we initialize $P$ by choosing $C=P/{\bf B}^2=10^{-6}$, ensuring a magnetically dominated plasma.

The simulations are evolved for $30$ ns, utilizing the unsplit solver and a Courant number of 0.8. We use 6 levels of refinement, corresponding to an equivalent resolution of $256\times 256$ zones. The reconstruction is piecewise-linear (second order) with characteristic limiting, whereas the Riemann solvers employed are the HLLD and Roe.

The resulting density profile is shown in Figure 34.50. The refinement closely follows the propagation of the discontinuity which reaches the analytically predicted location ($r=0.3$ cm) at $t=30$ ns. The lineouts for HLLD (green dots) and Roe (blue dots) shown on the right, display good agreement with the analytic solution (red line) and the discontinuity is sharply captured (resolved on two points). Our Figure 34.50 can be directly compared to Figures 2a and 3a of Velikovich et al. (2012).

Figure: Magnetized Noh problem in cylindrical coordinates. Left: density snapshot along with AMR levels after 30 ns for the HLLD Riemann solver. Right: Lineouts of density close to the origin, $r=[0, 0.6]$ cm, for HLLD and Roe, superposed on the analytic solution.

34.2.5 MHD Rotor

The two-dimensional MHD rotor problem (Balsara and Spicer, 1999), magnetoHD/Rotor, is designed to study the onset and propagation of strong torsional Alfvén waves, which is thereby relevant for star formation. The computational domain is a unit square $[0,1]\times[0,1]$ with non-reflecting boundary conditions on all four sides. The initial conditions are given by

$$\rho(x,y)=\left\{ \begin{array}{l@{\quad}l} 10 & r \leq r_0 1+9f(r) & r_0 < r < r_1 1 & r \geq r_1 \end{array} \right.$$ (34.20)

$$u(x,y)=\left\{ \begin{array}{l@{\quad}l} -f(r)u_0(y-0.5)/r_0 & r ... ...r_0 -f(r)u_0(y-0.5)/r & r_0 < r < r_1 0 & r \geq r_1 \end{array} \right.$$ (34.21)

$$v(x,y)=\left\{ \begin{array}{l@{\quad}l} f(r)u_0(x-0.5)/r_0 & r \leq r_0 f(r)u_0(x-0.5)/r & r_0 < r < r_1 0 & r \geq r_1 \end{array} \right.$$ (34.22)

$$p(x,y) = 1$$ (34.23)

$$B_x(x,y) = \frac{5}{\sqrt{4\pi}}$$ (34.24)

$$B_y(x,y) = 0,$$ (34.25)

where $r_0=0.1,r_1=0.115,r=\sqrt{(x-0.5)^2+(y-0.5)^2},w=B_z=0$ and a taper function $f(r)=\bigl(r_1-r\bigr)/\bigl(r-r_0\bigr)$. The value $\gamma=1.4$ is used. The initial set-up is occupied by a dense rotating disk at the center of the domain, surrounded by the ambient flow at rest with uniform density and pressure. The rapidly spinning rotor is not in an equilibrium state due to the centrifugal forces. As the rotor spins with the given initial rotating velocity, the initially uniform magnetic field in $x$-direction will wind up the rotor. The rotor will be wrapped around by the magnetic field, and hence start launching torsional Alfvén waves into the ambient fluid. The angular momentum of the rotor will be diminished in later times as the evolution time increases. The circular rotor will be progressively compressed into an oval shape by the build-up of the magnetic pressure around the rotor. The results shown in Figure 34.51 were obtained using the StaggeredMesh MHD solver using 6 refinement levels. The divergence free evolution of the magnetic fields are well preserved as illustrated in Figure 34.52.

Figure: The Rotor problem at $t=0.15$ (a) Density (b) Pressure (c) Mach number (d) Magnetic pressure.

[a] [b] [c] [d]

Figure: Divergence of magnetic fields using the StaggeredMesh solver at $t=0.15$ for the Rotor problem.

34.2.6 MHD Current Sheet

The two-dimensional current sheet problem, magnetoHD/CurrentSheet, has recently been studied by Gardiner and Stone (2005) in ideal MHD regime. The two current sheets are initialized and therefore magnetic reconnections are inevitably driven. In the regions the magnetic reconnection takes place the magnetic flux approaches vanishingly small values, and the loss in the magnetic energy is converted into heat (thermal energy). This phenomenon changes the overall topology of the magnetic fields and hence affects the global magnetic configuration.

The square computational domain is given as $[0,2]\times[0,2]$ with periodic boundary conditions on all four sides. We initialize two current sheets in the following:

$$B_y = \left\{ \begin{array}{l@{\quad}l} B_0 & 0.0 \leq x <0.5 -B_0 & 0.5 \leq x <1.5 B_0 & 1.5 \leq x \leq 2.0 \end{array} \right. ,$$ (34.26)

where $B_0=1$. The other magnetic field components $B_x, B_z$ are set to be zeros. The $x$ component of the velocity is given by $u=u_0\sin2\pi y$ with $u_0=0.1$, and all the other velocity components are initialized with zeros. The density is unity and the gas pressure $p=0.1$.

The changes of the magnetic fields seed the magnetic reconnection and develop formations of magnetic islands along the two current sheets. The small islands are then merged into the bigger islands by continuously shifting up and down along the current sheets until there is one big island left in each current sheet.

The temporal evolution of the magnetic field lines from $t=0.0$ to $t=5.0$ are shown in Figure 34.53(a) - Figure 34.53(f) on a $256\times 256$ uniform grid. In Figure 34.54 the same problem is resolved on an AMR grid with 6 refinement levels, showing the current density $j_z$ along with the AMR block structures at $t=4.0$. The StaggeredMesh solver was used for this problem.

Figure: The temporal evolutions of field lines for the MHD CurrentSheet problem. Equally spaced 60 contour lines are shown at time (a) $t=0.0$ (b) $t=1.0$ (c) $t=2.0$ (d) $t=3.0$ (e) $t=4.0$ (f) $t=5.0$.

[a] [b] [c] [d] [e] [f]

Figure: Current density at $t=4.0$ using the StaggeredMesh solver for the MHD CurrentSheet problem.

34.2.7 Field Loop

The 2D and 3D field loop advection problems (magnetoHD/FieldLoop) are known to be stringent test cases in multidimensional MHD. In this test problem we consider a 2D advection of a weakly magnetized field loop traversing the computational domain diagonally. Details of the problem has been described in Gardiner and Stone (2005).

The computational domain is $[-1,1]\times[-0.5,0.5]$, with a grid resolution $256\times148$, and doubly-periodic boundary conditions. With this rectangular grid cell, the flow is not symmetric in $x$ and $y$ directions because the field loop does not advect across each grid cell diagonally and hence the resulting fluxes are different in $x$ and $y$ directions. The density and pressure are unity everywhere and $\gamma=5/3$. The velocity fields are defined as,

$$\mathbf U=u_0( \theta, \theta,1)$$ (34.27)

with the advection angle $\theta$, given by $\theta=\tan^{-1}(0.5)\approx 26.57^{\circ}$. For the choice of the initial advection velocity we set $u_0=\sqrt 5$. The size of domain and other parameters were chosen such that the weakly magnetized field loop makes one complete cycle by $t=1$. It is important to initialize the magnetic fields to satisfy $\nabla \cdot\mathbf B=0$ numerically in order to avoid any initial nonzero error in $\nabla \cdot\mathbf B$. As suggested in Gardiner and Stone (2005), the magnetic field components are initialized by taking the numerical curl of the $z$-component of the magnetic vector potential $A_z$,

$$B_x=\frac{\partial A_z}{\partial y}, \;\;\;\; B_y=-\frac{\partial A_z}{\partial x},$$ (34.28)

where

$$A_z=\left\{ \begin{array}{l@{\quad}l} A_0(R-r) & r \leq R \\ 0 & \mbox{otherwise.} \end{array} \right. ,$$ (34.29)

By using this initialization process, divergence-free magnetic fields are constructed with a maximum value of $\nabla \cdot\mathbf B$ in the order of $10^{-16}$ at the chosen resolution. The parameters in (34.29) are $A_0=10^{-3}$ and a field loop radius $R=0.3$. This initial condition results in a very high plasma beta $\beta=p/B_p=2\times10^{6}$ for the inner region of the field loop. Inside the loop the magnetic field strength is very weak and the flow dynamics is dominated by the gas pressure.

The field loop advection is integrated to a final time $t=2$. The advection test is found to truly require the full multidimensional MHD approach (Gardiner and Stone, 2005, 2008; Lee and Deane, 2008). Since the field loop is advected at an oblique angle to the $x$-axis of the computational domain, the values of $\partial B_x / \partial x$ and $\partial B_y / \partial y$ are non-zero in general and their roles are crucial in multidimensional MHD flows. These terms, together with the multidimensional MHD terms $\mathbf{A}_{B_x}$ and $\mathbf{A}_{B_y}$, are explicitly included in the data reconstruction-evolution algorithm in the USM scheme (see Lee and Deane, 2008). During the advection a good numerical scheme should maintain: (a) the circular symmetry of the loop at all time: a numerical scheme that lacks proper numerical dissipation results in spurious oscillations at the loop, breaking the circular symmetry; (b) $B_z=0$ during the simulation: $B_z$ will grow proportional to $w\nabla \cdot \mathbf{B}\Delta t$ if a numerical scheme does not properly include multidimensional MHD terms.

From the results in Figure 34.55, the USM scheme maintains the circular shape of the loop extremely well to the final time step. The scheme successfully retains the initial circular symmetry and does not develop severe oscillations.

Figure: The field loop advection problem using the StaggeredMesh solver at time $t=2$ with the Roe Riemann solver. (a)$B_p$ with the MEC at $t=2$. The color scheme between $2.32\times 10^{-25}$ and $7.16\times 10^{-7}$ was used. (b)Magnetic field lines with the MEC at $t=2$. 20 contour lines of $A_z$ between $-2.16\times 10^{-6}$ and $2.7\times 10^{-4}$ are shown.

[a] [b]

A variant 3D version of this problem (Gardiner and Stone, 2008) is also available, and is illustrated in in Fig. 34.56. This problem is considered to be a particularly challenging test because the correct solution requires inclusion of the multidimensional MHD terms to preserve the in-plane dynamics. Otherwise, the failure in preserving the in-plane dynamics (i.e., growth in the out-of-plane component) results in erroneous behavior of the field loop. As shown in Fig. 34.56, the USM solver successfully maintains the circular shape of the field loop and maintains the out-of-plane component of the magnetic field to very small values over the domain. This figure compares very well with Fig. 2 in Gardiner and Stone, 2008.

Figure: Thresholded images of the field loop advection problem at times $t=0.0, 1.5,$ and $2.0$ using a uniform grid size $128\times 128\times 256$.

[$B_p$ at $t=0.0$] [$B_p$ at $t=1.5$] [$B_p$ at $t=2.0$]

34.2.8 3D MHD Blast

A 2D version of the MHD blast problem was studied by Zachary et al. (Zachary, Malagoli, and Colella, 1994) and we consider a variant 3D version of the MHD spherical blast wave problem here. This problem leads to the formation and propagation of strong MHD discontinuities, relevant to astrophysical phenomena where the magnetic field energy has strong dynamical effects. With a numerical scheme that fails to preserve the divergence-free constraint, unphysical states could be obtained involving negative gas pressure because the background magnetic pressure increases the strength of magnetic monopoles.

This problem can be computed in various magnetized flow regimes by considering different magnetic field strengths. The computational domain is a square $[-0.5, 0.5]\times[-0.5,0.5]\times[-0.5,0.5]$ with a maximum refinement level 4. The explosion is driven by an over-pressurized circular region at the center of the domain with a radius $r=0.1$. The initial density is unity everywhere, and the pressure of the ambient gas is $0.1$, whereas the pressure of the inner region is $1000$. The strength of a uniform magnetic field in the $x$-direction are 0, $50/\sqrt{4\pi}$ and $100/\sqrt{4\pi}$. This initial condition results in a very low-$\beta$ ambient plasma state, $\beta=2.513\times10^{-4}$. Through this low-$\beta$ ambient state, the explosion emits almost spherical fast magneto-sonic shocks that propagate with the fastest wave speed. The flow has $\gamma=1.4$.

With this strong magnetic field strength, $B_x=100/\sqrt {4\pi }$, shown in Figure 34.57, the explosion now becomes highly anisotropic as shown in the pressure plot in Figure 34.57. The Figure shows that the displacement of gas in the transverse $y$-direction is increasingly inhibited and hydrodynamical shocks propagate in both positive and negative $x$-directions parallel to $B_x$. This process continues until total pressure equilibrium is reached in the central region. This problem is also available in 2D setup.

Figure: The MHD blast test at time $t=0.01$ using the unsplit staggered mesh MHD solver. Density (top half) and magnetic pressure (bottom half) plots for three different strengths of $B_x$. (a) $B_x = 0$ (b) $B_x=50/\sqrt {4\pi }$ (c) $B_x=100/\sqrt {4\pi }$.

[a] [b] [c]