Draftversion February2,2008 PreprinttypesetusingLATEXstyleemulateapjv.6/22/04 RAM: A RELATIVISTIC ADAPTIVE MESH REFINEMENT HYDRODYNAMICS CODE Weiqun Zhang1 KavliInstitute forParticleAstrophysicsandCosmology,StanfordUniversity,P.O.Box20450, MS29,Stanford,CA94309 and Andrew I. MacFadyen InstituteforAdvancedStudy, Princeton,NJ08540 Draft versionFebruary 2, 2008 ABSTRACT 6 We have developed a new computer code, RAM, to solve the conservative equations of special rel- 0 ativistic hydrodynamics (SRHD) using adaptive mesh refinement (AMR) on parallel computers. We 0 haveimplementedacharacteristic-wise,finitedifference,weightedessentiallynon-oscillatory(WENO) 2 scheme using the full characteristic decomposition of the SRHD equations to achieve fifth-order ac- curacy in space. For time integration we use the method of lines with a third-order total variation n diminishing (TVD) Runge-Kutta scheme. We have also implemented fourth and fifth order Runge- a J Kutta time integration schemes for comparison. The implementation of AMR and parallelization is basedontheFLASHcode. RAMismodularandincludesthecapabilitytoeasilyswaphydrodynamics 0 1 solvers, reconstruction methods and physics modules. In addition to WENO we have implemented a finitevolumemodulewiththepiecewiseparabolicmethod(PPM)forreconstructionandthemodified 2 Marquina approximate Riemann solver to work with TVD Runge-Kutta time integration. We exam- v ine the difficulty of accurately simulating shear flows in numerical relativistic hydrodynamics codes. 1 We show that under-resolved simulations of simple test problems with transverse velocity compo- 8 nentsproduceincorrectresultsanddemonstratetheabilityofRAMtocorrectlysolvetheseproblems. 4 RAM has been tested in one, two and three dimensions and in Cartesian, cylindrical and spherical 5 coordinates. We have demonstrated fifth-order accuracy for WENO in one and two dimensions and 0 performed detailed comparisonwith other schemes for which we show significantly lowerconvergence 5 rates. Extensive testing is presented demonstrating the ability of RAM to address challenging open 0 questions in relativistic astrophysics. / h Subject headings: hydrodynamics – methods: numerical – relativity p - o 1. INTRODUCTION tivisticnumericalsimulationsduringthepastdecadedue r t Many astrophysical phenomena involve gas moving at to the development of modern numerical methods and s the increasing power of computers. We refer the reader a relativisticspeeds. Classicalsourcesincludeactivegalac- : ticnuclei(AGN),microquasars,pulsarwindnebulaeand to a comprehensive review of numerical relativistic hy- v drodynamics by Mart´ı & Mu¨ller (1999) and references gamma-raybursts(GRBs). Morerecently,conclusiveev- i X idence has mounted that a subset ofcore collapse super- therein. The so-called high-resolution shock-capturing (HRSC) methods, which are based on the fact that r novae are associated with long duration (τ &5 s) GRBs a (e.g.,Hjorth et al.2003; Stanek et al.2003). Astrophys- special relativistic hydrodynamics (SRHD) with causal equation of state is a hyperbolic system of conservation ical processes at the endpoint of stellar evolution thus laws (Anile 1989), are particularly promising in modern appear capable of accelerating flows to ultra-relativistic numerical simulations. They can achieve very high ac- speed (W &100, where W denotes Lorentz factor.) curacy and handle ultra-relativistic flows, strong shocks Additionally, the fading afterglows of cosmological and contact discontinuities extremely well. GENESIS GRBs are now observed with an increasing rate by the (Aloy et al. 1999) is a very efficient scheme based on SWIFT satellite (Gehrels et al. 2004). GRB afterglows HRSC techniques. are produced after the gamma-ray producing relativis- The high-order essentially non-oscillatory (ENO) tic flow transfers its energy to the circum-burst medium shock capturing schemes (Harten et al. 1987) have been in the form of a strong relativistic shock. Of particular verysuccessfulinnumericallysolvinghyperbolicsystems interest is evidence that the ejecta producing GRBs and e.g., the Eulerian gas dynamics equations. ENO-based theirafterglowsarebeamedintojets. Thequalityofcur- methodshavepreviouslybeenimplementedforcomputa- rentafterglowobservationsrequirehigh-resolutionmulti- tionalrelativistichydrodynamics(Dolezal & Wong1995; dimensional simulations of jetted relativistic shocks for Donat et al. 1998; Del Zanna & Bucciantini 2002). The interpretation. weighted essentially non-oscillatory (WENO) schemes Relativistic hydrodynamics simulations are more dif- (e.g., Liu, Osher & Chan 1994; Jiang & Shu 1996) are ficult than Newtonian simulations (Norman & Winkler recent developments following the philosophy of ENO 1986). However, there has been much progress in rela- schemes. Electronicaddress: [email protected] Adaptive Mesh Refinement (AMR) has played an in- Electronicaddress: [email protected] creasinglyimportantrole inmanybranchesofnumerical 1 ChandraFellow astrophysics including e.g., cosmology, star formation, 2 Zhang & MacFadyen jets, interacting binaries, stellar wind collisions, and su- and stress-energy of a fluid: pernova explosion and remnant evolution (see Norman (ρuµ) =0, (1) 2004, and references therein). Examples of recent work ;µ in numerical astrophysics include using AMR to solve and the equations of hydrodynamics (Plewa & Mu¨ller 2001), (Tµν) =0 (2) ;ν magnetohydrodynamics (MHD) (Balsara 2001), special where ρ is the rest mass density measured in the fluid relativistichydrodynamics(SRHD)(Hughes et al.2002), frame, uµ = W(c,u) is the fluid four-velocity, W is the general relativistic hydrodynamics (D¨onmez 2004) and Lorentz factor, c is the speed of light, u is the classical general relativistic MHD (Anninos et al. 2005). three-velocity,Tµν isthestress-energytensorofthefluid In this paper, we describe a new modular, highly andthesubscript denotesthecovariantderivative. For accurate, special relativistic hydrodynamics code with ;µ adaptive mesh refinement, RAM (Relativistic Adaptive a perfect fluid the stress-energy tensor is Mesh). The modular design of the code allows us to Tµν =ρhuµuν +pgµν, (3) easily change between various algorithms. We have im- whereh 1+ǫ+p/ρis therelativisticspecificenthalpy, plementedthefifth-orderWENOschemeofJiang & Shu ≡ ǫ is the specific internal energy, p is the pressure and (1996) in SRHD for the first time as one ofthe modules. gµν is the inverse metric which here we take to be the We have also implemented an HRSC module using the Minkowskimetricthoughextensiontocurvedspacetimes modified Marquina flux formula (Marquina et al. 1992; is evident. Aloy et al. 1999) and piecewise parabolic reconstruction SRHDcanbewrittenasasetofconservationlaws(see, (Colella & Woodward 1984) similar to the GENESIS e.g. Mart´ı & Mu¨ller 1999), code (Aloy et al. 1999). We also use piecewise linear re- construction (PLM) with both WENO and HRSC and ∂U 3 ∂Fj perform comparisons between all methods for several + =0, (4) standardtestcases. Wepresentdetailedtestsofthecode ∂t ∂xj j=1 demonstrating its ability to handle the extreme resolu- X where the conserved variable U is given by tion requirements of simulating ultra-relativistic flows. A primary motivation in writing the RAM code is to U=(D,S1,S2,S3,τ)T, (5) study the relativistic explosions producing cosmological GRBs. We plan to use RAM to simulate the relativistic and the fluxes are given by flows produced at the endpoint of stellar evolution. In Fj =(Dvj,S1vj+pδj ,S2vj+pδj ,S3vj+pδj ,Sj Dvj)T, particular, numerical simulations of the collapsar model 1 2 3 − (6) for GRBs (Woosley 1993; MacFadyen & Woosley 1999) here δµ is the Kronecker symbol and vj is the velocity. are very challenging because of the ultra-relativistic ν The conservedvariables U include rest mass density, D, speed and detailed micro-physics involved in the prob- the momentum density, Sj, and energy density, τ. They lem. Relativistic hydrodynamical simulations of jets in are measured in the laboratory frame, and are given by collapsars have been done successfully by several au- (assuming the speed of light c=1), thors (Aloy et al. 2000; Zhang et al. 2003, 2004) using the GENESIS method. We plan to use RAM to extend D=ρW (7) these calculations to higher resolution over a larger dy- Sj=ρhW2vj (8) namic range in lengthscale made possible by AMR. In 2 we present the fundamental equations of SRHD τ=ρhW2 p ρW, (9) § − − anddescribetheconservedvariablesevolvedbyourcode. where j = 1,2,3, W is the Lorentz factor. The equa- In 3 we describe the flux differenced semi-discrete § tions 4 are closed by an equation of state (EOS) given form of the equations we solve, the algorithms used for by p=p(ρ,ǫ). For an ideal gas, the EOS reads, construction of numerical fluxes and how the equations are integrated in time. In 4 we present results from p=(Γ 1)ρǫ, (10) § − standard test problems on a fixed uniform mesh includ- where Γ is the adiabatic index of the ideal gas. ing convergence tests. In 5 we describe the adap- § tive mesh refinement algorithm employed in RAM by 3. NUMERICAL SCHEMES utilizing components of the FLASH code (Fryxell et al. To numerically solve Eq. 4, each spatial dimension is 2000), which in turn adapted the PARAMESH package discretizedintocells. Usingthemethodoflines,thetime (MacNeice et al. 2000). In 6 we present test problems § dependent evolution of Eq. 4 can be expressed in the runwithAMRinone,two,andthreedimensionsinclud- semi-discrete form ingRiemannproblemswithtransversevelocity,a2Dax- isymmetric jet in cylindrical coordinates and a 3D blast dU Fx Fx i,j,k = i+1/2,j,k− i−1/2,j,k (11) wave. In 7 we summarize the results of the paper. In dt − ∆x § Appendix A we providedetails ofour implementationof Fy Fy SRHD in curvilinear coordinates. i,j+1/2,k− i,j−1/2,k − ∆y Fz Fz 2. GOVERNINGEQUATIONSOFSPECIAL i,j,k+1/2− i,j,k−1/2, RELATIVISTICHYDRODYNAMICS − ∆z The governing equations of special relativistic hydro- where Fx , Fy and Fz are the fluxes i±1/2,j,k i,j±1/2,k i,j,k±1/2 dynamics(SRHD)describetheconservationofrestmass at the cell interface. RAM Code 3 In order to achieve high-order accuracy in time, the approximateRiemann solver. The cell unknowns Ui are timeintegrationisdoneusingahigh-ordertotalvariation considered as cell averages. Thus, the second class can diminishing (TVD) Runge-Kutta scheme (Shu & Osher be considered as finite volume schemes. We use U-X to 1988),whichcombinesthefirst-orderforwardEulersteps denote the second class of schemes, where X stands for andinvolvespredictionandcorrection. Forexample,the the reconstruction scheme. third-order accuracy can be achieved via 3.1. F-X Schemes: Reconstruction Of Fluxes U(1)=Un+∆tL(Un) (12) IntheF-Xschemes,wehaveimplementedtheweighted 3 1 1 U(2)= Un+ U(1)+ ∆tL(U(1)) (13) essentially non-oscillatory (WENO) scheme and piece- 4 4 4 wise linear method (PLM) for the reconstruction of Un+1=1Un+ 2U(2)+ 2∆tL(U(2)), (14) fluxes. 3 3 3 Essentially non-oscillatory (ENO) finite difference where L(U) is the right hand side of Eq. 12, Un+1 is schemes (Harten et al. 1987) for equations of hyperbolic the final value after advancing one time step from Un. conservationlawsuseadaptivestencilsincalculatingthe Besides the third-order Runge-Kutta method (RK3), fluxes across the cell interfaces so that high-order ac- thestandardfourth-orderRunge-Kutta(RK4)andfifth- curacy can be achieved and numerical oscillations near order Runge-Kutta method (RK5) (Lambert 1991) can discontinuities can be significantly reduced. Later to- also be used for time integration. We usually use TVD tal variation diminishing (TVD) Runge-Kutta methods RK3inourcalculationsunlessotherwisestated. Wehave were applied to ENO schemes to make the schemes alsoimplementedandtestedRK4andRK5andusethem morecomputationallyefficient,andtheENOreconstruc- in some calculations. tions were carried out on fluxes instead of cell averages Note the method of lines, treats information fromcor- (Shu & Osher1988, 1989). ENO schemeshavebeen fur- ner zones differently than methods using time-averaged ther improved by many researchers (see Shu 1997, for fluxes which make some use of corner information a review). One improved scheme is the fourth-order when implemented with Strang splitting. Recently, weighted essentially non-oscillatory (WENO) scheme dimensionally unsplit methods have been developed whichwasfirstintroducedbyLiu, Osher & Chan(1994). (Mignone, Plewa & Bodo 2005; Gardiner & Stone 2005) In the WENO schemes, a weighted combination of sev- whichusethecorner-transportupwind(CTU)methodof eralpossible stencils are used instead of just one stencil. (Colella 1990) to include corner information completely. This would improve the accuracy while keeping the es- RAM, however, uses the third-order Runge-Kutta al- sentiallynon-oscillatorypropertyneardiscontinuities. A gorithmfortimeintegrationconsistingofseveralforward new algorithm for computing the weights has been used Eulersubsteps whicharecorrectedto achievehighaccu- in a modified fifth-order WENO scheme developed by racy. Since the final results after one full time step are Jiang & Shu(1996). TheWENOschemeofJiang & Shu computedthroughseveralsubsteps,thecornercells(e.g., (1996) has been applied by Feng et al. (2004) to cosmo- (x ,y ,z )) of a zone at (x ,y ,z ) affect its evo- logical simulations. In this paper, we have implemented i+1 j+1 k+1 i j k lution in a full step, although this is not the case for an the fifth-order modified WENO scheme of Jiang & Shu individual substep ( 3.1 & 3.2). In this sense, RAM, (1996)inspecialrelativistichydrodynamics. TheWENO is dimensionally uns§plit. Co§mparison with codes using schemes for hyperbolic conservation laws have various CTUmethodsonmulti-dimensionaltestproblemswould forms(Shu1997). The particularWENO schemewe use beofvaluetoperforminthefuture. Wenotethatmulti- is the characteristic-wise flux splitting finite difference dimensional test problems we have run with RAM pre- scheme. serve symmetries present in the initial conditions (see, IntheWENOscheme,thepartialdifferentialequations e.g., Fig. 8). are discretized into cells in spatial dimensions (Eq. 12). To solve the fluxes across the cell interfaces, we have The evolution of the partial differential equations is implemented two schemes. Both schemes involve a re- solved by using a TVD Runge-Kutta scheme. The key construction step to compute the variables at the cell componentoftheWENOschemeistocomputethefluxes interfaces. In the first class of schemes ( 3.1), the re- at the cell interfaces. As an example, we shall consider construction is carried out on fluxes. The§interface flux the x-direction flux across the cell interface between xi is obtained by, and xi+1. We shall assume that states at xi−2, xi−1, x , and x are either in the computational domain F =F(F ,...,F ), (15) i+2 i+3 i+1/2 i−r i+s or can be supplied by boundary conditions. where the stencil (i r,i r + 1,...,i + s) depends First,weconstructaRoetypeaveragestateatthecell − − on the choice of the reconstruction scheme. This class interface x from the left state at x and right state i+1/2 i can be considered finite difference schemes. We use at x . The eigenvectors of the average state will be i+1 F-X to denote this class of schemes, where X stands used as the basis for the decomposition of fluxes, F , m for the reconstruction scheme. In the other class of where m = i 2,i 1,i,i + 1,i + 2,i + 3. We use − − schemes ( 3.2), the reconstruction is carried out on √ρh as weight (Eulderink & Mellema 1995) to compute § the unknown variables. Then the interface unknowns, the weighted averagedof pressure, density, and velocity. U−,+ =U(U ,...,U ),areusedtoobtaininterface Strictly speaking, this does not satisfy all of the Roe’s i+1/2 i−r i+s conditions (Roe 1981). However, we do not use the av- fluxes, F =R(U− ,U+ ), (16) eragestate tocompute the fluxdirectly. Thereforestrict i+1/2 i+1/2 i+1/2 Roe’sconditionsarenotrequired. Indeed, simply taking where and + denote the left and right side of the in- arithmetic averagesof the primitive variables also works − terface i+1/2, respectively. The flux function R is an very well. 4 Zhang & MacFadyen For special relativistic hydrodynamics equations, the In computing the WENO weights, a parameter, ǫ, is complete characteristic structure of these hyperbolic introduced to avoid denominator being zero. This is the equations have been derived by Donat et al. (1998). only free parameter in the WENO scheme. Moreover, Thus it is feasible to extend the characteristic-wise the results are insensitive to the value of ǫ as long as it WENO to special relativistic hydrodynamics knowing is a small number about 10−6. the eigenvalues and eigenvectors of the Jacobian matri- We use F-WENO to denote the above method, which ces of the equations. The eigenvalues of the Jacobian uses the WENO scheme for the characteristic-wise re- matrices, constructionofsplittedfluxes. Thereconstructionofthe ∂F(U) coefficientsc fromc canalsobecomputedusing B= , (17) i+1/2,n mn ∂U otherhigh-orderreconstructionalgorithmsinsteadofthe WENO reconstructionalgorithmof Jiang & Shu(1996). namely the speed of characteristic waves, can have dif- Besidesthe WENO,we havealsoimplemented the PLM ferent signs. In other words, these waves could prop- for the reconstruction with a generalized minmod slope agate towards different directions. The second step of limiter (Kurganov & Tadmor 2000). Given c , c , and the WENOschemeis tosplittheflux intoleft-goingand i−1 i c , the left-biased interface value reads, right-going fluxes so that we can treat them separately. i+1 This flux splitting approach can avoid entropy violation c =c +0.5minmod(θ(c c ), (23) i+1/2 i i i−1 and make the scheme more robust. We use the local − 0.5(c c ),θ(c c )), Lax-Freiderichs splitting given by, i+1− i−1 i+1− i F+ =F+αU (18) where 1≤θ ≤2, and the minmod function reads, F− =F αU, (19) minmod(x,y,z)= 1(sgn(x)+sgn(y)) (24) − 4 where α is the local maximum of the absolute values of (sgn(x)+sgn(z))min(x, y , z ), wavespeedfor states atx . The eigen- | | | | | | i−2,i−1,i,i+1,i+2,i+3 valuesareallpositivefortheright-goingfluxF+,whereas here the sgn function returns the sign of the number. they are negative for the left-going flux F−. This becomes the more diffusive normal minmod lim- Because of the obvious symmetry between the left- iter when θ = 1, and becomes the monotonized central- goingfluxF− andright-goingfluxF+,weshallonlycon- difference limiter of van Leer (1977) when θ = 2. We siderheretheright-goingwaves. Thestencilfortheright- usually use θ = 1.5 unless otherwise stated. We call goingwavesati+1/2,whichconsistsofx this method F-PLM. Our numerical tests show that the i−2,i−1,i,i+1,i+2 in the fifth-order WENO scheme, is slightly upwind- results of F-PLM are comparable to those of F-WENO biased. We can expand F+, where m = i 2,i ( 4). m − − § 1,i,i+ 1,i + 2, in terms of the five right eigenvectors (Donat et al. 1998) of the average state at i+1/2, R , 3.2. U-X Schemes: Reconstruction Of Unknowns n where n=1,2,3,4,5. The expansion is given by, Besides the F-X schemes ( 3.1), we have also im- § plemented another method of computing the flux at 5 F+ = c R , (20) the interface to work with the Runge-Kutta scheme m mn n for time integration of Eq. 12. In this class of n=1 X schemes, U-X, instead of reconstructing flux directly where,m=i 2,i 1,i,i+1,i+2. Itisverystraightfor- as in the F-X schemes, a left state and a right state − − ward to compute cmn since we also know the left eigen- at the interface are reconstructed, and then the flux vectors, Ln, of the Jacobian matrix. cmn is given by, at the interface is calculated using these two states. cmn =Ln·F+m. (21) Areslooltutoiofnrescheonctkm-ceatphtoudrisn,gin(cHluRdSinCg)thmeetshoo-cdasllelidkehigthhe- Given cmn, we can construct the coefficients ci+1/2,n GENESIS method (Aloy et al. 1999) and some high- forthe cellinterfaceatxi+1/2 bygivingdifferentweights resolution central schemes (Lucas-Serrano et al. 2004; to different cmn. The weights depend upon the smooth- Del Zanna & Bucciantini 2002), canbe consideredto be ness of the stencils, and smooth stencils are given more inthiscategory. Theirmaindifferenceishowtocompute weight. Thisresultsinthefifth-orderaccuracyinsmooth the flux at the interface given the left and right states region and ENO near discontinuities. For the details of and how to reconstruct interface states. In the HRSC howtheweightsarechosenintheWENOreconstruction methods, an approximate Riemann solver utilizing the scheme, we refer the readers to the work of Jiang & Shu information of the characteristic waves is used, whereas (1996). thecentralschemesusesimplifiedexpressionsfortheflux Using the coefficients ci+1/2,n, where n denotes char- without using the characteristic information beyond the acteristic waves,the right-goingflux at the cell interface fastestwavespeed,whichisneverthelessrequiredforthe at xi+1/2 is given by, Courant-Friedrich-Levy(CFL) condition. The reconstruction of states at the cell interfaces 5 F+ = c R . (22) is usually performed by a high-order interpolation i+1/2 i+1/2,n n method like the piecewise parabolic method (PPM) n=1 X (Colella & Woodward 1984) or piecewise linear method Because of the obvious symmetry, the left-going flux, (PLM).Wehaveimplementedbothreconstructionmeth- F− , can be computed in the similar procedure. Note ods in the RAM code. For the parameters in the i+1/2 that the left-going flux is also upwind-biased, namely PPM algorithm, we choose the values suggested by right-biased. Mart´ı & Mu¨ller(1996)foralmostallthetestswithPPM. RAM Code 5 We considerit undesirable to fine-tune these parameters is employed and a very small time step is used, the nu- to achieve better agreement with analytic solutions for merical simulation will stop. One possible remedy when specific tests. Doing so, we believe, can create a false the code runs out of methods is to set the pressure of senseofhowagivencodeperformsingeneralsimulations bad cells to a floor value. Fortunately this has never for which the exact solution is not previously known. happened in our simulations. For the PLM algorithm, we use a generalized minmod In each time step, conserved variables U = slope limiter as describedin 3.1. Again, the parameter (D,S1,S2,S3,τ)T are updated directly. In order to cal- § θ = 1.5 is used by default. The reconstruction is per- culatethe fluxesatthecellinterfaces,physicalvariables: formedontheso-calledprimitivevariablesinsteadofcon- pressure, proper rest mass density and velocity need to served variables because unphysical conserved variables be computed. The processes of recovering physical vari- may arise otherwise. The pressure and the proper den- ablesfromconservedvariablesinvolveaquarticequation sity are reconstructed directly, whereas velocities are re- for velocity (e.g., Duncan & Hughes 1994). Though the constructedusingacombinationofreconstructingthree- quartic can be solved analytically, computing the ana- velocityandLorentzfactor. Sinceboththevelocitiesand lytic solution is actually more expensive than a simple Lorentz factor are reconstructed, they are unlikely to be numerical root finder, such as Newton-Raphson itera- self consistent. To make them consistent, four-velocities tion. Before the Newton-Raphson iteration, a certain are derived by multiplying three-velocities with Lorentz condition, (τ +D)2 >D2+S2, is checked to make sure factor. Using these four-velocities, new consistent veloc- that physical results can be produced from given con- itiesandLorentzfactorcouldbe derived. Wefoundthat served variables. Sometimes, unphysical results, such as this procedure is usually more robust than using only negative pressure and large velocity v > 1, will be pro- three-velocities or only four-velocities for the construc- duced. Then the fall-back mechanism will be used until tion of velocity and Lorentz factor. physical results are produced or the maximum allowed To compute the flux across the cell interface given the number of fall-back is reached. reconstructed right and left states at the interface, we 4. TEST PROBLEMS haveimplemented severalRiemannsolversincluding the modified Marquina’s flux (Aloy et al. 1999), local Lax- For any numerical code, it is very important to do Friedrichs flux and relativistic HLLE (Schneider et al. substantial tests. We have done a series of tests with 1993) in the RAM code. A comparison of several our code. Some of the tests are the so-called Riemann schemes for computing the flux has been performed by problems, which consist of computing the decay of an Lucas-Serrano et al. (2004). They have shown that the initial discontinuity of two constant states. It is possible results are relatively independent of the flux schemes. to get exact solutions for relativistic Riemann problems For the results shown in this paper, the modified Mar- (Pons et al. 2000). We have also done some extensively quina’salgorithmisusedtocomputetheinterfacefluxes. studied problems which have no analytic solutions. In allcases,ourresultsarecomparabletopublishedresults. 3.3. Failsafe Time Integration and Root Finder In this section, we will present our numerical results for Since ourcode is explicit,the time stepfor integration somestandardtests. Wewillcomparefourschemes( 3): § is subject to the Courant-Friedrich-Levy (CFL) condi- the finite difference characteristic-wise WENO, the fi- tion. We usually choose a CFL number to be 0.2-0.5. nitedifferencecharacteristic-wisePLM,thefinitevolume Unfortunately, it is not always failsafe. Occasionally component-wisePPMandthe finite volume component- unphysical results can be produced for ultra-relativistic wise PLM, denoted by F-WENO, F-PLM, U-PPM, and flows unless a very small CFL number or a diffusive U-PLM, respectively. A CFL number of 0.5 is used for scheme is used. But a very small CFL number is some- these tests unless otherwise stated. timesunacceptablebecauseofthetimeitcosts. Tomake 4.1. One-Dimensional Riemann Problem 1 the code robust, a fall-back mechanism is employed in our code, though we found that such failures rarely oc- In this test, the one-dimensional numerical region cur. That is a normal CFL number is used for the first (0 x 1) initially consists of two constant states: try of time integrations. If it fails, the code will return pL ≤= 13≤.33, ρL = 10.0, vL = 0.0 and pR = 10−8, to the beginning of the failed step and use more diffu- ρR =1.0,vR =0.0,whereLstandsfortheleftstate,and siveschemesforreconstructionand/oruseasmallertime R the right state. The fluid is assumed to be an ideal step. The recalculation with more diffusive schemes is gas with an adiabatic index Γ=5/3. The initial discon- performed on the whole grid, when AMR is not being tinuity is at x = 0.5. The results are shown in Fig. 1. used. The fall-back method on AMR grids will be fur- In this test problem, the evolution of the initial discon- ther discussed in 5. The fall-back requires us to keep tinuity gives rise to a shock, a rarefaction wave, and a a copy of the orig§inal states of conserved variables for contactdiscontinuity. Thisisafairlyeasytest. Allmod- a period of the whole Runge-Kutta step in order to be ernspecialrelativistic hydrodynamicscodes cancapture able to fall back. Fortunately, the Runge-Kutta scheme the expected features, acquire correct positions of the weuse(Shu & Osher1988)requirestheoriginalstatesto shock front, contact discontinuity and rarefaction wave. be saved anyway (Eqns. 12, 13, & 14). Thus it does not To quantitatively measure the errors, we calculated the increase the memory usage of the computation. In the L1 norm errors, L1 = j∆xj|uj −u(xj)|, where u(xj) subsequent steps, the normal reconstructionscheme will is the exact solution at x and u is the numerical re- j j P be used and the CFL number will gradually increase to sult. The L norm errors of density for four schemes 1 its normal value if it was previously decreased. The fall- with various grid resolutions at t=0.4 are shown in Ta- back mechanism makes the code almost failsafe. If un- ble1. Theaccuracyofourresultsiscomparabletothatof physicalresultspersistevenwhenthe first-ordermethod Lucas-Serranoet al.(2004);Mart´ı & Mu¨ller(1996). The 6 Zhang & MacFadyen 1.0 (a) (b) 1.0 p/20 v (a) p/1000 v (b) rho/10 rho/10 0.8 0.8 0.6 0.6 p/20 p/20 0.4 0.4 rho/10 rho/10 0.2 0.2 v v 0.0 0.0 1.0 (c) (d) 1.0 p/1000 v (c) p/1000 v (d) rho/10 rho/10 0.8 0.8 0.6 0.6 p/20 p/20 0.4 0.4 rho/10 rho/10 0.2 0.2 v v 0.0 0.0 0. 0 0.2 0.4 0.6 0.8 1. 0 0.2 0.4 0.6 0.8 1. 0 0 .0 0.2 0.4 0.6 0.8 1 .0 0.2 0.4 0.6 0.8 1 .0 Fig. 1.— One-dimensional Riemann problem 1 at t = 0.4. Fig. 2.— One-dimensional Riemann problem 2 at t = 0.4. Results for four schemes: (a) F-WENO; (b) F-PLM; (c) U-PPM; Results for four schemes: (a) F-WENO; (b) F-PLM; (c) U-PPM; (d) U-PLM are shown. The computational grid consists of 400 (d) U-PLM are shown. The computational grid consists of 400 zones. Numericalresultsareshowninsymbols,whereastheexact zones. Numericalresultsareshowninsymbols,whereastheexact solution is shown in solid lines. We show proper mass density solution is shown in solid lines. We show proper mass density (square),pressure(triangle)andvelocity(plussign). (square),pressure(triangle)andvelocity(plussign). TABLE 1 Fig. 2. In this test problem, the evolution of the ini- L1 errorsofthedensity forthe1DRiemann tial discontinuity gives rise to a right-moving shock, a Problem 1. Four schemeswith various left-moving rarefaction wave, and a contact discontinu- resolutions with uniformspacingareshown att=0.4. ity in between. Behind the shock, there is an extremely thin dense shell. The width of the shell is only 0.01056 Scheme Na L1 Error Convergence Rate at t = 0.4. With 400 uniform zones for 0 x 1, ≤ ≤ the thin shell is only covered by 4.2 zones. Due to the F-WENO 100 1.31e-1 200 7.25e-2 0.85 smearing at the contact discontinuity and the shock, it 400 3.32e-2 1.1 is not surprising that the thin shell is not well resolved 800 2.08e-2 0.67 in our results with only 400 zones. At this resolution 1600 1.00e-2 1.1 themaximumdensityinthe shellforF-WENO,F-PLM, 3200 5.07e-3 0.98 F-PLM 100 1.47e-1 U-PPM and U-PLM is 79%, 69%, 72% and 63% of the 200 8.50e-2 0.79 analytic value. However, there is no difficulty in resolv- 400 4.06e-2 1.1 ing the thin shell with increased resolution. We have 800 2.33e-2 0.80 calculated the L norm errors of density for various nu- 1600 1.22e-2 0.93 1 3200 7.48e-3 0.71 mericalresolutions(seeTable2)andachieveconvergence U-PPM 100 1.27e-1 asexpectedforproblemswithsharpdiscontinuities. The 200 7.30e-2 0.80 results are consistent with other published results (e.g., 400 3.47e-2 1.1 Lucas-Serranoet al. 2004). 800 1.97e-2 0.82 1600 9.77e-3 1.0 3200 5.10e-3 0.94 4.3. One-Dimensional Riemann Problem 3 U-PLM 100 1.32e-1 200 8.57e-2 0.62 Inthistest,theone-dimensionalnumericalregion(0 ≤ 400 3.86e-2 1.2 x 1)initiallyconsistsoftwoconstantstates: p =1.0, L 800 2.27e-2 0.77 ρ ≤= 1.0, v = 0.9 and p = 10.0, ρ = 1.0, v = 0.0, 1600 1.15e-2 0.98 L L R R R where L stands for the left state, and R the right state. 3200 6.48e-3 0.83 Thefluidisassumedtobeanidealgaswithanadiabatic aNumberofgridpoints index Γ = 4/3. The initial discontinuity is at x = 0.5. TheresultsareshowninFig.3. Inthisproblemastrong reverse shock forms in which post-shock oscillations are order of the convergence rate is about 1. This is consis- visible for the U-PPM and U-PLM simulations, espe- tent with the fact that there are discontinuities in the ciallyinthepressureprofile. IntheF-WENOsimulation, solution. the post-shock pressure oscillations is smaller those in theU-Xsimulations. Inthe F-PLMsimulationthepost- 4.2. One-Dimensional Riemann Problem 2 shock pressure oscillation is nearly invisible to the eye, In this test, the one-dimensional numerical region thoughtheyarestillpresentatthe0.1%level. Reducing (0 x 1) initially consists of two constant states: theCFLnumberfrom0.5decreasesthepost-shockoscil- p ≤= 10≤00.0, ρ = 1.0, v = 0.0 and p = 10−2, lations but does not eliminate them completely. Table 3 L L L R ρ = 1.0, v = 0.0, where L stands for the left state, presents the L norm errors and L order convergence R R 1 1 and R the right state. The fluid is assumed to be an rates for this problem. Note that the U-PPM scheme ideal gas with an adiabatic index Γ = 5/3. The initial is more accurate and converges faster for this problem discontinuity is at x = 0.5. The results are shown in thanother schemes. This is due to the factthat U-PPM RAM Code 7 TABLE 2 TABLE 3 L1 errorsofthedensity forthe1DRiemann L1 errorsofthedensity forthe1DRiemann Problem 2. Four schemeswith various Problem3. Four schemeswith various resolutions with uniformspacingareshown resolutions areshown att=0.4. att=0.4. Scheme Na L1 Error Convergence Rate Scheme Na L1 Error Convergence Rate F-WENO 100 9.97e-2 F-WENO 100 2.10e-1 200 6.29e-2 0.67 200 1.42e-1 0.56 400 3.01e-2 1.1 400 9.29e-2 0.61 800 1.69e-2 0.83 800 5.54e-2 0.75 1600 9.48e-3 0.83 1600 2.54e-2 1.1 3200 5.24e-3 0.86 3200 1.51e-2 0.75 F-PLM 100 1.12e-1 F-PLM 100 1.96e-1 200 6.98e-2 0.68 200 1.42e-1 0.46 400 3.45e-2 1.0 400 1.06e-1 0.42 800 1.94e-2 0.83 800 7.21e-2 0.56 1600 1.13e-2 0.78 1600 3.92e-2 0.88 3200 6.54e-3 0.79 3200 2.44e-2 0.68 U-PPM 100 9.72e-2 U-PPM 100 2.18e-1 200 5.60e-2 0.80 200 1.52e-1 0.52 400 2.49e-2 1.2 400 9.52e-2 0.68 800 1.30e-2 0.94 800 5.42e-2 0.81 1600 6.06e-3 1.1 1600 2.67e-2 1.0 3200 3.11e-3 0.96 3200 1.67e-2 0.68 U-PLM 100 9.53e-2 U-PLM 100 2.13e-1 200 6.32e-2 0.59 200 1.65e-1 0.37 400 2.99e-2 1.1 400 1.25e-1 0.41 800 1.78e-2 0.75 800 8.68e-2 0.53 1600 1.04e-2 0.78 1600 4.49e-2 0.95 3200 6.10e-3 0.77 3200 2.71e-2 0.73 aNumberofgridpoints aNumberofgridpoints 1.0 v (a) v (b) lemswithvelocitycomponentstransversetothedirection 0.8 p/25 p/25 of propagation of the main flow. In Newtonian hydro- 0.6 dynamics the transverse momentum is simply advected 0.4 with the flow and is not coupled directly to the equa- 0.2 rho/10 rho/10 tion of motion in the longitudinal direction. No serious difficulty is introduced by the presence of transverse ve- 0.0 locity components though transverse kinetic energy dis- 1.0 v (c) v (d) sipatedthroughviscousdissipationcanalterthelongitu- 0.8 p/25 p/25 dinalmotionby affecting the pressure. However,inrela- 0.6 tivistic flow, transverse velocities are directly coupled to 0.4 the dynamics along all directions by the Lorentz factor 0.2 rho/10 rho/10 which depends on all velocity components. This cou- pling makes relativistic Riemann problems with trans- 0.0 verse velocity much more difficult to solve correctly in 0. 0 0.2 0.4 0.6 0.8 1. 0 0.2 0.4 0.6 0.8 1. 0 a numerical code. In particular, higher resolution is Fig. 3.— One-dimensional Riemannproblem3at t=0.4. Re- needed. Under-resolved simulations produce incorrect sultsforfourschemes: (a)F-WENO;(b)F-PLM;(c)U-PPM;(d) shock positions along the normal direction. Here, we U-PLMareshown. Thecomputational gridconsists of400zones. Numerical results are shown in symbols, whereas the exact solu- presentrelativelyeasy1Dtests ofRiemannproblemsin- tionisshowninsolidlines. Weshowpropermassdensity(square), cluding transverse velocity which can be resolved with pressure(triangle)andvelocity(plussign). Thisproblemcontains moderate resolution. In 6.1 we show tests requiring astrongreverseshockinwhichpost-shockpressureoscillationsare § high resolution and demonstrate the ability of adaptive visibleinpanels(c)and(d). mesh refinement to accurately simulate very relativistic flows with significant transverse velocity components. In this test, the one-dimensional numerical region schemecapturesthe contactdiscontinuityin fewerzones (0 x 1) initially consists of two constant states: (3). p ≤= 10≤00.0, ρ = 1.0, v = 0.0, v = 0.0 and L L xL yL p = 10−2, ρ = 1.0, v = 0.0, v = 0.99, where 4.4. One-Dimensional Riemann Problem With R R xR yR L stands for the left state, and R the right state. The Non-Zero Transverse Velocity: Easy Test fluid is assumed to be an ideal gas with an adiabatic in- Many problems of interest in hydrodynamics involve dex Γ = 5/3. The initial discontinuity is at x = 0.5. strong shear flows. Astrophysical jets include shearing The results are shown in Fig. 4. In this test, the pres- layerswhere mixing of ambientmaterialinto the fast jet ence of transverse velocity alters the structure of the flow may be important. It is therefore important to test Riemann problem (Pons et al. 2000; Rezzolla & Zanotti the ability of numerical codes to handle Riemann prob- 2002). Thisproblemisrelativelyeasybecausethetrans- 8 Zhang & MacFadyen 1.0 p/1000 (a) p/1000 (b) 1.0 (a) (b) v v v v 0.8 0.8 p/8e9 p/8e9 0.6 0.6 0.4 rho/25 rho/25 0.4 rho/5e5 rho/5e5 0.2 0.2 0.0 0.0 1.0 p/1000 (c) p/1000 (d) 1.0 (c) (d) v v v v 0.8 0.8 p/8e9 p/8e9 0.6 0.6 0.4 rho/25 rho/25 0.4 rho/5e5 rho/5e5 0.2 0.2 0.0 0.0 0. 0 0.2 0.4 0.6 0.8 1. 0 0.2 0.4 0.6 0.8 1. 0 0 .0 0.2 0.4 0.6 0.8 1 .0 0.2 0.4 0.6 0.8 1 .0 Fig. 4.— “Easy” one-dimensional Riemann problem with non- Fig. 5.— One-dimensional shock heating problem in planar zero transverse velocity (vy = 0.99c) at t = 0.4. Results for four geometry at t=2.0. Results for four schemes: (a) F-WENO; (b) schemes: (a) F-WENO; (b) F-PLM; (c) U-PPM; (d) U-PLM are F-PLM; (c) U-PPM; (d) U-PLM are shown. The parameter θ in shown. The computational grid consists of 400 zones. Numerical the minmod slope limiter for U-PLM is set to be 1.0 in this test. resultsareshowninsymbols,whereas theexactsolutionisshown The computational grid consists of 100 zones. Numerical results in solid lines. We show proper mass density (square), pressure areshowninsymbols,whereastheexactsolutionisshowninsolid (triangle)andvelocityinx-direction(plussign). lines. We show proper mass density (square), pressure (triangle) andvelocity(plussign). TABLE 4 L1 errorsof thedensity forthe“easy”1D The shock heatingproblemis astandardtest to study Riemannproblemwith non-zerotransverse the ability of a code to handle very strong shocks with velocity att=0.4. Four schemeswith sufficiently few zones and without excessive post shock variousresolutions using auniformgridare shown. oscillations. Coldgasflowstowardareflectingboundary at x = 1.0 and a reverse strong shock forms and propa- Scheme Na L1 Error Convergence Rate gates to the left decelerating the gas to zero speed. The gashasaninitialspeedofv =0.9999999999=1.0 10−10 F-WENO 100 7.58e-1 − 200 3.92e-1 0.95 (corresponding to a Lorentz factor of W = 70710.675) 400 2.31e-1 0.76 andinitialdensityofρ=1.0. Due toconservationofen- 800 1.18e-1 0.97 ergyandthefactthegasinitiallyhasnearlyzerointernal 1600 6.58e-2 0.84 energycomparedto kinetic energy (we use a smallvalue 3200 3.44e-2 0.94 F-PLM 100 8.26e-1 ǫ0 = 0.003 for numerical reasons), the specific internal 200 4.59e-1 0.85 energy after the shock is simply ǫ = W 1. The com- 400 2.77e-1 0.73 pressionratioacrosstheshockcanbearb−itrarilylargein 800 1.49e-1 0.89 the relativistic case and is given by 1600 8.00e-2 0.90 3200 4.63e-2 0.79 Γ+1 Γ U-PPM 100 8.48e-1 σ = + ǫ (25) 200 4.25e-1 1.0 Γ 1 Γ 1 400 2.41e-1 0.82 − − 800 1.27e-1 0.92 where Γ is the adiabatic index of the gas. In this case 1600 6.43e-2 0.99 Γ=4/3soσ =282845.70andtheshockvelocityisgiven 3200 3.34e-2 0.95 by v = (Γ−1)W|v| = 0.33332862. In Fig. 5 we show our U-PLM 100 9.00e-1 s W+1 200 4.72e-1 0.93 results for a uniform mesh of 100 zones compared with 400 2.88e-1 0.71 the analytic solution at t=2.0. 800 1.52e-1 0.92 We note that for this problem with a constant state 1600 8.86e-2 0.78 3200 4.95e-2 0.84 behindtheshock,asopposedtoathinshellwhichmight aNumberofgridpoints more naturally occur in a realistic flow, the maximum Lorentz factor that can be achieved numerically is just limited by floating point precision. The parameter θ in the minmod slope limiter for U- verse velocity is in the cold gas of the right state, not in PLMissettobe1.0inthistesttoeliminatestrongoscil- the relativistically hotleft state or in the rarefactionfan lationswhichwouldappearifthedefaultvalueθ =1.5is which subsequently propagates into it. The simulation usedinU-PLM.Inthis testwithU-PPM,strongerdissi- agrees well with the analytic solution even at the rela- pationis alsoneededto avoidcrashes. Inparticular,one tively modest resolution of 400 zones. Table 4 presents of the PPM parameters, ǫ(2), is set to be 10−4 instead the L norm errors and L order convergence rates for 1 1 of the default value ǫ(2) = 1.0 (see Colella & Woodward this problem. 1984, for the meaning of this parameter). The reflecting wall at x = 1.0 poses difficulties for numerical simula- 4.5. One-Dimensional Shock Heating Problem In tions and gives rise to visible errors for zones near the Planar Geometry wall. The errors of density at the nearest zone to the RAM Code 9 reflecting wall are 3.9%, 2.4%, 7.4%, and 4.3%, for F- 350 WENO,F-PLM,U-PPM,andU-PLM,respectively. For 300 this particular problem with very strong shock but sim- 250 plestructure,morediffusiveschemesF-PLMandU-PLM P 200 perform better than F-WENO and U-PPM. 150 4.6. Isentropic Smooth Flows 100 2.2 Most tests presented in this paper involve strong 2.0 shocks. In addition to tests with discontinuities, it is 1.8 alsoveryimportantto test the capability ofthe schemes o 1.6 h to handle smooth flows. In this section, we present two r 1.4 tests involving the isentropic evolution of smooth flows 1.2 in 1D and 2D. 1.0 0.6 In the first test, a one-dimensional isentropic smooth 0.5 pulse is set up in a uniform reference state, and the run 0.4 stopssometimebeforeashockisformed. Thetestissim- 0.3 ilar to the convergence tests performed by Colella et al. v 0.2 (2006) for their Newtonian hydrodynamics code. The 0.1 initial density structure at t=0 is given by 0.0 ρ (x)=ρ (1+αf(x)), (26) -0.2 0.0 0.2 0.4 0.6 0.8 1 .0 0 ref x where ρ is the density of the reference state and ref Fig. 6.— One-dimensional isentropic flow. The initial pulse is shownontheleft,andthepulseatt=0.8ontheright. Pressure, ((x/L)2 1)4 : x <L f(x)= − | | (27) density and velocity are shown in the top, middle and bottom 0 : otherwise, panelsrespectively. Thesolidlinesareexact solutions. Numerical (cid:26) resultsofF-WENOatt=0(plussigns)andt=0.8(triangle)are here α is the amplitude of the pulse and L the width shown. of the pulse. The pressure is given by the isentropic relation, p = KρΓ, where K is a constant. The initial velocity of the reference state is vref = 0. The initial that the F-WENO scheme, which is fifth-order accurate velocity inside the pulse at t = 0 is set up by assuming in space and third-order accurate in time, has an L or- 1 one of the two Riemann invariants, der convergence rate of 3 4 for this one-dimensional ∼ − test with smooth flow. For the F-WENO-RK4 and F- 1 1+v 1 √Γ 1+c J = ln( ) ln( − s), (28) WENO-RK5 schemes, which are fourth and fifth-order − 2 1 v − √Γ 1 √Γ 1 c accurate in time integration respectively, order of con- − − − − s vergence rates up to 5 can be achieved. Note that isconstantacrossthewholeregion,wherecs isthesound the difference between∼the results of RK4 and RK5 is speed. We note that the other Riemann invariant, very small. For schemes other than WENO, including 1 1+v 1 √Γ 1+c schemes using RK4 or RK5, the order of convergence s J+ = ln( )+ ln( − ) (29) rate is about 2 (Table 5). These results clearly show 2 1 v √Γ 1 √Γ 1 c − − − − s that the fifth-order WENO scheme is very accurate and is not constant. The exact solution of the test can be converges very quickly for smooth flows. obtained by using standard characteristic analysis. The PPM is generally 3rd-order accurate in space, but the pulse will have a smooth shape before a shock eventu- flattening procedure and the requirement of monotonic ally forms. The width and height of the pulse does not profilesdegradestheaccuracyoftheschemeatplaceslike change before the shock forms. But the pulse will be- localextremaorwherethesecondderivativesofvariables come increasingly asymmetric as the shape of the front are large. Thus its accuracy can be as low as first-order ofthepulsebecomessteeperduringthepropagation(see in some places. Moreover, in the U-PPM scheme, the Fig. 6). Behind the moving pulse, the fluid goes back to reconstruction is carried out on primitive variables, not the reference state. conservative variables. The conservative variables are Our computational region for this test is 0.35 x cellaverages. Buttheprimitivevariablesconvertedfrom − ≤ ≤ 1. The reference state is, p = 100,ρ = 1,v = 0, conservativevariablesarenotexactlycellaverages. How- ref ref ref the amplitude of the pulse is α = 1.0, and the width is ever,the PPMreconstructionalgorithmregardsthemas L= 0.3. The adiabatic index in the equation of state is cell averages. Therefore the interface values of primitive Γ=5/3. variables may not have third-order accuracy. Thus it is We have run this test with different schemes and var- reasonable that the U-PPM scheme does not achieve a ious numerical resolution. Because the flow is very third-order convergence rate. smooth, all methods perform very well for this test. Inthesecondtest,wehaveperformedtwo-dimensional Fig. 6 shows the exact solution and the numerical re- calculationstoassesstheconvergencerateoftheWENO sults of F-WENO with 80 uniform zones. The results schemeinmulti-dimension. Thecomputationalregionin of the convergence rates are shown in Table 5. Various this test consists of a two-dimensional box in Cartesian Runge-Kutta methods, including the third-order TVD coordinates with 0.0 x 3.75 and 0.0 y 5.0. The ≤ ≤ ≤ ≤ scheme, standard fourth-order scheme (RK4) and fifth- boundary conditions are periodic for allfour sides of the orderscheme(RK5)havebeenusedinthistest. Wefind box. Like the one-dimensional test of isentropic flows, 10 Zhang & MacFadyen TABLE 5 TABLE 6 L1 errorsof thedensity forthe 1Disentropic L1 errors ofthe density forthe2Disentropicflow flowproblematt=0.8. Results with various problematt=2.4. Results with variousresolutions resolutions using a uniformgridareshown. usinga uniformgridareshown. Schemea Nb L1 Error Convergence Rate Schemea Nb L1 Error Convergence Rate F-WENO 80 2.07e-3 F-WENO 48×64 7.35e-2 160 1.10e-4 4.2 96×128 4.43e-3 4.1 320 1.70e-5 2.7 192×256 8.04e-4 2.5 640 1.47e-6 3.5 384×512 9.62e-5 3.1 1280 1.58e-7 3.2 768×1024 1.12e-5 3.1 2560 1.91e-8 3.1 F-WENO-RK4 48×64 7.24e-2 5120 2.37e-9 3.0 96×128 4.75e-3 3.9 F-WENO-RK4 80 1.87e-3 192×256 4.70e-4 3.3 160 1.18e-4 4.0 384×512 3.18e-5 3.9 320 1.31e-5 3.2 768×1024 1.24e-6 4.7 640 6.80e-7 4.3 F-WENO-RK5 48×64 7.19e-2 1280 2.54e-8 4.7 96×128 4.67e-3 3.9 2560 7.91e-10 5.0 192×256 4.61e-4 3.3 5120 2.38e-11 5.1 384×512 3.13e-5 3.9 F-WENO-RK5 80 1.87e-3 768×1024 1.22e-6 4.7 160 1.17e-4 4.0 aThethird,fourth(RK4),andfifth-order(RK5)Runge-Kutta 320 1.30e-5 3.2 methodsareusedwiththefifth-orderWENOscheme. 640 6.82e-7 4.3 1280 2.54e-8 4.7 bNumberofgridpointsinxandy-direction 2560 8.01e-10 5.0 5120 2.40e-11 5.1 F-PLM 80 8.79e-3 there is a static uniform reference state, which is set to 160 4.05e-3 1.1 320 1.22e-3 1.7 pref = 100,ρref = 1,vref = 0. Pulses which are periodic 640 3.10e-4 2.0 in space are added to the system. Along the direction 1280 7.83e-5 2.0 of k = (4/5,3/5), the profile is periodic with a spatial 2560 1.96e-5 2.0 5120 4.92e-6 2.0 period of S = 3.0, and the profile is constant along the F-PLM-RK4 80 8.85e-3 directionperpendiculartothevectork. Thusthespatial 160 4.06e-3 1.1 periods along the x and y-direction are 3.75 and 5.0, 320 1.21e-3 1.7 respectively. Notethatthese areconsistentwiththe size 640 3.11e-4 2.0 1280 7.84e-5 2.0 ofthe computationalbox withperiodic boundaries. The 2560 1.97e-5 2.0 pulses move along the direction of the vector k. The 5120 4.93e-6 2.0 initial density profile is given by ρ (d) (Eqs. 26 & 27), 0 U-PPM 80 1.11e-2 where d, the distance to the center of the nearest pulse, 160 2.47e-3 2.2 320 7.02e-4 1.8 is givenby d=mod(k r+S/2,S) S/2,here mod(a,b) 640 1.38e-4 2.3 returns the reminder o·f the division−a/b, and r=(x,y). 1280 2.92e-5 2.2 The amplitude of the pulse is α = 1.0, and the width is 2560 6.48e-6 2.2 L = 0.9. The adiabatic index in the equation of state 5120 1.52e-6 2.1 U-PPM-RK4 80 1.10e-2 is Γ = 5/3. Similar to the one-dimensional test, the 160 2.56e-3 2.1 initial pressure is given by the isentropic relation, and 320 5.74e-4 2.2 the initial velocity by assuming the Riemann invariant, 640 1.34e-4 2.1 J is constant. 1280 3.10e-5 2.1 − 2560 7.33e-6 2.1 Wehaverunthetwo-dimensionaltestusingtheWENO 5120 1.82e-6 2.1 scheme with three Runge-Kutta methods: RK3, RK4 U-PLM 80 1.12e-2 andRK5. ResultsoftheL normerrorsandconvergence 160 3.56e-3 1.7 1 rateareshowninTable6. Similartotheone-dimensional 320 1.03e-3 1.8 640 2.61e-4 2.0 test, both F-WENO-RK4 and F-WENO-RK5 performs 1280 6.50e-5 2.0 better than F-WENO, which uses RK3. But the dif- 2560 1.62e-5 2.0 ference of errorsbetween F-WENO-RK4and F-WENO- 5120 4.03e-6 2.0 RK5isverysmall. Therefore,forthistestitisnotworth U-PLM-RK4 80 1.12e-2 160 3.56e-3 1.7 using the more expensive RK5 for time integration. 320 1.02e-3 1.8 640 2.60e-4 2.0 4.7. Two-Dimensional Tests: Wind Tunnel With Step 1280 6.49e-5 2.0 2560 1.62e-5 2.0 Inordertotesttheabilityofourcodetohandlestrong 5120 4.04e-6 2.0 shocks in multiple dimensions, we have performed stan- aRK4 and RK5 denote the fourth and fifth-order dard tests published in the literature. Runge-Kutta methods, respectively. The third-order The Emery step (Emery 1968; Woodward & Colella Runge-Kutta method (RK3) is used unless otherwise 1984) consists of a horizontal wind flowing into a step, stated. bNumberofgridpoints representedas a reflectingboundarycondition. The cor- ner of the step represents a singular point of the rar- efactionfan. Asthewindcollideswiththestepareverse shockpropagatesbackintothewindformingabowshock