A MUSTA-FORCE algorithm for solving partial differential equations of relativistic hydrodynamics
J. Porter-Sobieraj, M. Słodkowski, D. Kikoła, J. Sikorski, P. Aszklar
AA MUSTA-FORCE algorithm for solvingpartial differential equations of relativistichydrodynamics
J. Porter-Sobieraj ,M. Słodkowski , D. Kikoła ,J. Sikorski P. Aszklar Warsaw University of Technology, Faculty of Mathematics and InformationScience, Koszykowa 75, 00-662 Warsaw, Poland Warsaw University of Technology, Faculty of Physics, Koszykowa 75, 00-662Warsaw, Poland University of Warsaw, Faculty of Physics, Hoża 69, 00-681 Warsaw, Poland,
Abstract
Understanding event-by-event correlations and fluctuations is cru-cial for the comprehension of the dynamics of heavy ion collisions. Rela-tivistic hydrodynamics is an elegant tool for modeling these phenomena;however, such simulations are time-consuming, and conventional CPUcalculations are not suitable for event-by-event calculations. This workpresents a feasibility study of a new hydrodynamic code that employsgraphics processing units together with a general MUSTA-FORCEalgorithm (Multi-Stage Riemann Algorithm - First Order Centeredscheme) to deliver a high-performance yet universal tool for event-by-event hydrodynamic simulations. We also investigate the performanceof selected slope limiters that reduce the amount of numeric oscillationsand diffusion in the presence of strong discontinuities and shock waves.The numerical results are compared to the exact solutions to assessthe code’s accuracy.
Keywords: relativistic hydrodynamic, simulation of heavy ion collisions,quark-gluon plasma, high energy nuclear physics, numerical algorithms,MUSTA-FORCE, parallel computing, CUDA/GPU1 a r X i v : . [ phy s i c s . c o m p - ph ] F e b Motivation
Relativistic hydrodynamics simulations are widely used in the modeling ofnuclear processes in high-energy nuclear physics when examining the propertiesof the quark gluon plasma (QGP). Detailed information regarding the reactionsthat take place at the microscopic level is not required. Hydrodynamic modelformalism treats QGP as a perfect fluid and assumes a single equation ofstate. The bulk and hot nuclear system can be described using hydrodynamicconservation laws and then solved numerically [1, 2]. Hydrodynamic modelsare extremely successful in describing experimental results for particles withlow transverse momentum which is a behavior of bulk nuclear matter. On theother hand, jets (narrow spays of hadrons and other particles produced bythe hadronization of a high energy quark or gluon) are widely used to probethe properties of the QGP. This is an approach analogous to tomography: anexternal, penetrating probe, whose properties (like a production mechanism)are under experimental and theoretical control, is shot through the medium.We can then infer the properties of the analyzed system from the modificationof the probe’s energy. Jets are such probes - external to the QGP. Becausetheir production requires a large momentum transfer, they are produced veryearly in the collision, in the initial hard interaction, before the QGP phase.Their production is described well by perturbative quantum chromodynamics(pQCD) calculations. Thus they are excellent tools for QGP research.The mechanism involved in energy loss due to interactions with nuclearmatter has been a topic of extensive theoretical and experimental studies overthe last two decades. However, the energy dissipated by jets can also alterthe properties of bulk nuclear matter (for instance, so-called elliptic flow) inthe intermediate transverse momentum range. There is little understandingof such effects. For such studies, we need to efficiently model the soft particleevolution with a high spatial resolution to capture the jet-induced modificationof the characteristics of the bulk nuclear matter. Moreover, the Cartesiancoordinate system is preferred to ensure a high spatial resolution that isconstant throughout the evolution of the system. Such calculations arenecessary to fully understand the proprieties of this unique state of nuclearmatter.The main motivation of our work is therefore to develop a new applica-tion to study the physics of heavy ion collisions, enabling the execution ofhydrodynamic simulations on high-resolution Cartesian grids. The goal is tocombine two areas of heavy ion physics that are usually treated separately:bulk physics, described by hydrodynamic models, and the physics of jets.Such an approach will allow us to investigate in detail how energy depositedby jets affects the behavior of soft particles. In turn, it will help to understand2etter the properties of nuclear matter under extreme conditions.Moreover, event-by-event correlations and fluctuations play a significantrole in the understanding of heavy ion collision dynamics. Unfortunately,hydrodynamic simulations are time- and computing-power-consuming, andevent-by-event calculations are challenging when conventional CPU computingis used. The achievements of hydrodynamic models in high-energy physicsmotivated us to work on a new high-performance program that facilitates event-by-event simulation with high numerical precision using graphics processingunits.Equations of ideal relativistic hydrodynamics in one spatial dimension canbe solved both analytically and numerically as shown e.g. in [3, 4, 5]. On athree dimensional grid, however, an analytical solution is difficult to achieveand the problem needs to be addressed via numerical algorithms. A number ofsuch algorithms already exist solving relativistic hydrodynamic equations forideal and viscous fluid models either directly [6, 7, 8, 9, 10] or retrieving thesolution from a particle-based kinematic model [11]. These solutions, however,incur a heavy computational cost in 3+1 dimensions. Hence, to ensurereasonable performance, the simulation is often done in reduced dimensionality,or one of the dimensions (usually rapidity) is represented in low resolution,by only a few layers. This approach is motivated by certain symmetries thatareapproximately held in heavy ion collisions. Full (3+1)-dimensional, highresolution simulations are currently very time-consuming computations fortraditional, single-threaded software.This problem can be solved by employing general purpose computing ongraphics processing units (GPGPU). The paper [12] has shown that a GPGPU-based implementation of the SHASTA algorithm [13] can provide up to twoorders of magnitude improvement in performance. In addition, our previouswork, based on the universal multi-stage (MUSTA-FORCE) algorithm [14,15] also proved that hydrodynamic computations can be performed withina reasonable timeframe. Although the implemented method has a highnumerical cost and complexity, it scales perfectly with parallel computations.Our implementation turned out to be over 200 times faster than a sequentialimplementation on a CPU [16]. It produces good statistics and high spacialand temporal resolutions at the same time.The speed-up gained by using a GPU made it possible to develop a valuablenew tool for high energy nuclear science. Our current GPU implementationallows a device to perform a massive number of simulations, and furthermore,the MUSTA-FORCE algorithm is a very universal tool. Its strength lies in thefact that it uses simple central schemes and does not require any knowledgeof the physical process’s details. However, in the case of large gradients it isnecessary to apply a slope limiter to reduce numerical oscillations. Therefore,3e focused here on how to improve the accuracy of the MUSTA-FORCEalgorithm. The hydrodynamic simulation results presented in this paper showthat the MUSTA-FORCE method is very sensitive to the choice of slope andslope limiter, and the way that it is applied. We used such a procedure foreach dimension separately, but in general it could be applied in other ways(e.g. a common slope limiting value could be chosen for all the variables).There is no general procedure in a multi-dimensional and non-scalar case,and always some experimentation is necessary for both a particular system ofequations, and perhaps even a particular problem. Under the conditions of oursimulation the schemes must be especially sensitive to problems containingboth strong discontinuities and smooth solution features.The paper is organized as follows. The MUSTA-FORCE algorithm, to-gether with the slope limiters, is described in detail in Section 2. Section 3presents the results of using such slope limiters for standard benchmarks innuclear physics and compares them with known analytical solutions. Wediscuss their ability to achieve high order accuracy in smooth regions whilemaintaining stable, non-oscillatory and sharp discontinuity transitions. Thelast section concludes our paper.
Relativistic hydrodynamics simulations are based on hyperbolic partial differ-ential equations in the form: ∂U∂t + ∂F ( U ) ∂x + ∂G ( U ) ∂y + ∂H ( U ) ∂z = 0 (1)and an equation of state: p = p ( e, n ) (2) U = ( E, M x , M y , M z , R ) is a vector of conserved quantities in the laboratoryrest frame ; E is the energy density, M x , M y and M z are the momentumdensities in the x , y and z Cartesian coordinates, respectively, and R isa conserved charge density. Vectors of fluxes F , G , H in the x , y , and z F ( U ) = ( E + p ) v x M x v x + pM y v x M z v x Rv x G ( U ) = ( E + p ) v y M x v y M y v y + pM z v y Rv y H ( U ) = ( E + p ) v z M x v z M y v z M z v z + pRv z (3)where v is the velocity, and p is pressure, defined by the energy e and chargedensity n in the fluid rest frame , where velocity v vanishes ( v = (0 , , E = ( ε + p ) γ − pM i = ( ε + p ) γ v i , i = x, y, z (4) R = nγ where ε is the energy density and γ = √ − v is the Lorentz factor. Eq. 4defines the transformation from rest frame variables to conserved variablesused in integration. To start hydrodynamic evolution, an initial state is required as input.The most basic is parametrizations based e.g. on Glauber-like modelsin the transverse plane (see [17] for a review), and Bjorken’s solution in thelongitudinal direction.Other approaches involve models based on color glass condensate (CGC),which describe a Lorentz contracted and slowed down, fast moving particle;pQCD+saturation model [18], or the string rope model [19].5hese models describe a smooth, averaged initial state. However, sincethe hydrodynamic equations are nonlinear, a solution with an averaged initialstate is not equivalent to the average of solutions with fluctuating initialconditions. Because of this, event-by-event calculations became a major pointof interest.Fluctuating initial conditions can be obtained using e.g. Monte–CarloGlauber [20, 21, 22] or CGC, SPheRIO, NeXus [23], NeXSPheRIO [24], andmodels like EPOS [25] or UrQMD [26].
For time propagation the standard Runge–Kutta methods are employed [29].For numerical stability only total variation diminishing methods are used.In general, a Runge–Kutta method for Eq. 1 can be written in the form: U n (0) = U n U n ( i ) = i − X k =0 ( α ik U n ( k ) + ∆ tβ ik L ( U n ( k ) )) ,i = 1 , . . . , m (5) U n +1 = U n ( m ) where the upper index without parentheses denotes the time step, the lowerindex denotes integration step, L is a numerical recipe to calculate the negativeflux gradient in Eq. 1 and α, β are constant coefficients given for a particularmethod.For second order accuracy the following method is used: U n (1) = U n + ∆ tL ( U n ) U n +1 = 12 ( U n + U n (1) + ∆ tL ( U n (1) )) (6)and for third order accuracy: U n (1) = U n + ∆ tL ( U n ) U n (2) = 34 U n + 14 U n (1) + 14 ∆ tL ( U n (1) ) (7) U n +1 = 13 U n + 23 U n (2) + 23 ∆ tL ( U n (2) )It is apparent that apart from the additional computational cost due to moreevaluations of L , these methods introduce the need for an additional storageregister for each of the conserved variables. As this can be an issue for largeresolution simulations, a low storage version of the third order method canbe used. 6 .4 Hybrid MUSTA-FORCE Algorithm To obtain a general and accurate solution for Eq. 1 and Eq. 2, we use ahybrid MUlti–STAge (MUSTA) approach [15, 27]. This utilizes a centeredflux in a predictor–corrector loop, solving the Riemann problem numerically,i.e. without using a priori information about waves.In order to calculate flux F i + the algorithm, in a one dimensional case,is as follows:1. Introduce auxiliary variables U ( l ) L and U ( l ) R and their fluxes F ( l ) L and F ( l ) R .2. Set U L = U i , U R = U i +1 .3. Calculate F ( l ) i + using a centered flux, U ( l ) L , U ( l ) R , F ( l ) L and F ( l ) R . If l reached a maximum number of iterations, stop.4. Solve Riemann problem locally: U ( l +1) L = U ( l ) L − ∆ t ∆ x (cid:18) F ( l ) i + − F ( l ) L (cid:19) ,U ( l +1) R = U ( l ) R − ∆ t ∆ x (cid:18) F ( l ) R − F ( l ) i + (cid:19) . (8)5. Go back to step 3.As a centered flux in step 3, we use the First ORder CEntered (FORCE)scheme: F force i + = 12 (cid:16) F lw i + + F lf i + (cid:17) (9)where F lw i + is the Lax-Wendroff type flux: F lw i + = F
12 ( U L + U R ) − α ∆ t x ( U R − U L ) ! (10)and F lf i + is the Lax-Friedrichs type flux: F lf i + = 12 ( F L + F R ) − ∆ x α ∆ t ( U R − U L ) (11)In a three-dimensional case α = 3, but other values may also be considered.To achieve second order accuracy in space and time, we extend ouralgorithm with the MUSCL-Hancock scheme. The basic idea of this schemeis to use more cells to interpolate inter-cell values and evolve them half a timestep: 7. Replace cell average values U ni by a piecewise linear function inside i -thcell: U i ( x ) = U ni + ( x − x i )∆ x ∆ i (12)where ∆ i is a slope vector and will be defined later.In the local coordinates the points x = 0 and x = ∆ x correspond toboundaries of the cell x i − and x i + . The values at these points are U Li = U ni − ∆ i / U Ri = U ni + ∆ i / U Li and U Ri by a time ∆ t :˜ U Li = U Li + 12 ∆ t ∆ x ( F ( U Li ) − F ( U Ri ))+ 12 ∆ t ∆ y ( G ( U Li ) − G ( U Ri ))+ 12 ∆ t ∆ z ( H ( U Li ) − H ( U Ri ))˜ U Ri = U Ri + 12 ∆ t ∆ x ( F ( U Li ) − F ( U Ri )) (13)+ 12 ∆ t ∆ y ( G ( U Li ) − G ( U Ri ))+ 12 ∆ t ∆ z ( H ( U Li ) − H ( U Ri ))3. Use ˜ U Li and ˜ U Ri as U L and U R in MUSTA.A simple choice for the slope ∆ i in Eq. 12 is:∆ i = 12 ( U ni +1 − U ni − ) (14)which indeed results in a second-order accurate algorithm. However, aspredicted by Godunov’s theorem, it has the unpleasant effect of producingspurious oscillations in the case of strong gradients. Complex study of MUSTAschemes can be found in [28]. To avoid such oscillations and to solve problems that appear in the presenceof shocks, discontinuities or sharp changes, flux limiting and slope limiting8ethods have been proposed [30, 31]. We employed a slope limiting method;instead of ∆ i as in Eq. 14 we use˜∆ i = ξ ( r i )( U i − U i − ) (15)in Eq. 12, where ξ is called the slope limiter and r i = U i +1 − U i U i − U i − . (16)Then one can calculate U Li and U Ri using the following relations U Li = U i − ξ (1 /r i )( U i +1 − U i ) ,U Ri = U i + ξ ( r i )( U i − U i − ) . (17)There are a number of possible choices for ξ , each with its own character-istics and features. The four different limiters investigated in this paper areMinbee (MB) [32], Superbee (SB) [33], van Albada (VA) [34] and van Leer(VL) [35]. They are expressed as: ξ MB ( r ) = max(0 , min(1 , r )) (18) ξ SB ( r ) = max(0 , min(2 r, , min( r, ξ VA ( r ) = r + rr + 1 (20) ξ VL ( r ) = r + | r | | r | (21)All the four limiters are symmetric, i.e. they meet the symmetry property: ξ ( 1 r ) = ξ ( r ) r (22)that provides that forward and backward gradients are treated in the samemanner. Otherwise, the results on the left would differ from those on theright despite initial symmetry in the system.The limiters are designed to reduce the scheme to first order accuracy nearshocks, and keep higher order in smooth areas. Introducing non–linearityin this way reduces spurious oscillations and retains good accuracy of thesolution.The two most extreme slope limiters are Minbee and Superbee. The firstone is the most dissipative and the second - the least. Between them lies theadmissible region for second order total variation diminishing (TVD) limiters.Here it should also be stressed, that results are very sensitive to the choiceof slope in Eq. 14 and the slope limiter—both the formula for ξ , and the waythat it is applied. 9 .6 Calculating the hydrodynamic flux To complete the chapter about numerical methods for relativistic hydrody-namics, one more detail must still be dealt with. It is easy to change fromrest frame variables and velocities to the conserved variables used in theintegration. The inverse transformation, however, is not trivial.It is needed, however, each time we compute the flux, since the equationof state (2) is defined as a function of e and n , and not E and R . As in [36],we can invert equations (4) and get e = E − M vn = R √ − v (23) v = ME + p ( e, n )= ME + p ( E − M v, R √ − v )The set of equations (23) can be solved numerically, starting from the lastone as the fixed point equation for v . c h a r g e d e n s i t y [ G e V / f m ] distance from center [fm]simulationanalytic 00.10.20.30.40.50.60.70.8 -3 -2 -1 0 1 2 3 4 v e l o c i t y [ c − ] distance from center [fm]simulationanalytic Figure 1: Sod shock tube, charge density in the local rest frame of the fluid(left) and velocity (right), MUSTA-FORCE with no limiter.
In the case of hydrodynamics simulations, fully (3+1)-dimensional simulationsin space and in time are much-desired, to describe the system’s evolution10ithout any assumptions regarding its symmetries and without decreasingthe dimension of the problem. Moreover, event-by-event calculations becomea major point of interest, with fluctuating initial conditions and with a largeamount of statistics. Such simulations are extremely expensive in terms ofcomputing power and require an efficient and fast computer code.Therefore, we implemented the hydrodynamics simulation algorithm on aGPU using an NVIDIA CUDA framework [37]. Due to the large numericalgrids in the simulations and the complexity of the computations, we usedsurface memory to store the simulation data. The state of the system issaved as 5 single precision floating-point numbers per lattice cell – the energydensity, conserved charge density, and 3 momenta density; all in the laboratoryreference frame . The maximum grid that fits then within the surface memorylimitations is 240 . To verify the simulation reliability, the MUSTA-FORCE algorithm itself wastested (along with the four slope limiters – Minbee, Superbee, van Albada andvan Leer) against three analytical solutions to relativistic hydrodynamics –the Sod shock tube [38, 39], the Hubble-like expansion [40] and ellipsoidal flow[41, 42]. For each of them, plots of chosen variables are presented togetherwith the theoretical curves.The number of stages and the order of Runge–Kutta method were fittedexperimentally. In all the presented cases four MUSTA stages guaranteedstability and good numerical accuracy in acceptable calculation time on GPU.The third order accurate Runge–Kutta method was used for time integration.Many numerical experiments for different parameters were carried out c h a r g e d e n s i t y [ G e V / f m ] distance from center [fm]simulationanalytic 00.10.20.30.40.50.60.70.8 -3 -2 -1 0 1 2 3 4 v e l o c i t y [ c − ] distance from center [fm]simulationanalytic Figure 2: Sod shock tube, charge density in the local rest frame of the fluid(left) and velocity (right), MUSTA-FORCE with with Minbee limiter.11o verify correctness, convergence and robustness of the MUSTA-FORCEmethod. We present here the results for one representative set of parametersper analytical solution. The initial conditions and parameters like time step∆ t and grid spacing ∆ x , ∆ y and ∆ z were given for each experiment separately.In all the cases the Courant number was less than 1. More detailed analysisof CFL numbers and stability limits in the MUSTA approach can be foundin [43]. The first test is a solution to the Riemann problem. This is a one dimensionalsolution, whose initial state comprises two regions of stationary fluid with acharge and pressure discontinuity in the middle.When the discontinuity is big enough, a relativistic shock wave appearsin the solution. The initial conditions (given in Table 1 together with otherparameters) were chosen to produce such a shock wave.parameter valuegrid size 500grid spacing 0 . . p L p R n L n R c h a r g e d e n s i t y [ G e V / f m ] distance from center [fm]simulationanalytic 00.10.20.30.40.50.60.70.8 -3 -2 -1 0 1 2 3 4 v e l o c i t y [ c − ] distance from center [fm]simulationanalytic Figure 3: Sod shock tube, charge density in the local rest frame of the fluid(left) and velocity (right), MUSTA-FORCE with Superbee limiter.12 c h a r g e d e n s i t y [ G e V / f m ] distance from center [fm]simulationanalytic 00.10.20.30.40.50.60.70.8 -3 -2 -1 0 1 2 3 4 v e l o c i t y [ c − ] distance from center [fm]simulationanalytic Figure 4: Sod shock tube, charge density in the local rest frame of the fluid(left) and velocity (right), MUSTA-FORCE with van Albada limiter. c h a r g e d e n s i t y [ G e V / f m ] distance from center [fm]simulationanalytic 00.10.20.30.40.50.60.70.8 -3 -2 -1 0 1 2 3 4 v e l o c i t y [ c − ] distance from center [fm]simulationanalytic Figure 5: Sod shock tube, charge density in the local rest frame of the fluid(left) and velocity (right), MUSTA-FORCE with van Leer limiter.The solution is divided into waves: the shock wave, the contact disconti-nuity, and the rarefaction wave.The results are presented in Figs. 1–5.For the no limiter case (Fig. 1) with the MUSTA-FORCE algorithm,oscillations at the contact points of waves become visible. The shock is alsosmeared out from the left side. For Minbee (Fig. 2) the overshoot is almostentirely gone, at the cost of some visible diffusion, especially in the shockwave region. Superbee in Fig. 3 is the most compressive limiter. The shocksare sharpened, but additional oscillations are introduced. Also the overshootclearly visible in the velocity plot is enhanced. The van Albada limiter inFig. 4 trades some sharpness for better rendition of the velocity profile; f andfor van Leer (Fig. 5) the shocks are sharp and the oscillations are gone, butthe overshoot is still significant.To sum up, van Albada and the Minbee limiter seem to be closest to the13nalytic solution, and those two will be presented in rest of the tests.
This is a three-dimensional, spherically symmetrical solution of matter thatexpands uniformly. The velocity is proportional to the distance from thecenter v = ~rt . The energy density is given by: e = e τ √ t − r ! c s ) (24)In our case the solution is well defined for r < t . For the test we set r < t − . e = v = 0) outside this region. This means that thesolution is exact only in the central area—on the periphery the matter willexpand into the vacuum, so a rarefaction wave is expected. The solution usesan ultra-relativistic equation of state p = c s e .Initial parameters are given in Table 2.parameter valuegrid size 120 grid spacing 0 . . e c s τ t y = z = 0 sections areshown through the three dimensional solution.The results are similar to those in the previous test. The schemes wereaccurate in the middle region, and the Minbee limiter has shown less diffusionthan van Albada at the sides. The last test uses the ellipsoidal solution, which is a generalized Hubble–likesolution (with velocity proportional to ~r ), a gaussian profile and vanishingpressure p = 0. The variables are given by the following equations:14 .1110 -6 -4 -2 0 2 4 6 e n e r g y d e n s i t y [ G e V / f m ] distance from center [fm]simulationanalytic 00.10.20.30.40.50.60.70.80.91 -6 -4 -2 0 2 4 6 v e l o c i t y [ c − ] distance from center [fm]simulationanalytic Figure 6: Hubble-like expansion, energy density in the local rest frame of thefluid (left) and velocity (right), MUSTA-FORCE with Minbee limiter. e n e r g y d e n s i t y [ G e V / f m ] distance from center [fm]simulationanalytic 00.10.20.30.40.50.60.70.80.91 -6 -4 -2 0 2 4 6 v e l o c i t y [ c − ] distance from center [fm]simulationanalytic Figure 7: Hubble-like expansion, energy density in the local rest frame of thefluid (left) and velocity (right), MUSTA-FORCE with van Albada limiter. e = C e Q i ( t + T i ) exp − b e t τ ! (25) n = C n Q i ( t + T i ) exp − b n t τ ! (26) ~v = a ( t ) xt , a ( t ) yt , a ( t ) zt ! (27)where τ = r t − P i a i x i , a i ≡ a i ( t ) = t/ ( t + T i ) and C e , C n , b e , b n , T i areconstants; i = 1 , , y = z = 0 sections are15arameter valuegrid size 120 grid spacing 0 . . t C e C n . b e b n T . T . T . e n e r g y d e n s i t y [ G e V / f m ] distance from center [fm]simulationanalytic 00.20.40.60.81 -6 -4 -2 0 2 4 6 v e l o c i t y [ c − ] distance from center [fm]simulationanalytic Figure 8: Ellipsoidal flow, energy density in the local rest frame of the fluid(left) and velocity (right), MUSTA-FORCE with Minbee limiter.shown through the three dimensional solution. This solution, despite beingpressureless, most resembles a physically relevant situation. The Minbeelimiter has shown again less oscillations than others at the sides.
As a result of this paper, we have implemented and tested a MUSTA-FORCEalgorithm dedicated to solving conservative field equations. The universalMUSTA-FORCE algorithm proved to be quite efficient and robust. Despitepromising initial results, it turned out that the MUSTA-FORCE producednumerical oscillations. Fig. 10 illustrates two dimensional energy densitycross section in the xy plane of the MUSTA-FORCE with Albada limiter.16 e n e r g y d e n s i t y [ G e V / f m ] distance from center [fm]simulationanalytic 00.20.40.60.81 -6 -4 -2 0 2 4 6 v e l o c i t y [ c − ] distance from center [fm]simulationanalytic Figure 9: Ellipsoidal flow, energy density in the local rest frame of the fluid(left) and velocity (right), MUSTA-FORCE with van Albada limiter.It confirms the clearly visible anisotropy of the numerical oscillations. Toaddress this issue, we performed a detailed study of universal slope limitersto reduce oscillations or asymmetric propagation. In all test cases the Minbeelimiter has shown the best numerical properties.This implementation of a relativistic hydrodynamic code is designedto run efficiently on contemporary graphics processing units, which havemany times more computing power compared to ordinary CPU processors.Benchmarks and comparison to an equivalent implementation of some of thesealgorithms in C show speedups of over 2 orders of magnitude. Thanks tosuch performance levels, event-by-event analyses can be conducted efficientlywith a relatively low cost, just a few modern workstations. Our work alsoprovides a framework for implementing other algorithms for high-statisticsevent-by-event hydrodynamic simulations.All the calculations were made as fully (3+1)-dimensional simulations inspace and in time, without any assumptions regarding its symmetries andwithout decreasing the dimension of the problem. Such (3+1)-dimensionalnumerical simulation on high-resolution Cartesian grids is a novel result. Itallows us to run hydrodynamic simulations, provide information about prop-erties of the relativistic bulk nuclear matter (QGP) by studying jet-mediuminteraction and jet-induced flow. Furthermore event-by-event simulationsallow us to investigate initial condition fluctuation of the heavy ion collisionin the pre-equlibrium stage and it might have an impact on the dynamic inthe nuclear medium. 176 -4 -2 0 2 4 6 x [fm]-6-4-20246 y [ f m ] Figure 10: Ellipsoidal flow in the xy plane, energy density in the local restframe of the fluid, MUSTA-FORCE with van Albada limiter. References [1] E. Gourgoulhon, "An introduction to relativistic hydrodynamics,"
EASPublications Series , vol. 21, 2008, pp. 43–79.[2] P. Huovinen and P. V. Ruuskanen, "Hydrodynamic models for heavy ioncollisions,"
Ann. Rev. Nucl. Part. Sci. , vol. 56, 2006, pp. 163–206.[3] V. Schneider, U. Katscher, D. H. Rischke, B. Waldhauser, J. A. Maruhnand C.-D. Munz, "New algorithm for ultra-relativistic numerical hydro-dynamics,"
Journal of Computational Physics , vol. 105, 1993, pp. 92-107.[4] D. S. Balsara, "Riemann solver for relativistic hydrodynamics,"
Journalof Computational Physics , vol. 114, 1994, pp. 284-297.[5] H. Cheng, H. Yang, and Y. Zhang, "Riemann problem for the ChaplyginEuler equations of compressible fluid flow."
International Journal of onlinear Sciences and Numerical Simulation , 11 no. 11, 2010, pp. 985-992.[6] Y. Akamatsu, S. Inutsuka, C. Nonaka and M. Takamoto, "A new schemeof casual viscous hydrodynamics for relativistic heavy-ion collisions:a Riemann solver for quark-gluon plasma," Journal of ComputationalPhysics , vol. 256, 2014, pp. 34-54.[7] I. A. Karpenko, P. Huovinen and M. Bleicher, "A 3+1 dimensionalviscous hydrodynamic code for relativistic heavy ion collisions,"
ComputerPhysics Communications , vol. 185, no. 11, pp. 3016-3027, 2014.[8] Y. Tachibana and T. Hirano, "Momentum transport away from a jet inan expanding nuclear medium,"
Phys. Rev. , vol. C90, no. 2, pp. 021902,2014.[9] P. Bożek, "Flow and interferometry in (3 + 1)-dimensional viscoushydrodynamics,"
Phys. Rev. C , vol. 85, no. 11, pp. 034901–034909, 2012.[10] S. V. Akkelin, Y. Hama, Iu. A. Karpenko and Yu. M. Sinyukov, "Hydro-kinetic approach to relativistic heavy ion collisions,"
Phys. Rev. C , vol. 78,no. 3, pp. 034906–034920, 2008.[11] I. Sagert, W. Bauer, D. Colbry, J. Howell, R. Pickett, A. Staber and T.Strother, "Hydrodynamic shock wave studies within a kinetic Monte Carloapproach,"
Journal of Computational Physics , vol. 266, 2014, pp. 191-213[12] J. Gerhard, V. Lindenstruth and M. Bleicher, "Relativistic hydrody-namics on graphic cards,"
Comput.Phys.Commun. , vol. 184, no. 2 2013,pp. 311-319.[13] J. P. Boris and D. L. Book, "Flux-corrected transport. I. SHASTA, afluid transport algorithm that works,"
Journal of Computational Physics ,vol. 11, no. 1, pp. 38–69, 1973.[14] E. F. Toro,
Multi-stage predictor-corrector fluxes for hyperbolic equa-tions , Isaac Newton Institute for Mathematical Sciences Preprint SeriesNI03037-NPA, University of Cambridge, UK, 2003.[15] E. F. Toro, "MUSTA: a multi-stage numerical flux,"
Appl. Numer. Math. ,vol. 56, no. 10, 2006, pp. 1464–1479.[16] S. Cygert, J. Porter-Sobieraj, D. Kikola, J. Sikorski and M. Slodkowski,"Towards an efficient multi-stage Riemann solver for nuclear physics19imulations," in Computer Science and Information Systems (FedCSIS),2013 Federated Conference on , 2013, pp. 441–446.[17] M. L. Miller, K. Reygers, S. J. Sanders and P. Steinberg, "Glaubermodeling in high energy nuclear collisions,"
Ann.Rev.Nucl.Part.Sci. ,vol. 57, 2007, pp. 205–243.[18] K. J. Eskola, H. Niemi and P. V. Ruuskanen, "Elliptic flow from pQCD+ saturation + hydro model,"
J. Phys. G35 , 2008.[19] V. K. Magas, L. P. Csernai and D. Strottman, "Effective string ropemodel for the initial stages of ultra-relativistic heavy ion collisions,"
Nuclear Physics A , vol. 712, no. 1-2, 2002, pp. 167–204.[20] B. Alver, M. Baker, C. Loizides and P. Steinberg, "The PHOBOS GlauberMonte Carlo", 2008, arXiv:0805.4411.[21] W. Broniowski, M. Rybczynski and P. Bozek, "GLISSANDO: Glauberinitial-state simulation and more,"
Comput.Phys.Commun. , vol. 180,no. 1, 2009, pp. 69–83.[22] M. Rybczynski, G. Stefanek, W. Broniowski and P. Bozek, "GLIS-SANDO 2: Glauber initial-state simulation and more..., ver. 2,"
Com-put.Phys.Commun. , vol. 185, no. 6, 2014, pp. 1759-1772.[23] H. J. Drescher, S. Ostapchenko, T. Pierog and K. Werner, "Initialcondition for QGP evolution from NEXUS,"
Phys.Rev. C , vol. 65:054902,2002.[24] R. Derradi de Souza, J. Takahashi and T. Kodama, "Effects of initialstate fluctuations in the final state elliptic flow measurements using theNeXSPheRIO model,"
Phys. Rev. , vol. C85, pp. 054909, 2012.[25] T. Pierog, Iu. Karpenko, S. Porteboeuf, S. and K. Werner, "New Devel-opments of EPOS 2," 2010, arXiv:1011.3748.[26] S. A. Bass, M. Belkacem, M. Bleicher, M. Brandstetter, L. Bravinaet al. "Microscopic models for ultrarelativistic heavy ion collisions,"
Prog.Part.Nucl.Phys. , vol. 41, 1998, pp. 255–369.[27] E. F. Toro,
Riemann solvers and numerical methods for fluid dynamics:A practical introduction , Springer, 1999.[28] E.F. Toro, and V.A. Titarev, "MUSTA fluxes for systems of conservationlaws"
Journal of Computational Physics , 216(2), pp.403-429., 2006.2029] E. Constantinescu and A. Sandu, "Explicit time stepping methods withhigh stage order and monotonicity properties," in Proc. of the 9th Inter-national Conference on Computational Science , 2009, pp. 293–301.[30] A. Harten and S. Osher, "Uniformly high-order accurate nonoscillatoryschemes,"
SIAM Journal on Numerical Analysis , 1987, pp. 279–309.[31] M. Berger, S. M. Murman and M. J. Aftosmis, "Analysis of slope limiterson irregular grids",
In 43rd AIAA Aerospace Sciences Meeting and Exhibit ,2005.[32] P. K. Sweby and M. J. Baines, "Convergence of Roe’s scheme for the gen-eral non-linear scalar wave equation,"
Reading Univ. Numerical AnalysisReport , 1981.[33] P. L. Roe, "Some contributions to the modelling of discontinuous flows,"
In Proc. AMS/SIAM Seminar , San Diego, 1983.[34] G. D. van Albada, B. van Leer and W. W. Roberts, "A comparativestudy of computational methods in cosmic gas dynamics,"
Astronomyand Astrophysics , vol. 108, 1982, pp. 76–84.[35] B. van Leer, "Towards the ultimate conservative difference scheme,Monotonicity and conservation combined in a second order scheme,"
J.Comp.Phys. , vol. 14, 1974, pp. 361–370.[36] D. H. Rischke, "Fluid dynamics for relativistic nuclear collisions,"
InHadrons in Dense Matter and Hadrosynthesis , Springer Verlag, 1999.[37] NVIDIA Corporation: CUDA C Programming Guide, 2014.[38] K. W. Thompson, "The special relativistic shock tube,"
Journal of FluidMechanics , vol. 171, 1986, pp.365–375.[39] J. M. Martí and E. Müller, "Numerical hydrodynamics in special relativ-ity,"
Living Rev. Relativity , vol. 6, 2003.[40] M. Chojnacki, W. Florkowski and T. Csörgö, "Formation of Hubble-likeflow in little bangs,",
Physical Review C , vol. 71, no. 4, pp. 044902, 2005.[41] Yu. M. Sinyukov and Iu. A. Karpenko, "Quasi-inertial ellipsoidal flowsin relativistic hydrodynamics," 2005, arXiv:nucl-th/0505041.[42] Y. M. Sinyukov and I. A. Karpenko, "Ellipsoidal flows in relativistichydrodynamics of finite systems"
Acta Physica Hungarica Series A, HeavyIon Physics , vol. 25, no. 1, pp. 141–147, 2006.2143] P.M. Blakely, N. Nikiforakis, and W.D. Henshaw, "Assessment of theMUSTA approach for numerical relativistic hydrodynamics"