RAM: A Relativistic Adaptive Mesh Refinement Hydrodynamics Code

1.6 MB·

by  Weiqun Zhang

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: [e-mail] Adaptive Mesh Refinement (AMR) has played an in-
Electronicaddress: [e-mail] 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