An MPI-CUDA implementation of an improved Roe method for two-layer shallow water systems

An MPI-CUDA implementation of an improved Roe method for two-layer shallow water systems Marc de la Asunci´ona, Jos´e M. Mantasa, Manuel J. Castrob, E. D. Fern´andez-Nietoc aDpto. Lenguajes y Sistemas Inform´aticos, Universidad de Granada bDpto. An´alisis Matem´atico, Universidad de M´alaga cDpto. Matem´atica Aplicada I, Universidad de Sevilla Abstract The numerical solution of two-layer shallow water systems is required to simulate accurately stratified fluids, which are ubiquitous in nature: they appear in atmospheric flows, ocean currents, oil spills, . . . Moreover, the implementation of the numerical schemes to solve these models in realistic scenarios imposes huge demands of computing power. In this paper, we tackle the acceleration of these simulations in triangular meshes by exploiting the combined power of several CUDA-enabled GPUs in a GPU cluster. For that purpose, an improvement of a path conservative Roe type finite volume scheme which is specially suitable for GPU implementation is presented, and a distributed implementation of this scheme which uses CUDA and MPI to exploit the potential of a GPU cluster is developed. This implementation overlaps MPI communication with CPU-GPU memory transfers and GPU computation to increase efficiency. Several numerical experiments performed on a cluster of modern CUDA-enabled GPUs show the efficiency of the distributed solver. Keywords: Shallow water simulation, GPU cluster computing, finite volume schemes, unstructured meshes, CUDA, MPI Preprint submitted to Journal of Parallel and Distributed Computing March 4, 2011 1. Introduction The two-layer shallow water system has been used as the numerical model to simulate several phenomena related to stratified geophysical flows such as atmospheric flows, ocean currents, oil spills or tsunamis generated by underwater landslides. The simulation of these phenomena gives to place to very long lasting simulations in big computational domains. Moreover, some of these phenomena (tsunami propagation or oil spills) could require real time calculation. Therefore, extremely efficient numerical schemes and implementations are needed to be able to analyze those problems in practical execution times. Since the numerical schemes to solve shallow water systems usually exhibit a high degree of potential parallelism, the development of parallel versions of these schemes for high performance platforms seems to be a suitable way of achieving the required performance in realistic applications. A cost effective way of obtaining a substantially higher performance in these applications consists in using Graphics Processor Units (GPUs). These architectures make it possible to obtain performances that are orders of magnitude faster than a standard CPU and are growing in popularity among the scientific and engineering community [1, 2]. Moreover, several GPU programming toolkits such as CUDA [3] have been developed to facilitate the programming of GPUs for general purpose applications. There are previous proposals to port finite volume one-layer shallow water solvers to a GPU by using a graphics-specific programming language [4, 5, 6], but currently most of the proposals to simulate shallow flows on a single GPU are based on the CUDA programming model. A CUDA solver for one2 layer system based on the first order finite volume scheme presented in [7] is described in [8] to deal with structured regular meshes. The extension of this CUDA solver for two-layer shallow water system is presented in [9]. There also exist proposals to implement, using CUDA-enabled GPUs, high order schemes to simulate one-layer systems [10, 11, 12] and to implement first-order schemes for one and two-layer systems on triangular meshes [13]. Although the use of single GPU systems makes it possible to satisfy the performance requirements of several applications which are not computationally expensive, this situation is not frequent in most realistic applications. Many applications require to handle huge meshes and large number of time steps. Moreover, in some applications real time accurate predictions (for instance, to approximate the effect of an unexpected oil spill) could be required. The characteristics of these applications suggest to combine the power of multiple GPUs to satisfy the performance requirements. One approach to use several GPUs in these problems is based on programming shared memory multi-GPU desktop systems. These platforms have been used to accelerate considerably fluid dynamic [14] and shallow water [15] simulations by combining shared memory programming primitives to manage threads in CPU and CUDA to program the GPU. Although this isa cost-effective approach, these platforms only offer a reduced number of GPUs and more flexible systems are desirable to answer to the growing performance requirements of many realistic applications. A more flexible approach to obtain the required performance involves to use clusters of GPU-enhanced computers where each node is equipped with a single GPU or with a multi-GPU system. The computation on GPU clusters 3 could make it possible to scale the reduction in execution time according to the number of GPUs (which can be easily increased). Therefore, this approach is more flexible than using a multi-GPU desktop system and the memory limitations of a GPU-enhanced node can be overcome by suitably distributing the data among the nodes, enabling us to simulate significantly larger realistic models and with greater precision. The use of GPU clusters to accelerate data intensive computations is gaining in popularity [16, 17, 18, 19]. In [20], a scheme to solve one-layer shallow water systems is implemented on a GPU cluster for real-time tsunami simulation. Most of the proposals to exploit GPU clusters in scientific computing use MPI [21] to implement the communication among the processes of the distributed system and CUDA [3] to program the GPU (or GPUs) associated to each node. A common way to reduce the remote communication overhead in these distributed implementations consists in using non-blocking communication MPI functions to overlap the data transfers between nodes of the cluster with the GPU computation and the CPU-GPU data transfers. This work deals with the acceleration of the numerical solution of twolayer shallow water systems by exploiting the parallelization of an improved finite volume scheme for unstructured meshes on GPU clusters. For that purpose, a distributed implementation of this scheme for a GPU cluster is developed by using MPI and CUDA. This implementation incorporates an efficient management of the distributed unstructured mesh and mechanisms to overlap computation with communication. The outline of the article is as follows: the next section describes the underlying mathematical model and presents an improvement of a first order 4 Roe type finite volume scheme, called IR-Roe scheme. Section 3 describes a data parallel version of the IR-Roe scheme. The efficiency of several implementations of the IR-Roe scheme and the classical Roe scheme [7, 22] is compared in Section 4. In the two next sections we describe a single and a multi-GPU distributed implementation, respectively, of the method for triangular meshes using the CUDA framework. Section 7 shows the experimental results obtained when the implementations are applied to solve an internal dam break problem on a cluster of 4 NVIDIA Fermi GPUs. Finally, Section 8 summarizes the main conclusions and presents the future work. 2. An efficient numerical scheme 2.1. The two-layer shallow water system Let us consider the system of equations governing the 2d flow of two superposed immiscible layers of shallow fluids in a subdomainΩ⊂R2: ∂W ∂t+∂F1 ∂x(W)+∂F2 ∂y(W)=B1(W)∂W ∂x+B2(W)∂W ∂y+S1(W)∂H ∂x+S2(W)∂H ∂y, (1) where W=!h1q1,x q1,y h2q2,x q2,y "T, F1(W)=#q1,x q2 1,x h1 +1 2gh2 1 q1,xq1,y h1 q2,x q2 2,x h2 +1 2gh2 2 q2,xq2,y h2$T F2(W)=#q1,y q1,xq1,y h1 q2 1,y h1 +1 2gh 2 1q2,y q2,xq2,y h2 q2 2,y h2 +1 2gh2 2$T , Sk(W)=!0gh1(2 −k)gh1(k−1) 0 gh2(2 −k)gh2(k−1) "T,k=1,2, Bk(W)=  0P1,k(W) rP2,k(W)0 Pl,k(W)=     0 0 0 −ghl(2 −k)00 −ghl(k−1) 0 0      l=1,2. 5 Index 1 in the unknowns makes reference to the upper layer and index 2 to the lower one; gis the gravity and H(x), the depth function measured from a fixed level of reference; r=ρ1/ρ2is the ratio of the constant densities of the layers (ρ1<ρ2) which, in realistic oceanographical applications, is close to 1. Finally, hi(x,t) and qi(x,t) are, respectively, the thickness and the mass-flow of the i-th layer at the point xat time t, and they are related to the velocities ui(x,t) = (ui,x(x,t),u i,y(x,t)), i=1,2 by the equalities: qi(x,t)=ui(x,t)hi(x,t),i=1,2. Let us define the matrices Ak(W)=Jk(W)−Bk(W), k=1,2 where Jk(W)=∂Fk ∂W(W) are the Jacobians of the fluxes Fk, and we assume that (1) is strictly hyperbolic. Let us also remark that the system (1) verifies the property of invariance by rotations. Effectively, let us define Tη=  Rη0 0Rη  ,R η=     1 0 0 0ηxηy 0−ηyηx      , and let us denote Fη(W)=F1(W)ηx+F2(W)ηy,B(W) = (B1(W),B 2(W)), and S(W)=(S1(W),S 2(W)). Then Fη(W)=T−1 ηF1(TηW),T ηB(W)·η=B1(TηW),T ηS(W)·η=S1(TηW) (2) Moreover, it is easy to check that TηWverifies the system ∂t(TηW)+∂ηF1(TηW)=B1(TηW)∂ηW+S1(TηW)∂ηH+Qη⊥,(3) where Qη⊥=Tη#−∂η⊥Fη⊥(W)+B(W)·η⊥∂η⊥W+S(W)·η⊥∂η⊥H$. 6 2.2. The IR-Roe Numerical Scheme To discretize (1) the computational domain Ωis decomposed into cells or finite volumes: Vi⊂R2. Here, it is assumed that the cells are triangles. Let us denote by Lthe number of triangles of the mesh. Given a finite volume Vi,|Vi|will represent its area; Ni∈R2its center; Nithe set of indexes jsuch that Vjis a neighbor of Vi;Γij the common edge of two neighboring triangles Viand Vj, and |Γij |its length; ηij =(ηij,x,ηij,y) the normal unit vector at the edge Γij pointing towards the triangle Vj; and Wn ithe constant approximation to the average of the solution in the triangle Viat time tnprovided by the numerical scheme. Let us now briefly describe the Roe type scheme for system (1) that can be defined taking into account the property of invariance by rotations [23]: Wn+1 i=Wn i−∆t |Vi|+ j∈Ni|Γij|FIR−ROE ij −(4) where FIR−ROE ij −is defined as follows: 1. Let us define Wη=[h1q1,ηh2q2,η]T=Tη(W)[1,2,4,5], and Wη⊥= [q1,η⊥q2,η⊥]T=Tη(W)[3,6], where W[i1,··· ,is]is the vector defined from vector W, using its i1-th, . . . , is-th components. 2. Let Φ− ηij be the 1D numerical Roe flux associated to the 1D two-layer shallow-water system defined by the 1-st, 2-nd, 4-th and 5-th equations of system (3) where the term Qη⊥ ij has been neglected: Φ− ηij =P− ij (F1(Wηij ,j)−F1(Wηij ,i)−Bij(Wηij ,j −Wηij ,i)−Sij(Hj−Hi)) +F1(Wηij ,i). 7 where F1(Wηij )=F1(Tηij W)[1,2,4,5],Bij (Wηij,j −Wηij,i )=,B1,ij(Tηij (Wj− Wi))-[1,2,4,5],Sij(Hj−Hi)=,S1,ij(Hj−Hi)-[1,2,4,5]. Finally, P− ij =1 2Kij(I−sgn(Dij))K−1 ij , where Iis the identity matrix, Kij is the matrix whose columns are the eigenvectors of the matrix Aij, and sgn(Dij) is the diagonal matrix whose coefficients are the signs of the eigenvalues of Aij, being Aij =         0 1 0 0 gh1,ij −u2 1,ηij 2u1,ηij gh1,ij 0 0 0 0 1 rgh2,ij 0gh2,ij −u2 2,ηij 2u2,ηij         with uk,ηij =uk,ij ·ηij,k=1,2 and hk,ij =hk,i+hk,j 2, where uk,l,ij = √hk,iuk,l,i+√hk,j uk,l,j √hk,i+√hk,j ,k=1,2, l=x, y. 3. Let us define Φ− η⊥ ij =4(Φ− ηij )[1]u∗ 1,η⊥ ij (Φ− ηij )[3]u∗ 2,η⊥ ij 5T , where u∗ k,η⊥ ij is defined as follows u∗ k,η⊥ ij =     qk,η⊥ ij ,i hk,i if (Φηij )[2k−1] >0 qk,η⊥ ij ,j hk,j otherwise k=1,2. Let us remark that Φ− η⊥ ij is the numerical flux associated to the 3-rd and 6-th equations of system (3) where, again, the term Qη⊥ ij has been neglected. Its derivation has been done following the main ideas of the HLLC method for the shallow water system introduced in [24] as qk,η⊥ ij , k=1,2 can be seen as a passive scalar that is convected by the flow. 4. Finally, the global numerical flux is defined by FIR−ROE ij −=T−1 ηij F− ij , where F− ij =4(Φ− ηij )[1] (Φ− ηij )[2] (Φ− η⊥ ij )[1] (Φ− ηij )[3] (Φ− ηij )[4] (Φ− η⊥ ij )[2] 5. 8 A CFL condition must be imposed to ensure stability of both schemes: 1 2 ∆t Vi+ j∈Ni|Γij%Dij%∞≤γ,1≤i≤L, with 0 <γ≤1.(5) As in the case of systems of conservation laws, when sonic rarefaction waves appear it is necessary to modify the numerical scheme to get entropysatisfying solutions. For instance, the Harten-Hyman entropy fix technique [25] can be easily adapted here. Let us also remark that the scheme is path-conservative in the sense introduced by Pares in [26] and [22]. It is well-balanced for stationary solutions corresponding to water at rest. More general results concerning the consistency and well-balanced properties of Roe schemes have been studied in [22] and [27]. 3. Parallelization of the scheme Figure 1a shows a graphical description of the parallel algorithm, obtained from the description of the IR-Roe numerical scheme given in Section 2. The main calculation phases are identified with circled numbers, and the main sources of data parallelism are represented with overlapping rectangles. Initially, the finite volume mesh is constructed from the input data. Then the time stepping process is repeated until the final simulation time is reached. At the (n+ 1)-th time step, Equation (4) must be evaluated to update the state of each cell. Each of the main calculation phases present a high degree of parallelism because the computation at each edge or volume is independent with respect to that performed at other edges or volumes: 9 the submesh so that all the communication volumes that are adjacent to a particular submesh appear consecutively in the array. For example, in Figure 2b, note that volumes 10-12 are sent to the lower submesh, while volumes 12 and 13 are sent to the right submesh, thus overlapping the sendings. In order to perform this arrangement, firstly we build a list of communication volumes for each adjacent submesh. For example, in Figure 2b we would have two lists: [10, 11, 12] and [12, 13]. Now, for each communication volume of the submesh that must be sent to two MPI processes, we build a pair (p1,p 2), meaning that the volume must be sent to processes p1and p2. Figure 2c shows an example centered on submesh 4, where all the pairs are specified. Once all the pairs have been built, we perform a reordering of them (and their elements if necessary) so that we get a list of consecutive processes. In Figure 2c, the pairs are reordered obtaining: (0,1), (1,2), (2,3) and (5,6). This gives the consecutive list of processes 0, 1, 2, 3, 5 and 6. Now we carry out the adequate swaps: 1. In the list storing the volumes that are adjacent to submesh 0, we put the volume shared with submesh 1 at the end. 2. In the list storing the volumes that are adjacent to submesh 1, we put the volume shared with submesh 0 at the start, and the volume shared with submesh 2 at the end. 3. We continue processing the list in the same way until it finishes. Finally we join all the lists of communication volumes adequately to get the definitive block of communication volumes. Note that a submesh must know the ordering of the communication volumes that receives from another submesh. Therefore, at this point all sub16 meshes must send this information to their adjacent submeshes. Note also that this algorithm does not work when a submesh has two volumes that must be sent to the same submeshes. This is reflected by a duplicated pair in the sequence of pairs, but we can always perform a domain decomposition where this does not occur. It neither works when a submesh is formed by a single volume, but this will never happen in a real problem. 6.3. Multi-GPU Code We have implemented two versions of the multi-GPU algorithm: one with blocking MPI sends and receives, and another one which overlaps MPI communication with CPU-GPU memory transfers and kernel computation. Algorithm 1 shows the general steps of the non-overlapping implementation. In lines 3-5 we send to each adjacent submesh the communication volumes that are adjacent to it. Then, in lines 6-8 we receive from the same submeshes their communication volumes that are adjacent to our submesh. We have used bufferized MPI send operations with a given buffer to avoid deadlocks with big meshes. Lines 9-10 copy the received communication volumes to GPU memory. In line 14 a MPI reduction is performed to obtain the global minimum ∆tin all the MPI processes. Lines 16-17 copy the new states of the communication volumes of our submesh from device to host. Algorithm 2 shows the general steps of the overlapping implementation. Lines 3-5 (the reception of the communication volumes from the adjacent submeshes) overlap with lines 6-7 (the copy of the new states of our communication volumes from device to host). Then, lines 8-10 (the sending to each adjacent submesh of the communication volumes that are adjacent to it) 17 Algorithm 1 Non-overlapping multi-GPU algorithm 1: n←number of adjacent submeshes 2: while (t < tend)do 3: for i= 1 to ndo 4: Send communication volumes to adjacent submesh i 5: end for 6: for i= 1 to ndo 7: Receive communication volumes from adjacent submesh i 8: end for 9: CudaMemcpy(Layer 1 of comm. volumes from host to device) 10: CudaMemcpy(Layer 2 of comm. volumes from host to device) 11: processEdges<<<grid, block>>>(...) 12: computeDeltaTVolumes<<<grid, block>>>(...) 13: ∆t←getMinimumDeltaT(...) 14: MPI Allreduce(∆t, min ∆t,...) 15: computeVolumeStates<<<grid, block>>>(...) 16: CudaMemcpy(Layer 1 of comm. volumes from device to host) 17: CudaMemcpy(Layer 2 of comm. volumes from device to host) 18: t←t+ min ∆t 19: end while overlap with line 11 (the processing of the non-communication edges, since these edges do not need external data to be processed). In line 12 we wait for the communication volumes of the adjacent submeshes to arrive. Once they have arrived, in lines 13-14 we copy them to GPU memory. In line 15 only the communication edges are processed. 7. Experimental Results In this section we will test the single and multi-GPU implementations described in Sections 5 and 6, respectively. The test problem and the parameters are the same that were used in Section 4. 18 Algorithm 2 Overlapping multi-GPU algorithm 1: n←number of adjacent submeshes 2: while (t < tend)do 3: for i= 1 to ndo 4: Receive comm. volumes from adjacent submesh i(Non-blocking) 5: end for 6: CudaMemcpy(Layer 1 of comm. volumes from device to host) 7: CudaMemcpy(Layer 2 of comm. volumes from device to host) 8: for i= 1 to ndo 9: Send comm. volumes to adjacent submesh i(Non-blocking) 10: end for 11: processEdges<<<grid, block>>>(Non-communication edges) 12: MPI Waitall (Comm. volumes from adjacent submeshes) 13: CudaMemcpy(Layer 1 of comm. volumes from host to device) 14: CudaMemcpy(Layer 2 of comm. volumes from host to device) 15: processEdges<<<grid, block>>>(Communication edges) 16: computeDeltaTVolumes<<<grid, block>>>(...) 17: ∆t←getMinimumDeltaT(...) 18: MPI Allreduce(∆t, min ∆t,...) 19: computeVolumeStates<<<grid, block>>>(...) 20: t←t+ min ∆t 21: end while We have used the Chaco software [31] to divide a mesh into equally sized submeshes, the OpenMPI implementation [32] and the GNU compiler. All the programs were executed in a cluster formed by four Intel Xeon servers with 8 GB RAM each one, connected with a Gigabit Ethernet switch. Graphics cards used were two Tesla C2050 and two GeForce GTX 570. Since the GTX 570 card provides better performance for our programs than the Tesla, it is suitable to use the four cards to measure strong scalability taking the Tesla as the reference card. Table 2 shows the execution times in seconds for all the meshes and number of GPUs. Figure 3a shows graphically the 19 Tesla 2 Tesla C2050 2 Tesla + 2 GTX 570 Volumes C2050 Non-Overlap Overlap Non-Overlap Overlap 4016 0.0090 0.012 0.014 0.011 0.015 16040 0.045 0.038 0.040 0.027 0.031 64052 0.29 0.19 0.19 0.12 0.12 256576 2.10 1.18 1.16 0.69 0.63 1001898 15.63 8.40 7.98 4.60 4.08 2000608 45.37 23.34 23.00 12.08 11.66 3000948 82.84 43.52 41.83 23.34 21.23 Table 2: GPU execution times in seconds for the IR-Roe method. speedups obtained with the single GPU program executed on a Tesla C2050 with respect to the CPU versions of the IR-Roe method used in Section 4. Figure 3b shows the speedups obtained with the multi-GPU implementations with respect to one Tesla. As it can be seen, using a Tesla C2050, for big meshes we have reached a speedup of 32 and 10 with respect to monocore and quadcore CPU versions, respectively. As expected, the overlapping multi-GPU implementation outperforms the non-overlapping version, and the weak and strong scaling reached by the overlapping version are close to perfect for up to four GPUs. 8. Conclusions and future work In this paper we have presented an improvement of a first order wellbalanced Roe type finite volume solver for two-layer shallow water system. This numerical scheme has proved to be computationally more efficient than the classical Roe scheme and is more suitable to be implemented in modern CUDA-enabled GPUs than the classical Roe scheme. A multi-GPU dis20 (a) Tesla C2050 speedup with respect to serial and quadcore CPU versions (b) Multi-GPU speedup with respect to one Tesla C2050 Figure 3: Speedups obtained for one and several GPUs. tributed implementation of this scheme that works on triangular meshes has been implemented using MPI and CUDA. Numerical experiments carried out on a GPU cluster have shown the efficiency of this solver, obtaining weak and strong scaling close to perfect for up to four GPUs by overlapping MPI communications with CPU-GPU memory transfers and GPU computation. As further work, we propose to extend the proposal to enable high order numerical schemes and to integrate a dynamic load balancing strategy (which is necessary in problems, such as flood simulations, where the computational load for each spatial subdomain could vary dramatically). Acknowledgements The authors acknowledges partial support from the DGI-MEC projects MTM2008-06349-C03-03, MTM2009-11923 and MTM2009-07719. 21 References [1] M. Rumpf, R. Strzodka, Graphics Processor Units: New Prospects for Parallel Computing, in: Numerical Solution of Partial Differential Equations on Parallel Computers, Vol. 51 of Lecture Notes in Computational Science and Engineering, Springer, 2005, pp. 89–134. [2] J. Owens, M. Houston, D. Luebke, S. Green, J. Stone, J. Phillips, GPU computing, Proceedings of the IEEE 96 (5) (2008) 879–899. [3] NVIDIA Corporation, NVIDIA CUDA C Programming Guide 3.2, 2010. [4] T. Hagen, J. Hjelmervik, K.-A. Lie, J. Natvig, M. O. Henriksen, Visual simulation of shallow-water waves, Simulation Modelling Practice and Theory 13 (8) (2005) 716–726, Programmable Graphics Hardware. [5] M. Lastra, J. M. Mantas, C. Ure˜na, M. J. Castro, J. A. Garc´ıaRodr´ıguez, Simulation of shallow-water systems using graphics processing units, Mathematics and Computers in Simulation 80 (3) (2009) 598– 618. [6] W.-Y. Liang, T.-J. Hsieh, M. T. Satria, Y.-L. Chang, J.-P. Fang, C.-C. Chen, C.-C. Han, A GPU-Based Simulation of Tsunami Propagation and Inundation, in: Proceedings of the 9th International Conference on Algorithms and Architectures for Parallel Processing, ICA3PP ’09, Springer-Verlag, Berlin, Heidelberg, 2009, pp. 593–603. [7] M. Castro, J. Garc´ıa-Rodr´ıguez, J. Gonz´alez-Vida, C. Par´es, A parallel 2d finite volume scheme for solving systems of balance laws with non22 conservative products: Application to shallow flows, Computer Methods in Applied Mechanics and Engineering 195 (19-22) (2006) 2788–2815. [8] M. de la Asunci´on, J. M. Mantas, M. J. Castro, Simulation of onelayer shallow water systems on multicore and CUDA architectures, The Journal of Supercomputing (2010) 1–9. [9] M. de la Asunci´on, J. M. Mantas, M. J. Castro, Programming CUDAbased GPUs to simulate two-layer shallow water flows, in: Euro-Par 2010 - Parallel Processing, Vol. 6272 of Lecture Notes in Computer Science, Springer Berlin / Heidelberg, 2010, pp. 353–364. [10] M. J. Castro, S. Ortega, M. de la Asunci´on, J. M. Mantas, GPU computing for shallow water flow simulation based on finite volume schemes, Bol. Soc. Esp. Mat. Apl. 50 (2010) 27–45. [11] A. Brodtkorb, T. Hagen, K.-A. Lie, J. Natvig, Simulation and visualization of the saint-venant system using GPUs, Computing and Visualization in Science (2011) 1–13. [12] J. Gallardo, S. Ortega, M. de la Asunci´on, J. M. Mantas, Twodimensional compact third-order polynomial reconstructions. Solving nonconservative hyperbolic systems using GPUs, Journal of Scientific Computing (2011) 1–23. [13] M. J. Castro, S. Ortega, M. de la Asunci´on, J. M. Mantas, J. M. Gallardo, GPU computing for shallow water flow simulation based on finite volume schemes, Comptes Rendus M´ecanique 339 (2-3) (2011) 165–184, High Performance Computing. 23 [14] J. Thibault, I. Senocak, Accelerating incompressible flow computations withaPthreads-CUDA implementation onsmall-footprint multiGPU platforms, The Journal of Supercomputing (2010) 1–27. [15] M. L. Saetra, A. R. Brodtkorb, Shallow water simulations on multiple GPUs, Proceedings of the Para 2010 Conference, Lecture Notes in Computer Science (2010) Accepted for publication. [16] D. Komatitschand, G. Erlebacher, D. G¨oddeke, D. Mich´ea, High-order finite-element seismic wave propagation modeling with MPI on a large GPU cluster, J. Comput. Phys. 229 (2010) 7692–7714. [17] Z. Fan, F. Qiu, A. Kaufman, S. Yoakum-Stover, GPU Cluster for High Performance Computing, in: Proceedings of the 2004 ACM/IEEE conference on Supercomputing, SC ’04, IEEE Computer Society, Washington, DC, USA, 2004, pp. 47–. [18] Y. Zhang, F. Mueller, X. Cui, T. Potok, Data-intensive document clustering on graphics processing unit (GPU) clusters, J. Parallel Distrib. Comput. 71 (2011) 211–224. [19] R. Abdelkhalek, H. Calendra, O. Coulaud, J. Roman, G. Latu, Fast Seismic Modeling and Reverse Time Migration on a GPU Cluster, in: The 2009 High Performance Computing & Simulation - HPCS’09, Leipzig Allemagne, 2009, Best Paper Award at HPCS’09 Total. [20] M. Acu˜na, T. Aoki, Real-time tsunami simulation on multi-node GPU cluster, Supercomputing (2009) [Poster]. 24 [21] Message Passing Interface Forum, MPI: A Message Passing Interface Standard, Univ. of Tennessee, Knoxville, Tennessee. [22] M. Castro, E. Fern´andez, A. Ferreiro, A. Garc´ıa, C. Par´es, High order extension of Roe schemes for two dimensional nonconservative hyperbolic systems, J. Sci. Comput. 39 (2009) 67–114. [23] E. D. Fern´andez-Nieto, Modelling an numerical simulation of submarine sediment shallow flows: transport and avalanches, Bol. Soc. Esp. Mat. Apl. 49. [24] E. D. Fern´andez-Nieto, D. Bresch, J. Monnier, A consistent intermediate wave speed for a well-balanced HLLC solver, Comptes Rendus Mathematique 346 (13-14) (2008) 795 – 800. [25] J. H. A. Harten, Self-adjusting grid methods for one-dimensional hyperbolic conservation laws, J. Comp. Phys. 50 (1983) 235–269. [26] C. Par´es, High order extension of Roe schemes for two dimensional nonconservative hyperbolic systems, SIAM J. Num. Anal. 44 (2006) 300– 321. [27] C. Par´es, M. Castro, On the well-balance property of Roe’s method for nonconservative hyperbolic systems. applications to shallow-water systems, M2AN 38 (5) (2004) 821–852. [28] B. Chapman, G. Jost, R. van der Pas, Using OpenMP: Portable Shared Memory Parallel Programming, 2007. 25