A robust numerical scheme for highly compressible magnetohydrodynamics: Nonlinear stability, implementation and tests

A R OBUST NUMERICAL SCHEME F OR HIGHL Y COMPRESSIBLE MA GNETOHYDR OD YNAMICS: NONLINEAR ST ABILITY, IMPLEMENT A TION AND TESTS K. W AAGAN, C. FEDERRA TH, AND C. KLINGENBER G Abstract. The ideal MHD equations are a central mo del in astrophysics, and their solution relies up on stable numerical sc hemes. W e present an implemen tation of a new metho d, whic h possesses excellent stabilit y properties. Numerical tests demonstrate that the theoretical stability prop erties are v alid in practice with negligible compromises to accuracy . The result is a highly robust scheme with state-of-the-art efficiency . The sc heme’s robustness is due to entrop y stabilit y , p ositivity and properly discretised P ow ell terms. The implementation tak es the form of a mo dification of the MHD mo dule in the FLASH co de, an adaptiv e mesh refinement code. W e compare the new sc heme with the standard FLASH implemen tation for MHD. Results show comparable accuracy to standard FLASH with the Ro e solver, but highly improv ed efficiency and stability , particularly for high Mach num ber flows and low plasma β . The tests include 1D sho ck tub es, 2D instabilities and highly supersonic, 3D turbulence. W e consider turbulent flows with RMS sonic Mach num b ers up to 10, t ypical of gas flo ws in the interstellar medium. W e in vestigate both strong initial magnetic fields and magnetic field amplification b y the turbulent dynamo from extremely high plasma β . The energy spectra sho w a reasonable decrease in dissipation with grid refinemen t, and at a resolution of 512 3 grid cells we iden tify a narro w inertial range with the exp ected p ow er-law scaling. The turbulent dynamo exhibits exponential growth of magnetic pressure, with the gro wth rate twice as high from solenoidal forcing than from compressive forcing. Two v ersions of the new scheme are presen ted, using relaxation-based 3-wa ve and 5-wa ve approximate Riemann solvers, resp ectively . The 5-wa ve solver is more accurate in some cases, and its computational cost is close to the 3-wa ve solv er. 1. Introduction In order to mo del complex nonlinear astroph ysical phenomena, numerical solutions of the equations of ideal MHD are cen tral. Examples include the dynamics of stellar atmospheres, star formation and accretion discs. Magnetised gas flo ws in astroph ysics are t ypically highly compressible, nonlinear and often sup ersonic. Hence, obtaining stable numerical results in such regimes is extremely challenging if the numerical sc heme m ust also be b oth accurate and efficien t. Here we presen t a new MHD solv er that preserv es stable, ph ysical solutions to the compressible MHD equations by construction, while at the same time sho wing impro v ed efficiency and comparable accuracy to standard schemes based on the v ery accurate, but less stable Roe appro ximate Riemann solver [59]. Date : October 15, 2018. K. W aagan would lik e to thank Dr. Mic hael Kn¨ olker at the High Altitude Observ atory in Colorado for advice and support with this w ork. K. W aagan is partially supported by NSF Gran ts DMS07-07949, DMS10-08397, DMS10-08058, and ONR Grant N000140910385. C. F ederrath is grateful for support from the Landesstiftung Baden-W ¨ urttem b erg via their program International Col- laboration I I under Grant P-LS-SPI I/18, and has received funding from the Europ ean Researc h Council under the European Communit y’s Seven th F ramework Programme (FP7/2007-2013 Grant Agreemen t No. 247060). C. F ederrath ac kno wledges com- putational resources from the HLRB II pro ject Grant pr32lo at the Leibniz Rec henzentrum Garc hing for running the turbulence simulations. The FLASH co de w as developed in part by the DOE-supported Alliances Center for Astroph ysical Thermonuclear FLASHes (ASC) at the Universit y of Chicago. 1 W e consider numerical solutions of the ideal MHD system (here written in conserv ation form, and in three dimensions, letting I 3 denote the 3 × 3 identit y matrix) (1.1) ρ t + ∇ · ( ρ u ) = 0 , ( ρ u ) t + ∇ · ( ρ u ⊗ u + ( p + 1 2 | B | 2 ) I 3 − B ⊗ B ) = 0 , E t + ∇ · [( E + p + 1 2 | B | 2 ) u − ( B · u ) B ] = 0 , B t + ∇ · ( B ⊗ u − u ⊗ B ) = 0 , ∇ · B = 0 , with the sp ecific internal energy e , such that the total energy density is given b y E = ρe + 1 2 ρ u 2 + 1 2 | B | 2 , and the pressure giv en b y the equation of state p = p ( ρ, e ). The Cartesian components of the v elo cit y are denoted b y u = ( u, v , w ). The system fits the generic form of a conserv ation law U t + ∇ · F ( U ) = 0, except for the restriction on ∇ · B . How ever, if this restriction is satisfied at the initial time t = 0, it automatically holds at later times t > 0 for the exact solution. Sho ck conditions are common in astrophysics, hence (1.1) should b e understo o d in a weak sense. The lack of regularity and large scale ranges in astrophysical flo ws make it particularly challenging to devise numerically stable schemes. Recent dev elopments in nonlinear stability analysis ha v e made suc h schemes p ossible, how ever. The central stability notions are entrop y stability and p ositivit y of mass density ρ and internal energy ρe . This pap er presen ts an implementation of the robust p ositiv e second order scheme from [72], whic h uses the en tropy stable approximate Riemann solvers of [4, 5]. Sev eral codes are a v ailable for high-p erformance compressible MHD simulations. W e men tion Athena [65], AstroBEAR [13], Chombo [60], Enzo [73], FLASH [23], Nirv ana [75], Pluto [49], RAMSES [22] and V A C [68]. The positive sc heme of [72] has so far b een unav ailable in suc h co des. W e hav e implemented it as an alternativ e MHD mo dule in the FLASH co de. W e fo cus here on a comparison with the original FLASH co de. Benc hmark tests in one and t w o spatial dimensions are presen ted. T ransition to unstable turbulence-lik e flo w is particularly emphasised. Additional tests of the algorithms ma y b e found in [5, 35] and [72]. As an example application we consider turbulence in molecular clouds, which is c haracterised by high Mac h n umbers and the presence of magnetic fields. Before presen ting the results, w e give a review of the numerical schemes and their theoretical prop erties. 2. Numerical schemes In this section we give a quick summary of the numerical techniques used, and their justification. Consider first a system of conserv ation laws in one spatial dimension (we consider the x -dimension) U t + F ( U ) x = 0. The numerical sc hemes w e will use are of the finite v olume type (2.1) S ∆ t U i = U i − ∆ t h  F i + 1 2 − F i − 1 2  , where U i are av erages o ver in terv als (or ’cells’) of length h indexed by i at a giv en time t , and the op erator S ∆ t up dates the cell av erages to time t + ∆ t . The numerical fluxes F are ev aluated at the cell in terfaces, hence U is conserved and w e call the sc heme conserv ative. 2.1. First order accuracy in one space dimension. First order accurate schemes can b e given b y (2.2) F i + 1 2 = F ( U i , U i +1 ) . T ypically , F ( · , · ) is giv en by an approximate Riemann solver [3, 67, 41]. The follo wing are considered here: The Ro e solver and the HLLE (Harten, Lax, v an Leer and Einfeldt) solv er as implemen ted in FLASH, and the multi-w av e HLL-type solvers of [4, 5] (named HLL3R and HLL5R, or HLLxR). Other HLL-type solv ers are giv en in [27, 50, 26, 42, 44, 25] among others. The numerical flux should b e consistent (i.e. F ( U, U ) = F ( U )), and satisfy appropriate stability criteria. F or gas and plasma dynamical systems a strong and physically meaningful stability criterion is implied by the second law of thermo dynamics, which reads in one dimension: (2.3) ( ρs ) t + ( ρus ) x ≥ 0 , s b eing the sp ecific en tropy of the system. This inequality holds, by definition, for an exact Riemann solver, and should also hold for an approximate Riemann solver. The HLLxR solv ers w ere pro v ed in [4, 5] to satisfy (2.3), and consequently the resulting conserv ativ e sc heme satisfies a discrete version of (2.3). W e refer to an appro ximate Riemann solver satisfying (2.3) as b eing entrop y stable. F or the Ro e and HLLE solv ers there 2 are no rigorous pro ofs of en tropy stability . F or HLLE this app ears not to b e a problem, while Ro e solvers yield non-physical sho cks at sonic p oints. Therefore, Ro e solvers are typically equipp ed with a so-called en tropy-fix, whic h adds extra dissipation at sonic p oints. The Roe solver also tends to pro duce negative mass density v alues, while HLLE and HLLxR prov ably preserv e p ositive mass density . Positivit y of mass densit y and en trop y stability imply that the internal energy remain p ositive. The p ositivity of density and in ternal energy ensures that the hyperb olicity of the e quations do es not break down (in whic h case the n umerical schemes no longer mak e sense), and implies that the computed conserved quan tities are stable in L 1 [55]. Entrop y stabilit y and p ositivity are referred to as nonlinear stability conditions since they apply directly to nonlinear systems. The ’R’ in HLLxR comes from the deriv ation of the solv ers from a relaxation appro ximation to the MHD equations, hence they are referred to as relaxation solvers. The underlying relaxation approximation has formally b een found to b e entrop y dissipative through a Chapman-Enskog analysis [4]. The corresp onding relaxation solv er for hydrodynamics [2], was tested in an astrophysical turbulence setting in [34]. Appro ximate Riemann solvers tend to v ary in accuracy according to the lev el of detail they take in to accoun t. T ypically , some w av e mo des are b etter resolved than others. F or the solvers considered here, the HLLE solver can only optimally resolve the fast w av es, while the Ro e solv er can, in principle, resolv e all isolated wa v es as w ell as the exact solv er. The relaxation solv ers lie somewhere in betw een: The HLL3R can optimally resolve material con tacts and tangential discon tinuities at v anishing Alfv´ en sp eed, and the HLL5R additionally resolv es velocity shears at v anishing Alfv ´ en sp eed. Hence HLL5R is a true generalisation of the HLLC solv er for the Euler equations. These differences in resolution most prominen tly app ear when the w av e mo des in question are close to stationary , i.e. muc h slow er than h/ ∆ t . 2.2. Second order accuracy in one space dimension. In order to obtain second-order accuracy , some in terp olation, i.e. reconstruction of the states from cell av erages, has to b e employ ed. W e consider here the conserv ative MUSCL-Hanco c k sc heme, introduced in [70] (see also [67]). W e refer to [72] for the non- conserv ative case. Let W denote the primitive state v ariables ( ρ, u , B , p ). Awa y from discontin uities the equations can b e rewritten as W t + A ( W ) W x = 0 for a matrix A ( W ). The MUSCL-Hancock algorithm go es as follows: (1) Reconstruction: ev aluate discrete differences D W i . F or oscillation con trol, w e use the MC-limiter (monotonised central limiter), so for each component of W i w e take (2.4) D W i = σ i min  2 | W i +1 − W i | , 1 2 | W i +1 − W i − 1 | , 2 | W i − W i − 1 |  with (2.5) σ i =      1 , W i +1 − W i > 0 , W i − W i − 1 > 0 − 1 , W i +1 − W i < 0 , W i − W i − 1 < 0 0 , otherwise . (2) Prediction step: ev aluate (2.6) W c i = W i − ∆ t 2 h A ( W i ) D W i (3) Ev aluate the cell edge v alues (2.7) W − i = W c i − 1 2 D W i , W + i = W c i + 1 2 D W i . (4) Use the cell edge v alues as input to the numerical flux in the conserv ative sc heme, (2.8) S ∆ t U i = U i − ∆ t h  F  U + i , U − i +1  − F  U + i − 1 , U − i  , where U denotes the conserved state v ariables. Using primitive v ariables as the basis for reconstructing the states w as recommended for example in [12]. It ensures that material contact discontin uities are repro duced exactly (and also shear wa v es at v anishing 3 Alfv ´ en sp eed). Alternativ ely , one can use that the primitive form W t + A ( W ) W x = 0 may be diagonalised as (2.9) ( R j i ) t + λ j ( R j i ) x = 0 , j = 1 , 2 , . . . , d with W i = X i R i . The matrix X i is given by the eigenv ectors X j i of A ( W i ). The reconstructed gradien ts ma y then b e ev aluated using the relation D W i = X i D R i , given gradien ts D R i , and we get (2.10) W ± i = X i R ± i = X i ( R i − ∆ t 2 h λR i ± 1 2 D R i ) = W i − ∆ t 2 h A ( W i ) D W i ± 1 2 D W i , where λ is the diagonal matrix having the eigen v alues λ j as entries. The original FLASH co de uses the c haracteristic v ariable-based reconstruction, while we use the primitive v ariables with the relaxation solvers. The limiter (2.4) forces the scheme to b ehav e like the first order scheme near a discontin uity , hence ensuring that the entrop y dissipation of the approximate Riemann solver kicks in. Even so, it turns out to b e imp ortan t to ensure the p ositivit y of mass density and internal energy . In [72], a limiting pro cedure was presen ted such that the MUSCL-Hanco ck sc heme is p ositive for any p ositive approximate Riemann solver. It consists of a modification that applies to any piecewise linear reconstruction of the W or R -v ariables. This limiter (more sp ecifically the ’ W -reconstruction’ of [72]) is used with the relaxation solvers. 2.3. Multidimensionalit y. Clearly our scheme is based on the treatmen t of one-dimensional Riemann problems. Multidimensional Riemann problems are muc h to o complicated to b e used as a basis of schemes, so v arious means hav e b een suggested to deriv e m ultidimensional schemes from one-dimensional schemes. W e only consider Cartesian grids here, where one can rely on dimensional splitting. This metho d consists of first solving the system ignoring the y and z - deriv atives, then ignoring the x - and z -deriv atives etc. The splitting adds an error, and in order to maintain second order accuracy , the order of the directions must b e rev ersed b etw een each time step (Strang splitting). Unsplit metho ds hav e b een dev elop ed to eliminate this error completely (e.g. in [39, 38]). W e limit our discussion to the split scheme, and fo cus on issues sp ecific to the multidimensional MHD equations. Note how ever that these issues, as well as the prop osed solutions, are not sp ecific to split schemes. Due to considering the directions separately , w e hav e to mak e sense of on e-dimensional problems where the longitudinal magnetic field comp onen t is non-constant, in other words we can no longer assume ∇ · B = 0 initially . In order to handle the resulting more general one-dimensional problems, w e use the follo wing approac h of [56], consisting of mo difying the ev olution equation for B to (2.11) B t + ∇ · ( B ⊗ u − u ⊗ B ) − u ∇ · B = 0 . Equation (2.11) implies that (2.12) ( ∇ · B ) t + ∇ · ( u ∇ · B ) = 0 , consequen tly ∇ · B remains zero analytically , and numerical errors in ∇ · B are transp orted by the flow. In the original FLASH co de, the momen tum and energy equations were also modified, whic h is what was actually recommended in [56]. W e found it sufficient for s tabilit y in [5] and [72] to only include the P o well term for the induction equation. The one-dimensional Po well system is only conserv ativ e when ( B n ) is constant, so the sc hemes must take the form (2.13) S ∆ t U i = U i − ∆ t h  F l ( U + i , U − i +1 ) − F r ( U + i − 1 , U − i )  − ∆ tS i ( U i − 1 , U i , U i +1 ) . In the original FLASH co de, the source contributions from the cell edges are ignored, and the Po w ell source term is simply ev aluated by a central discretisation. In [5] and [72], we found that it w as essential for stabilit y , b oth practically and theoretically , to discretise the P ow ell term in a prop er upwind manner. This w as carried out in [4] for the first order sc heme b y extending the relaxation system to the P o well system, see also [25]. The second order extension can be found in [72], and n umerical tests demonstrating the imp ortance of including a prop erly discretised Po well source are in [25, 35] and [72]. W e include div ergence cleaning in the form of the parab olic cleaning method of [47]; see also [14]. W e use it as already implemen ted in FLASH. It has the adv an tage of not expanding the numerical stencil of the sc heme. FLASH also has the non-local pro jection metho d of [6] implemen ted, whic h could in principle be used 4 Sc heme HLL3R / HLL5R FLASH-Ro e FLASH-HLLE Riemann solver HLL3R/HLL5R Ro e HLLE Reconstruction p ositiv e, primitive c haracteristic characteristic P ow ell term B only full full P ow ell term discretisation up wind cen tral cen tral T able 1. Summary of the sc hemes. also with the new schemes. An alternative to using (2.11) is the staggered mesh (or constrained transp ort) approac h (review ed in [69]), whic h in one dimension essentially means ev aluating B n at the cell in terface instead of as a cell av erage. The numerical stability results w e rely on (entrop y stability and positivity) are not av ailable for staggered metho ds, but staggered schemes ha v e the adv antage of guaranteeing that ∇ · B = 0 to approximation order in smo oth regions. Note that all approximate Riemann solvers w e consider ma y also b e used as comp onents of staggered metho ds. 2.4. Isothermal gas. As a simple mo del of cooling (e.g., due to radiation), it is common to assume that the gas temperature is constant. This means that p = c 2 s ρ , where the constan t c s is the sound sp eed. A practical wa y to implemen t isothermality consists of (i) solving the full equations for ideal gas with a γ close to 1 within eac h time step, and (ii) pro jecting to the isothermal state betw een each time step. W e use γ = 1 . 00001. This pro jection metho d is preferred due its simplicity , accuracy , and generalisability . By generalisabilit y w e mean that th e pro jection to p = c 2 s ρ is easily replaced with more elaborate co oling models. The pro jection metho d requires robust underlying numerics, since the full energy equation has to b e solved prop erly . 2.5. The sc hemes. The changes made to the original FLASH co de in order to implement the new sc heme are as describ ed ab o ve. The time step size is chosen from the maximum w av e sp eed of the previous time step and the maximum of q c 2 s + B 2 ρ o ver all cell a verages ( c s denoting sound sp eed). W e use a CFL-num b er of 0.8 throughout, unless otherwise noted. W e used version 2.5 of FLASH, but it is straigh tforward to implement the new scheme in other finite volume co des. F or n umerical testing, we compare the new schemes HLL3R and HLL5R with the basic FLASH schemes for MHD. FLASH offers b oth a Ro e and an HLLE approximate Riemann solv er. The relev ant schemes are summarised in T able 1. All sc hemes are finite volume sc hemes of MUSCL-Hanco c k t yp e, with differences in the (approximate) Riemann solv er, reconstruction metho d and P o w ell terms as describ ed abov e. The names of the schemes are taken from the Riemann solvers for simplicit y , although other algorithmic differences may also b e imp ortan t. T o distinguish the schemes from the Riemann solvers, we use the ‘sans serif ’ font to indicate the scheme and normal fon t type for the Riemann solv ers (see T ab. 1). The thermo dynamics are either adiabatic with an ideal gas equation of state or isothermal. 3. Numerical studies 1D In one spatial dimension, the differences b et ween the sc hemes are giv en mainly b y the accuracy of the appro ximate Riemann solver and the stability prop erties of b oth the appro ximate Riemann solver and the reconstruction. 3.1. Brio-W u sho c k tub e. W e first consider initial sho ck tub e data from [8] that hav e become a standard test case. It contains a v ariety of n umerically challenging w a v e types o ccurring in MHD flows, in particular a fast rarefaction, a comp ound wa ve, a material con tact, and a slow sho c k. The initial data are given by U = U l for x < 0 . 5, and U = U r for x > 0 . 5, with γ = 5 / 3 and ρ l = 1 , u l = 0 , B l = (0 . 75 , 1 , 0) , p l = 1 ρ r = 0 . 125 , u r = 0 , B r = (0 . 75 , − 1 , 0) , p r = 0 . 1 . Figure 3.1 sho ws small differences b et ween HLL3R and FLASH-Ro e on this standard test with FLASH-Roe sligh tly sharper on the compound wa ve. FLASH-HLLE is similar to HLL3R , except that the material contact w as more smeared out than in HLL3R (see Fig. 3.2). W e conclude that all the sc hemes p erform reasonably 5 Figure 3.1. Brio-W u sho ck tub e. Resolution h = 256 − 1 , time t = 0 . 1. FLASH-Ro e and HLL3R are compared. The w a v es seen are from left to right, a fast rarefaction, a slow mo de comp ound w a ve, a material contact discon tin uity and a slow sho c k. The smaller amplitude righ t-going Alfv´ en and fast wa v es are omitted for clarity . Figure 3.2. The same exp eriment as Figure 3.1. FLASH-HLLE and HLL3R are compared. 6 Figure 3.3. Low β expansion wa ve. Resolution h = 256 − 1 . Adiabatic case, v elo cit y ˆ u = 3 . 1 at time t = 0 . 05. FLASH-Ro e and HLL3R are compared. Figure 3.4. Low β expansion wa ve, isothermal gas. Resolution h = 256 − 1 . V elocity ˆ u = 3 . 1 at time t = 0 . 05. FLASH-HLLE and HLL3R are compared. w ell on this test case, with some differences in accuracy dep ending on the level of detail in the approximate Riemann solver. 3.2. Lo w β expansion w av es. This is a test from [5]. It consists of rarefactions into a region of lo w plasma β (defined as β = 2 p B 2 ). The initial data are given by U = U l for x < 0 . 5, and U = U r for x > 0 . 5, with γ = 5 / 3 and ρ l = 1 , u l = ( − ˆ u, 0 , 0) , B l = (1 , 0 . 5 , 0) , p l = 0 . 45 ρ r = 1 , u r = ( ˆ u, 0 , 0) , B r = (1 , 0 . 5 , 0) , p r = 0 . 45 . Hence the sound sp eed is initially √ 3 / 2 ≈ 0 . 87. A t ˆ u = 3 . 1 there are only small differences b et ween FLASH-Ro e and HLL3R (Fig. 3.3). W e also p erformed this test for the isothermal case with c 2 s = 0 . 45. FLASH-Ro e could only handle ˆ u - v alues up to ab out 2.9. F or higher v alues it failed to run due to negative pressure and density . In con trast, FLASH-HLLE and HLL3R prov ed quite stable also in the isothermal case. They are compared in Figure 3.4. HLL3R is able to represen t low er densities than FLASH-HLLE during the expansion. 4. Numerical studies 2D In tw o dimensions w e consider tw o standard test cases from the numerical literature and a Kelvin- Helmholtz instability . The tests exhibit strong sho c k wa ves and complicated wa ve interactions. Regarding n umerical stability , in addition to the approximate Riemann solver and the reconstruction, the influence of differen t discretisations of the Po well term are inv estigated. 7 Figure 4.1. Orszag-T ang test. Pressure at time t = 0 . 5, with resolution h = 256 − 1 . Left: FLASH-HLLE , righ t: HLL3R . The result from FLASH-Ro e was visually hard to discern from HLL3R . Figure 4.2. B x along a slice in the y -direction for the Orszag-T ang test. HLL3R (dashed) and FLASH-HLLE (solid). The result from FLASH-Ro e was indistinguishable from HLL3R . 4.1. Orszag-T ang. This is a standard test case for ideal MHD schemes. It consists of smo oth velocity and magnetic field profiles going through sho c k formations, complicated w av e in teractions and ev en tually forming instabilities. The initial data are given b y γ = 5 / 3, (4.1) ( ρ, u , B , p ) = (1 , sin( π y ) , sin( π x ) , 0 , − sin( π y ) , sin(2 π x ) , 0 , 0 . 6) . A t t = 0 . 5 w e observed only minor differences b etw een HLL3R and FLASH-Ro e , while FLASH-HLLE gav e a smo other pressure p eak along the central curren t sheet (see Fig. 4.1). Plotting B x along the cut at x = 0 . 428 (Fig. 4.2) sho ws that FLASH-HLLE gives a slightly more smeared-out current sheet. HLL3R and FLASH-Ro e ga ve v ery similar results at t = 0 . 5. A t t = 1 a magnetic island can b e observed at the domain centre with FLASH-Ro e and HLL3R , while with FLASH-HLLE this instability is suppressed by n umerical viscosity (Fig. 4.3). The HLL5R was also tested and ga v e results very similar to HLL3R and FLASH-Ro e . W e consider computational times in T able 2. They are from a computation on t wo pro cessors, in the form of total CPU time in seconds and in the p ercentage of that computation time sp ent in the sweep routine (whic h is the only differing one). This should be taken with a grain of salt, as these n umbers v ary a bit from time to time. W e to ok the av erage of tw o runs. It app ears that HLL3R giv es ab out a 10% sp eed-up compared to FLASH-Ro e . FLASH-HLLE is also a little slow er than HLL3R , which is due to the different reconstruction pro cedures. All in all we would exp ect FLASH-HLLE to b e faster than FLASH-Roe , since only 8 Figure 4.3. Orszag-T ang test. Pressure at time t = 1 with resolution h = 256 − 1 . T op to b ottom: FLASH-Ro e , HLL3R , FLASH-HLLE . 9 Sc heme HLL3R HLL5R FLASH-HLLE FLASH-Ro e N time steps 477 539 476 482 Sw eep time fraction 43.6% 47.6% 47.2% 50.0% Av erage time / step 0.555 s 0.520 s 0.571 s 0.604 s T able 2. Computational times for the Orszag-T ang test. the computation of the numerical flux differs. The time to compute the numerical flux in HLL3R should b e somewhere in b et ween the time for FLASH-Ro e and FLASH-HLLE , and in addition HLL3R has a simpler reconstruction and a slightly more elaborate P ow ell term than FLASH-Ro e and FLASH-HLLE . The HLL5R sc heme is the most efficient, how ever it requires slightly smaller time steps. The shorter time steps are not surprising, since the wa ve speeds of the HLL5R Riemann solver may be up to 10% larger than with the corresp onding 3-wa ve solver in some configurations; see [5]. Consequently , there would b e a roughly 10% reduction in time step size. 4.2. Cloud-sho c k interaction. The initial conditions are the same as in [25]. The computational domain is ( x, y ) ∈ [0 , 2] × [0 , 1] with artificial Neumann type b oundary conditions. The initial conditions consist of a sho c k moving to the righ t, lo cated at x = 0 . 05, and a circular cloud of density ρ = 10 and radius r = 0 . 15 cen tred at ( x, y ) = (0 . 25 , 0 . 5). W e to ok γ = 5 / 3 and ρ =      3 . 86859 if x < 0 . 05, 10 . 0 if ( x − 0 . 25) 2 + ( y − 0 . 5) 2 < 0 . 15 2 , 1 . 0 otherwise, (4.2) u = ( (11 . 2536 , 0 , 0) if x < 0 . 05, (1 . 0 , 0 , 0) otherwise, (4.3) p = ( 167 . 345 if x < 0 . 05, 1 . 0 otherwise, (4.4) B = ( (0 , 2 . 18261820 , − 2 . 18261820) if x < 0 . 05, (0 , 0 . 56418958 , 0 . 56418958) otherwise. (4.5) Hence, the sho ck has a sonic Mach n um b er around 11.8 and an Alfv´ enic Mac h num b er around 19.0. The solution b ecomes very complicated as the sho ck interacts with the cloud, and instabilities are observed along filaments of the c loud. W e used these initial data to test the ability of the new scheme to handle adaptiv e mesh refinement. Eigh t levels of refinemen t w ere used with the highest resolution corresponding to a resolution of 1024 × 2048 grid cells. The refinemen t criterion w as based on the mass density ρ with the FLASH co de parameters refine cutoff = 0 . 75 and derefine cutoff = 0 . 25. Mass density profiles from HLL3R and FLASH-Ro e are shown in Figure 4.4. The fron t of the cloud looks more rounded with HLL3R , and the structures of the turbulent wak e differ somewhat. Both FLASH-Ro e and HLL3R agree qualitatively well, but since there is no analytic solution it is unclear how to make a quan titativ e assessmen t of the differences. HLL5R gav e v ery similar results to HLL3R in this test. FLASH-HLLE roughly repro duced the shap e of the cloud front, but ga v e a different structure in the turbulen t wak e. 4.3. Kelvin-Helmholtz instabilit y. The Kelvin-Helmholtz instability is imp ortant as a mo del for the onset of turbulence. F rom a numerical point of view, we found it illustrativ e of some differences b et w een the Riemann solvers. W e consider a case where b oth the v elo cit y and magnetic fields switch p olarity along a straight line. The initial data are given b y p = 3 / 5, γ = 5 / 3, ( v , w , B y , B z ) = 0 and ρ = 2 , u = 0 . 5 , B x = 0 . 5 , (4.6) for y > 0, and ρ = 1 , u = − 0 . 5 , B x = − 0 . 5 , (4.7) 10 Figure 4.4. Cloud-sho ck interaction with 8 levels of refinement, corresp onding to 1024 × 2048 cells. Grid and mass densit y . T op: FLASH-Ro e , b ottom: HLL3R for y < 0. W e p erturb the velocity in the y -direction b y (4.8) v = 0 . 0125 sin(2 π x ) exp  − 100( y − 0 . 5) 2  . A t y = − 0 . 5 and y = 0 . 5 we imp ose reflecting b oundary conditions. Hence, we consider a s mall p erturbation to a tangen tial discontin uity in b oth u and B along the line y = 0. Initially the relativ e sonic Mac h num b er is 1, the magnetosonic Mac h num b er is 2 / √ 5, while plasma β = 8 γ . W e used a resolution of 512 2 cells. In n umerical exp erimen ts of this type, instabilit y growth tends to o ccur at the small scales, and therefore dep end on numerical viscosit y . In [34], this was reflected in a strong dependence b etw een grid resolution and instabilit y growth rate. Mass densit y profiles at time t = 1 are sho wn in Figure 4.5. The higher n umerical viscosit y of the HLLE solv er results in FLASH-HLLE suppressing all instabilities at this time. HLL3R yields a sharper shear-discontin uit y than FLASH-HLLE , but the roll-ups still did not dev elop. In contrast, HLL5R and FLASH-Ro e pro duce well-dev eloped Kelvin-Helmholtz roll-ups. The different b ehaviour of the sc hemes can b e explained b y the appro ximate Riemann solvers. By the remarks at the end of section 2.1, Ro e and HLL5R exactly resolv e the unperturb ed initial configuration, while HLL3R smears out the velocity shear slightly . HLLE smears out the material and magnetic contact discon tinuities of the initial configuration, with the result that the instability is en tirely suppressed at early times. The presence of small perturbations to grid-aligned stationary discon tinuities makes this test case particularly sensitive to these prop erties of the Riemann solvers. In particular, we see the effect of exactly 11 Figure 4.5. Kelvin-Helmholtz instability at time t = 1. T op to b ottom: FLASH-Ro e , HLL5R , HLL3R , FLASH-HLLE . 12 Figure 4.6. Kelvin-Helmholtz at time t = 9, FLASH-Ro e (left) and HLL5R (right) . resolving the discontin uities, i.e., having a sc heme that is ’well-balanced’ with resp ect to the unperturb ed state. In non-stationary or non-grid-aligned cases, none of the sc hemes ha v e this well-balancing prop ert y , hence we expect differences to b e less pronounced. A t late times of the instabilit y dev elopment, it is harder to assess the performance of the sc hemes. Results from FLASH-Ro e and HLL5R are shown in Figure 4.6. Both schemes give a similar structure of tw o larger v ortices. FLASH-Ro e pro duces small-scale w a v es in the surrounding medium, whic h are p ossibly a numerical artefact. One p ossible source of spurious features in FLASH-Ro e is the central discretisation of the Po well source term. Such errors ha v e b een p ointed out in [25], [24], [35] and [72]. Finally , in order to demonstrate the robustness of HLL3R and HLL5R , w e tested an isothermal case with a higher sonic Mac h num b er. The sound sp eed was set to 0 . 1, and the initial data were ( v , w , B y , B z ) = 0 ev erywhere, ρ = 1 , u = 0 . 5 , B x = 2 , for y > 0, and ρ = 1 , u = − 0 . 5 , B x = − 2 , for y < 0. The p erturbation was again given by (4.8). Hence, the relative sonic Mach num b er was 10, and the relative Alfv´ enic Mac h num b er was 1 / √ 2. Resolution was uniform with 512 2 grid cells. The standard sc hemes, FLASH-Ro e and FLASH-HLLE failed on this case due to unphysical states, but the new sc hemes, HLL3R and HLL5R pro duced physical results. Density profiles from HLL3R and HLL5R are sho wn in Fig. 4.7. The resulting instability con tains vortex-lik e stru ctures with complicated internal features. Due to the c haotic nature of these solutions, w e cannot ev aluate the differences seen b et w een the schemes. 13 Figure 4.7. High sonic Mac h num b er Kelvin-Helmholtz test for isothermal gas at time t = 1. T op: HLL3R . Bottom: HLL5R . 14 5. Supersonic MHD turbulence In this section we inv estigate the p erformance of three of the MHD sc hemes, FLASH-Ro e , HLL3R and HLL5R on the mo delling of sup ersonic MHD turbulence. T esting the schemes in the regime of highly com- pressible turbulen t flo ws is of particular in terest for astroph ysical applications ranging from molecular cloud formation on galactic scales do wn to scales of individual star-forming accretion disks. These disks can launch p o w erful sup ersonic jets, driven by magneto-cen trifugal forces [1] and magnetic tow ers [45]. Here how ev er, w e concentrate on intermediate scales of molecular cloud ev olution, which is largely determined b y the inter- pla y of supersonic turbulence and magnetic fields. Understanding the nature and origin of sup ersonic MHD turbulence is considered to b e the key to understanding star formation in molecular clouds [46, 15, 61, 48]. 5.1. Numerical setup and p ost-pro cessing. W e start with a short description of the type of setup often adopted in n umerical mo dels of molecular cloud turbulence. It is called turbulence-in-a-box, because a three- dimensional b ox with p erio dic b oundary conditions and random forcing is used to excite turbulence. W e use an isothermal equation of state throughout all turbulence-in-a-b ox runs, i.e., the pressure, p th , is related to the densit y b y p th = c 2 s ρ with the constant sound sp eed, c s . F or these tests we used an initially uniform density , h ρ i = 1 . 93 × 10 − 21 g cm − 3 , corresp onding to a mean particle num b er density of h n i = h ρ i / ( µm p ) ≈ 500 cm − 3 , with the mean molecular w eight, µ = 2 . 3 and the proton mass, m p , and a b o x size of L 3 = (4 p c) 3 . The sound sp eed for the t ypical gas temp eratures in molecular clouds (ab out 11 K) is roughly c s = 0 . 2 km s − 1 . F or driving the turbulence we use the same forcing mo dule as used in [17, 18, 64, 58], whic h is based on the sto c hastic Ornstein-Uhlenbeck pro cess [16]. The sp ectrum of the forcing amplitude is a paraboloid in F ourier space in a small wa ven um b er range 1 < | k | < 3, p eaking at | k | = 2. This corresp onds to injecting most of the kinetic energy on scales of half of the b o x size. The turbulence forcing pro cedure is explained in detail in [19]. W e use a solenoidal (divergence-free) forcing throughout this study unless noted otherwise. F or the comparison of the differen t MHD sc hemes the forcing amplitude w as set suc h that the turbulence reac hes a ro ot-mean-squared (RMS) Mac h num ber, M ≈ 2 in the fully dev elop ed, statistical steady regime. The resolution was fixed at 256 3 grid cells. W e also consider a highly sup ersonic case with M ≈ 10 in section 5.6, for which we inv estigate the resolution dep endence b y comparing runs with 128 3 , 256 3 and 512 3 grid zones. This study , how ever, could only b e follo w ed with the new t yp es of solv ers presented here (HLL3R and HLL5R), b ecause the Ro e solver is numerically unstable for highly sup ersonic turbulence. As shown in previous studies [18, 64, 19] turbulence is fully developed after ab out tw o large-scale eddy turno ver times, 2 T , where one large-scale eddy turno ver time is defined as T = L/ (2 M c s ). The RMS velocity is thus defined as V = L/ (2 T ) = M c s . In order to test the MHD sc hemes in a strongly magnetised case, we set the initial uniform magnetic field strength to B 0 = 8 . 8 × 10 6 Gauss in the z -direction, which gives an initial plasma beta, β = p th /p m , 0 ≈ 0 . 25, with the magnetic pressure, p m , 0 = B 2 0 / (8 π ). No initial fluctuations of the magnetic field were added, suc h that the initial RMS magnetic field is equal to B 0 . Ho wev er, the turbulence twists, stretches, and folds the initially uniform magnetic field, such that the RMS field strength increases due to turbulent dynamo amplification until saturation (see the review by Brandenburg and Subramanian 2005 [7]; and section 5.7). The Alfv´ enic Mach num ber is M A = V /v A = M p β / 2 ≈ 0 . 7, i.e., we consider slightly sub- Alfv ´ enic turbulence. The HLL3R scheme, how ever, was used in b oth the highly sub-Alfv´ enic ( M A  1) and the highly sup er-Alfv ´ enic ( M A  1) turbulent regimes b y [10, 9], who mo delled turbulent flows with β = 0 . 01 . . . 10, M A = 0 . 4 . . . 50 and M = 2 . . . 20. All results w ere av eraged ov er one eddy turno ver time, T , in the regime of fully dev elop ed turbulence, t > 2 T [19]. Since we pro duced output files ev ery 0 . 1 T , we computed F ourier sp ectra and probability distribution functions for each individual snapshot and av eraged afterwards ov er 11 snapshots for 2 . 3 ≤ t/T ≤ 3 . 3. The a veraging procedure allows us to study the temp oral fluctuations of all statistical measures and to improv e the statistical significance of the similarities and the differences seen b et w een the different sc hemes. This is an imp ortant adv antage ov er comparing instan taneous snapshots, because it must b e k ept in mind that individual snapshots can be differen t due to intermitten t turbulen t fluctuations. These fluctuations can either make the results lo ok very similar or v ery different b etw een the different schemes, if one compares instantaneous snapshots only . Th us, a meaningful comparison can only b e made by time a veraging the results and considering the statistically av erage b ehaviour of the schemes (see also, [58]). 15 Figure 5.1. Time evolution of the RMS Mach num ber, M for FLASH-Ro e , HLL3R and HLL5R , sho wn as solid, dotted and dashed lines, resp ectiv ely . T urbulence is fully dev elop ed for t & 2 T , where T is the large-scale eddy turnov er time. 5.2. Time evolution. Figure 5.1 shows the time evolution of the RMS Mach num b er, M = p h u 2 i /c s . The magnetised fluid is accelerated from rest by the random sto chastic forcing as discussed in the previous section. It reaches a statistical steady state at around t ≈ 2 T , such that turbulence is fully developed for t & 2 T (see also, [18, 19]). W e can th us safely use time-a veraging for 2 . 3 ≤ t/T ≤ 3 . 3 of the F ourier sp ectra discussed in the next section. 5.3. F ourier sp ectra. As a typical measure used in turbulence analyses (e.g., [21, 34, 32]), we sho w F ourier sp ectra of the v elo cit y and magnetic fields in Figure 5.2, decomp osed into their rotational part, E ⊥ and their longitudinal part, E k . F or the velocity sp ectra (top panel), E ⊥ u ∝ | u ⊥ | 2 is a measure for the rotational motions of the turbulent flo w, while E k u ∝ | u k | 2 is a measure for the compressions induced by shocks. Since here we drive the turbulence with a solenoidal random force [19], and since the RMS Mach n umber is only sligh tly sup ersonic ( M ≈ 2) the bulk of the turbulent motions is in the rotational part, E ⊥ u . F rom the v elo cit y spectra we conclude that FLASH-Ro e is the least dissipativ e sc heme of the three sc hemes tested. The HLL3R is slightly more dissipativ e, which is seen from the faster loss of kinetic energy at wa v en umbers k & 20. The HLL5R is in b etw een the FLASH-Roe and the HLL3R sp ectrum. Given that the limited n umerical resolution does not allo w for an accurate determination of the turbulence sp ectrum in the inertial range, the differences for k & 20 are acceptable (see also section 5.6 b elow). [32] and [19] find that about 30 grid cells are necessary to resolve the turbulence on the smallest scales. This means that the turbulence on all w av enum b ers k & N / 30 ≈ 8 may b e affected by the limited resolution (here N = 256), which means that the differences b etw een the different sc hemes for k & 20 can b e safely ignored. The magnetic field sp ectra (Figure 5.2, b ottom panel) display small differences for k & 10. As for the velocity sp ectra the differences occur on scales smaller than the scales that are affected b y numerical resolution, which means that these differences are negligible in applications of sup ersonic MHD turbulence. The rotational part of the magnetic field, E ⊥ m , shows slightly less dissipation for the HLL3R and HLL5R compared FLASH-Ro e . The longitudinal part, E k m , is a measure for ( ∇ · B ) 2 . The FLASH-Ro e scheme k eeps the div ergence of the magnetic field smaller by an order of magnitude compared to HLL3R and the HLL5R on all scales. Ho wev er, all schemes main tain the ∇ · B constraint within acceptable v alues. The ratio R E k m d k / R ( E k m + E ⊥ m ) d k < 5 . 5 × 10 − 7 for all times, hence n umerical ∇ · B effects did not hav e any significan t dynamical influence. 5.4. Probabilit y distribution functions. T o further inv estigate the dissipative prop erties of the differen t sc hemes w e sho w probabilit y distribution functions (PDFs) of the v orticity , |∇ × u | , in Figure 5.3 (top left panel). The vorticit y is a measure of rotation in the turbulen t flow. Stronger n umerical dissipation leads to 16 Figure 5.2. Comparison of the F ourier sp ectra of the turbulent velocity ( top p anel ) and the turbulent magnetic field ( b ottom p anel ) for the FLASH-Ro e , HLL3R and HLL5R , shown as solid, dotted and dashed lines, resp ectively . The numerical resolution was 256 3 grid cells. Both the velocity and the magnetic field spectra w ere decomposed into their rotational and compressible parts, E ⊥ and E k , respectively . The temp oral v ariations of the sp ectra are on the order of three times the line width in this plot. a faster decay of small-scale eddies. Figure 5.3 (top left panel) shows that the most dissipative among the three sc hemes tested, HLL3R pro duces smaller v alues of the mean vorticit y by ab out 10% compared to the FLASH-Ro e . The mean vorticit y achiev ed with HLL5R is equiv alent to FLASH-Ro e . The RMS v alues of the v orticity for HLL3R and HLL5R are 11% and 1%, resp ectively smaller than for FLASH-Ro e . This confirms the exp ected ranking of the schemes in terms of their numerical dissipation from the previous section: the HLL3R is the most dissipativ e, while the HLL5R is only marginally more dissipative than FLASH-Ro e . The PDF of the divergence of the v elo cit y field, ∇ · u is shown in the top righ t panel of Figure 5.3. The asymmetry betw een rarefaction ( ∇ · u > 0) and compression ( ∇ · u < 0) is t ypical for compressible, supersonic turbulence (see, e.g., [54, 64]). FLASH-Ro e pro duces sligh tly more zones with higher compression, which can b e attributed to its lo wer numerical dissipation compared to HLL3R and HLL5R (cf., the F ourier sp ectra discussed in the previous subsection). How ever, all schemes agree to within the error bars, indicating strong temp oral v ariations of compressed regions due to the intermitten t nature of the sho cks. The PDF of the gas densit y is an essential ingredien t for analytic mo dels of star formation (e.g., [52, 37, 29]). All those mo dels need the densit y PDF to estimate the mass fraction ab ov e a given density threshold b y integrating the PDF. Th us, man y studies hav e fo cused on inv estigating the density PDF in a sup ersonic turbulen t medium (e.g., [71, 53, 54, 33, 51, 43, 31, 40, 17]). Figure 5.3 (b ottom left panel) shows the PDF of the logarithmic density , s = ln( ρ/ h ρ i ), where h ρ i is the mean volume densit y . In this transformation of the density , the PDF follows closely a Gaussian distribution, i.e., a log-normal distribution in ρ with a standard deviation of σ 2 s = ln  1 + b 2 M 2  (see, e.g., [19]). F rom this expression we derive a turbulence 17 Figure 5.3. Probability distribution functions of the vorticit y , |∇ × u | ( top left ), the di- v ergence of the v elo cit y , ∇ · u ( top right ), the logarithmic density , s = ln( ρ/ h ρ i ) ( b ottom left ), and the magnetic pressure, p m /p m , 0 = B 2 /B 2 0 ( b ottom right ), for FLASH-Ro e , HLL3R and HLL5R , shown as solid, dotted and dashed lines, resp ectiv ely . The numerical resolution w as 256 3 grid cells. The grey hatched regions indicate the 1- σ temp oral fluctuations of the distributions. forcing parameter of b = 0 . 35 ± 0 . 05 for the giv en RMS Mac h n umber, M = 2 . 3 ± 0 . 2 and σ s = 0 . 71 ± 0 . 04 for all three sc hemes. This is in v ery go o d agreemen t with the expected forcing parameter for purely solenoidal forcing as applied in this study , b = 1 / 3 (see [17, 19]). Finally , w e show the PDF of the magnetic pressure, p m = B 2 / (8 π ) in units of the initial magnetic pressure, p m , 0 , in Figure 5.3 (b ottom right panel). No significant differences are seen b etw een the three differen t schemes. Thus, the magnetic pressure distribution is quite robust. 5.5. Computational p erformance of the MHD schemes. All runs with the three different schemes, FLASH-Ro e , HLL3R and HLL5R w ere p erformed on the same machine, the HLRB I I at the Leibniz Sup er- computing Center in Munich in a mo de of parallel computation (MPI) with 64 CPUs. In T able 3 w e list the w all clo c k time p er step, the total n um b er of steps, the total amoun t of CPU hours spent, and the CFL n umber used for each sc heme in our driven MHD turbulence exp eriments with a numerical resolution of 256 3 grid cells. Only for FLASH-Ro e w e had to use a CFL num b er of 0.2 instead of 0.8 (b oth HLL3R and HLL5R ), b ecause it b ecame unstable for CFL n umbers > 0 . 2 and crashed due to the pro duction of cells with negativ e densities. The time p er step w as estimated from av eraging ov er 0 . 1 T within the last turnov er time, i.e., b etw een 3 . 2 ≤ t/T ≤ 3 . 3 b y av eraging ov er 309, 70 and 71 steps for the FLASH-Ro e , HLL3R and HLL5R , resp ectively to a v oid v ariations pro duced b y I/O pro cessing. Considering the time p er step, HLL3R and HLL5R are ab out 4% and 3%, resp ectively faster than FLASH-Ro e . This sp eed-up of HLL3R is slightly less than the sp eed-up measured for the Orszag-T ang test (cf. T ab. 2). The difference b et w een HLL3R and HLL5R is almost negligible. Indeed, the t w o corresponding solv ers ha ve a very similar structure, except that 18 Sc heme Time p er step [s] T otal steps T otal CPU hours CFL FLASH-Ro e 69 . 8 ± 0 . 5 16,502 20,480 0.2 HLL3R 67 . 3 ± 0 . 5 1,920 2,297 0.8 HLL5R 67 . 9 ± 0 . 5 2,065 2,493 0.8 T able 3. Performance of the MHD schemes in driven MHD turbulence with M = 2 for a numerical resolution of 256 3 grid cells. All runs were performed in a mo de of parallel computation with 64 CPUs on the SGI Altix 4700 (HLRB I I) at the Leibniz Sup ercomputing Cen ter in Munic h. Note that the run with the FLASH-Ro e scheme was unstable for CFL n umbers > 0 . 2 in this test. the HLL5R solver typically goes through one more if-test compared to the HLL3R solver. HLL5R go es through a few more time steps, whic h is either b ecause of its slightly more restrictiv e CFL constrain t, or flo w features such as lo cal differences Alfv´ en sp eed. The total n um b er of steps and the total CPU time sp ent by t end = 3 . 3 T for each run are shown in the third and fourth column of T able 3. Given the factor of 4 b etw een the CFL num b er used with FLASH- Ro e compared to HLL3R and HLL5R it is not surprising that the total n um b er of steps is muc h higher for FLASH-Ro e . How ever, the num b er of steps for FLASH-Ro e is ab out 8.6 and 8.0 times higher compared to HLL3R and HLL5R resp ectively , which is more than the difference in the CFL num b ers. The reason for this is that FLASH-Ro e is close to b eing unstable even with CFL = 0.2. This pro duces timesteps that are often significan tly smaller than necessary . T aken altogether, for the problem setup discussed in this section, HLL3R and HLL5R are factors of 8.9 and 8.2 resp ectively faster than FLASH-Roe , which is reflected by the total amoun t of CPU time sp ent for a successful MHD turbulence run with M ≈ 2. F or higher Mach n um b er runs, how ever, FLASH-Ro e was extremely unstable, whic h did not allow us to make this comparison of the sc hemes at higher Mac h num b er. 5.6. Resolution study. As outlined in the introduction, interstellar turbulence is highly supersonic with t ypical Mach num b ers of the order of M ≈ 10. As discussed in the comparison of the different MHD sc hemes in the previous subsections, the FLASH-Ro e scheme turned out to b e highly unstable even for low Mac h n umber turbulence. This limitation is ov ercome with the HLL3R and HLL5R sc hemes. In this section, w e apply our new HLL3R scheme to MHD turbulence with M = 10 as an example of the high Mac h num b er turbulence typical of the in terstellar medium and discuss numerical resolution issues. The three runs analysed in this section w ere run at n umerical resolutions of 128 3 , 256 3 and 512 3 grid cells with a slightly smaller initial magnetic field, B 0 = 4 . 4 × 10 − 6 Gauss, which giv es an initial plasma beta of β ≈ 1. First, w e show slices through the three-dimensional b ox at y = L/ 2 of the density and magnetic pressure fields in Figure 5.4. With increasing resolution more small-scale structure is resolved. The density and magnetic pressure are correlated due to the compression of magnetic field lines in sho cks. V ector fields sho wing the v elo cit y and the magnetic field structure are ov erlaid on the slices of the densit y and magnetic pressure, resp ectively . Numerical dissipation b ecomes smaller with higher resolution. Since the inertial range of the turbulence is defined as all the length scales muc h smaller than the forcing scale and muc h larger than the dissipation scale, a minimum numerical resolution of ab out 512 3 grid cells is required to separate the inertial range from the forcing and dissipativ e scales. T o see how dissipation depends on numerical resolution w e plot the v elo cit y and magnetic field sp ectra in Figure 5.5 for the three different resolutions in analogy to Figure 5.2. The turbulence is driv en on large scales corresponding to k ≈ 2 (injection scale, L inj ), while the sp ectra b ecome steep er and b egin dropping to zero at w av en umbers k & 10, 20, and 40 for 128 3 , 256 3 , and 512 3 , indicating the onset of dissipation (on scales ` dis ), resp ectiv ely . Th us, the range of scales separating the injection scale from the dissipation scale increases with increasing resolution. How ever, the inertial range is defined as all scales, ` with L inj  `  ` dis . Purely hydrodynamical simulations of driven turbulence with up to 1024 3 grid cells show that the inertial range just b ecomes apparent for resolutions & 512 3 , whic h means that we cannot see a conv ergence of the slop e in the sp ectra shown in Figure 5.5. How ever, measuring 19 Figure 5.4. Slices of the density ( left ) and magnetic pressure ( right ) at y = L/ 2 and t = 2 T for n umerical resolutions of 128 3 ( top ), 256 3 ( midd le ), and 512 3 ( b ottom ). 20 Figure 5.5. Comparison of the F ourier sp ectra of the turbulent velocity ( top p anel ) and the turbulen t magnetic field ( b ottom p anel ) for HLL3R in M = 10, β = 1 MHD turbulence. The influence of the grid resolution is shown: 128 3 (dotted), 256 3 (dashed) and 512 3 (solid). Both the velocity and the magnetic field spectra w ere decomposed into their rotational and compressible parts, E ⊥ and E k , resp ectively , as in Figure 5.2. the p ow er-law slop e, E ∝ k β of the v elo cit y sp ectra in the w av en umber range 5 . k . 10 at 512 3 yields slop es of β ≈ − 1 . 6 ± 0 . 1, and β ≈ − 1 . 9 ± 0 . 1, for E ⊥ u and E k u , resp ectively . The solenoidal part, E ⊥ u is in go o d agreement with the Kolmogorov sp ectrum ( β = − 5 / 3) [36, 21], while the compressible part, E k u is significan tly steep er and closer to Burgers turbulence ( β = 2) [11], applicable to a sho c k-dominated medium. Note that E k m clearly decreases with increasing resolution. 5.7. T urbulen t dynamo. Magnetic fields are ubiquitous in molecular clouds, but it remains contro v ersial whether these fields hav e an influence on the cloud dynamics (see, e.g., [46]). How ever, it is widely accepted that magnetic fields pla y a significan t role on the scales of protostellar cores, where they lead to the generation of sp ectacular jets and outflows, launched from the protostellar disks, a pro cess for which a wound-up magnetic field seems to b e the key [1, 45]. Thus, w e test in this section whether our new MHD schemes can represen t w ound-up magnetic field configurations with the same turbulence-in-a-b o x approac h as in the previous subsections. The magnetic pressure can b ecome comparable to the thermal pressure in dense cores due to the am- plification of the magnetic field through first, compression of magnetic field lines, and second, due to the winding, twisting and folding of the field lines b y v orticity , a pro cess called turbulent dynamo (see [7] for a comprehensive review of turbulent dynamo action in astrophysical systems). Magnetic field amplification in the early univ erse during the formation of the first stars and galaxies was discussed analytically in [63], sho wing that the magnetic pressure can reach levels of ab out 20% of the thermal pressure in primordial mini-halos, thus potentially influencing the fragmen tation of the gas. Numerical simulations also show that 21 Figure 5.6. Shows the magnetic pressure as a function of time (for 200 eddy turnov er times, T ) for solenoidal forcing (solid line) and compressiv e forcing (dashed line). The dash-dotted horizon tal line shows the thermal pressure. The turbulen t dynamo works with an exp onential gro wth rate of ab out 0 . 60 /T for solenoidal forcing and 0 . 28 /T for compressiv e forcing (dotted lines). F or comparison, typical amplification rates found in subsonic, solenoidally driven turbulence are ab out 0 . 5 /T (e.g., [62]). the turbulent dynamo is an efficient pro cess to amplify small seed magnetic fields during the formation of the first stars and galaxies ([74, 66]). In present-da y star formation the influence of magnetic fields on the scales of dense cores w ere inv estigated in numerical work ([57, 30]), concluding that magnetic fields strongly affect fragmen tation of dense gas. Understanding the magnetic field gro wth due to the turbulent dynamo is thus crucial for future studies of star formation. Many studies hav e fo cused on subsonic turbulence (e.g., [62]), with only very few contributions on the supersonic regime. F or instance, [28] only studied the turbulen t dynamo for mildly supersonic Mach n umbers, M . 2 . 5. F or molecular clouds ho wev er, the highly supersonic regime is more relev an t. Figure 5.6 sho ws a comparison study of the turbulent dynamo operating in the sup ersonic regime ( M = 5). W e used the HLL3R for this test with a numerical resolution of 128 3 grid cells. Unlike the previous turbulence runs discussed abov e w e started from an extremely small initial magnetic field, B 0 = 4 . 4 × 10 − 16 Gauss, whic h corresp onds to an initial plasma beta of β ≈ 10 20 . W e used t wo limiting cases to drive the turbulence: purely solenoidal (div ergence-free) forcing as abov e, and purely compressive (curl-free) forcing (see [19] for a detailed analysis of the differences of solenoidal and compressive turbulence forcings and their role in the context of molecular cloud dynamics and star formation). Figure 5.6 shows that the turbulent dynamo leads to an exp onential growth of the magnetic pressure ov er more than 15 orders of magnitude within about 100 large-scale eddy turnov er times, T . The dynamo growth rate is ab out twice as large for solenoidal forcing ( ≈ 0 . 60 /T ) compared to compressive forcing ( ≈ 0 . 28 /T ), due to the higher a v erage vorticit y generated by solenoidal forcing (see [19]), whic h mak es the small-scale dynamo more efficient. The actual gro wth rates, ho wev er, dep end on the kinematic and magnetic Reynolds num b ers. Since w e did not add physical viscosity and magnetic diffusivit y , these num b ers are controlled by numerical viscosit y and diffusivit y , and thus b y n umerical resolution. A resolution study of the dynamo in recent FLASH simulations is presented in [66] and in [20], the kinematic and magnetic Reynolds num b ers found in dynamo simulations similar to the ones presen ted here are around 200 for a n umerical resolution of 128 3 grid cells. A more detailed analysis of the Mac h num b er dep endence of the turbulent dynamo amplification in solenoidal and compressive forcings is in preparation. The dynamo saturates after ab out 70 T and 140 T for solenoidal and compressiv e forcing, resp ectiv ely . The saturation level is close to the thermal pressure in b oth cases. F or compressiv e forcing the saturated magnetic pressure is about 5% of the thermal pressure, while it is 40% for solenoidal forcing. This n umerical test shows that the turbulent dynamo works with the new HLL3R scheme. This is an imp ortan t 22 test as it shows that the scheme repro duces the exp ected amplification of the magnetic pressure due to the winding-up of magnetic field lines. 6. Summar y W e presented an implementation of an accurate, efficient and highly stable numerical metho d for MHD problems. The metho d is reviewed here, and presen ted in detail in [72]. It is implemented as a mo dification of the FLASH co de [23], which enables large-scale, m ulti-pro cessor sim ulations and adaptive mesh refinement. In [72] it w as found that our metho d could handle significan tly larger ranges of the sonic Mac h n um b er and plasma β than a standard MHD sc heme. This w as confirmed in this paper by comparisons with the standard FLASH code. The algorithmic changes underlying the increased stabilit y can b e brok en down in to three parts: (1) An en tropy stable approximate Riemann solver that preserves p ositivity of density and internal energy ([4, 5]). (2) F or second order accuracy , a reconstruction metho d that ensures p ositivit y ([72]). (3) In multidimensions, a stable discretisation of the P ow ell system ([72]). All these ingredien ts w ere essential in obtaining the desired stability and efficiency of the o verall scheme. The differen t elements of the new scheme hav e b een studied separately in previous pap ers. While [72] fo cused on the p ositive second-order algorithm and m ultidimensionality , only a single appro ximate Riemann solv er, HLL3R, w as considered. The present study contrasts a standard scheme to the combination of these three new ingredients. The new sc heme was implemented in tw o versions featuring the 3-wa ve (HLL3R) and 5-wa v e (HLL5R) appro ximate Riemann solvers of [5] resp ectively , while there w ere tw o standard implementations of the FLASH co de (version 2.5), using the Ro e and the HLLE appro ximate Riemann solvers. W e observed some increase of numerical dissipation compared to the Ro e solv er of FLASH, but it was minor, and due to the replacemen t of the Ro e solver with the robust and efficient HLL-type solvers. The HLLE solver was found to be the most dissipativ e, while HLL5R show ed almost iden tical dissipation properties and accuracy to the Ro e solver. HLL3R was rank ed b et w een HLLE and Ro e in terms of accuracy . As a physical application, we ha v e considered forced MHD turbulence at high sonic Mach num b er. W e w ere able to compare the new and old sc hemes at RMS sonic Mach n umber 2 with an initial plasma β = 0 . 25. The schemes were all found to give similar and reasonable results, but the new sc hemes HLL3R and HLL5R w ere altogether ab out eigh t times more efficien t in this test. The Ro e solver-based scheme in the FLASH co de was sligh tly less dissipative, but had to be run at a four times low er CFL num ber to b e stable. At RMS sonic Mach n um b er 10 only the new schemes yielded physical results. W e found reasonable dependence of dissipation on numerical resolution at this Mach num b er, and w ere able to infer a small inertial range from the velocity pow er sp ectra at 512 3 resolution. Finally , we studied the turbulent dynamo action at RMS Mac h n um b er 5 with the new scheme. So far, there ha ve b een very few studies of turbulent dynamos in the supersonic regime. W e found dynamo-generated exponential gro wth rates of the magnetic pressure that differed according to the type of forcing mechanism, i.e., solenoidal versus compressive forcing [17, 18, 19]. T urbulence simulations with the new scheme across a wide range of mo dest to large sonic and Alfv´ enic Mach n umbers hav e been presented in [10, 9]. The tw o relaxation-based approximate Riemann solv ers HLL3R and HLL5R hav e previously not b een compared at high order and in higher spatial dimensions than 1D. In man y cases we found them to give similar results, with HLL3R b eing sligh tly more efficient. How ever, w e presen ted one case, a 2D Kelvin- Helmholtz instability , where the more detailed HLL5R w as significan tly less viscous. This was b ecause of v elo cit y shears parallel to b oth the grid and the magnetic field lines. In a three-dimensional turbulence run at Mach 2, w e also found HLL5R to b e somewhat less dissipative than HLL3R, giving results that were very close to the Ro e solv er-based FLASH version. The new FLASH MHD mo dule is freely av ailable up on contact with the corresp onding author. References [1] R. D. Blandford and D. G. Pa yne. Hydromagnetic flows from accretion discs and the pro duction of radio jets. MNRAS , 199:883–903, 1982. [2] F ran¸ cois. Bouch ut. Entropy satisfying flux vector splittings and kinetic BGK models. Numer. Math. , 94(4):623–672, 2003. 23 [3] F ran¸ cois Bouch ut. Nonlinear stability of finite volume metho ds for hyp erb olic c onservation laws and wel l-b alanc e d schemes for sourc es. F rontiers in Mathematics. Basel: Birkh¨ auser. viii, 135 p., 2004. [4] F ran¸ cois Bouch ut, Christian Klingen b erg, and Knut W aagan. A multiw a ve appro ximate Riemann solver for ideal MHD based on relaxation I - theoretical framework. Numerische Mathematik , 108(1):7–41, 2007. [5] F ran¸ cois Bouch ut, Christian Klingen b erg, and Knut W aagan. A multiw a ve appro ximate Riemann solver for ideal MHD based on relaxation II - numerical implementation with 3 and 5 wav es. Numerische Mathematik , 115(4):647–679, 2010. [6] J. U. Brackbill and D. C. Barnes. The effect of nonzero product of magnetic gradient and B on the n umerical solution of the magnetohydrodynamic equations. Journal of Computational Physics , 35:426–430, 1980. [7] A. Brandenburg and K. Subramanian. Astroph ysical magnetic fields and nonlinear dynamo theory. Phys. Rep. , 417:1–209, 2005. [8] M. Brio and C.C. W u. An up wind differencing sc heme for the equations of ideal magnetohydrodynamics. J. Comput. Phys. , 75(2):400–422, 1988. [9] C. M. Brun t, C. F ederrath, and D. J. Price. A metho d for reconstructing the PDF of a 3D turbulent density field from 2D observ ations. MNRAS , 405:L56–L60, 2010. [10] C. M. Brunt, C. F ederrath, and D. J. Price. A metho d for reconstructing the v ariance of a 3D physical field from 2D observ ations: application to turbulence in the interstellar medium. MNRAS , 403:1507–1515, 2010. [11] J. M. Burgers. A mathematical mo del illustrating the theory of turbulence. A dv. Appl. Mech. , 1:171–199, 1948. [12] P . Colella and P .R. W oo dward. The piecewise parabolic metho d (PPM) for gas-dynamical simulations. J. Comput. Phys. , 54:174–201, 1984. [13] A. J. Cunningham, A. F rank, P . V arni ` ere, S. Mitran, and T. W. Jones. Sim ulating Magnetohydrodynamical Flow with Constrained T ransp ort and Adaptive Mesh Refinement: Algorithms and T ests of the AstroBEAR Code. ApJS , 182:519–542, 2009. [14] A. Dedner, F. Kemm, D. Kr¨ oner, C.-D. Munz, T. Schnitzer, and M. W esen b erg. Hyperb olic divergence cleaning for the MHD equations. J. Comput. Phys. , 175(2):645–673, 2002. [15] B. G. Elmegreen and J. Scalo. Interstellar T urbulence I: Observ ations and Pro cesses. ARA&A , 42:211–273, 2004. [16] V. Eswaran and S. B. Pope. An examination of forcing in direct numerical simulations of turbulence. Computers and Fluids , 16:257–278, 1988. [17] C. F ederrath, R. S. Klessen, and W. Schmidt. The Densit y Probabilit y Distribution in Compressible Isothermal T urbulence: Solenoidal versus Compressive F orcing. ApJ , 688:L79–L82, 2008. [18] C. F ederrath, R. S. Klessen, and W. Schmidt. The F ractal Density Structure in Sup ersonic Isothermal T urbulence: Solenoidal V ersus Compressive Energy Injection. ApJ , 692:364–374, 2009. [19] C. F ederrath, J. Roman-Duv al, R. S. Klessen, W. Schmidt, and M.-M. Mac Lo w. Comparing the statistics of in terstellar turbulence in simulations and observations. Solenoidal versus compressive turbulence forcing. A&A , 512:A81, 2010. [20] C. F ederrath, S. Sur, D. R. G. Schleic her, R. Banerjee, and R. S. Klessen. A new Jeans resolution criterion for (M)HD simulations of self-gravitating gas: Application to magnetic field amplification by gravit y-driven turbulence. ApJ , 2011. submitted [21] Uriel F risc h. T urbulenc e . Cambridge Univ. Press, 1995. [22] S. F romang, P . Hennebelle, and R. T eyssier. A high order Godunov scheme with constrained transp ort and adaptive mesh refinement for astrophysical magnetohydrodynamics. A&A , 457:371–384, 2006. [23] B. F ryxell, K. Olson, P . Ric ker, F. X. Timmes, M. Zingale, D. Q. Lam b, P . MacNeice, R. Rosner, J. W. T ruran, and H. T ufo. Flash: An adaptive mesh h ydro dynamics code for modeling astroph ysical thermon uclear flashes. The Astrophysic al Journal Supplement Series , 131(1):273–334, 2000. [24] F. F uchs, K. H. Karlsen, S. Mishra, and N.H. Risebro. Stable upwind schemes for the magnetic induction equation. Math. Mo del. Num. Anal , 43:825–852, 2009. [25] F. F uchs, A McMurry, S. Mishra, Risebro N.H., and K. W aagan. Approximate Riemann solvers and stable high-order finite volume schemes for multi-dimensional ideal MHD. Communication in Computational Physics , 9:324–362, 2011. [26] K.F. Gurski. An HLLC-type approximate Riemann solver for ideal magnetohydrodynamics. SIAM J. Sci. Comput. , 25(6):2165–2187, 2004. [27] Amiram Harten, Peter D. Lax, and Bram v an Leer. On upstream differencing and Go dunov-t yp e sc hemes for hyperb olic conserv ation laws. SIAM R ev. , 25:35–61, 1983. [28] N. E. L. Haugen, A. Brandenburg, and A. J. Mee. Mach num b er dep endence of the onset of dynamo action. MNRAS , 353:947–952, 2004. [29] P . Henneb elle and G. Chabrier. Analytical Theory for the Initial Mass F unction: CO Clumps and Prestellar Cores. ApJ , 684:395–410, 2008. [30] P . Hennebelle and R. T eyssier. Magnetic pro cesses in a collapsing dense core. II. F ragmentation. Is there a fragmentation crisis? A&A , 477:25–34, 2008. [31] A.-K. Jappsen, R. S. Klessen, R. B. Larson, Y. Li, and M.-M. Mac Low. The stellar mass sp ectrum from non-isothermal grav oturbulent fragmentation. A&A , 435:611–623, 2005. [32] S. Kitsionas, C. F ederrath, R. S. Klessen, W. Schmidt, D. J. Price, L. J. Dursi, M. Gritschneder, S. W alch, R. Pion tek, J. Kim, A.-K. Jappsen, P . Ciecielag, and M.-M. Mac Low. Algorithmic comparisons of decaying, isothermal, supersonic turbulence. A&A , 508:541–560, 2009. [33] R. S. Klessen. One-Poin t Probability Distribution F unctions of Sup ersonic T urbulen t Flows in Self-gra vitating Media. ApJ , 535:869–886, 2000. 24 [34] C. Klingenberg, W. Schmidt, and K. W aagan. Numerical comparison of Riemann solv ers for astroph ysical hydrodynamics. Journal of Computational Physics , 227:12–35, 2007. [35] Christian Klingen berg and Kn ut W aagan. Relaxation solvers for ideal MHD equations -a review. A cta Mathematic a Sci- entia , 30(2):621 – 632, 2010. [36] A. N. Kolmogorov. Dissipation of energy in lo cally isotropic turbulence. Dokl. A kad. Nauk SSSR , 32:16–18, 1941. [37] M. R. Krumholz and C. F. McKee. A General Theory of T urbulence-regulated Star F ormation, from Spirals to Ultralumi- nous Infrared Galaxies. ApJ , 630:250–268, 2005. [38] D. Lee and A. E. Deane. An unsplit staggered mesh scheme for multidimensional magnetohydrodynamics. Journal of Computational Physics , 228:952–975, 2009. [39] D. Lee, A. E. Deane, and C. F ederrath. A New Multidimensional Unsplit MHD Solver in FLASH3. In N. V . P ogorelov, E. Audit, P . Colella, & G. P . Zank, editor, Astr onomic al Society of the Pacific Confer enc e Series , volume 406 of Astr o- nomic al So ciety of the Pacific Confer enc e Series , page 243, 2009. [40] M. N. Lemaster and J. M. Stone. Density Probability Distribution F unctions in Sup ersonic Hydro dynamic and MHD T urbulence. ApJ , 682:L97–L100, 2008. [41] R.J. LeV eque. Numerical Metho ds for Conservation L aws . Birkh¨ auser, Basel, Switzerland; Boston, U.S.A., 2nd edition, 1992. [42] Shengtai Li. An HLLC Riemann solver for magneto-hydrodynamics. J. Comput. Phys. , 203(1):344–357, 2005. [43] Y. Li, R. S. Klessen, and M.-M. Mac Lo w. The F ormation of Stellar Clusters in T urbulent Molecular Clouds: Effects of the Equation of State. ApJ , 592:975–985, 2003. [44] T. J. Linde. A thr ee-dimensional adaptive multifluid MHD model of the heliospher e . PhD thesis, UNIVERSITY OF MICHIGAN, 1998. [45] D. Lynden-Bell. On why discs generate magnetic tow ers and collimate jets. MNRAS , 341:1360–1372, 2003. [46] M.-M. Mac Lo w and R. S. Klessen. Control of star formation by sup ersonic turbulence. R eviews of Modern Physics , 76:125–194, 2004. [47] Barry Marder. A method for incorporating Gauss’ la w in to electromagnetic pic co des. J. Comput. Phys. , 68(1):48–55, 1987. [48] C. F. McKee and E. C. Ostriker. Theory of Star F ormation. ARA&A , 45:565–687, 2007. [49] A. Mignone, G. Bo do, S. Massaglia, T. Matsakos, O. T esileanu, C. Zanni, and A. F errari. PLUTO: a Numerical Code for Computational Astrophysics. In JENAM-2007, ”Our non-stable Universe”, held 20-25 August 2007 in Y er evan, Armenia. Abstr act bo ok, p. 96-96 , pages 96–96, 2007. [50] T ak ahiro Miy oshi and Kany a Kusano. A multi-state HLL appro ximate Riemann solver for ideal magnetoh ydrodynamics. J. Comput. Phys. , 208(1):315–344, 2005. [51] E. C. Ostriker, J. M. Stone, and C. F. Gammie. Density , V elo cit y , and Magnetic Field Structure in T urbulent Molecular Cloud Models. ApJ , 546:980–1005, 2001. [52] P . P adoan and ˚ A. Nordlund. The Stellar Initial Mass F unction from T urbulent F ragmentation. ApJ , 576:870–879, 2002. [53] P . Padoan, ˚ A. Nordlund, and B. J. T. Jones. The univ ersality of the stellar initial mass function. MNRAS , 288:145–152, 1997. [54] T. Passot and E. V´ azquez-Semadeni. Densit y probability distribution in one-dimensional p olytropic gas dynamics. Phys. Rev. E , 58:4501–4510, 1998. [55] Benoit Perthame and Chi-W ang Sh u. On p ositivity preserving finite v olume schemes for Euler equations. Numer. Math. , 73(1):119–130, 1996. [56] Kenneth G. Pow ell. An approximate Riemann solver for magnetohydrodynamics (that works in more than one dimension). T echnical report, Institute for Computer Applications in Science and Engineering (ICASE), 1994. [57] D. J. Price and M. R. Bate. The impact of magnetic fields on single and binary star formation. MNRAS , 377:77–90, 2007. [58] D. J. Price and C. F ederrath. A comparison betw een grid and particle metho ds on the statistics of driv en, sup ersonic, isothermal turbulence. MNRAS , 406:1659–1674, 2010. [59] P . L. Ro e. Approximate Riemann Solv ers, Parameter V ectors, and Difference Schemes. Journal of Computational Physics , 43:357, 1981. [60] R Samtaney , P Colella, T J Ligo cki, D F Martin, and S C Jardin. An adaptive mesh semi-implicit conservativ e unsplit method for resistive mhd. Journal of Physics: Confer enc e Series , 16(1):40, 2005. [61] J. Scalo and B. G. Elmegreen. Interstellar T urbulence I I: Implications and Effects. ARA&A , 42:275–316, 2004. [62] A. A. Sc hekochihin, S. C. Cowley, S. F. T aylor, J. L. Maron, and J. C. McWilliams. Sim ulations of the Small-Scale T urbulent Dynamo. ApJ , 612:276–307, 2004. [63] D. R. G. Schleic her, R. Banerjee, S. Sur, T. G. Arshakian, R. S. Klessen, R. Beck, and M. Spaans. Small-scale dynamo action during the formation of the first stars and galaxies. I. The ideal MHD limit. A&A , 522:A115, 2010. [64] W. Schmidt, C. F ederrath, M. Hupp, S. Kern, and J. C. Niemeyer. Numerical sim ulations of compressiv ely driven interstellar turbulence: I. Isothermal gas. A&A , 494:127, 2009. [65] J. M. Stone, T. A. Gardiner, P . T eub en, J. F. Hawley, and J. B. Simon. Athena: A New Code for Astrophysical MHD. ApJS , 178:137–177, 2008. [66] S. Sur, D. R. G. Schleic her, R. Banerjee, C. F ederrath, and R. S. Klessen. The Generation of Strong Magnetic Fields During the F ormation of the First Stars. ApJ , 721:L134–L138, 2010. [67] Eleuterio F. T oro. R iemann solvers and numeric al methods for fluid dynamics. A pr actic al intr oduction. 2nd e d. Berlin: Springer. xix, 624 p. , 1999. 25 [68] G. T oth. A General Co de for Mo deling MilD Flows on Parallel Computers: V ersatileAdvection Co de. In Y. Uc hida, T. Kosugi, & H. S. Hudson, editor, IAU Col lo q. 153: Magnetodynamic Phenomena in the Solar A tmospher e - Prototyp es of Stel lar Magnetic A ctivity , page 471, 1996. [69] G. T´ oth. The ∇· B = 0 Constrain t in Shock-Capturing Magnetohydrodynamics Co des. Journal of Computational Physics , 161:605–652, 2000. [70] Bram v an Leer. On the relation b etw een the upwind-differencing sc hemes of Go dunov, Engquist–Osher and Ro e. SIAM Journal on Scientific and Statistic al Computing , 5(1):1–20, 1984. [71] E. V´ azquez-Semadeni. Hierarchical Structure in Nearly Pressureless Flows as a Consequence of Self-similar Statistics. ApJ , 423:681, 1994. [72] K. W aagan. A positive MUSCL-Hancock scheme for ideal magnetohydrodynamics. J. Comput. Phys. , 228(23):8609–8626, 2009. [73] H. Xu, D. C. Collins, M. L. Norman, S. Li, and H. Li. A Cosmological AMR MHD Module for Enzo. In B. W. O’Shea & A. Heger, editor, First Stars III , volume 990 of Americ an Institute of Physics Conferenc e Series , pages 36–38, 2008. [74] H. Xu, H. Li, D. C. Collins, S. Li, and M. L. Norman. T urbulence and Dynamo in Galaxy Cluster Medium: Implications on the Origin of Cluster Magnetic Fields. ApJ , 698:L14–L17, 2009. [75] U. Ziegler. Self-gravitational adaptive mesh magnetohydrodynamics with the NIR V ANA co de. A&A , 435:385–395, 2005. (Knut W aagan) Center for Scientific Comput a tion and Ma thema tical Modeling (CSCAMM) University of Mar yland, CSIC Building 406, College P ark, MD 20742-3289, USA E-mail address : kwaagan@cscamm.umd.edu (Christoph F ederrath) Zentrum f ¨ ur Astronomie der Universit ¨ at Heidelberg, Institut f ¨ ur Theoretische Astrophysik, Alber t-Ueberle-Str. 2, D-69120 Heidelberg, Germany E-mail address : chfeder@ita.uni-heidelberg.de (Christoph F ederrath, current address) Ecole Normale Sup ´ erieure de L yon, CRAL, 69364 L yon, France (Christian Klingenberg) Dep ar tment of Ma thema tics, W ¨ urzburg University, Am Hubland 97074 W ¨ urzburg, Germany E-mail address : klingenberg@mathematik.uni-wuerzburg.de 26