14.3 Magnetohydrodynamics (MHD)

14.3.1 Description

The FLASH4 code includes two magnetohydrodynamic (MHD) units that represent two different algorithms. The first is the eight-wave model (8Wave) by Powell et al. (1999) that is already present in FLASH2. The second is a newly implemented unsplit staggered mesh algorithm (USM or StaggeredMesh). It should be noted that there are several major differences between the two MHD units. The first difference is how each algorithm enforces the solenoidal constraint of magnetic fields. The eight-wave model basically uses the truncation-error method, which effectively removes the effects of unphysical magnetic monopoles if they are generated during simulations. It does not, however, completely eliminate monopoles that are spurious in a strict physical law. To improve such unphysical effects in simulations, the unsplit staggered mesh algorithm uses the constrained transport method (Evans and Hawley, 1988) to enforce divergence-free constraints of magnetic fields. This method is shown to maintain magnitudes of $\nabla\cdot\textbf{B}\xspace$ substantially low, e.g., to the orders of $10^{-12}$ or below, in most simulations. The second major difference is that the unsplit staggered mesh algorithm uses a directionally unsplit scheme to evolve the MHD governing equations, whereas the eight-wave method uses a directionally splitting method as in FLASH2. In general, the splitting method is shown to be robust, relatively straightforward to implement, and generally faster than the unsplit method. The splitting method, however, does generally introduce splitting errors when solving one-dimensional subproblems in each sweep direction for multidimensional MHD equations. This error gets introduced in simulations because (i) the linearized Jacobian flux matrices do not commute in most of the nonlinear multidimensional problems (LeVeque, 1992; LeVeque, 1998), and (ii) in MHD, dimensional-splitting based codes are not able to evolve the normal (to the sweep direction) magnetic field during each sweep direction (Gardiner and Stone, 2005).

Note that the eight-wave solver uses the same directionally splitting driver unit Driver/DriverMain/split as the PPM and RHD units do, while the unsplit staggered mesh solver (StaggeredMesh) has its own independent unsplit driver unit Driver/DriverMain/unsplit.

Both MHD units solve the equations of compressible ideal and non-ideal magnetohydrodynamics in one, two and three dimensions on a Cartesian system. Written in non-dimensional (hence without $4\pi$ or $\mu_0$ coefficients) conservation form, these equations are

$$\frac{\partial\rho}{\partial t} + \nabla\cdot(\rho\textbf{v}\xspace ) = 0$$ (14.19)

$$\frac{\partial\rho\textbf{v}\xspace }{\partial t} + \nabla\cdot (\rho\textbf{v}\xspace \textbf{v}\xspace -\textbf{B}\... ...xtbf{B}\xspace )+\nabla p_{\ast} = \rho\textbf{g}\xspace + \nabla\cdot{\bf\tau}$$ (14.20)

$$\frac{\partial \rho E}{\partial t} + \nabla\cdot(\textbf{v}\xspace (\rho E+ p_{\ast})- \textbf{B}\xspa... ... T) + \nabla\cdot(\textbf{B}\xspace \times(\eta\nabla\times\textbf{B}\xspace ))$$ (14.21)

$$\frac{\partial \textbf{B}\xspace }{\partial t} + \nabla\cdot(\textbf{v}\xspace \textbf{B}\xspace -\textbf{B}\xspace \textbf{v}\xspace ) = -\nabla\times(\eta\nabla\times\textbf{B}\xspace )$$ (14.22)

where

$$p_{\ast} = p+\frac{B^2}{2},$$ (14.23)

$$E = \frac{1}{2}v^2+\epsilon+\frac{1}{2}\frac{B^2}{\rho},$$ (14.24)

$${\bf\tau} = \mu \left( (\nabla\textbf{v}\xspace )+(\nabla\textbf{v}\xspace )^{\mbox{\tiny T}}- \frac{2}{3}(\nabla\cdot\textbf{v}\xspace )\mathbf{I} \right)$$ (14.25)

are total pressure, specific total energy and viscous stress respectively. Also, $\rho$ is the density of a magnetized fluid, $\textbf{v}\xspace$ is the fluid velocity, $p$ is the fluid thermal pressure, $T$ is the temperature, $\epsilon$ is the specific internal energy, $\textbf{B}\xspace$ is the magnetic field, $\textbf{g}\xspace$ is the body force per unit mass, for example, due to gravity. $\tau$ is the viscosity tensor, $\mu$ is the coefficient of viscosity (dynamic viscosity), $\mathbf{I}$ is the unit (identity) tensor, $\sigma$ is the heat conductivity, and $\eta$ is the resistivity. The thermal pressure is a scalar quantity, so that the code is suitable for simulations of ideal plasmas in which magnetic fields are not so strong that they cause temperature anisotropies. As in regular hydrodynamics, the pressure is obtained from the internal energy and density using the equation of state. The two MHD units support general equations of state and multi-species fluids. Also, in order to prevent negative pressures and temperatures, a separate equation for internal energy is solved in a fashion described earlier in the hydrodynamics chapter.

The APIs of the MHD units are fairly minimal. The units honor all of hydrodynamics unit variables, interface functions and runtime parameters described in the above hydrodynamics unit chapter (see Chp:Hydrodynamics Unit). In addition, both the eight-wave and the unsplit staggered mesh units declare additional plasma variables and runtime parameters, which are listed in Table 14.7 and Table 14.8.

14.3.2 Usage

In the current release, the eight-wave unit serves as a default MHD solver. In order to choose the unsplit staggered mesh unit for MHD simulations, users need to include +usm in a setup script. The default eight-wave unit will be automatically chosen if there is no such specification included.

A word of caution: The eight-wave solver is only compatible with native grid interpolation in AMR simulations. This is because the solver only uses two layers of guard cells in each coordinate direction. The choice -gridinterpolation=native is automatically adopted if +8wave is specified in setup, otherwise, -gridinterpolation=native should be explicitly included in order to use the eight-wave solver without specifying +8wave. For instance, running a script ./setup magnetoHD/BrioWu -1d -auto +8wave will properly setup the BrioWu problem for the 8Wave solver, and ./setup magnetoHD/BrioWu -1d -auto +usm for the StaggeredMesh solver.

Supported configurations: Both MHD units currently support the uniform grid with FIXEDBLOCKSIZE and NONFIXEDBLOCKSIZE modes, and the adaptive grid with PARAMESH on Cartesian geometries, as well as 2D cylindrical (R-Z). When using AMR grids, the eight-wave unit supports both PARAMESH 2 and PARAMESH 4, while only PARAMESH 4 is supported in the unsplit staggered mesh solver because face-centered variables are only fully supported in PARAMESH 4.

Table 14.7: Additional solution variables used in the MHD units.

Variable	Type	Description
magx	PER_VOLUME	$x$-component of magnetic field
magy	PER_VOLUME	$y$-component of magnetic field
magz	PER_VOLUME	$z$-component of magnetic field
magp	(GENERIC)	magnetic pressure
divb	(GENERIC)	divergence of magnetic field

Table 14.8: Additional runtime parameters used in the MHD units.

Variable	Type	Default	Description
UnitSystem	string	“none”	System of units in which MHD calculations are to be performed. Acceptable values are “none” “CGS” and “SI”.
killdivb	logical	.true.	Enable/disable divergence cleaning.
flux_correct	logical	.true.	Enable/disable flux correction on AMR grid.

14.3.3 Algorithm: The Unsplit Staggered Mesh Solver

A directionally unsplit staggered mesh algorithm (USM), which solves ideal and non-ideal MHD governing equations (14.19) $\sim$ (14.22) in multiple dimensions, is a new MHD solver. Since FLASH4-beta, a full 3D implementation has been included as a new official release (Lee, 2013). The unsplit staggered mesh unit is based on a finite-volume, high-order Godunov method combined with a constrained transport (CT) type of scheme which ensures the solenoidal constraint of the magnetic fields on a staggered mesh geometry. In this approach, the cell-centered variables such as the plasma mass density $\rho$, plasma momentum density $\rho\textbf{v}\xspace$ and total plasma energy $\rho E$ are updated via a second-order MUSCL-Hancock unsplit space-time integrator using the high-order Godunov fluxes. The rest of the cell face-centered (staggered) magnetic fields are updated using Stokes' Theorem as applied to a set of induction equations, enforcing the divergence-free constraint of the magnetic fields. Notice that this divergence-free constraint is automatically guaranteed and satisfied in pure one-dimensional MHD simulations, but special care must be taken in multidimensional problems.

The overall procedure of the unsplit staggered mesh scheme can be broken up into the following steps (Lee, 2006; Lee and Deane, 2009; Lee, 2013):

• Quasi-linearization: This step replaces the nonlinear system (14.19) $\sim$ (14.22) with an approximate, quasi-linearized system of equations.

• Data Reconstruction-evolution: This routine calculates and evolves cell interface values by half time step using a second-order MUSCL-Hancock TVD algorithm (Toro, 1999). The approach makes use of a new method of 'multidimensional characteristic analysis' that can be achieved in one single step, incorporating all flux contributions from both normal and transverse directions without requiring any need of solving a set of Riemann problems (that is usually adopted in transverse flux updates). In this step the USM scheme includes the multidimensional MHD terms in both normal and transverse directions, satisfying a perfect balance law for the terms proportional to $\nabla \cdot \mathbf{B}$ in the induction equations.

• An intermediate Riemann problem: An intermediate set of high-order Godunov fluxes is calculated using the cell interface values obtained from the data reconstruction-evolution step. The resulting fluxes are then used to evolve the normal fields by a half time step in the next procedure.

• A half time step update for the normal fields: The normal magnetic fields are evolved by a half time step using the flux-CT method at cell interfaces, ensuring the divergence-free property on a staggered grid. This intermediate update for the normal fields and the half time step data from the data reconstruction-evolution step together provide a second-order accurate MHD states at cell interfaces.

• Riemann problem: Using the second-order MHD states calculated from the above procedures, the scheme proceeds to solve the Riemann problem to obtain high-order Godunov fluxes at cell interfaces.

• Unsplit update of cell-centered variables: The unsplit time integrations are performed using the high-order Godunov fluxes to update the cell-centered variables for the next time step.

• Construction of electric fields: Using the high-order Godunov fluxes, the cell-cornered (edged in 3D) electric fields are constructed. The unsplit staggered mesh scheme computes a new modified electric field construction (MEC) scheme that includes first and second multidimensional derivative terms in Taylor expansions for high-order interpolations. This modified electric field construction provides enhanced accuracy by explicitly adding proper amounts of dissipation as well as spatial gradients in its interpolation scheme.

• Flux-CT scheme: The electric fields from the MEC scheme are used to evolve the cell face-centered magnetic fields by solving a set of discrete induction equations. The resulting magnetic fields satisfy the divergence-free constraint up to the accuracy of machine round-off errors.

• Reconstruct cell-centered magnetic fields: The cell-centered magnetic fields are reconstructed from the divergence free cell face-centered magnetic fields by taking arithmetic averages of the cell face-centered fields variables.

Note that the procedure required in solving one-dimensional MHD equations is much simpler than solving the multidimensional ones and only involves the first through third and the fifth steps in the above oulined scheme. The choices of TVD slope limiters available in the unsplit staggered mesh scheme (see Table 14.9) includes the minmod limiter as well as the compressible limiters such as vanLeer or mc limiter. Another choice, called hybrid limiter, can be used to provide a mixed type of limiters as described in Balsara (2004). In this choice, one uses a compressible limiter to produce a crisp representation for linearly degenerate waves (e.g., an entropy wave and left- and right-going Alfvén waves). To this end, a compressible limiter can be applied to the density and the magnetic fields variables, where these variables contribute much of the variations in such linearly degenerate waves. Other variables, the velocity field components and pressure, constitute four genuinely nonlinear wave families (i.e., left- and right-going fast/slow magneto-sonic waves) in MHD. These genuinely nonlinear wave families inherently behave according to their self steepening mechanism and one can simply use a diffusive but robust minmod limiter. Another limiter, called limited, is also available (see details in Toro, 1999, 2nd Ed., section 13.8.4), and users need to specify a runtime parameter $\beta$ (LimitedSlopeBeta in flash.par) if this limiter is chosen for a simulation.

The unsplit staggered mesh unit solves a set of discrete induction equations in multi-dimensional problems to proceed temporal evolutions of the staggered magnetic fields using electric fields. For instance, in a two-dimensional staggered grid, the unsplit staggered mesh unit solves a two-dimensional pair of discrete induction equations that were found originally by Yee (1966):

$${b}^{n+1}_{x,i+{1}/{2},j} = {b}^{n}_{x,i+{1}/{2},j}- \frac{\Delta t}{\Delta y} \Bigr\{ {E}^{n+1/2}_{z,i+{1}/{2},j+{1}/{2}} -{E}^{n+1/2}_{z,i+{1}/{2},j-{1}/{2}} \Bigl\},$$ (14.26)

$${b}^{n+1}_{y,i,j+{1}/{2}} = {b}^{n}_{y,i,j+{1}/{2}}- \frac{\Delta t}{\Delta x} \Bigr\{ -{E}^{n+1/2}_{z,i+{1}/{2},j+{1}/{2}} +{E}^{n+1/2}_{z,i-{1}/{2},j+{1}/{2}} \Bigl\}.$$ (14.27)

The superindex $n+1/2$ in the above equations simply indicates an intermediate timestep right after the temporal update of the cell-centered variables.

A three-dimensional schematic figure of the staggered grid geometry with collocations of edge-based values (electric fields $E$) and face based values (magnetic fields $b$) is shown in Figure 14.2.

One of the main advantages of using the CT-type of scheme is that the cell face-centered magnetic fields ${b}^{n+1}_{x,i+{1}/{2},j}$ and ${b}^{n+1}_{y,i,j+{1}/{2}}$, which are updated via the above induction equations, satisfy the divergence-free constraint locally. The numerical divergence of the magnetic fields is defined as

$$(\nabla\cdot\textbf{B}\xspace )^{n+1}_{i,j}= \frac{{b}^{n+1}_{x,i... ...Delta x}+ \frac{{b}^{n+1}_{y,i,j+{1}/{2}}-{b}^{n+1}_{y,i,j-{1}/{2}} }{\Delta y}$$ (14.28)

and it remains zero to the accuracy of machine round-off errors, provided that $nabla\cdot\textbf{B}\xspace )^{n}_{i,j}=0$.

On an AMR grid, the unsplit staggered mesh scheme uses a direct injection method as a default to preserve divergence-free prolongation to the cell face-centered fields variables. This method is one of the simplest approaches that is offered by PARAMESH 4 to maintain the divergence-free constraint in prolongation. This simple method ensures the solenoidal constraint well enough where the fields are varying smoothly, but can introduce oscillations in regions of steep field gradient. In such cases Balsara's prolongation algorithm can be useful. Both prolongation algorithms are supported and enabled using runtime parameters in the unsplit staggered mesh solver (see Table 14.9 below).

To solve the above induction equations (14.26) and (14.27) in a flux-CT type scheme, it is required to construct cell edge-based electric fields. The simplest choice is to use the cell face-centered high-order Godunov fluxes and take an arithmetic average to construct cell-cornered (edge-based in 3D) electric fields:

$${E}^{n+1/2}_{z,i+{1}/{2},j+{1}/{2}} = \frac{1}{4}\Bigr\{-F^{n+{1}/{2}}_{B_y,i+{1}/{2},j}- F^{n+{1}/{2}}... ...1} + G^{n+{1}/{2}}_{B_x,i,j+{1}/{2}}+ G^{n+{1}/{2}}_{B_x,i+1,j+{1}/{2}} \Bigl\}$$

$$= \frac{1}{4}\Bigr\{E^{n+{1}/{2}}_{z,i+{1}/{2},j}+ E^{n+{1}/{2}}_{z... ...},j+1}+ E^{n+{1}/{2}}_{z,i,j+{1}/{2}}+ E^{n+{1}/{2}}_{z,i+1,j+{1}/{2}} \Bigl\},$$ (14.29)

where $F_{B_y}$ and $G_{B_x}$ represent the $x$ and $y$ high-order Godunov flux components corresponding to the magnetic fields $B_y$ and $B_x$, respectively (see details in Balsara and Spicer, 1999).

A high-order accurate version is also available by turning on a logical switch E_modification in the unsplit staggered mesh scheme, which takes Taylor series expansions of the cell-cornered electric field ${E}^{n+1/2}_{z,i+{1}/{2},j+{1}/{2}}$ in all directions, followed by taking an arithmetic average of them (Lee, 2006; Lee and Deane, 2009).

The last step in the unsplit staggered mesh scheme is to reconstruct the cell-centered magnetic fields $B_{x,i,j}$ and $B_{y,i,j}$ from the divergence-free face-centered magnetic fields. The unsplit staggered mesh scheme takes arithmetic averages of the face-centered fields variables to obtain the cell-centered magnetic fields, which is sufficient for second order accuracy. After obtaining the new cell-centered magnetic fields, the total plasma energy may need to be corrected in order to preserve the positivity of the thermal temperature and pressure (Balsara and Spicer, 1999; Tóth, 2000). This energy correction is very useful especially in problems involving very low $\beta$ plasma flows.

There are several choices available for calculating high order Godunov fluxes in the unsplit staggered mesh scheme. The default solver is Roe's linearized approximate solver, which takes into account all seven waves that arise in the MHD equations. The Roe solver can adopt one of the two entropy fix routines (Harten, 1983; Harten and Hyman, 1983) in ordetblrefr to avoid unphysical states near strong rarefaction regions. As all seven waves are considered in Roe's solver, high numerical resolutions can be achieved in most cases. However, Roe's solver still can fail to maintain positive states near very low densities even with the entropy fix. In this case, computationally efficient and positively conservative Riemann solvers such as HLL (Einfeldt et al., 1991), HLLC (S. Li, 2005), or HLLD (Miyoshi and Kusano, 2005) can be used to maintain positive states in simulations. A hybrid type of Riemann solver which combines using the Roe solver for high accuracy and HLLD for stability is also available.

The USM solver has been recently extended also for 2D and 2.5D cylindrical (R-Z) geometries, both for uniform grids and AMR. In the cylindrical implementation, we followed the guidelines of Mignone et al. (2007) and Skinner & Ostriker (2010). Special care was also taken to ensure a divergence free interpolation of the staggered magnetic field components, when grid movements occur in AMR. This novel prolongation scheme is based on the methods described in Balsara (2001, 2004) and Li & Li (2004). More information regarding the cylindrical implementation can be found in Tzeferacos et al. (2012, in print) whereas new test problems, provided with this release, are available at Sections 34.2.3 and 34.2.4. Handling of different geometries will be available in future releases.

14.3.3.1 Slowly moving shock handling in PPM

A new dissipative mechanism is a hybridized slope limiter for PPM that combines a new upwind biased slope limiter with a conventional TVD slope limiter (Lee, D., Astronum Proc. 2010). This hybrid upwind limiter reduces spurious numerical oscillations near discontinuities, and therefore can compute sharp, monotonized profiles in compressible flows when using PPM, especially in Magnetohydrodynamics (MHD) slowly moving shock regions. (See more in Chapter 34.2.1.)

By the nature of very small numerical dissipations in PPM, unphysical oscillations in discontinuous MHD solutions can appear in a specific flow region, referred to as a slowly moving shock (SMS) relative to the grid. This new approach handles numerical non-oscillatory MHD solutions using PPM in SMS regions. The SMS should not be confused with so-called ”slow MHD shock” which corresponds to two slow waves (i.e., $u \pm c_s$) in MHD, where $c_s$ is the slow magneto-acoustic velocity.

The method first detects a local, slowly moving shock, and considers an upwind direction to compute a monotonicity-preserving slope limiter. This new approach, in addition to improving the numerical solutions in MHD to levels that reduce (or eliminate) such oscillatory behaviors while preserving sharp discontinuities in MHD, is also simple to implement. The method has been verified against the results from other high-resolution shock-capturing (HRSC) methods such as MUSCL and WENO schemes.

In order to enable the SMS treatment for PPM, users set a runtime parameter upwindTVD=.true. in flash.par. The SMS method does require to have an extended stencil, and users should specify +supportPPMUpwind in setup. See more in Section 14.1.3.

Figure 14.2: A 3D control volume on the staggered grid with the cell center at $(i,j,k)$. The magnetic fields are collocated at the cell face centers and the electric fields at the cell edge centers. The line integral of the electric fields $\int_{{\partial \mathcal{F}_n}} \mathbf{E}\cdot\mathbf{T}dl$ along the four edges of the face $\mathcal{F}_{x,i+1/2,j,k}$ gives rise to the negative of the rate of change of the magnetic field flux in $x$-direction through the area enclosed by the four edges (i.e., the area of $\mathcal{F}_{x,i+1/2,j,k}$).

Table 14.9: Runtime parameters used in the unsplit staggered mesh MHD (physics/Hydro/HydroMain/unsplit/MHD_StaggeredMesh) solver additional to those described for the unsplit hydro solver (physics/Hydro/HydroMain/unsplit/Hydro_Unsplit).

Variable	Type	Default	Description
killdivb	logical	.true.	On/off $\nabla \cdot \mathbf{B}=0$ handling on the staggered grid
E_modification	logical	.true.	Enable/disable high-order electric field construction
E_upwind	logical	.false.	Enable/disable an upwind update for induction equations
energyFix	logical	.true.	Enable/disable energy correction
facevar2ndOrder	logical	.true.	Turn on/off a second-order facevar update
ForceHydroLimit	logical	.false.	On/off pure Hydro mode
prolMethod	string	"INJECTION_PROL"	Use either direct injection method ("INJECTION_PROL") or Balsara's method ("BALSARA_PROL") in prolonging divergence-free magnetic fields stored in face-centered variables
RiemannSolver	string	"ROE"	"HLLD" is additionally available for MHD, "Hybrid" is also available for MHD.

Stability limit: As described in the unsplit hydro solver unit (physics/Hydro/HydroMain/unsplitHydro_Unsplit), the USM MHD solver can take a wide range of CFL limits in all three dimensions (i.e., CFL $<$ 1). However, in some circumstances where there are strong shocks and rarefactions, shockLowerCFL=.true. could be useful to gain more numerical stability by using. It is also helpful to use (1) artificial viscosity and flattening, or (2) lower order reconstruction scheme (e.g., MH), or (2) diffusive Riemann solver such as HLL-type, or LLF solvers, or (3) a reduced CFL accordingly.

Divergence-free prolongation of magnetic fields on AMR in the unsplit staggered mesh solver: It is of importance to preserve divergence-free evolutions of magnetic fields in MHD simulations. Moreover, some special cares are required in prolonging divergence-free magnetic fields on AMR grids. One simple straightforward way in this aspect is to prolong divergence-free fields to newly created children blocks using direct injection. This injection method therefore inherently preserves divergence-free properties on AMR block structures and works well in most cases. This method is a default in the unsplit staggered mesh solver and can also be enabled by setting a runtime parameter prolMethod = "INJECTION_PROL". Another way, proposed by Balsara (2001), is also available in the unsplit staggered mesh solver and can be chosen by setting prolMethod = "BALSARA_PROL". Both prolongation methods are supported in MHD's 2.5D and 3D simulations. In 2 and 2.5D cylindrical geometry however, since neither method takes into account geometrical factors, we use a modified prolongation algorithm based on Balsara (2004) and Li&Li (2004). This is the default option and is activated by choosing prolMethod = "BALSARA_PROL". The need for this special refinement requires to have an MHD's own customized implementation of Simulation_customizeProlong.F90 placed in the source/Simulation/SimulationMain/magnetoHD/.

14.3.4 Algorithm: The Eight-wave Solver

The eight-wave magnetohydrodynamic (MHD) unit is based on a finite-volume, cell-centered method that was proposed by Powell et al. (1999). The unit uses directional splitting to evolve the magnetohydrodynamics equations. Like the PPM and RHD units, this MHD unit makes one sweep in each spatial direction to advance physical variables from one time level to the next. In each sweep, the unit uses AMR functionality to fill in guard cells and impose boundary conditions. Then it reconstructs characteristic variables and uses these variables to compute time-averaged interface fluxes of conserved quantities. In order to enforce conservation at jumps in refinement, the unit makes flux conservation calls to AMR, which redistributes affected fluxes using the appropriate geometric area factors. Finally, the unit updates the solution and calls the EOS unit to ensure thermodynamical consistency.

After all sweeps are completed, the MHD unit enforces magnetic field divergence cleaning. In the current release only a diffusive divergence cleaning method (truncation-error method), which was the default method in FLASH2, is supported and the later release of the code will incorporate an elliptic projection cleaning method.

The ideal part of the MHD equations are solved using a high-resolution, finite-volume numerical scheme with MUSCL-type (van Leer, 1979) limited gradient reconstruction. In order to maximize the accuracy of the solver, the reconstruction procedure is applied to characteristic variables. Since this may cause certain variables, such as density and pressure, to fall outside of physically meaningful bounds, extra care is taken in the limiting step to prevent this from happening. All other variables are calculated in the unit from the interpolated characteristic variables.

In order to resolve discontinuous Riemann problems that occur at computational cell interfaces, the code by default employs a Roe-type solver derived in Powell et al. (1999). This solver provides full characteristic decomposition of the ideal MHD equations and is, therefore, particularly useful for plasma flow simulations that feature complex wave interaction patterns. An alternative Riemann solver of HLLE type, provided in FLASH2, has not been incorporated into FLASH4 yet.

Time integration in the MHD unit is done using a second-order, one-step method due to Hancock (Toro, 1999). For linear systems with unlimited gradient reconstruction, this method can be shown to coincide with the classic Lax-Wendroff scheme.

A difficulty particularly associated with solving the MHD equations numerically lies in the solenoidality of the magnetic field. The notorious $\nabla\cdot\textbf{B}\xspace =0$ condition, a strict physical law, is very hard to satisfy in discrete computations. Being only an initial condition of the MHD equations, it enters the equations indirectly and is not, therefore, guaranteed to be generally satisfied unless special algorithmic provisions are made. Without discussing this issue in much detail, which goes well beyond the scope of this user's guide (for example, see Tóth (2000) and references therein), we will remind the user that there are three commonly accepted methods to enforce the $\nabla\cdot\textbf{B}\xspace$ condition: the elliptic projection method (Brackbill and Barnes, 1980), the constrained transport method (Evans and Hawley, 1988), and the truncation-level error method (Powell et al., 1999). In the current release, the truncation-error cleaning methods is provided for the eight-wave MHD unit, and the constrained transport method is implemented for the unsplit staggered mesh MHD units (see Sec:usm_algorithm for details).

In the truncation-error method, the solenoidality of the magnetic field is enforced by including several terms proportional to $\nabla\cdot\textbf{B}\xspace$. This removes the effects of unphysical magnetic tension forces parallel to the field and stimulates passive advection of magnetic monopoles, if they are spuriously created. In many applications, this method has been shown to be an efficient and sufficient way to generate solutions of high physical quality. However, it has also been shown (Tóth, 2000) that this method can sometimes, for example in strongly discontinuous and stagnated flows, lead to accumulation of magnetic monopoles, whose strength is sufficient to corrupt the physical quality of computed solutions. In order to eliminate this deficiency, the eight-wave MHD code also uses a simple yet very effective method originally due to Marder (1987) to destroy the magnetic monopoles on the scale on which they are generated. In this method, a diffusive operator proportional to $\nabla\nabla\cdot\textbf{B}\xspace$ is added to the induction equation, so that the equation becomes

$$\frac{\partial \textbf{B}\xspace }{\partial t}+\nabla\cdot(\textb... ...space \nabla\cdot\textbf{B}\xspace +\eta_a\nabla\nabla\cdot\textbf{B}\xspace ,$$ (14.30)

with the artificial diffusion coefficient $\eta_a$ chosen to match that of grid numerical diffusion. In the FLASH code, $\eta_a=\frac{\lambda}{2} (\frac{1}{\Delta x}+\frac{1}{\Delta y}+\frac{1}{\Delta z})^{-1}$, where $\lambda$ is the largest characteristic speed in the flow. Since the grid magnetic diffusion Reynolds number is always on the order of unity, this operator locally destroys magnetic monopoles at the rate at which they are created. Recent numerical experiments (Powell et al., 2001) indicate that this approach can very effectively bring the strength of spurious magnetic monopoles to levels that are sufficiently low, so that generated solutions remain physically consistent. The entire $\nabla\cdot\textbf{B}\xspace$ control process is local and very inexpensive compared to other methods. Moreover, one can show that this process is asymptotically convergent (Munz et al., 2000), and each of its applications is equivalent to one Jacobi iteration in solving the Poisson equation in the elliptic projection method. The caveat is that this method only suppresses but does not completely eliminate magnetic monopoles. Whether this is acceptable depends on the particular physical problem.

As an alternative way to eliminate magnetic monopoles completely, the FLASH2 code includes an elliptic projection method, in which the unphysical divergence of the magnetic field can be removed to any desired level down to machine precision. As yet, this method has not been made available in FLASH4 and will be supported in a later release.

14.3.5 Extended MHD

FLASH includes a complete generalized Ohm's law that incorporates all extended MHD terms of the Braginskii formulation. In other words, Hall, Biermann battery, anisotropic magnetic resistivity, Seebeck, Nernst, and cross-field terms are all included as shown in the generalized Ohm's Law below (Eq. 14.31):

$$\mathbf{E} = - \mathbf{u} \times \mathbf{B} +\frac{\mathbf{J}}{n_e e} \times \mathbf{B} - \frac{\nabla P_e}{n_e e}$$ (14.31)

$$+ \frac{\eta_\parallel}{\epsilon_0 \omega^2_{pe} \tau_{ei}}\mathbf{... ...{\eta_\wedge}{\epsilon_0 \omega^2_{pe} \tau_{ei}}(\mathbf{b} \times \mathbf{J})$$

$$- \frac{k\beta_\parallel}{e} \mathbf{b}(\mathbf{b}\cdot\nabla T_e)-... ... T_e)\times \mathbf{b}) - \frac{k\beta_\wedge}{e}(\mathbf{b} \times \nabla T_e)$$

where $\mathbf{E}$ is the electric field, $\mathbf{B}$ the magnetic field, $\mathbf{u}$ the plasma velocity, $\mathbf{J}$ the current density, $n_e$ the electron number density, $e$ the electron charge, $P_e$ the electron pressure, $k$ is Boltzmann's constant, $\mathbf{\hat{b}}$ the unit vector in the direction of the magnetic field, $\omega_{pe}$ the electron plasma frequency, $\epsilon_0$ the permittivity of vacuum, $\tau_{ei}$ the electron-ion collision time, $T_e$ the electron temperature, and $\eta$ and $\beta$ are respectively the resistivity and the thermoelectric coefficients in dimensionless form.

Each of these terms can be enabled or disabled in FLASH at run time via the flash.par file. For example, a typical flash.par file for extended MHD runtime parameters should look more or less like this:

# Hall effect -----------------------------------------
use_Hall = .true. 

# Biermann battery (Source implementation)-------------
use_Biermann = .true.  
hy_biermannSource = .true.

# Biermann battery (Flux implementation)-------------
use_Biermann = .true.  
use_Biermann3T = .true.

# Magnetic Resistivity --------------------------------
useMagneticResistivity  = .true.
resistivitySolver = "explicit"
resistivityForm = "fullaniso"
useCrossMagRes = .true.  # the eta_wedge term


# Thermoelectric Terms --------------------------------
use_Seebeck  = .true.
use_CrossField  = .true.
use_Nernst  = .true.

One way of simulating a diffusive process is to add those physics units in a simulation Config file. For example, a snippet of Config file can look like this:

# Diffusive time step calculation and implicit solver(s)
REQUIRES physics/Diffuse/DiffuseMain/Unsplit

# Magnetic Resistivity
REQUIRES physics/materialProperties/MagneticResistivity/MagneticResistivityMain/DaviesWen

# Thermoelectric
REQUIRES physics/materialProperties/Thermoelectric/ThermoelectricMain/DaviesWen

For Hall effect, Biermann battery, and Thermoelectric transport, there is no specific units to add, the only requirement is to require "source/physics/Hydro/HydroMain/unsplit/MHD_StaggeredMesh" (+usm) and for Biermann battery the MultiTemp implementation (+3t). Typically, for HEDP simulations, the setup will be done using "+usm3t". Also, for realistic Thermoelectric transport, a coefficients implementation should be chosen in the Config file as shown above. Below we detail the Hall, Biermann, Resistivity, and Thermoelectric implementations.

Including diffusive terms: New treatments has been made in including the non-ideal diffusive terms in the USM-MHD solver in the FLASH3.2 release and remain unchanged since then. In the previous releases before FLASH3.2, all the non-ideal diffusive terms had to be grouped together in the "resistive MHD" part of the unit by turning the flag variable useMagneticResistivity on, and the non-ideal terms were included all together. In the FLASH3.2 release, each individual term can be separately included by turning each corresponding logical variables in run time. The Diffusion time step is computed using the Diffuse_computeDt.F90 routine in Driver_computeDt.F90 to provide a more consistent way of computing a non-ideal time step.

14.3.5.1 Hall effect

14.3.5.1.1 How to use it ?

As indicated above (14.3.5), the user needs to include at least the Unsplit Staggered Mesh MHD solver ("+usm" during setup) in order to call the Hall implementation. Then the runtime parameter use_Hall has to be set to = .true.. Adding the Hall effect in the MHD model introduces two new kind of waves (see 14.3.5.1.5) that need to be taken into account for stability. It is why the cfl runtime parameter will generally have to be modified by the user for simulations including Hall physics. It is common for Hall simulations to use cfl values smaler than 0.1. The Hall effect is available in FLASH in 1D/2D/3D Cartesian geometries as well as in 2D Cylindrical geometry. We also require the user to set the runtime parameter conserveAngField to =.true. when using Hall effect in 2D cylindrical geometry.

14.3.5.1.2 The Hall effect in the MHD model

Implementing the Hall effect in a MHD code consists to add one term in the induction equation:

$$\frac{\partial \mathbf{B}}{\partial t} = - \nabla \times \mathbf{... ...}= \mathbf{E_{ind}} + \mathbf{E_{Fric}} + \mathbf{E_{Bier}} + \mathbf{E_{Hall}}$$ (14.32)

where each electric field is obtained from the generalized Ohm's law. $\mathbf{E_{ind}}=-\mathbf{u}\times \mathbf{B}$ is the inductive electric field, $\mathbf{E_{Fric}}$ is the electric field associated to the existence of friction between electrons and ions (one component of this field, and the one implemented in the FLASH code, is the ohmic electric field $\mathbf{E_{ohm}}=\eta \mathbf{j}$), $\mathbf{E_{Bier}}=- \nabla p_e /en_e$ is the Biermann electric field (discussed in the next section) and finally:

$$\mathbf{E_{Hall}} = \frac{\mathbf{j} \times \mathbf{B}}{en_e}$$ (14.33)

is the Hall electric field. We note that other fields exist in the generalized ohm's law (inertial...) but are not described here.

Because the electrical current $\mathbf{j}$ is defined as $\mathbf{j}=e n_e (\mathbf{u_i} - \mathbf{u_e})$, the sum of the inductive and Hall terms can be written as $\mathbf{E_{ind}} + \mathbf{E_{Hall}} = - \mathbf{u_e}\times \mathbf{B}$. This indicates that the "frozen-in" behavior is actually valid only for the electron population. However, in many cases (see next section), it is a good approximation to consider the field to be frozen in the bulk (ion) population. Here we are interested in cases where this approximation fails and thus where the presence of the Hall term is necessary.

In FLASH, the Hall effect is implemented by updating the face-centered fluxes for magnetic field and total energy equations. By firt writing the magnetic energy equation (in cgs):

$$\frac{\partial (B^2/8\pi)}{\partial t} + \nabla \cdot \left( \frac{c}{4 \pi} \mathbf{E}\times \mathbf{B}\right)= - \mathbf{j}\cdot \mathbf{E}$$ (14.34)

we notice that the Hall electric field is not responsible for energy generation source $\mathbf{j}\cdot \mathbf{E_{Hall}} = 0$ and only the Poynting vector ( $c\mathbf{E}\times \mathbf{B}/4\pi$) is affected by the presence of this additional electric field. This specificity of the Hall effect means that there is no need to modify other equations in the MHD model in order to satisfy total energy conservation (something we will have to do for Biermann). The Poynting contribution is added in FLASH in the face-centered total energy fluxes. Below we detail the flux formulation for the magnetic field equation, both in cartesian and cylindrical coordinates.

14.3.5.1.3 Formulation for a cartesian implementation

In cartesian coordinates, the induction equation 14.32 can be written as:

$$\frac{\partial\mathbf{B}}{\partial t} + \frac{\partial}{\partial ... ...e{ \begin{pmatrix} -E_y \\ E_x \\ 0 \end{pmatrix}}_\text{z-flux} = \mathbf{0}$$ (14.35)

or equivalently, using matrix notation:

$$\centering \frac{\partial \mathbf{B}}{\partial t} + \nabla \cdot \bar{\bar{\mathbf{F}}} = \mathbf{0}$$ (14.36)

where the flux matrix is expressed in terms of the electric field components:

$$\centering \bar{\bar{\mathbf{F}}} = \begin{bmatrix}0 & -E_z & E_y E_z & 0 & -E_x -E_y & E_x & 0 \end{bmatrix}$$ (14.37)

with the tensor divergence defined as $\nabla \cdot \bar{\bar{\mathbf{F}}}\vert_i = \partial_j F_{ji}$.

Expressing the Hall electric field components (from 14.33):

$$\centering \mathbf{E_{Hall}} = \begin{pmatrix}(j_yB_z - j_zB_y)/en_e (j_zB_x - j_xB_z)/en_e (j_xB_y - j_yB_x)/en_e \end{pmatrix}$$ (14.38)

the Hall flux matrix is explicitly given by:

$$\centering \bar{\bar{\mathbf{F}}}_{Hall} = \begin{bmatrix}0 & -(j... ...B_y)/en_e -(j_zB_x - j_xB_z)/en_e & (j_yB_z - j_zB_y)/en_e & 0 \end{bmatrix}$$ (14.39)

14.3.5.1.4 Formulation for a cylindrical implementation

In cylindrical (r,$\theta$,z) coordinates, the induction equations 14.32 writes as:

$$\centering \frac{\partial \mathbf{B}}{\partial t} = - \begin{pmat... ...rE_{\theta})}{\partial r} -\frac{1}{r}\frac{E_r}{\partial \theta} \end{pmatrix}$$ (14.40)

which can be written as:

$$\frac{\partial\mathbf{B}}{\partial t} + \frac{1}{r}\frac{\partial... ...1}{r} \underbrace{ \begin{pmatrix} 0 \\ E_z \\ 0 \end{pmatrix}}_\text{source}$$ (14.41)

the "source term" on the right side here (often called "geometrical" source term) comes from the fact that we can write $\partial E_z/\partial r = 1/r \partial (rE_z)/\partial r - E_z/r$. Eq.14.41 can also be written without including the source term and rather having everything written in terms of fluxes. It's actually what is being done currently in FLASH and that's why it requires the user to set the runtime parameter conserveAngField to = .true. for the 2D cylindrical case.

The Hall electric field in 2D axisymmetric cylindrical is explicitely given by:

$$\centering \mathbf{E_{Hall}} = \begin{pmatrix}(j_{\theta} B_z - j... ... (j_zB_r - j_rB_z)/en_e (j_xB_{\theta} - j_{\theta} B_x)/en_e \end{pmatrix}$$ (14.42)

14.3.5.1.5 Time step limitation due to the Hall effect

As mentionned before, the addition of the Hall term in the MHD model is responsible for the possibility for the plasma to sustain two additional propagation modes. The phase velocity of these modes are computed in Hydro_computeDt.F90 and compared with other present speeds (sound, Alfven, fast and slow magneto-acoustic...) to deterine the temporal time step dt. To study the possibility for the plasma to generate Hall-induced waves, we keep only the Hall term in the induction equation and express it in this form:

$$\frac{\partial \mathbf{B}}{\partial t} = \underbrace{ - \frac{1}{... ...e^{2}}\nabla n_e \times (\mathbf{j}\times \mathbf{B}) }_\text{Hall drift waves}$$ (14.43)

the first right side term is responsible for the generation of the so called "Whistler waves" whereas the second term is responsible for the generation of the "Hall drift waves".

It is easy to show that (see Huba, 2003), applying the plane wave perturbation on equation 14.43 as well as using the Maxwell-Ampere law in the non-relativistic limit ( $\mu_0 \mathbf{j}=\nabla \times \mathbf{B}$), the dispersion relation (in parallel propagation) for the whistler waves is given by:

$$\omega_{wh} = k_{\parallel} v_{A} \frac{k_{\parallel} c}{\omega_{pi}} = k_{\parallel} v_{A} (k_{\parallel}L_i)$$ (14.44)

It can be clearly seen from this relation that the Whislter waves (Hall-induced) will become faster than the Alfven waves (inductively-induced) if $k_{\parallel}L_i>1$ or expressed in an other way if the spatial scale of the perturbation ( $\sim k_{\parallel}^{-1}$) is smaller than the ion inertial length $L_i$. This demonstrate why hall physics becomes really relevant for a physical problem on scales smaller than $L_i$.

By looking at the second right side term in equation 14.43, it can be shown that the Hall drift waves propagates in inhomogeneous plasmas and in the $\mathbf{B} \times \nabla n_e$ direction (assumed to be the direction of the $\mathbf{k}$ vector). A simple dispersion relation can be obtained in the limit $kL_{n_e}«1$:

$$\omega_{hd} = k v_{A} \frac{ c}{L_{n_e}\omega_{pi}} = k v_{A} \frac{L_i}{L_{n_e}}$$ (14.45)

with $L_{n_e}=n_e/\vert\nabla n_e\vert$ being the spatial scale of the density gradient. Here again, the Hall-induced Hall drift waves will propagates faster than the Alfven speed only if $L_{n_e}<L_i$, that is to say if the spatial scale of the density gradient is smaller than the ion inertial length.

14.3.5.2 Biermann battery mechanism

14.3.5.2.1 How to use it ?

As indicated above (14.3.5), the user needs to include the Unsplit Staggered Mesh MHD solver ("+usm" during setup) in order to call the Biermann implementation. We currently offer support only for the 3T implementation ("+3t" during setup). Users wanting to use Biermann should thus setup their simulation with "+usm3t". Then the runtime parameter hy_useBiermann has to be set to = .true.. Two different versions of the Biermann battery mechanism are currently available in FLASH: a source version which should be considered as the most robust and reliable implementation as well as a flux-based version which should be used for specific problems with user-specified boundaries conditions. To use the source version, the user needs to set the runtime parameter hy_biermannSource to = .true. whereas for the flux-based version the runtime parameter use_Biermann3T has to be set to = .true.. If only hy_useBiermann = .true. is specified in the par file then the source version is selectionned by default. A strong difference between the two versions lies in the fact that the flux-based one introduces a new kind of wave and thus new stability considerations that have to be taken into account by the user through the runtime parameter cfl. The Biermann battery mechanism, for both versions, is available in FLASH in 2D/3D Cartesian geometries as well as in 2D Cylindrical geometry. We also require the user to set the runtime parameter conserveAngField to =.true. when using Biermann battery mechanism in 2D cylindrical geometry. Finally, we recommand to the user to always set the runtime parameter shockDetect to = .true. with the consequence of removing Biermann B-field generation in cells tagged as shock.

14.3.5.2.2 The Biermann battery mechanism in the MHD model

Implementing the Biermann battery mechanism in a MHD code consists of adding terms in the induction equation as well as in the internal energy equations. Regarding only the Biermann term, the induction equation writes (see eq.14.32) as:

$$\frac{\partial \mathbf{B}}{\partial t} = - \nabla \times \mathbf{E_{Bier}}$$ (14.46)

with the Biermann electric field given by:

$$\mathbf{E_{Bier}} = - \frac{\nabla p_e}{e n_e}$$ (14.47)

This mechanism describes the possibility for a plasma to self-generate magnetic fields if there exist some crossed gradients of density and temperature. It can be seen by rewriting eq.14.46 (and inserting 14.47) as:

$$\frac{\partial \mathbf{B}}{\partial t} = \frac{k_B}{e n_e}\nabla T_e \times \nabla n_e$$ (14.48)

From this last equation it's also clear that Biermann-generated B-fields will be possible only in 2D and 3D geometries. As mentionned above, including Biermann in a MHD model should always be accompanied by the addition of one or several terms in the internal energy equations because of energy conservation considerations. This can be understand by looking at the magnetic energy equation (eq.14.34) and noting that the Biermann electric field is responsible for a source term $-\mathbf{j}\cdot\mathbf{E_{Bier}}\char93 0$. To maintain the conservation of total energy, one can then understand that this source terms must be compensated somewhere in another form of energy sink (where does the energy necessary to generate the B-fields comes from ?). This point is addressed in two different ways in FLASH:

• In the source version, we just aim to keep track of the energy by considering that the increase of magnetic energy due to the Biermann mechanism in a given cell is compensated directly by a decrease of electron internal energy.

• In the flux version, we follow a more physical argument by observing that the presence of the Biermann electric field in the Ohm's law is a direct consequence of keeping track of the electron pressure work on the electron population (which is, in MHD, assumed to be equilibrated by the work of an neutralizing electric field: $E_{Bier}$). Therefore, it seems logical that one needs to take into account the full electron enthalpy transport in this case and this transport happens with the velocity of the electron population $\mathbf{v_e}=\mathbf{v} - \mathbf{j}/en_e$. In conclusion, the presence of the Biermann term in the induction equation should be accompanied by the following electron internal energy equation (without any source, collision or conduction process, for simplicity):

  $$\frac{\partial \epsilon_{e}}{\partial t}+ \nabla(\epsilon_{e} \ma... ...(\mathbf{v}-\mathbf{j}/en_e)+ p_{e}\nabla\cdot (\mathbf{v}-\mathbf{j}/en_e) = 0$$ (14.49)

  
  

  It is then easy to show that with such an equation for electron internal energy and in the presence of $E_{Bier}$ in the Ohm's law, the total energy equation (magnetic + kinetic + internal) can be written in a fully conservative form. As we will demonstrate in 14.3.5.2.5, the addition of these new terms introduces a new kind of wave: the Thermal-magnetic waves and a corresponding phase velocity that needs to be taken into account for stability.

14.3.5.2.3 Formulation for a Cartesian implementation

Using the formalism described in 14.3.5.1.3, we write explicitely the Biermann electric field in cartesian coordinates:

$$\centering \mathbf{E_{Bier}} = \begin{pmatrix}- \partial_x p_e/en_e - \partial_y p_e/en_e - \partial_z p_e/en_e \end{pmatrix}$$ (14.50)

and thus the Biermann flux matrix is explicitly given by:

$$\centering \bar{\bar{\mathbf{F}}}_{Biermann} = \begin{bmatrix}0 &... ...rtial_x p_e/en_e \partial_y p_e/en_e & \partial_x p_e/en_e & 0 \end{bmatrix}$$ (14.51)

14.3.5.2.4 Formulation for a Cylindrical implementation

Assuming the electric field being composed only of the Biermann one, from 14.40 and in the case of an axisymmetric configuration ( $\partial /\partial \theta$ = 0), eq.14.41 can be written as (here $E_{\theta}=0$):

$$\frac{\partial\mathbf{B}}{\partial t} + \frac{\partial}{\partial ... ...race{ \begin{pmatrix} 0 \\ E_r \\ 0 \end{pmatrix}}_\text{z-flux} = \mathbf{0}$$ (14.52)

In this case also we require the user to set the runtime parameter conserveAngField to = .true. for the 2D cylindrical case. It can be seen that in this geometry only a $\theta-$magnetic field can be generated. The Biermann electric field in 2D axisymmetric cylindrical is explicitely given by:

$$\centering \mathbf{E_{Biermann}} = - \begin{pmatrix}\partial_r p_e /en_e 0 \partial_z p_e /en_e \end{pmatrix}$$ (14.53)

14.3.5.2.5 Time step limitation due to the Biermann battery mechanism

As mentionned previously, the flux version of Biermann (see 14.3.5.2.2) allows the propagation of a new kind of waves: the Thermal-magnetic waves (Pert, 1977). To derive the dispersion relation of these waves, suppose a cyldrincal axisymmetric situation where the density depends only on the coordinate r whereas the background temperature is homogeneous. Then consider this set of equations:

$$\frac{\partial \mathbf{B}}{\partial t} = \nabla \times \left(\frac{\nabla p_e}{e n_e} \right)$$ (14.54)

$$\frac{\partial p_e}{\partial t} = \nabla \left( \frac{p_e}{en_e} ... ...f{j}\right) + (\gamma - 1)p_e \nabla\cdot\left( \frac{\mathbf{j}}{en_e} \right)$$ (14.55)

where we assumed $\epsilon_e = p_e/(\gamma - 1)$. Then we apply perturbation theory in cylindrical coordinates (with $T_{e0}=cst$ and $\beta_0(r) = ln(n_{e0}(r))$):

$$\frac{\partial B_{1\theta}}{\partial t} = - \frac{1}{e n_{e0}}\frac{\partial \beta_0}{\partial r}\frac{\partial p_{e1}}{\partial z}$$ (14.56)

$$\frac{\partial p_{e1}}{\partial t} = -(\gamma - 1)\frac{k_B T_{e0... ..._0 e}\frac{\partial \beta_0}{\partial r}\frac{\partial B_{1\theta}}{\partial z}$$ (14.57)

Considering plane waves $\propto exp(i(\omega t - k z))$, we get:

$$\omega B_{1\theta} = \frac{1}{e n_{e0}}\frac{\partial \beta_0}{\partial r} k p_{e1}$$ (14.58)

$$\omega p_{e1} = (\gamma - 1) \frac{k_B T_{e0}}{\mu_0 e}\frac{\partial \beta_0}{\partial r} k B_{1\theta}$$ (14.59)

and we finally get the dispersion relation for the Thermal-magnetic waves:

$$\omega^{2} = (\gamma - 1) \frac{k_B T_{e0}}{\mu_0 e^2 n_{e0}}\left(\frac{\partial \beta_0}{\partial r} \right)^2 k^2$$ (14.60)

from which we can identify the phase velocity $v_{TM}$ of Thermal-magnetic waves:

$$v_{TM} = \sqrt{\frac{\gamma - 1}{\gamma}} \frac{L_{i}}{L_{c,n_{e0}}}c_s$$ (14.61)

where $c_s=\sqrt{\gamma Z k_B T_{e0}/m_i}$ is the ion speed of sound, $L_i$ is the ion inertial length and $L_{c,n_{e0}}=\left\vert \partial \beta_{0} / \partial r \right\vert^{-1}$ is the characteristic length scale of the density gradient. It can clearly be seen that the Thermal-magnetic waves will be come faster than sound waves on scales smaller than the ion inertial length, in a similar way whistler amd Hall drift waves become faster than Alfven waves on this same scale. These non-dispersive waves propagate in inhomogeneous plasmas on iso-density surfaces. In FLASH, the time step related to this new caracteristic speed is computed in Hydro_computeDt.F90 and compared with others speeds in order to get the stable correct global time step. The cfl condition will generally have to be lowered in the presence of the flux-based version of Biermann.

14.3.5.3 Resistive MHD

14.3.5.3.1 How to use it ?

As indicated above, in order to use resistive MHD in FLASH, the user needs to set the runtime parameter useMagneticResistivity to = .true.. Currently, three implementations of resistive MHD are available in FLASH: the "Constant" version, "SpitzerHighZ", as well as the "DaviesWen" one. Depending on which one the user wants to use, the proper path needs to be required in the Config file:

# For Constant resistivity
REQUIRES physics/materialProperties/MagneticResistivity/MagneticResistivityMain/Constant

# For SpitzerHighZ resistivity
REQUIRES physics/materialProperties/MagneticResistivity/MagneticResistivityMain/SpitzerHighZ

# For DaviesWen resistivity
REQUIRES physics/materialProperties/MagneticResistivity/MagneticResistivityMain/DaviesWen

See 22.2 for more details on each coefficients implementation. It is important for the user to know that for all implementions the quantity actually computed and used is the magnetic diffusivity (in units of $l^2 t^{-1}$). Therefore, in the standard FLASH unit system (CGS), the value of the resistivity parameter should be in $cm^2 s^{-1}$. In the case of user-specific implementations of the “MagneticResistivity.F90" and “MagneticResistivity_fullState.F90" files, the values returned to the code should also always be in units of $l^2 t^{-1}$.

There are three parameters that control how resistive MHD is solved in the code: resistivityForm, useCrossMagRes, and resistivitySolver. The resistivityForm parameter accepts values of “parallel", “perpendicular" (default), or “fullaniso", and this controls the form of the resistivity terms in the generalized Ohm's law (eq. 14.31. A value of “parallel" solves $\mathbf{E} = \frac{\eta_\parallel}{\epsilon_0 \omega^2_{pe} \tau_{ei}}\mathbf{J}$, and a value of “perpendicular" solves $\mathbf{E} = \frac{\eta_\perp}{\epsilon_0 \omega^2_{pe} \tau_{ei}}\mathbf{J}$. A value of “fullaniso" will solve $\mathbf{E} = \frac{\eta_\parallel}{\epsilon_0 \omega^2_{pe} \tau_{ei}}\mathbf{J... ...lon_0 \omega^2_{pe} \tau_{ei}} (\mathbf{b} \times \mathbf{J}) \times \mathbf{b}$. Note that this has been rewritten as compared to eq. 14.31. When useCrossMagRes is set to .true., the $\eta_\wedge$ term is included and solved for by adding it to the Hall term.

The resistivitySolver parameter accepts values of “explicit" (default), “implicit", or “hybrid". This controls which type of solver(s) is being used to solve the resistive terms. The explicit solver is a flux-based approach, and requires that the time step be limited for numerical stability. Thus, when using this explict solver, the user should always set useDiffuse and useDiffuseComputeDtMagnetic to .true. (as well as give a reasonable value to dt_diff_factor). See 18.1.6 for details on the implicit solver and how to use it, but note that the implicit solver is currently limited to only work generally on the uniform grid implementation (setup shortcut +ug).

The “hybrid" solver splits the resistive terms and solves $\frac{\eta_\parallel}{\epsilon_0 \omega^2_{pe} \tau_{ei}}\mathbf{J}$ implicitly, while solving $\frac{(\eta_\perp-\eta_\parallel)}{\epsilon_0 \omega^2_{pe} \tau_{ei}} (\mathbf{b} \times \mathbf{J}) \times \mathbf{b}$ explicitly. Hence it only makes sense to use the hybrid solver with resistivityForm = "fullaniso".

14.3.5.4 Thermoelectric transport

14.3.5.4.1 How to use it ?

Currently, three implementations of thermoelectric coefficients are available in FLASH: the "Constant" version, "EppHain", as well as the "DaviesWen" one. More details on these coefficients can be found in 22.4. As indicated above, in order to use thermoelectric transport in FLASH, the user needs to set various runtime parameters to = .true. to turn on different terms. For convenience, the thermoelectric terms in eq. 14.31 are rewritten as:

$$\mathbf{E} = - \frac{k\beta_\parallel}{e} \nabla T_e - \frac{k(\b... ...T_e)\times \mathbf{b} + \frac{k\beta_\wedge}{e}(\nabla T_e \times \mathbf{b}) .$$ (14.62)

These terms are each turned on/off, respectively, with runtime parameters use_Seebeck, use_CrossField, and use_Nernst.

The same runtime parameters also turn on/off analogous terms in the heat flux equation:

$$\mathbf{q} = - \frac{k T_e \beta_\parallel}{e} \mathbf{J} - \frac... ...imes \mathbf{J}) - \frac{k T_e \beta_\wedge}{e}(\mathbf{b} \times \mathbf{J}) .$$ (14.63)

Hence, in the current version of FLASH, the user cannot run a simulation with a term on in the induction equation (eq. 14.62) while leaving it off in the heat flux equation (eq. 14.63).