Meshless MHD

Accurate, Meshless Methods for Magnetohydrodynamics

Philip F. Hopkins & Matthias J. Raives
TAPIR & The Walter Burke Institute for Theoretical Physics, Mailcode 350-17, California Institute of Technology, Pasadena, CA 91125, USA
E-mail:phopkins@caltech.edu
Submitted to MNRAS, April, 2015
Abstract

Recently, we explored new, meshless finite-volume Lagrangian methods for hydrodynamics: the “meshless finite mass” (MFM) and “meshless finite volume” (MFV) methods; these capture advantages of both smoothed-particle hydrodynamics (SPH) and adaptive mesh-refinement (AMR) schemes. We extend these to include ideal magnetohydrodynamics (MHD). The MHD equations are second-order consistent and conservative. We augment these with a divergence-cleaning scheme, which maintains . We implement these in the code GIZMO, together with state-of-the-art SPH MHD. We consider a large test suite, and show that on all problems the new methods are competitive with AMR using constrained transport (CT) to ensure . They correctly capture the growth/structure of the magnetorotational instability (MRI), MHD turbulence, and launching of magnetic jets, in some cases converging more rapidly than state-of-the-art AMR. Compared to SPH, the MFM/MFV methods exhibit convergence at fixed neighbor number, sharp shock-capturing, and dramatically reduced noise, divergence errors, & diffusion. Still, “modern” SPH can handle most test problems, at the cost of larger kernels and “by hand” adjustment of artificial diffusion. Compared to non-moving meshes, the new methods exhibit enhanced “grid noise” but reduced advection errors and diffusion, easily include self-gravity, and feature velocity-independent errors and superior angular momentum conservation. They converge more slowly on some problems (smooth, slow-moving flows), but more rapidly on others (involving advection/rotation). In all cases, we show divergence-control beyond the Powell 8-wave approach is necessary, or all methods can converge to unphysical answers even at high resolution.

keywords:
methods: numerical — hydrodynamics — instabilities — turbulence — cosmology: theory

1 Introduction: The Challenge of Existing Numerical Methods

Magnetic fields are an essential component in astrophysical hydrodynamics, and for many astrophysical problems can be reasonably approximated by ideal (infinite conductivity) magnetohydrodynamics (MHD). The MHD equations are inherently non-linear, however, and most problems require numerical simulations. But this poses unique challenges, especially for numerical methods which are Lagrangian (i.e. the mesh elements follow the fluid), rather than Eulerian (solved on a fixed grid).

In most discretizations (although see Kawai, 2013), evolving the MHD equations in time will lead to a violation of the “divergence constraint” (the requirement that ). Unfortunately, this cannot simply be ignored to treated as a “standard” numerical error term which should converge away with increasing resolution, because certain errors introduced by a non-zero are numerically unstable: they will eventually destroy the correct solution (even at infinite resolution) and/or produce unphysical results (e.g. negative pressures). Arguably the most elegant solution is the so-called “constrained transport” (CT) method of Evans & Hawley (1988), which maintains to machine precision; however, while there is no obvious barrier in principle to implementing this in meshless and unstructured mesh methods (see recent developments by Mocz et al., 2014a), it has thus far only been practical to implement for real problems in regular, Cartesian grid (or adaptive-mesh refinement; AMR) codes. But for many problems in astrophysics, Lagrangian, mesh-free codes have other advantages: they minimize numerical diffusion and over-mixing, move with the fluid so automatically provide enhanced resolution with the mass (in a continuous manner, which avoids low-order errors necessarily introduced by AMR refinement boundaries), couple simply and accurately to cosmological expansion and -body gravity codes, easily handle high Mach numbers, conserve angular momentum and naturally handle orbiting disks without prior knowledge of the disk geometry, avoid “grid alignment” and carbuncle instabilities (where the grid imprints preferred directions on the gas), and feature errors which are independent of the fluid bulk velocity (so can converge more rapidly when the fluid moves).

A variety of approaches have been developed to deal with these errors. The simplest commonly-used method, the so-called “Powell 8-wave cleaning,” simply subtracts the unstable error terms resulting from a non-zero from the equation of motion (Powell et al., 1999). This removes the more catastrophic numerical instabilities, but does not solve the convergence problem – many studies have shown that certain types of problems, treated with only this method, will simply converge to the wrong solution (Tóth, 2000; Mignone & Tzeferacos, 2010; Mocz et al., 2014a). And the subtraction necessarily violates momentum conservation, so one would ideally like the subtracted terms (the values) to remain as small as possible. Therefore more sophisticated “cleaning” schemes have been developed, the most popular of which have been variants of the Dedner et al. (2002) method: this adds source terms which transport the divergence away (in waves) and then damp it. This has proven remarkably robust and stable.

However, applications of these techniques in Lagrangian codes in astrophysics have remained limited. The most popular Lagrangian method, smoothed-particle hydrodynamics (SPH), suffers from several well-known errors that make MHD uniquely challenging. The SPH equations are not consistent at any order (meaning they contain zeroth-order errors; Morris 1996; Dilts 1999; Read et al. 2010); this introduces errors which converges away very slowly and causes particular problems for divergence-cleaning. Also, naive implementations of the equations are vulnerable to the tensile and particle pairing instabilities. And artificial diffusion terms, with ad-hoc parameters, are required in SPH to deal with discontinuous fluid quantities. As such, many previous implementations of MHD in SPH were unable to reproduce non-trivial field configurations, were extremely diffusive, or were simply unable to numerically converge; in turn key qualitative phenomena such as the magneto-rotational instability (MRI) and launching of magnetic jets could not be treated (see Swegle et al., 1995; Monaghan, 2000; Børve et al., 2001; Maron & Howes, 2003; Price & Rosswog, 2006; Rosswog & Price, 2007; Price & Bate, 2008; Dolag & Stasyszyn, 2009).

Recently, however, a number of breakthroughs have been made in Lagrangian hydrodynamics, with the popularization of moving-mesh and mesh-free finite-volume Godunov methods. Springel (2010); Duffell & MacFadyen (2011); Gaburov et al. (2012) have developed moving-mesh MHD codes, which capture many of the advantages of both AMR and SPH, using the Dedner et al. (2002) cleaning method. Meanwhile, Lanson & Vila (2008a, b); Gaburov & Nitadori (2011); Hopkins (2015) have developed a class of new, mesh-free finite volume methods which are both high-order consistent (convergent) and fully conservative. These are very similar to moving-mesh codes (in fact, Voronoi moving-meshes are technically a special case of the method). In Hopkins 2015, these are developed for hydrodynamics in the multi-method, hydrodynamics+gravity+cosmology code GIZMO, which is an extension of the -body gravity and domain decomposition algorithms from GADGET-3 (Springel, 2005) to include a variety of new hydrodynamic methods.111A public version of this code, including the full MHD implementation used in this paper, is available at http://www.tapir.caltech.edu/~phopkins/Site/GIZMO.html. Users are encouraged to modify and extend the capabilities of this code; the development version of the code is available upon request from the author. In Hopkins (2015), a broad range of test problems are considered, and it is shown that these also capture most of the advantages of AMR and SPH, while avoiding many of their disadvantages. Particularly important, these eliminate the low-order errors, numerical instabilities, and artificial diffusion terms which have plagued SPH. Gaburov & Nitadori (2011) considered a range of MHD test problems and found very encouraging preliminary results; they showed that these mesh-free methods could handle complicated non-linear problems like the MRI with accuracy comparable to state-of-the-art grid codes. Meanwhile, tremendous improvements have also been made in SPH (Ritchie & Thomas, 2001; Price, 2008; Wadsley et al., 2008; Cullen & Dehnen, 2010; Read & Hayfield, 2012; Saitoh & Makino, 2013; Hopkins, 2013; Tricco & Price, 2012, 2013; Tricco, 2015).

Therefore, in this paper, we extend the mesh-free MFM and MFV Lagrangian hydrodynamics in GIZMO to include MHD, and consider a systematic survey of a wide range of test problems, and compare state-of-the-art grid-based (AMR) codes, MFM, MFV, and SPH MHD methods, implemented in the same code. This includes problems such as the MRI and MHD jets which have been historically challenging. We show in all cases that the new meshfree methods exhibit good convergence and stable behavior, and are able to capture all of the important behaviors, even at low resolution. On some problem classes, we show they converge faster than state-of-the-art AMR codes using CT.

2 Numerical Methodology

2.1 Review of the New Meshless Methods

Paper I derives and describes the pure-hydrodynamic version of the numerical methods here in detail, including self-gravity and cosmological integration. This is almost entirely identical in MHD; therefore we will not repeat it. However we will very briefly review the new numerical methods.

The equations we solve are the standard finite-volume Godunov-type equations: the fundamental difference between our meshless methods and a moving-mesh is simply that the definition of the volume partition (how the volume is divided among different mesh-generating points or “cells/particles”) is distinct. The further difference between this and a fixed-grid code is, of course, that the mesh-generating points/cells move, and that their arrangement can be irregular (as opposed to a Cartesian grid).

In a frame moving with velocity , the homogeneous Euler equations in ideal MHD (and pure hydrodynamics, which forms the special case ) can be written as a set of hyperbolic partial differential conservation equations of the form

(1)

where refers to the inner product between the gradient operator and tensor , is the outer product, is the “state vector” of conserved (in the absence of sources) variables, the tensor is the flux of conserved variables, and is the vector of source terms

(2)

where is mass density, is the total specific energy ( the internal energy), is the sum of thermal and magnetic pressures, and is a scalar field defined below.

The meshless equations of motion are derived in Paper I in standard Galerkin fashion beginning from the integral form of the conservation laws, after multiplying Eq. 1 by a test function

(3)

Here the domain is such that , where is the number of spatial dimensions, is the co-moving derivative of any function , and is the normal vector to the surface , and the test function is an arbitrary differentiable Lagrangian function.

To transform this into a discrete set of equations, must chose how to partition the volume (for the “averaging/integration” step). If we choose a uniform Cartesian grid between uniformly spaced points, then we will recover the standard Godunov-type finite-volume grid-based equations of motion (like that in ATHENA and many popular AMR codes). If we choose a Voronoi tesselation between moving mesh-generating points, we recover a moving-mesh method similar to AREPO. For the new methods in Paper I, we partition the volume according to a continuous weighting function , such that the fraction of the differential volume at the point associated with the mesh-generating point at is given by

(4)

Where is any kernel/weight function and is a “kernel length.” Note that this guarantees a “partition of unity” (the volume is perfectly divided into cell-volumes ), and leads to a Voronoi-like partition, but with slightly smoothed boundaries between cells (which leads to some advantages, and some disadvantages, compared to moving meshes, where the boundaries are strict step functions). In the limit where goes to a delta function, the method becomes exactly a Voronoi-type moving-mesh method.222Here and throughout this paper, we will define the kernel size at the location of cell/particle as the effective cell side-length, based on the cell volume . In 1D/2D/3D, this is , , , respectively. The volume is calculated directly from the neighbor positions defining the volume partition (see Paper I). This exactly reproduces the grid spacing if the particles are arranged in a Cartesian grid. Note that, in principle, can be non-zero at . In MFM/MFV methods this has nothing to do with the effective cell/particle volume/size (conserved quantities are not “smoothed” over the kernel, but only averaged inside the a single cell of volume just like in a grid code), but instead reflects the size of the stencil (number of neighbor cells) between which fluxes are computed. As in grid codes, increasing the stencil size can increase diffusion, but does not directly alter the resolution scale.

This choice is combined with a second-order accurate moving-least squares matrix-based gradient operator, which has been utilized in many other methods (including grid-based codes; see Maron et al. 2012; Tiwari & Kuhnert 2003; Liu et al. 2005; Luo et al. 2008; Lanson & Vila 2008a, b; Mocz et al. 2014b). Eq. 3 can then be expanded and analytically integrated to yield a second-order accurate set of discrete evolution equations:

(5)

This is identical to the standard Godunov-type finite-volume equations. The term is simply the cell-volume integrated value of the conserved quantity to be carried with cell/particle (e.g. the total mass , momentum, or energy associated with the cell ); its time rate of change is given by the sum of the fluxes into/out of an “effective face area” , plus the volume-integrated source terms. The full mathematical derivation and expression for the is given in Paper I (§ 2.1).

We then use a standard MUSCL-Hancock type scheme for finite-volume Godunov methods to solve Eq. 5. This is commonly used in grid and moving-mesh codes (van Leer, 1984; Toro, 1997; Teyssier, 2002; Fromang et al., 2006; Mignone et al., 2007; Cunningham et al., 2009; Springel, 2010); it involves a slope-limited linear reconstruction of face-centered quantities from each mesh generating point (cell “location”), with a first-order drift/predict step for evolution over half a timestep, and then the application of a Riemann solver to estimate the time-averaged inter-cell fluxes for the timestep. See Paper I(§ 2 and Appendices A & B) for details. The points then move with the center-of-mass gas velocity.

In Paper I, we derive two variants of this method and implement them in GIZMO. First, the meshless finite-volume (MFV) method. This solves the Riemann problem between cells assuming the effective faces move with the mean cell velocity; this is analogous to a moving-mesh code, and includes mass fluxes between cells. Second, the meshless finite-mass (MFM) method. This solves the Riemann problem assuming the face deforms in a fully Lagrangian fashion; in this case there are no mass fluxes. The two are formally identical up to a difference in the non-linear (second-order) error terms in the fluxes, provided the cells move with the gas velocity. In practice, each has some advantages and disadvantages, discussed below.

2.2 Code Modifications for MHD

Everything described above is identical in hydodynamics and MHD; and all details of the code (except those specifically described below) are unchanged from Paper I.

As usual for finite-volume Godunov schemes, we explicitly evolve the conservative variables (integrated magnetic field over the volume partition corresponding to a mesh-generating point) and ; primitive variables and gradients are then constructed from these (e.g.  as in Paper I).

2.2.1 The Riemann Solver

As in Paper I, we solve the 1D, un-split Riemann problem in the rest-frame of the effective face between the cells. However we require a Riemann solver that allows . Since in the hydro case we use the HLLC solver, here we adopt the widely-used HLLD solver (Miyoshi & Kusano, 2005). This is accurate at the order required, and extremely well-tested (see e.g. Miyoshi & Kusano 2005; Stone et al. 2008 and modern versions of e.g. RAMSES and ENZO).

The frame motion is calculated for both the MFM and MFV methods as in Paper I. We emphasize that for our MFM method, it is required that we use a solver which explicitly includes the contact wave (i.e. the contact discontinuity, on either side of which mass is conserved). This is because in MFM we always solve the Riemann problem with frame moving exactly with the contact wave; attempting to simply find a frame where the mass flux vanishes in a simpler (such as the HLL or Rusanov) approximation leads to incorrect, and in many cases unphysical solutions.

2.2.2 Signal Velocities and Time-Stepping

As in Paper I, we limit timesteps with a standard local Courant-Fridrisch-Levy (CFL) timestep criterion; in MHD we must replace the sound speed with the fast magnetosonic speed:

(6)
(7)
(8)

Here is the separation between two points, is the unit vector , is the sound speed, is the Alfven speed, is the effective cell size defined above, refers to the maximum over all interacting neighbors of , and is the signal velocity (Whitehurst, 1995; Monaghan, 1997).

2.2.3 Divergence-Cleaning

Ideally, this would complete the method; but the above method cannot ensure the divergence constraint. To do this, we must add the following source terms:

(9)
(10)

The first term () represent the Powell et al. (1999) or “8-wave” cleaning, and subtracts the numerically unstable terms from non-zero . This is necessary to ensure numerical stability and Galilean invariance – most problems will crash or converge to incorrect solutions without this. The second term () follows the method of Dedner et al. (2002), who introduce a conservative scalar field which transports divergence away from the source and damps it. This is necessary to keep low, minimizing the resulting errors.

Following Dedner et al. (2002) and Gaburov & Nitadori (2011), it is straightforward to show that this leads to the following form for the discrete terms in our equation of motion (Eq. 5):333Eq. 9 is the continuum equation; we stress that we are not free to choose how we discretize it. To actually ensure numerical stability, the form of the terms must exactly match those terms from our Riemann solver solution which are unstable (e.g. the tensile terms). Likewise, the divergence cleaning must act specifically to reduce defined in the same manner, or it does not serve any useful purpose. It might be tempting, for example, to use the value of calculated from our second-order accurate matrix gradient estimator for the Powell terms, or to construct a pair-wise symmetrized version of Eq. 11 which manifestly maintains momentum conservation. However, these will not actually eliminate the unstable terms, and we have confirmed that they lead to catastrophic errors in our tests.

(11)
(12)
(13)

where the and terms are defined below.

We have some freedom to choose and . Following Tricco & Price (2012), we take , where is the maximum signal velocity as described above ( is simply a convenience parameter which re-normalizes the characteristic speed). The value is close to the fast magnetosonic wave speed, but also accounts for super-sonic cell approach velocities, which is critical for good behavior in highly supersonic compressions. We have experimented with variations in the dimensionless parameter , and find the best results for (values produce ineffective cleaning, values lead to numerical instability). We take (where is the effective cell length defined above).444In timestepping, we update for the term with the implicit solution ; this allows us to take larger timesteps without numerical instability. For (the damping speed), we have considered several choices, all of which give very similar results. These are detailed in Appendix D; our default choice is Eq. 50, which is closely related to the local fastest possible signal velocity. We have also experimented extensively with , and find a best compromise between stability and diffusivity for values ; we adopt as our default in all problems here.

Note that, unlike the original Dedner et al. (2002) formulation (and implementations in codes like AREPO and PLUTO), this means and are spatially variable. This allows us to maintain hierarchical timestepping (while forcing a constant imposes a global maximum timestep, a severe CPU cost penalty), and has many other advantages (see Appendix B). But it is not a priori obvious that this will maintain stability. However, in Appendix B we discuss this in detail and derive a rigorous stability criterion, which should be satisfied by our choices above, designed so that and are locally smooth (on the kernel scale). We confirm this stability in our numerical tests.

We caution that the discrete source terms, particularly the Powell terms which subtract -centered quantities, are not manifestly antisymmetric between cell pairs . This means that momentum and energy conservation in MHD are only accurate up to integration accuracy, times a term proportional to (unlike in hydrodynamics, where conservation can be ensured at machine accuracy).555Some of the Dedner terms can be made anti-symmetric without destroying the numerical stability of the scheme; for example the -flux correction described below, and the term (by using a single wavespeed for the whole problem). However in all our tests the conservation errors from those terms are always sub-dominant to the error from the (inescapable) Powell term (and all these errors vanish when ). So we find the best overall conservation properties result from using the most accurate possible cleaning scheme, rather than a partially-conservative (but less accurate) cleaning. Controlling is critical to minimize these errors.

These source terms also modify the Riemann problem. When we perform the reconstruction to obtain the left and right states (-side) and (-side), we can define a convenient coordinate system where (i.e. the -axis is normal to the effective face between cells and ). In this coordinate system, the normal-component of the -fields , will in general not be equal. But equal values (i.e. non-zero in the 1D problem) are required for a physical solution. Without divergence-cleaning (Powell-only), this is handled by simply replacing and with the mean value . Dedner et al. (2002) showed that with the source terms of Eq. 9, the infinitely-sharp discontinuity leads to a physical solution , , in infinitesimally small time:

(14)
(15)
(16)
(17)

here is the fastest wave speed in the local 1D problem, which can be computed only from the and values (it does not necessarily correspond to in Eq. 9 above).666We have also explored an alternative, two-wavespeed formulation of the and terms, discussed in Appendix F. For all tests here, the difference is small. This is separable from the full Riemann solution. So, in the Riemann problem, we first update and according to the above, then compute the full Riemann solution using the updated values (and usual and ). The flux of is then , where is the normal face velocity in the boosted frame (in which we solve the Riemann problem). Because is advected with the fluid, we follow Gaburov & Nitadori (2011) and simply take the flux to be for , and for , where is the mass flux (so this vanishes for our MFM method).777We have experimented with using for the -flux. However, this yields no improvement on any test problem, and the term from can introduce numerical instability under some circumstances. Therefore we use the simpler -flux.

The HLLD solver requires an initial guess for the left and right wavespeeds, to compute a solution. If we define as the normal-component of the reconstructed velocities, then we use , where is the updated fast magnetosonic wavespeed using the updated normal :

(18)

and .888The HLLD solver can fail in some very rare circumstances if “bad” guesses are used. We therefore check whether our first estimate produces a solution where the pressure is everywhere positive. If this fails, then we instead compute the Roe-averaged velocity and fast magnetosonic speed , and try and . We then check again; if this fails (which does not occur in any test problem here) we test a Lax-Friedrich estimate (, with from our first guess). If this somehow fails still, we go back and re-compute the interface using a piecewise-constant (first-order) approximation, then check the series of wavespeeds again. If this fails, the code exits with an error (this only occurs when unphysical values are input into the solver).

Finally, as in Paper I, we solve the Riemann problem in the boosted frame corresponding to the mean motion of the quadrature point between mesh-generating points. We must therefore boost back to the simulation frame. The de-boosted fluxes for cells follow Paper I, for the hydro terms, but with the additional terms for (for the fluxes to cell-):

(19)
(20)

The second equation just accounts for the energy flux associated with the corrected -flux.

Figure 1: Linear magnetosonic wave test problem (see § 3.1). Here a traveling, one-dimensional fast magnetosonic wave is propagated one wavelength; we then define the norm as the mean absolute error relative to the known analytic solution in density (top) or magnetic field (bottom; error here is equivalent to the numerical errors). We compare our new, meshless Lagrangian finite-volume Godunov methods (“meshless finite-mass” or MFM, and “meshless finite-volume” or MFV) from Paper I (see § 2.1), to the best current implementation of SPH MHD (see § 2.3), and to state of the art grid codes, here ATHENA run as a third-order PPM code using constrained transport (CT) to ensure to machine precision. Dotted line shows second-order convergence (); MFM/MFV and grid/PPM methods converge at this rate, as expected. Convergence is also good () in MFM/MFV for the divergence errors (, so these are fractionally very small). SPH shows some (slower) convergence until its known zeroth-order errors dominate; then errors flatten with resolution. This is reduced in SPH by increasing the kernel size. “SPH-lo” (standard ) uses the equivalent of in 3D (our default choice in all MFM/MFV runs). “SPH-hi” uses a 3D-equivalent .

2.3 The SPH MHD Implementation in GIZMO

As described in Paper I, GIZMO is a multi-method code: users can run with the MFM or MFV hydrodynamic methods, or SPH, if desired. We therefore update our SPH implementation to include MHD. The exact SPH equations are given in Appendix A.

Briefly, the non-magnetic implementation of SPH follows the “modern” P-SPH method developed in Hopkins (2013) and extended in Paper I. This includes state-of-art re-formulations of the SPH hydrodynamics equations to eliminate the known “surface tension” errors (Saitoh & Makino, 2013; Hopkins, 2013), Lagrangian-derived terms to account for variable smoothing lengths (Springel & Hernquist, 2002), addition of artificial diffusion terms for thermal energy (Price, 2008; Wadsley et al., 2008), higher-order switches for artificial diffusion terms to minimize unnecessary dissipation (Cullen & Dehnen, 2010), the use of higher-order kernel functions to allow larger SPH neighbor numbers and reduce the zeroth-order SPH errors (Dehnen & Aly, 2012), switches to prevent disparate time-stepping between neighbor particles (Saitoh & Makino, 2009), introduction of more accurate matrix-based gradient estimators (García-Senz et al., 2012), and conservative, more accurate coupling of SPH to gravity (Price & Monaghan, 2007; Barnes, 2012).

The SPH MHD implementation combines these improvements to SPH with the MHD algorithms from the series of papers by Tricco & Price (2012, 2013); Tricco (2015). These introduce artificial diffusion for magnetic fields (artificial resistivity), with a similar “switch” to reduce unnecessary dissipation, and re-discretize the MHD equations from the particle Lagrangian included the Dedner et al. (2002) and Powell et al. (1999) terms, so that the divergence cleaning actually acts on the tensile-unstable terms and these terms are properly subtracted (unlike most previous SPH-MHD implementations). We make some further improvements: following Price et al. (2012); Bate et al. (2014) and directly evolving the conserved quantities and , and slightly modifying the artificial resistivity terms to allow the method to capture cosmological field growth and MHD fluid-mixing instabilities.

3 Test Problems

We now consider a series of test problems. To make comparison as fair as possible, we will consider MFM/MFV/SPH implementations in the same code (GIZMO). Unless otherwise stated, we compare to fixed-grid results from ATHENA (Stone et al., 2008); this is representative of the state-of-the-art in finite-volume, non-moving mesh codes. Because we are interested in methods which can be applied to complicated, multi-physics systems, we use the same numerical implementation and identical values of purely numerical parameters (e.g.  and ), within each method, for all problems. It is of course possible to improve performance in any simple test problem by customizing/tweaking the method, but this usually entails a (sometimes serious) loss of accuracy on other problems. The implementations used here are therefore our attempts at a “best compromise” across all problems.

3.1 Linear Magnetosonic Waves: Testing Convergence

We begin by considering a simple linear one-dimensional magnetosonic wave.999See http://www.astro.princeton.edu/~jstone/Athena/tests/linear-waves/linear-waves.html The problem is trivial, but since virtually all numerical schemes become first-order at discontinuities, smooth linear problems with known analytic solutions are necessary to measure formal convergence. Following Stone et al. (2008), we initialize a unit-length (in the -direction), periodic domain, with polytropic gas, density , pressure , and magnetic field . We add to this a traveling fast magnetosonic wave101010For the grid and MFM/MFV methods, the perturbation is initialized by keeping the particles equidistant and modifying the conserved quantities. In SPH however, significantly better performance is obtained by initializing the density perturbation by using exactly equal particle masses with perturbed locations (this allows apparent convergence in SPH to extend further). We therefore show this case for SPH (as it corresponds to the more likely case in real problems). We note that this initial condition also improves our MFM/MFV performance, but by a smaller amount. The SPH performance is also substantially improved if we “turn off” artificial viscosity (or set its minimum to zero), but this is numerically unstable and produces catastrophic errors in most of our tests (see e.g. Hopkins, 2013, 2015; Hu et al., 2014; Rosswog, 2014). with amplitude , and allow the wave to propagate one wavelength; we then define the L1 error norm for the density (or any other variable) . We consider .

All methods we consider are able to evolve the wave. Fig. 1 plots the L1 norm for density (the velocity variables look similar), and . In , we see that both our MFM and MFV methods converge – as expected – with second-order accuracy. We compare to a state-of-the-art grid code, here ATHENA, run in the most accurate possible mode: PPM (formally a third-order reconstruction method), with the CTU integrator, and CT used to ensure . Despite the higher order of ATHENA and the fact that it uses CT to ensure errors remain at the machine-error level, the errors are nearly identical to our MFM/MFV results. The norm for directly measures the divergence errors; since we do not use CT, these are non-zero. However they are (1) very small, and (2) converge away appropriately – in fact, we see super-convergence () for our MFM/MFV methods. As in Paper I, the convergence rate in all variables is independent of the kernel neighbor number in MFM/MFV: the choice only controls the normalization of the errors (larger neighbor numbers reduce noise, but increase diffusion). Based on our experiments in Paper I, we find roughly optimal results using a 3D-equivalent neighbor number ( in 1D).

For SPH, we see slower convergence in , but it is reasonable at low resolution; but at high-resolution, the SPH zeroth-order errors (here, the “kernel bias” error in density estimation, and the systematic gradient error which results from imperfect particle order in the kernel) begin to dominate, and the errors flatten with resolution. This is more severe if we use a lower neighbor number; going to higher-order kernels and higher neighbor numbers suppresses the errors, although they still eventually appear. For true convergence, must increase with , as is well-known (see Zhu et al., 2014, and references therein). Here we compare 3D-equivalent (our default MFM/MFV choice used in all runs in this paper), henceforth referred to as “SPH-lo,” to 3D (henceforth “SPH-hi”).111111Increasing the neighbor number with particle number must be done carefully in SPH, as discussed in e.g. Price (2012); Dehnen & Aly (2012); Zhu et al. (2014), since one must avoid the pairing instability and also account for the fact that simply increasing neighbor number in SPH changes the resolution length. For SPH-lo we use the popular Schoenberg (1946) cubic spline kernel, so in 3D corresponds to an “effective” resolution scale of in units of the mean inter-particle spacing (following Dehnen & Aly 2012, , where is the standard deviation of the kernel). For SPH-hi we use the Schoenberg (1946) quintic spline, so in 3D corresponds to . Therefore we caution that the “effective resolution” of SPH-hi is slightly larger () than SPH-lo at fixed ; however the difference is much smaller than might naively be expected based on alone.

Figure 2: Brio-Wu shocktube (§ 3.2), at time ; we compare the exact solution to that computed at finite resolution (plotted region contains elements across the -direction) with different methods. High-order grid methods have converged well at this resolution, except for post-shock ringing in . MFM/MFV methods also show good convergence; but at this resolution, MFM still shows some small “overshoot” in the jumps at (more sensitive to our slope-limiter than the method itself), and both show some small (percent-level) errors in owing to the errors; however the fractional magnitude of is controlled well by our cleaning scheme (typical errors at this resolution; still below at jumps). Discontinuities are well-captured across cells/particles in the linear direction. SPH (with high ) captures all the key features, but at this resolution shows larger noise, -level overshoots in at rarefactions (), some suppression of internal energy around (owing to more smeared-out dissipation from divergence-cleaning), and significantly larger errors (reaching ); these do converge away but more slowly.
Figure 3: Brio-Wu shocktube, as Fig. 2, for additional methods. If we consider MFM/MFV using only the Powell et al. (1999) source terms to stabilize MHD, but no Dedner et al. (2002) divergence-control (as has often been done in the literature), we obtain incorrect shock jumps; most noticeably in . This error does not converge away. This is despite the fact that the formal errors are still small; the key is the terms in the Dedner et al. (2002) scheme that enter the Riemann problem and act specifically at discontinuities. We also compare SPH run with the same (lower) neighbor number as our MFM/MFV methods; here the noise is larger (as expected).
Figure 4: The Toth super-sonic shocktube (§ 3.3). The resolution is similar to the Brio-Wu test. MFM and grid/CT methods converge most rapidly to the exact solution (dotted), followed by MFV, which exhibits some residual noise in around at this resolution. MFM/MFV both control errors well; the large discontinuity introduces a small (percents-level) offset in which converges away (nearly ideal), this is especially challenging for divergence-cleaning methods. SPH shows larger errors, and noise; with some systematic offset in the jump around (as divergence cleaning is spread over several smoothing lengths), which converges slowly. The noise is reduced with larger neighbor number, but not eliminated. In all methods, the Powell scheme alone leads to systematically incorrect shock jumps in and ; as in Fig. 2, these do not converge.

3.2 Brio-Wu Shocktube: Capturing MHD Discontinuities & Controlling Noise

Next we consider the Brio & Wu (1988) shocktube; this tests whether the code can accurately represent uniquely MHD shocks, rarefactions, and contact discontinuities. We initialize a 2D periodic box (size , , with cells/particles) with left-state , and right-state , with . Figs. 2-3 compare the results at time . Note that our box is intentionally large and extends well beyond the “active” domain; the dynamically active region is only elements across.121212We have run this test with a number of different element configurations in our initial conditions, including: square, triangular, and hexagonal lattices, SPH and gravitational glasses, and random (Poisson) positions. We have also considered the case of equal particle masses (different particle spacing across the initial discontinuity) or unequal particle masses (same spacing). As expected from the derivation in Gaburov & Nitadori (2011), our MFM/MFV methods are only very weakly sensitive to these choices (slightly more noise appears in the “worst case” situation of random positions). Perhaps more surprising, our SPH results are also only weakly sensitive to this. We show results for a relaxed SPH glass but they are qualitatively identical to any of the other configurations. And we show explicitly in Paper I that the choice of equal/different particle masses makes a difference only in the magnitude of the errors just at the contact discontinuity.

At this resolution, the CT-based grid code is well-converged, except for some post-shock “ringing” most visible in . In MFM & MFV methods, the agreement with the exact solution is good at this resolution. At higher-resolution, we find that the MFM/MFV results are nearly indistinguishable from the exact solution. At lower resolution, there is some “overshoot” at the density discontinuity, with MFM – this is discussed in detail in Paper I; it is mostly sensitive to the choice of slope-limiter (not the basic numerical method). The effects are much smaller in MFV, owing to mass fluxes allowing more sharply-captured density discontinuity. This also causes a small pressure “blip” in MFM at the contact discontinuity, but this converges away. Shock jumps and discontinuities are captured across cells/particles in each direction; comparable to high-order grid methods.

As an indicator of the errors, we plot the dimensionless magnitude of ; our divergence cleaning keeps this generally at low values (), except at the magnetic shocks (the large discontinuities in ); but even there, the maximum value is still . The good divergence-cleaning is also manifest in ; at the same shocks, there are some jumps in generated, but this is returned to the correct value across particles, and the magnitude of the deviations from the analytic is typically at the sub-percent level, at this resolution.

With no divergence corrections whatsoever, the problem crashes. However, Fig. 3 shows that if we run with only the Powell terms in Eq. 9 (i.e. do not include the Dedner et al. (2002) divergence-cleaning and damping terms), this particular problem is not badly corrupted in most respects. As expected the divergence and deviation in are larger. More seriously, though, a systematically incorrect jump in appears, which does not converge away without divergence cleaning, even at higher resolution in the -direction. This problem occurs in MFM, MFV, and SPH.

SPH captures all the qualitative features. However, if we run with the same neighbor number (kernel size) as in our MFM/MFV methods, the noise is larger. This is shown in Fig. 3. If we instead use the equivalent of a 3D neighbor number of (as opposed to the we use for MFM/MFV), the noise is reduced, as expected. However, there is still much larger noise and post-shock ringing (compared to MFM/MFV), and significant overshoot in the velocities and at the rarefactions.131313In extensive experiments, we found that the magnitude of these errors in SPH, at fixed resolution, is determined by the artificial viscosity and resistivity schemes. If we simply assume a constant (large) artificial viscosity (i.e. disable our normal “shock detection” switch) we are able to eliminate most of the noise, and reduce the overshoot, in SPH in Figs. 2-3. However, such a choice would severely degrade the performance of our SPH implementation on almost every other problem we consider. As an attempted compromise, in Appendix A, we discuss modifications to the default Cullen & Dehnen (2010) viscosity scheme and Tricco & Price (2013) resistivity scheme, which we have used here. Using instead the default version of the viscosity scheme makes only small differences for large neighbor number (“SPH-hi”), but produces much larger (order-unity) noise levels at low-neighbor number (“SPH-lo”) on this test. In highly super-sonic tests, the difference is negligible. Divergence-cleaning works in SPH, but is much less effective, especially around the contact discontinuity, where the divergence errors reach . We stress though that these errors are resolution-dependent, and do converge away eventually. Better results (at a given resolution) can also be obtained on shocktube problems in SPH by increasing the strength of the artificial dissipation terms (viscosity/resistivity/conductivity); however this significantly degrades performance on other tests we consider. For an example of this test which demonstrates good agreement between SPH and the exact solution (by combining higher resolution and stronger dissipation), see Tricco & Price (2013).

We have also compared the Ryu & Jones (1995) MHD shocktube, which exhibits all seven MHD waves simultaneously. The qualitative results and differences between methods are the same as in the Brio-Wu test, so we do not show it here.

3.3 Toth Shocktube: The Critical Need for Divergence Cleaning Beyond Powell

Next we consider the Tóth (2000) shocktube; this tests super-sonic MHD shocks. It is particularly important because Mignone & Tzeferacos (2010) showed that a hyperbolic divergence-cleaning scheme such as the Dedner et al. (2002) method (or CT) is necessary to get the correct shock jump conditions at any resolution, in grid-based methods. We initialize a 2D periodic box with the same initial element configuration as the Brio & Wu (1988),141414Again, we have considered a variety of initial element configurations. In the Tóth (2000) shocktube we find these produce negligible differences. with left-state , and right-state and .

Fig. 4 compares the results at time (for clarity, we plot only a randomly-chosen subset of cells/particles for each method). Our MFM method does extremely well at this resolution; the only difference between it and high-resolution grid runs is a small ( level) deviation in introduced by errors at this resolution, and some noise at the shock in the small . The errors are generally extremely small, at super-sonic shocks and elsewhere. MFV is similar, although there is substantially more noise at the same resolution in and in the post-shock region; this is seen in pure hydro in Paper I and in Gaburov & Nitadori (2011); the larger noise translates to slightly larger errors at shocks, hence slower convergence in . Both MFM & MFV are indistinguishable from the exact solution at larger resolution in the -direction.

However, with the Powell-only (no divergence-cleaning) mode, we see that the shock is in the wrong place! The shock jump is systematically wrong in , and , and this leads to it being in the wrong place over time. Again, this appears even at infinite resolution. The noise, especially in , is also much larger, and errors, as expected, are factor larger.

In SPH, the noise is much larger, even with much larger 3D-equivalent . This is especially noticeable in , , and . In , the shock jump is also systematically over-estimated by , owing to the much larger () errors at the shock jump. As discussed below (§ 5), in SPH, divergence-cleaning cannot act over smaller scales than a few smoothing lengths, so the method has difficulty controlling the errors seen in the Powell-only case. The divergence errors are systematically larger by factors . With as in our MFM/MFV methods, the noise is yet larger (order-unity fractional noise in , ).

Figure 5: Field loop advection test (§ 3.4). An equilibrium, 2D field loop is advected uniformly across the domain; we plot (in units of the initial loop value, as labeled) at time , for tests with resolution. Here we use in initially Cartesian particle lattice with unequal-mass particles (the density within the loop is a factor higher than the background). This is the natural configuration for fixed-mesh codes but represents a near “worst case” for Lagrangian codes. In all cases, the initial conditions (top left) should be reproduced identically. In non-moving grid methods, even at arbitrarily high order, advection errors diffuse the loop, while numerical resistivity reduces the central field strength. In MFM/MFV/SPH, the advection errors are eliminated, and numerical resistivity at the center is reduced; however divergence-cleaning and “grid noise” from different particle masses around the loop edge produce some diffusion and noise (the peak value of in any cell at any time remains , however). The Powell-only scheme exhibits much more severe noise, because the errors are transported but not damped; this leads to non-linear corruption of the solution.
Figure 6: Quantitative decay of the box-averaged magnetic energy in the field loop test (Fig. 5), owing to numerical diffusion/resistivity. At infinite resolution, methods preserve the initial . Because they are Lagrangian, MFM/MFV/SPH methods show much less dissipation than high-order grid methods at the same resolution (default runs here are , but we compare MFM/grid runs for reference). SPH shows some spurious initial growth of as divergence-cleaning acts on zeroth-order kernel errors, but the subsequent decay rate is close to MFM/MFV. Otherwise on this test SPH is not as sensitive to . Powell-only methods are unstable on this problem, and lead to artificial field amplification.
      
Figure 7: As Figs. 5-6, but with a different initial particle configuration: here a triangular lattice with constant particle mass (with fixed total particle number, so the resolution is a factor of higher within the magnetized loop). In this case, the errors are reduced in all our Lagrangian methods. In particular, the noise in SPH, which is sensitive to both the local particle arrangement and differences in particle masses, is greatly reduced. Qualitatively, the features in all cases are identical, however.

3.4 Advection of a Field Loop: Minimizing Numerical Diffusion

The next test is a standard test of advection errors and numerical dissipation. We initialize a periodic 2D domain: inside a circle of about the origin, we set with . Outside the circle , and everywhere ; this is an equilibrium configuration that should simply be advected. In Fig. 5, we plot images of the magnetic energy density. In Fig. 6, we plot the total magnetic energy in the box as a function of time (which should remain constant at its initial value). For the sake of direct comparison between grid and mesh-free codes, in these plots we take the initial element configuration to be a Cartesian grid, with unequal-mass elements.

Advection of any configuration not perfectly aligned with the grid is challenging in grid codes; here the loop is continuously diffused away, at a rate that increases rapidly at lower resolution. In Lagrangian methods, on the other hand, stable configurations with bulk advection should be advected perfectly. In Paper I, we demonstrate that our MFM & MFV methods can advect arbitrary pressure-equilibrium hydrodynamic configurations (including arbitrary scalar quantities) to within machine accuracy. However, here the introduction of the divergence-cleaning source terms leads to some initial diffusion of in MFM/MFV. This is enhanced by the fact that the particle masses change discontinuously at the “edge” of the loop. But still, we clearly see the benefit of a Lagrangian method: convergence is much faster than in fixed-grid codes (even using CT); the dissipation in our simulation with MFM/MFV is approximately equivalent to that in ATHENA at . In the image, we see slightly more noise; this is the expected “grid noise” which is higher in meshless methods, but the diffusion is less (in particular, the “hole” which appears at the center owing to numerical resistivity is minimized).151515Of course, in any of our Lagrangian methods, it is trivial to obtain perfect evolution of the field loop test for arbitrarily long times, if we simply disable any dissipation terms. In MFM/MFV this amounts to invoking the “energy-entropy” switch described in Paper I (setting it always to “entropy,” i.e. using purely adiabatic fluxes), and in SPH it amounts to disabling the artificial dissipation terms (see examples in e.g. Rosswog & Price, 2007; Price, 2012). Then, because the system is in uniform motion there is zero advection and the evolution is trivial. The same is true for moving-mesh codes. But these changes make the methods numerically unstable (and would produce disastrous errors at any shocks or discontinuities). We therefore will not consider such modifications further.

If we use Powell-only divergence subtraction, the conservation errors associated with non-zero can actually lead to non-linear growth of , such that the total magnetic energy increases! The growth in this case is nearly resolution-independent, since it is sourced around the sharp discontinuity in the field at the edge of the loop. By the end of the simulation, the total magnetic energy in the Powell-only run has increased . The noise and asymmetry in the image have grown severely.

In SPH, there is an initial, brief but unphysical growth in ; this comes from the divergence-cleaning being less effective than in MFM/MFV (so it behaves like the Powell case); however once enough diffusion and particle re-arrangement has occurred, the divergence-cleaning operator can work effectively, and the energy decays at approximately the same rate as our MFM/MFV calculations.

As noted above, the noise and dissipation in the mesh-free methods is enhanced by the artificial jump in element masses at the edge of the loop (this is not a “natural” configuration for our mesh-free methods). We therefore consider in Fig. 7 the same test, with the same total number of elements in the box, with an initial triangular lattice configuration and equal-mass elements. The combination of reduced noise at the loop edge, and a factor higher resolution within the loop (since it is higher-density), does reduce the noise and errors significantly. However, the qualitative behavior in every case is identical.

Figure 8: Hawley-Stone current sheet (§ 3.5). We plot magnetic field lines (arrows indicate local field strength and direction) at time with , , in MFM. MFV, SPH, and grid-methods produce very similar results for these parameters. Reconnection along the current sheets leads to magnetic “islands” which grow and merge. Our new mesh-free methods are able to stably evolve the current sheet indefinitely, with for all cells at all times.
Figure 9: Stability limits of the current sheet problem in Fig. 8. Given an initial pressure , and velocity perturbation with amplitude , we consider the minimum and maximum for which the problem can be evolved stably to time (further to the top-right is more-stable). In the grid-based CT method of ATHENA, the total-energy formulation of the code, coupled with high-accuracy subtraction needed for accurate CT, and advection errors when the fluid moves over the grid, mean that the method will crash (negative pressures result) for or . Combining the duel-energy formalism from Paper I, with a Lagrangian method that moves with the fluid, and using divergence-cleaning instead of CT, we are able to stably evolve the system until reaches machine-error levels , and arbitrarily large (we have not considered larger only because the simulations become too expensive, not because they crash).

3.5 Hawley-Stone Current Sheet: Numerical Stability

This test follows Hawley & Stone (1995). In a 2D periodic domain with , , we initialize and , with for and otherwise. This is not a good test of algorithm accuracy, since the non-linear solution depends sensitively on the numerical dissipation in different methods. However, it is a powerful test of code robustness. Qualitatively, the solution should exhibit rapid reconnection along the initial current sheet, which will launch nonlinear polarized Alfven waves, that generate magnetosonic waves, while magnetic islands form, grow, and merge. For smaller and larger , it becomes more difficult for algorithms to evolve without crashing or returning unphysical solutions (e.g. negative pressures in the Riemann problem).

Fig. 8 shows the magnetic topology at time in a run with , . For these parameters, MFM, MFV, SPH, and grid methods (here, ATHENA) all look very similar.161616For extensive description of this test problem in ATHENA, see http://www.astro.virginia.edu/VITA/ATHENA/cs.html The real test arises when we vary and ; in Fig. 9, we plot the maximum and minimum which we are able to use in each algorithm before the code crashes or returns an unphysical result. ATHENA crashes after some small early-time evolution for or . The low- problem most likely owes to the fact that the method evolves total energy: when the magnetic energy dominates, this causes serious difficulty recovering the correct internal energy (since we must subtract two large numbers), eventually producing negative temperatures. Here we use a dual-energy formalism described in Paper I, which does not conserve total energy to machine error, but can handle essentially arbitrary ratios. So for our SPH, MFM, MFV implementations, is limited by essentially machine error (, depending on the formulation). For , we find similar results; the increased stability owes both to the same dual-energy formalism above, but also to the Lagrangian method, which eliminates the advection errors that, in grid-based codes, become larger with the local fluid velocity. In non-moving grid codes, eventually, at any resolution, there is some bulk velocity which will wipe out the correct physical solution completely (necessitating still-higher resolution); this is avoided in Lagrangian methods. We explore only values of up to in this test because the timestep becomes so small that it is impractical to evolve the system to late non-linear times, but we suspect the robustness of the algorithms should hold to similar machine error levels (i.e. allowing ).

Figure 10: The Orszag-Tang vortex (§ 3.6). We show images of density (in code units, as labeled) at time , in runs with elements (particles/cells). All methods develop the major qualitative features, though there is some additional smoothing in SPH. Note that the contact discontinuities and shocks are captured sharply in MFM/MFV methods.
Figure 11: Comparison of the Orszag-Tang problem from Fig. 10. We plot density , magnetic field components , , and , in horizontal slices at . With cells, MFM, MFV, and grid methods have converged well to the exact solution (except for some small smoothing of the sharpest features). Here, the “exact” line is a ATHENA PPM CT result; at resolution , MFM/MFV/unboosted-ATHENA results are indistinguishable. SPH performs well but shows further smoothing which converges more slowly. We also consider a grid simulation where the fluid is given an additional boost (a uniform ); the Lagrangian (MFM/MFV/SPH) methods are invariant to these boosts. But in grid codes the boost produces a smoothing of the features at ; these require grid resolution to converge away. Using only Powell divergence subtraction leads to a systematic error in which is small, but does not converge away.
Figure 12: Comparison of the Orszag-Tang problem as Fig. 10, but at time , when parts of the flow have broken up into chaotic turbulence. Unsurprisingly the differences between methods have grown, but the results are still qualitatively similar. Note that the additional SPH smoothing in Fig. 10 has here suppressed the formation of the central compact vortex; this feature is most sharply resolved in our MFM/MFV runs.
Figure 13: Quantitative comparison of the Orszag-Tang problem as Fig. 11, but at time . Here we take the slice through the center (), where differences are maximized. The density plot demonstrates the presence of the central vortex in particular; it is suppressed by smoothing in the grid code and disappears in the boosted grid, SPH, and Powell results at this resolution.
Figure 14: Resolution study of the Orszag-Tang vortex in our MFM method (MFV is essentially identical). We show both images as Fig. 10 and slices as Fig. 11, for several resolutions. Most major features are present even at low resolution; convergence with higher resolution is clear for all features. The formal and convergence accuracy is close to ideal scaling for this problem (given that it includes shocks and discontinuities).
Figure 15: Magnetic rotor (§ 3.7). We show images as Fig. 10 of the magnetic pressure (in code units, as labeled), for runs with resolution. Most features appear identical; however note some difficulty in SPH capturing the pressure extrema at , and additional noise in SPH (especially with low neighbor number). With Powell-only cleaning, similar errors and noise appear.

3.6 Orszag-Tang Vortex: Shock-Capturing & Super-Sonic MHD Turbulence

The Orszag-Tang vortex is a standard MHD test which captures a variety of MHD discontinuities, and develops super-sonic MHD turbulence, which is particularly challenging for many methods. In a periodic 2D domain of unit size, we take and set .

Fig. 10 compares images of the resulting density field at ; Fig. 11 quantitatively compares these by plotting the density and , components of in a slice through at the same time. At this (early) time, the flow contains a complicated set of shocks, but has not yet broken up into turbulence. Figs 12-13 consider the same at , when parts of the flow begin to evolve chaotically (amplifying the differences between methods). All runs here have resolution. In all the methods here, all of the key qualitative features are captured at , and most at , including several complicated shocks, discontinuities, and sharp features. At , non-moving grid, MFM, and MFV solutions are essentially indistinguishable they have converged very well to the exact solution already, with only a slight smoothing of the sharpest features (as expected). Their “effective resolution” appears identical. At , qualitative features are similar but quantitative differences appear; the MFM/MFV methods better resolve the central vortex/density peak and sharpest features in (not surprising, since their Lagrangian nature automatically provides higher resolution in this region), but there are stronger oscillations in the low-density regions. However, if we boost the system by a constant velocity, at fixed resolution the stationary-grid solution degrades; for a boost it is more comparable to a run at both times. Obviously the Lagrangian methods avoid this source of error. Divergence errors are well-controlled at both times even in super-sonic shocks.

SPH does nearly as well, although the features in are more smoothed (similar to a MFM/MFV/grid result) at , and this suppresses the appearance of the central density peak at . On this problem, there is not a dramatic difference, interestingly, between SPH with normal versus large neighbor numbers. As we saw in the shocktube tests, the divergence errors are larger by a factor of several in SPH compared to MFM/MFV.

If we use a Powell-only cleaning scheme, the results are not too badly corrupted by errors; however we do see some systematic offsets, particularly in around , and the position of the jumps around at , and the central density peak and discontinuity at , that do not appear to converge away.

Fig. 14 demonstrates the convergence in our MFM method (MFV results are nearly identical here). As we increase the resolution, we clearly see good convergence towards the exact solution in all features here. Quantitatively, the convergence in the and norms of the plotted quantities is first-order, as expected due to the presence of shocks (the same is true in ATHENA). Compare this to older SPH MHD implementations, which only saw convergence in some features, while others converged slowly or not at all (a well-known issue in SPH; Stasyszyn et al. 2013).

We have also considered the pure driven-turbulent box problems from Paper I (rms Mach numbers and ), as well as the ABC dynamo (Arnold et al., 1981) with small initial seed fields; we confirm that the turbulent power spectra and growth rate of magnetic energy in the box agree well between MFM, MFV, and grid methods. This echoes our conclusions for the pure hydrodynamic case in Paper I. With Powell-only cleaning, however, the growth rate of the magnetic energy in the highly super-sonic case is artificially high (growing faster than the flow-crossing time, indicating a clear numerical artifact). In SPH, the results depend on neighbor number, in a manner that we demonstrate in detail for the MRI test problem below.

Figure 16: Density and magnetic fields in a slice through in the MHD rotor problem in Fig. 15. The qualitative differences between methods are the same as in Fig. 11 (the “grid+boost” run uses the same bulk ). MFM/MFV/grid methods agree well. Note that boosted grid methods smooth the pressure spikes significantly, and Powell-only cleaning leads to systematically incorrect shock positions.

3.7 The Magnetic Rotor: Torsional MHD Waves

Next we consider the MHD rotor, a standard problem from Balsara & Spicer (1999) used to test strong torsional Alfven waves. A 2D domain of unit size and , is initialized with uniform , . Inside a circle of size , and ; this is surrounded by a ring-shaped transition region from with , with . At , and .

As with the Orszag-Tang vortex, Figs. 15-16 plot images of the magnetic energy density and field values in horizontal slices, at time . Again, all methods capture the key qualitative features. MFM/MFV and grid methods are very similar. Grid methods converge slightly more rapidly on the sharp field “spikes” at if the bulk velocity is nil, but if the fluid is boosted appreciably, the density spikes are noticeably more smoothed. In either case, MFM, MFV, and grid methods exhibit convergence at the same order. Divergence errors are again well-controlled, even around large discontinuities in .

Here, for SPH, using the same neighbor number as MFM/MFV (3D-equivalent 32 neighbors) leads to significant noise, and some systematic errors in around . These are resolved if a large neighbor number is used. However in both cases significant smoothing of the local extremum in the magnetic pressure (the brown patches at in Fig. 15) is evident, and is less accurately suppressed.

With only Powell-cleaning, significantly more noise is evident in the image; moreover, the systematic offset of the discontinuity positions in is clear, and this does not converge away.

Figure 17: Blastwave in a strongly magnetized medium (§ 3.8). We plot density (as labeled) in simulations. The blastwave expands asymmetrically along field lines (initial ), as expected. All methods capture the key behaviors; Lagrangian methods capture slightly more detail in the high-density regions, but with enhanced grid noise. Stationary-grid methods converge more slowly when the fluid is boosted (here, ); the errors break the symmetry of the solution noticeably at this resolution, but they do eventually converge away. Using only Powell/8-wave methods (no -damping) leads to physically incorrect shapes of the high-density features in the top-right and bottom-left corners; these do not converge away with resolution.
Figure 18: Slices through the MHD blastwave in Fig. 17. Most methods agree well; the convergence is good at this resolution except for SPH without a large neighbor number, or grids if the fluid is moving. The errors from the Powell-only (no-cleaning) method are dramatic: in , we see factor systematically incorrect results.
Figure 19: Magnetic Rayleigh-Taylor instability (§ 3.9). We show the density (as labeled), in both early (top, ) and late-time (bottom, ) results at resolution. MFM, MFV, and grid methods (with no bulk motion) all converge to the correct solution. In the early stages they are nearly identical; even in late stages there is little difference (small differences appear at late times in MFV owing to growth of the slightly-larger grid noise at early times). As in the pure-hydro case, in non-moving grid methods, “boosting” the fluid grid by a uniform velocity slows down mode growth and breaks symmetry, unless significantly higher resolution is used. In SPH, reasonable results can be obtained if large neighbor numbers are used (although some large-scale asymmetry appears in the non-linear solution because of imperfect -cleaning); with smaller neighbor number (the same as used in MFM/MFV methods), low-order errors in SPH dominate and no mode grows. However, the SPH results are also very sensitive to the form of the artificial dissipation (artificial resistivity or AR) employed: we compare the results if we use the fast magnetosonic velocity in AR, as advocated by Tricco & Price (2013) (which reduces noise in super-sonic problems), instead of our default choice (the Alfven velocity). For details, see Appendix A. This produces too much dissipation around the contact discontinuity, which damps the -field to sub-critical values; so the instability behaves (incorrectly) like the pure-hydro case.
Figure 20: Divergence errors in the MHD RT instability from Fig. 19. We plot (as labeled), for the same times as Fig. 19. The grid methods here use CT so to machine precision. MFM/MFV methods maintain (for every particle/cell) even into the late non-linear evolution (most of the values at are boundary effects, in fact). In SPH, the zeroth-order errors around a contact discontinuity lead to less-accurate cleaning in these regions. In the artificial resistivity (AR) case, around these discontinuities; this is not because is large, but because the AR over-damps at the discontinuity. In the Powell case, order-unity errors appear at late times.
Figure 21: Magnetic Kelvin-Helmholtz (KH) instability (§ 3.9). We show density (as labeled) in the non-linear phases and (KH growth timescale ), at resolution. Differences between methods resemble the RT test (Fig. 19). MFM, MFV, and non-boosted grid results agree well into the non-linear growth phase (where different grid noise effects lead to small departures). Boosting the fluid in the grid case or using typical neighbor numbers for SPH break symmetry and suppress the instability at this resolution. Divergence errors in each case closely resemble those in Fig. 20 around the contact discontinuities.
Figure 22: Magnetic “blob” test (§ 3.9). We show density (arbitrary units) at and , with initial vertical field (top) and (bottom) in code units. The weak-field case is essentially identical to the pure hydrodynamic case studied in Paper I, as expected; a combination of non-linear Rayleigh-Taylor and Kelvin-Helmholtz instabilities “shred” the cloud. The strong-field case suppresses cloud destruction. MFV is shown, but differences between methods are qualitatively identical to the RT and KH tests above.

3.8 Magnetized Blastwave: Strong Shocks & Grid-Alignment Effects

Next we consider a strong blastwave in a magnetic field; this is another standard problem, which tests the ability of the code to handle strong shocks and preserve symmetry. In the hydrodynamic version of this test, it is also very challenging for grid-based codes to avoid strong grid-alignment and preferential propagation of shocks along the grid (see Paper I); however these effects are reduced by the asymmetric forces in the MHD problem.

We initialize a 2D periodic domain of unit size with , zero velocity, , , and pressure within a circle at the center of initial size and outside. Figs. 17-18 show images and slices through the solution at time , for runs. Visually, the MFM/MFV, SPH (with high neighbor number), and grid (with zero boost to the fluid) solutions look similar. The blast expands correctly, with the cavity growing more rapidly in the direction along the field lines. There is some additional detail visible in the high-density regions in the Lagrangian solutions (the “dimpling” in the upper-right and lower-left); this is real and appears at slightly higher resolution in the grid calculation as well. The SPH solution is slightly more smoothed along the shock fronts. Here we also show the Powell-only and boosted-grid results; unlike the Orszag-Tang and rotor problems, here the differences are plainly visible by-eye. The Powell case develops incorrect features; the boosted-grid case loses symmetry.

Quantitatively, we see these effects in Fig. 18. Note that MFM/MFV methods exhibit very small values; SPH exhibits larger , but still quite small. In general, different methods agree fairly well quantitatively. However, for the Powell-only case, the value of , in particular, is seriously wrong in the upper-right portion of the solution. We emphasize that the error is worse if we re-simulate this with Powell-only cleaning but a resolution of . The failure of this method to guarantee the correct shock jumps corrupts the entire late-time solution. Note that, in SPH (using the full Dedner et al. (2002) cleaning), the magnitude of the errors is not much smaller than our Powell-only run, but this incorrect behavior does not appear. This emphasizes that the magnitude of is not, by itself, the whole story; non-linear errors are damped by the Dedner et al. (2002) approach (at the cost of some additional local diffusion) which would otherwise build up coherently at later times. Even quite large values of can be tolerated, if the algorithm includes a proper damping formulation (as our MFM/MFV/SPH implementations do).

Figure 23: Magneto-rotational instability (MRI; § 3.10) test problem. We show the results from a 2D axisymmetric shearing box (horizontal/vertical axes correspond to radial/vertical coordinates), with vertical set so the fastest-growing modenumber is . We plot the azimuthal/toroidal component of the magnetic field (the behavior of the radial component is similar, but with reversed sign), in an MFM calculation at resolution (units are scaled to the maximum/minimum values in each frame as labeled, since the absolute value of grows exponentially). The initial seed noise is amplified on the correct timescale, and from times , the dominant mode pattern is clearly visible. At late times, the non-linear modes break up into turbulence, as expected. MFM, MFV and high-order CT-based grid methods (see Guan & Gammie 2008) produce similar results.
Figure 24: MRI, as Fig. 23 (same colorscale; just showing the half of the box), at times when the linear mode dominates, for runs with initial set so that the fastest-growing mode corresponds to , as labeled (MFM shown, MFV is nearly identical). In each case, the correct mode is clearly visible. We find approximately cells/particles are required across each “node” in the linear direction to resolve the correct mode growth, the same as the number of cells in PPM grid methods.
Figure 25: MRI (with mode number ), as Fig. 23, as a function of resolution and numerical method. MFM/MFV methods are very similar, and resolve the MRI with with as few as elements. Modern SPH MHD is, in fact, also able to capture the MRI modes, if a larger neighbor number is used. If we use a smaller neighbor number in SPH, the low-order errors completely swamp the correct solution.

3.9 MHD Rayleigh-Taylor & Kelvin-Helmholtz Instabilities: Fluid Mixing in MHD

Fluid mixing instabilities are astrophysically important and historically challenging for SPH methods; the hydrodynamic forms of these are discussed at length in Paper I. In MHD, non-zero suppresses the growth of small-scale modes. If magnetic fields are too strong, no interesting modes grow. If fields are too weak, the case is essentially hydrodynamic. But there is an interesting, MHD-specific regime when the fields strengths are near-marginal; we consider this in a Rayleigh-Taylor (RT) problem.

We take initial conditions from Paper I and Abel (2011): in a 2D domain with (periodic boundaries) and (reflecting boundaries), we take , , and with , , . Initial pressures are assigned to produce a gradient in hydrostatic equilibrium with a uniform gravitational acceleration in the direction (at the interface, so ). An initial -velocity perturbation with is applied in the range (otherwise the velocities are zero).171717Following the method in Paper I, we use the routines generously provided by R. O’Leary to construct the mesh-free initial conditions, then re-interpolate this onto the ATHENA grid. This is critical to ensure that the same “seed” grid noise is present in both methods, which in turn is necessary to see similar behavior in the late-time, non-linear phase of evolution.

Fig. 19 shows the resulting density field at intermediate and late times, in a simulation. In MFM, MFV, grid/AMR with non-moving fluid, and SPH (with sufficiently large neighbor number), the linear growth of the field is essentially identical; this is consistent with our pure-hydro results in Paper I. Even the non-linear, late-time results agree reasonably well (although there is some symmetry-breaking in SPH which owes to less-accurate -cleaning, even at large neighbor number). There is slightly more “grid noise” in MFV at this resolution (which leads to small differences in the late-time evolution). Just like in the pure hydrodynamic case, in SPH, the results are very sensitive to neighbor number: if a smaller neighbor number is used (say, the same as we use in our MFM/MFV methods), then the initial seed mode is too weak – it is overwhelmed by the zeroth-order errors in the method, and no modes grow. Similarly, if we boost the fluid by a constant velocity in stationary-grid methods, advection errors break symmetry and dramatically suppress the growth of the mode at this resolution (much higher resolution is required to match the accuracy of the Lagrangian methods). All of this is consistent with our results from the hydrodynamic case in Paper I.

If we consider the Powell-only case, the linear mode evolution appears reasonable – the growth rates are only slightly suppressed. However, in the non-linear phase, the solution is totally corrupted! The non-linear errors accumulate, if only Powell-type schemes are used (once again, because in this method, divergence errors are only transported, not suppressed), until they overwhelm the real solution. Clearly, tests restricted to the linear regime are not sufficient to validate divergence-cleaning schemes.

In SPH, we also find another difficulty unique to MHD. Numerical stability in SPH requires somewhat ad-hoc artificial dissipation terms for , the “artificial resistivity.” As discussed in Appendix A, and Rosswog & Price (2007); Tricco & Price (2013), this carries ambiguities that are not present in hydrodynamics: the appropriate “signal velocity” for the resistivity could be the sound speed, Alfven velocity, or magnetosonic speed, and in some cases resistivity should be applied in rarefactions. The correct answer depends on the type of MHD discontinuity. In Godunov methods such as MFM/MFV/grid codes, the correct form of the dissipation is provided by the appropriate solution to the Riemann problem. But in SPH, even with the state-of-the-art “switches” used here, this is difficult to assign correctly. By default, we use the Alfven velocity in this switch. But Tricco & Price (2013) show that this can lead to too-low a resistivity in super-sonic MHD turbulence, in turn producing shock breakup and serious noise (see their Fig. 7). They recommend increasing the dissipation by using the fast magnetosonic speed instead. However, if we do that for this problem, it produces excessive dissipation of the magnetic field around the contact discontinuity.181818For the test in Fig. 19 labeled “ AR,” we keep everything about our SPH MHD method identical (including the “switch” and maximum viscosity ), except replace the Alfven velocity in Eq. 29 with the fast magnetosonic speed. This makes the problem behave (incorrectly) like the pure-hydrodynamic case.

In Fig. 20, we show full 2D maps of the divergence . in MFM/MFV, these are extremely well-controlled, with median values and maxima . In our default SPH implementation, they are larger by a factor of a few. In the “ AR” SPH run, we see much larger values along the contact discontinuity. These are not, however, caused by poor -control; rather, excessive dissipation of around the discontinuity leads to a local sharp depression of (the denominator), as opposed to . In the Powell runs the divergence errors are (as expected) much larger.

We have also compared the magnetized Kelvin-Helmholtz instability, shown in Fig. 21. The initial conditions follow McNally et al. (2012) (see also Paper I), a 2D setup with resolution elements following two streams with , , and and , with a amplitude initial seed mode and small interface region between the streams. We add a uniform , about a factor below the critical value which suppresses the instability, and show results at and (where the KH growth time is ). Here, we obtain qualitatively identical conclusions to our RT test. In the linear and even non-linear stages, MFM/MFV and non-boosted grid results agree well, with more small-scale structure in the non-linear stages in the meshless methods (owing to increased grid noise). Quantitatively, the total magnetic energy in the box grows, with excellent agreement between these methods and converged solutions until (well into the non-linear phase). However, the boost completely suppresses mode growth in stationary grid methods at this resolution. In SPH, a reasonable, answer can be obtained with sufficiently high neighbor number and an appropriate choice for the artificial resistivity, but the instability is totally suppressed if we use typical neighbor numbers and behaves as if the field were much weaker in the “ AR” case. Comparing , we see the same behavior around the contact discontinuity as in the RT test; the relative performance of different methods is essentially identical, although in all cases the errors are systematically smaller by a factor of a couple.

Finally, in Fig. 22, we also compare an MHD version of the “blob” test from Agertz et al. (2007). The setup is identical to our detailed study of different methods on this problem in Paper I (featuring a cold cloud in pressure equilibrium in a hot wind tunnel), but with an initially constant field in the vertical direction. As expected the field strongly suppresses the non-linear RT and KH instabilities that tend to disrupt the cloud in the hydrodynamic case (for a detailed study, see Shin et al., 2008). Our qualitative conclusions (comparing different methods) are identical to the tests above.

Figure 26: Divergence errors in the MRI test in Fig. 25 at resolution. We plot (as labeled), for two times, when the fastest-growing mode dominates, and when the system has broken up into turbulence (Fig. 23). In both cases, the errors are well-controlled in MFM/MFV; in SPH some regions reach larger in the turbulent phase; however these almost entirely correspond to regions of nearly-vanishing (i.e. ), so the errors are not dynamically significant.
Figure 27: Growth of the MRI modes in Figs. 23-25. We measure the Fourier amplitude in the azimuthal (; left) and radial (; center) magnetic fields, as well as the average magnetic energy ( in the box (right). We compare each to the expectations of linear theory (dotted black line). We compare three different resolutions (, , , in progressively thicker lines). MFM/MFV methods give similar results; in both cases the simulations converge to the correct linear growth rate quickly at resolutions above