Efficient multilayer shallow-water simulation system based on GPUs

Efficient multilayer shallow-water simulation system based on GPUs Miguel Lastraa,∗, Manuel J. Castro D´ıazb, Carlos Ure˜naa, Marc de la Asunci´ona aDepto. Lenguajes y Sistemas Inform´aticos. E.T.S Ingen´ıeria Inform´atica. Periodista Daniel Saucedo Aranda s/n. Universidad de Granada. bDepto. An´alisis Matem´atico. Facultad de Ciencias.Campus de Teatinos s/n. Universidad de M´alaga. Abstract The computational simulation of shallow stratified fluids is a very active research topic because these types of systems are very common in a variety of natural environments. The simulation of such systems can be modelled using multilayer shallow-water equations but do impose important computational requisites especially when applied to large domains. General Purpose Computing on Graphics Processing Units (GPGPU) has become a vivid research field due to the arrival of massively parallel hardware platforms (based on graphics cards) and adequate programming frameworks which have allowed important speed-up factors with respect to not only sequential but also parallel CPU based simulation systems. In this work we present a proposal for the simulation of shallow stratified fluids with an arbitrary number of layers using GPUs. The designed system does fully adapt to the many-core architecture of modern GPUs and several experiments have been carried out to illustrate its scalability and behavior on different GPU models. We propose a new elaborated 3D computational scheme for an underlying 2D mathematical model. This scheme allowed implementing a system capable of handling an arbitrary number of layers. The system adds no overhead when used for two-layer scenarios, compared to an existing 2D system specifically designed for just two layers. Our proposal is aimed at creating a GPU-based computational scheme suitable for the simulation of multilayer large-scale real-world scenarios. Keywords: Shallow-water, finite volumes, simulation, GPU, CUDA 1. Introduction The simulation of free surface or internal waves in shallow stratified fluids are commonly modelled by the multilayer shallow-water equations formulated as a conservation law with non-conservative products and source terms. Stratified fluids are ubiquitous in nature: they appear in atmospheric flows, ocean currents, estuarine systems,... This is the situation, for instance, in the Strait of Gibraltar, where surface water from the Atlantic inflows over saltier westwards-flowing Mediterranean water. Simulating those phenomena requires very long lasting simulations in big computational domains which is why very efficient implementations are needed to be able to analyze those problems in affordable computational times. A 2D multilayer shallow-water system can be discretized by the natural extension to 2D domains of a first order PVM path-conservative type finite volume scheme introduced in [8]. PVM schemes have been introduced in [8] in the framework of balance laws and non-conservative hyperbolic systems. They are defined in terms of viscosity matrices computed by a suitable polynomial evaluation of a Roe Matrix. These ∗Corresponding author Email addresses: [email protected] (Miguel Lastra), [email protected] (Manuel J. Castro D´ıaz), [email protected] (Carlos Ure˜na), [email protected] (Marc de la Asunci´on) 1Corresponding author. Tel: +34958246144 Preprint methods have the advantage that they only need some information about the eigenvalues of the system and no spectral decomposition of the Roe Matrix is required. As consequence, they are much faster than Roe schemes for systems with an increasing number of unknowns. In this work the PVM-2U method has been chosen among the PVM schemes introduced in [8] as it provides the best results, concerning computational time and accuracy. This method can be seen as a the natural extension of the scheme introduced by Degond et al. in [7] for non-conservative systems. Since the appearance of computing frameworks like Cuda [14] and openCL [12] the huge computational capabilities of modern Graphics Processing Units (GPUs) have been unveiled for computationally intensive physically based simulations. Nevertheless, elaborated software designs are required to take full advance of the processing power offered by GPUs. The Cuda computational model, provided by Nvidia [15], requires the computational tasks to be mapped to a set of threads that will be run in parallel. Each of these threads will run the same program (kernel) where the number of threads running at a point in time will depend on the hardware requirements of each thread. At the lower level, threads are run in groups of 32 threads called warps and it is necessary to avoid code divergences inside a warp to avoid serialization. On the other hand, coalesced memory accesses are also essential to allow several threads access the required data stored in memory by only issuing one physical memory access for a group of threads. When a warp requires a memory access operation, the whole warp is stopped until that memory operation has finished and another warp is run. Modern GPUs provide a large amount of registers to avoid costly context switches by keeping a large amount of threads in a ready to run state. Taking into account GPUs currently provide in the order of thousands of computing cores, a large amount of threads is required to avoid computing cores not receiving enough workload and idling. This means any simulation scheme must be designed to provide a high degree of thread level parallelism. The higher the number of threads that can be kept active, the higher the value of a metric called occupancy. On the other hand, as in many architectures, computing cores have pipelines which should be kept as full as possible to maximize the performance. This requires an adequate level of instruction level parallelism which means having certain degree of independence between the instructions run by each thread. This fact requires that each thread must also receive an adequate level of workload and as much independence as possible from memory accesses (arithmetic intensity). The computing cores in Nvidia branded hardware are grouped into multiprocessors (SMXs) and threads are also grouped into thread blocks. Each thread block is run on the same SMX and this allows using an intra-block synchronization mechanism and also intra-block collaboration through the use of a high-speed shared memory (shared among the threads of each block). Dividing the computations that make up the simulation program into independent threads and grouping them into blocks requires also to consider the resources required by each thread and several hardware level limits like the maximum resident threads and blocks per SMX to keep the occupancy level high (although a higher occupancy does not automatically translate into a higher performance). There are many examples in literature of successfully using GPUs for shallow-water simulation [5, 4, 2, 10, 3], even before frameworks like Cuda were available, or widely used, and physical simulations were run using the standard graphics pipeline [11, 13]. In this work we concentrate on multilayer environments with an arbitrary number of layers simulated using Cuda on Nvidia GPUs. This paper is organized as follows: in Section 2 the mathematical model that has been used is described, Section 3 describes de numerical model, Section 4 addresses the multilayer computational model proposed in this work and Section 5 presents the set of experiments that have been performed. Finally, Section 6 presents the conclusions that were reached. 2 2. Mathematical model Let us consider the system of equations governing the 2D flow of msuperposed immiscible layers of shallow homogeneous fluids with constant densities in a subdomain D⊂R2:                        ∂thl+∂xqx,l +∂yqy,l = 0, ∂tqx,l +∂x q2 x,l hl +1 2gh2 l!+∂yqx,lqy,l hl+ghl∂x X k>l hk+X k<l ρk ρl hk−H!= 0. ∂tqy,l +∂xqx,lqy,l hl+∂y q2 y,l hl +1 2gh2 l!+ghl∂y X k>l hk+X k<l ρk ρl hk−H!= 0 (1) where l= 1,· · · , m being mis the number of layers, hl(x, t) and ql(x, t) = (qx,l(x, t), qy,l(x, t)) are, respectively, the thickness and the mass-flow of the l-th layer at point xat time t, and they are related to the mean velocities ul(x, t)=(ux,l(x, t), uy,l(x, t)), by the equalities: ql(x, t) = ul(x, t)hl(x, t); H(x) is the basin depth measured from a fixed reference level, gis the gravity constant and ρjthe densities of each layer verifying 0< ρ1<· · · < ρm. Notice that h1is the height of the layer of fluid on the top and hmis the height of the layer of fluid over the bottom (See Figure 1). Figure 1: Multilayer stratified flow This system can be written under the structure of a system of balance laws with non-conservative products ∂tw+∂xFx(w) + ∂yFy(w) + Bx(w)∂xw+By(w)∂yw=Gx(w)∂xH+Gy(w)∂yH(2) where w=   w1 . . . wm   , wl=  hl qx,l qy,l  =  hl ux,lhl uy,lhl , 3 Fx(w) =    Fx,1 . . . Fx,m   , Fx,l =     ux,lhl u2 x,lhl+1 2gh2 l ux,luy,lhl     , Fy(w) =    Fy,1 . . . Fy,m   , Fy,l =      uy,lhl ux,luy,lhl u2 y,lhl+1 2gh2 l       Gx(w) =    Gx,1 . . . Gx,m   , Gx,l =  0 ghl 0 , Gy(w) =    Gy,m . . . Gy,m   , Gy,l =  0 0 ghl . Bx(w) is the 3m×3mmatrix defined by Bx(w) =      0B1,2 xB1,3 x· · · B1,m x B2,1 x0B2,3 x· · · B2,m x . . .. . .0· · · . . . Bm,1 xBm,2 xBm,3 x· · · 0      where Bl,k x,k, l = 1,· · · , m,k=lare 3 ×3 matrices defined by Bl,k x=                 0 0 0 ρk/ρlhl0 0 0 0 0  if k < l   0 0 0 hl0 0 0 0 0  otherwise. By(w) is the 3m×3mmatrix defined by By(w) =      0B1,2 yB1,3 y· · · B1,m y B2,1 y0B2,3 y· · · B2,m y . . .. . .0· · · . . . Bm,1 yBm,2 yBm,3 y· · · 0      where Bl,k y,k, l = 1,· · · , m,k=lare 3 ×3 matrices defined by Bl,k y=                  0 0 0 0 0 0 ρk/ρlhl0 0  if k < l  0 0 0 0 0 0 hl0 0  otherwise. Let us define the matrices Aα(w) = Jα(w) + Bα(w), α=x, y where Jα(w) = ∂Fα ∂w (w) are the Jacobians of the fluxes Fα, and we assume that w∈Ω⊂R3m, so that (1) is strictly hyperbolic. 4 Let us also remark that the system (1) verifies the property of invariance by rotations. Effectively, given η∈R2with ∥η∥= 1, let us define Tη=     R1,1 η0· · · 0 0R2,2 η· · · 0 . . .. . .. . .. . . 0 0 · · · Rm,m η     , Rl,l η=  1 0 0 0ηxηy 0−ηyηx , l = 1,· · · , m and let us denote Fη(w) = Fx(w)ηx+Fy(w)ηy,B(w)=(Bx(w), By(w)), and G(w) = (Gx(w), Gy(w)) then Fη(w) = T−1 ηFx(Tηw), TηB(w)·η=Bx(Tηw), TηG(w)·η=Gx(Tηw).(3) Moreover, it is easy to check that Tηwverifies the system ∂t(Tηw) + ∂ηFx(Tηw) + Bx(Tηw)∂ηw=Gx(Tηw)∂ηH−Rη⊥,(4) where Rη⊥=Tη∂η⊥Fη⊥(w) + B(w)·η⊥∂η⊥w−G(w)·η⊥∂η⊥H. Finally, let us remark that there is no explicit formula for the eigenvalues and eigenvectors of matrix Aη(w) = Ax(w)ηx+Ay(w)ηy, but it is well known in oceanography that the fastest and the slowest wave of system (1) could be approximated by the expressions λη L= ¯uη−v u u tg m X l=1 hl, λη R= ¯uη+v u u tg m X l=1 hl(5) where ¯uηis defined by ¯uη=Pm l=1 ul·η hl Pm l=1 hl .(6) 3. Numerical Scheme To discretize (1) the computational domain Dis decomposed into subsets with a simple geometry, called cells or finite volumes: Vi⊂R2. Here, it is assumed that the cells are rectangular with edges parallel to the Cartesian axes. Let us denote by Tthe mesh and by NV the number of cells. Given a finite volume Vi,|Vi|will represent its area; Ni∈R2its center; Nithe set of indexes jsuch that Vjis a neighbor of Vi;Eij the common edge of two neighboring cells Viand Vj, and |Eij |its length (∆x (respectively ∆y) the length of the horizontal (respectively vertical) edges); ηij = (ηij,x, ηij,y) the normal unit vector at edge Eij pointing towards the cell Vjand wn ithe approximation to the average of the solution in the cell Viat time tnprovided by the numerical scheme: wn i∼ =1 |Vi|ZVi w(x, tn)dx. Let us now briefly describe the numerical scheme that we use to approximate the solution of the system (1). Here, we use the natural extension to 2D problems of the PVM-2U scheme introduced in [8] : wn+1 i=wn i−∆t |Vi|X j∈Ni |Eij|F− ij (7) where F− ij is defined taking into account the property of invariance by rotations of system (1). The following steps are performed to define F− ij : 5 1. We define wηij =h1, qηij 1, h2, qηij 2,· · · , hm, qηij m, T=Tηij (w)[1,2,4,5··· ,3m−2,3m−1], and wη⊥ ij =qη⊥ ij 1, qη⊥ ij 2,· · · , qη⊥ ij mT =Tηij (w)[3,6,··· ,3m], where w[i1,··· ,is]is the vector defined from vector w, using its i1-th, . . . , is-th components, qηij l=ql·ηij and qη⊥ ij l=ql·η⊥ ij,l= 1,· · · , m. 2. Let Φ− ηij be the 1D numerical PVM-2U flux associated to the 1D multilayer shallow-water system defined using the 1-st, 2-nd, 4-th and 5-th, · · · , (3m−2)-th and (3m−1)-th equations of system (4) where the term Rη⊥ ij has been neglected: Φ− ηij =1 2(F(wηij j)− F(wηij i) + Bij(wηij j−wηij i)− Gij(Hj−Hi) (8) −Qij((wηij j−wηij i)−(A∗ ij)−1Gij(Hj−Hi)) + FC(wηij i) (9) where F(wηij ) =                   uηij 1h1 (uηij 1)2h1+g 2h2 1 uηij 2h2 (uηij 2)2h2+g 2h2 2 . . . uηij mhm (uηij m)2hm+g 2h2 m                   ,Gij =                 0 ghij 1 0 ghij 2 . . . 0 ghij m                 ,FC(wηij ) =                  uηij 1h1 (uηij 1)2h1 uηij 2h2 (uηij 2)2h2 . . . uηij mhm (uηij m)2hm                  , Bij(wηij j−wηij i) =                     0 ghij 1 m X k=2 (hk,j −hk,i)! 0 ghij 2 m X k=3 (hk,j −hk,i) + 1 X k=1 ρk ρ2 (hk,j −hk,i)! . . . 0 ghij m m−1 X k=1 ρk ρm (hk,j −hk,i)!                     and hij l=hl,j +hl,i 2, l = 1,· · · , m. Matrix Qij is defined by Qij =αij 0Id +αij 1Aij +αij 2A2 ij (10) 6 with αij 0=(SM)2Sm(sign(Sm)−sign(SM)) (Sm−SM)2, αij 1=SM(|SM|−|Sm|) + Sm(sign(SM)Sm−SMsign(Sm)) (Sm−SM)2, αij 2=Sm(sign(Sm)−sign(SM)) (Sm−SM)2, (11) where SM=   Sij Lif |Sij L|≥|Sij R|, Sij Rotherwise (12) and Sm=   Sij Rif |Sij L|≥|Sij R|, Sij Lotherwise. (13) Sij Rand Sij Lare two estimations of the fastest and slowest waves, respectively, related to the Riemann problem associated to edge Eij. Here, the following estimations are used: Sij R= max ληij R,j , λij R, Sij L= min ληij L,i , λij L, being ληij R,j =uηij j+v u u tg m X l=1 hl,j, λij R=uηij ij +v u u tg m X l=1 hij l ληij L,i =uηij i−v u u tg m X l=1 hl,i, λij L=uηij ij −v u u tg m X l=1 hij l where uηij k=Pm l=1 hl,k ul,k ·ηij Pm l=1 hl,k , k =i, j and uηij ij =Pm l=1 uij lhij l Pm l=1 hij l , with uij l=phl,juηij l,j +phl,iuηij l,i phl,j +phl,i , l = 1,· · · , m and uηij l,k =ul,k ·ηij , k =i, j. 7 Aij is the 2m×2mmatrix defined by Aij =                   0 1 0 0 · · · 0 0 ghij 1−(uij 1)22uij 1ghij 10· · · ghij 10 0 0 0 1 · · · 0 0 ρ1 ρ2 ghij 20ghij 2−(uij 2)22uij 2· · · ghij 20 . . .. . .. . .. . .· · · . . .. . . 0 0 0 0 · · · 0 1 ρ1 ρm ghij m0ρ2 ρm ghij m0· · · ghij m−(uij m)22uij m                   (14) and A∗ ij is defined from Aij by setting uij l= 0, l= 1,· · · , m. Taking into account the expression of Qij (10) and A∗ ij, Φ− ηij can be written as Φ− ηij =1 2(Eij −(αij 0e Iij +αij 1Eij +αij 2AijEij)) + FC(wηij i) (15) where Eij =F(wηij j)− F(wηij i) + Bij(wηij j−wηij i)− Gij(Hj−Hi) and f Iij =                  f1 j−f1 i−(f2 j−f2 i) qηij 1,j −qηij 1,i f2 j−f2 i−(f3 j−f3 i) qηij 2,j −qηij 2,i . . . fm j−fm i qηij m,j −qηij m,i                  with fl α= m X k=l hk,α −Hα, α =i, j. Finally, let us remark that Eij can be written in an equivalent way as Eij =FC(wηij j)− FC(wηij i) + Pij 8 where Pij is a discretization of the pressure terms given by Pij =                      0 ghij 1f1 j−f1 i 0 ghij 2ρ1 ρ2(f1 j−f1 i)−(f2 j−f2 i)+ (f2 j−f2 i) . . . 0 ghij m m−1 X k=1 ρk ρm(fk j−fk i)−(fk+1 j−fk+1 i)+ (fm j−fm i)!                      . Equivalently the flux corresponding to cell wjcan be defined as: Φ+ ηij =1 2(Eij + (αij 0e Iij +αij 1Eij +αij 2AijEij)) − FC(wηij j).(16) 3. Let us define Φ− η⊥ ij =           (Φ− ηij )[1] uη⊥ ij ,∗ 1 (Φ− ηij )[3] uη⊥ ij ,∗ 2 . . . (Φ− ηij )[2m−1] uη⊥ ij ,∗ m           ,(17) where uη⊥ ij ,∗ lis defined as follows uη⊥ ij ,∗ l=           qη⊥ ij l,i hl,i if (Φ− ηij )[2l−1] >0 qη⊥ ij l,j hl,j otherwise l= 1,· · · , m. Let us remark that Φ− η⊥ ij is the numerical flux associated to the 1-st, 6-th, · · · , 3m-th equations of system (4) where, again, the term Rη⊥ ij has been neglected. Its derivation has been done following the main ideas of the HLLC method for the shallow-water system introduced in [9] as qη⊥ ij l,l= 1,· · · , m can be seen as a passive scalar that is convected by the flow. Equivalently, the flux corresponding to cell wjcan be defined as Φ+ η⊥ ij =−Φ− η⊥ ij . 9 to an arbitrary number of layers. In the multilayer system the bigger kernel is the one associated to the computation of F− l,ij and λij max which requires 42 registers. It can therefore be concluded that the strategy of a finer-grained design of the multilayer system offers a better performance for this type of problem as it increases the thread level parallelism while achieving a sufficient instruction level parallelism for each kernel instance. This increased performance comes at the cost of a larger memory footprint because of the intermediate data storage requirements. 5.1.1. Comparison with a CPU two-layer system The multilayer code was also compared to a native CPU two-layer implementation. The same 10 seconds of simulation of the dam break mentioned in the previous section with 256 ×256 volumes, requires 1489 seconds on the CPU (Intel Core i7-2600 CPU @ 3.40GHz). Using four OpenMP threads, the run time is 412 seconds. This results in a speed-up of over 250×with respect to the single-threaded CPU version and over 70×with respect to the multithreaded CPU version. This CPU version includes no special optimizations apart from the OpenMP parallelization and the compiler optimizations produced by the -O3 compiler flag. 5.2. Multilayer experiments Two experiments were performed to test the performance of our work. The first experiment used a domain consisting in a 10×10 square. The basin depth is defined by the function: H(x)=1−0.5 exp (−10 (x−5)2) and all layers have a thickness value of 1, except for three layers: •The two layers affected by the presence of the simulated dam. The dam height is 4.5 and the distance to the top neighboring layer is 0.5. •The height of layer mis defined by H(x). The function that defines the thickness of each layer is shown in equation 21. hl(x) =          H(x) if l=m 1 if l= 1,...,(d−2),(d+ 1),...,(m−1) (5.5DA(x)) + (1 −DA(x)) if l=d (0.5DA(x)) + 5 (1 −DA(x)) if l= (d−1) (21) where d=m 2is the layer on which the dam is defined. In this experiment the center of the circular dam is located at c= (2,5). Regarding the initial conditions, ql(x,0) = 0, l = 1,...m. The CFL value used was 0.9 and ρl= 2 ρl−1. Wall boundary conditions were imposed. The test scenario scheme with six layers and basin is shown in Figure 7. There are also video captures of different simulation scenarios available on the complementary material website: http://dicits.ugr.es/ software/MultiLayerSW/ In Table 2 the results for different number of layers and 2D sizes are shown. The number of volumes (millions of) processed per second is also shown considering the number of simulation iterations performed, and the total number of computational volumes. In the first test, the total thickness of the water increases as layers are added, and so does the density difference between the first and the last layer. This results in a decreasing ∆tvalue and therefore an increasing number of iterations as the number of layers increases. For the second experiment, the total thickness of the water was set to 10. This means that as layers are added, the thickness of each layer decreases. Additionally the density values were set so that ρ1= 0.1 and ρm= 1.0 and all the intermediate ρivalues are interpolated lineally. The rest of the experiment parameters are the same as for the previous one. This test scenario produces a fixed ∆tvalue at each iteration step and a constant number of iterations. In this case, as layers as added, the simulation scenario remains fixed but more resolution is added to the dimension represented by the number of layers. The simulation times for the second experiment are shown for a 512 ×512 domain together with the results of the previous experiment in Figure 8. This plot shows that when keeping the number of iterations 16 #Volumes per layer #Layers Time #Vols/s (millions) 256 ×256 4 12.7s 110 256 ×256 8 34.1s 100 256 ×256 16 112.7s 79 512 ×512 4 93.7s 121 512 ×512 8 254.5s 109 512 ×512 16 865.6s 84 1024 ×1024 4 732.7s 125 1024 ×1024 8 1981.6s 113 1024 ×1024 16 6848.1s 86 Table 2: GTX680 - Multilayer experiments. Time expressed in seconds Figure 7: Test scenario with six layers constant, doubling the number of layers does only increase the run time by a factor between 2 and 3. This is a very good result taking into account the inter-layer dependencies become more important as the number of layers increases. A comparison using a Tesla K20 GPU (2495 Cuda cores @ 706Mhz) was also performed. This GPU offers a greater amount of Cuda cores (60% more) than the GTX 680 (1536 Cuda cores @ 1058Mhz) but at a lower clock speed (33% less). Table 3 shows the results obtained. Finally, a single experiment was performed using a GTX Titan GPU (2688 Cuda cores @ 876Mhz). This GPU provides 75% more cores than the GTX680 at a 15% lower clock rate and the results obtained are shown in table 4 #Volumes per layer #Layers GTX 680 time Tesla K20 time 256 ×256 8 34,1s 29s 512 ×512 16 865.6s 802s Table 3: GTX 680-Tesla K20 multilayer experiment. Time expressed in seconds Taking into account the 1.5×speed-up factor obtained using the GTX Titan and the up to 1.11×factor obtained with the Tesla K20 with respect to the results of the GTX 680 model, we believe the system shows a good scalability with respect to the number of available Cuda cores. The important clock rate difference of the K20 model explains that the results obtained with this GPU are worse than the ones that could be expected by only considering its number of Cuda cores. 17 Figure 8: 512 ×512 simulation times #Volumes per layer #Layers GTX 680 time GTX Titan time 256 ×256 6 23s 15,2s Table 4: GTX 680-GTX Titan multilayer experiments. Time expressed in seconds 5.3. Volume reordering We have also analyzed the benefits of the volume reordering scheme with respect to the ad-hoc volume/data structure mapping using the first of the two test scenarios introduced in section 5.2. #Layers Reordered volumes Ad-hoc order 2 6.1s 5.8s 6 25.3s 26.1s 12 75.6s 82.37s Table 5: Volume reordering benefits. 65 535 (256 ×256) volumes per layer. Time expressed in seconds Table 5 shows that the reordering scheme does produce some overhead because it requires a higher amount of operations to compute volume and edge indices and this fact produces a slightly higher run time for a low number of layers. As the number of edges increases, which also increases the number of memory accesses, the benefit of the reordering is clear. 6. Conclusions An efficient GPU based multilayer simulation scheme has been presented which allows an arbitrary number of layers as it is fully parameterized. The fact that it is based on a 2D mathematical model that 18 has been implemented by using a 3D approach, a higher level of thread level parallelism is achieved which offers a good scalability as the number of layers increases. The efficiency of the proposed simulation scheme has been shown in terms of the number of volumes which are processed per second and the speed-up values with respect to previous two-layer systems, using both their GPU and their single and multi-threaded CPU version. The proposed general scheme introduces no overhead when compared to existing GPU-based systems which were specifically designed for two layers. Moreover, it achieves a run time reduction from one to two orders of magnitude when compared to multi and single-threaded CPU implementations. Our simulation software processes around 100 million of 3D volumes per second for scenarios with between 4 and 16 layers on a single GPU. The results show that the simulation scheme that has been presented can be applied to large-scale problems to simulate real-world phenomena like lake or oceanic currents. Acknowledgments The work of M. J. Castro was supported in part by the Spanish Government and FEDER through the Research project MTM2012-38383-C02-01, and by the Andalusian Government through the projects P11-FQM8179 and P11-RNM7069 References [1] M. de la Asunci´on, J.M. Mantas, M.J. Castro, Programming CUDA-based GPUs to simulate two-layer shallow water flows, in: P. D’ambra, M. Guarracino, D. Talia (Eds.), Euro-Par 2010, volume 6272 of Lecture Notes in Computer Science, Springer, 2010, pp. 353–364. [2] M. de la Asunci´on, J.M. Mantas, M.J. Castro, Simulation of one-layer shallow water systems on multicore and CUDA architectures, Journal of Supercomputing 58 (2011) 206–214. [3] M. de la Asunci´on, J.M. Mantas, M.J. Castro, E.D. Fern´andez-Nieto, An MPI–CUDA implementation of an improved roe method for two-layer shallow water systems, Journal of Parallel and Distributed Computing. Special Issue on Accelerators for High–Performance Computing 72 (2012) 1065–1072. [4] A.R. Brodtkorb, T.R. Hagen, K.A. Lie, J.R. Natvig, Simulation and visualization of the Saint-Venant system using GPUs, Computing and Visualization in Science 13 (2010) 341–353. [5] A.R. Brodtkorb, M.L. Sætra, Shallow water simulations on multiple GPUs, in: PARA 2010: State of the Art in Scientific and Parallel Computing, Reykjavik (Iceland), pp. 56–66. [6] M.J. Castro, E.D. Fern´andez-Nieto, A.M. Ferreiro, J.A. Garc´ıa, C. Par´es, High order extensions of Roe schemes for two-dimensional nonconservative hyperbolic systems, Journal of Scientific Computing 39 (2009) 67–114. [7] P. Degond, P.F. Peyrard, G. Russo, P. Villedieu, Polynomial upwind schemes for hyperbolic systems, Comptes Rendus de l’Acad´emie des Sciences - Series I - Mathematics 328 (1999) 479 – 483. doi:http://dx.doi.org/10.1016/S0764-4442(99) 80194-3. [8] M.J.C. D´ıaz, E.D. Fern´andez-Nieto, A class of computationally fast first order finite volume solvers: PVM methods., SIAM J. Scientific Computing 34 (2012). [9] E.D. Fern´andez-Nieto, D. Bresch, J. Monnier, A consistent intermediate wave speed for a well-balanced HLLC solver, Comptes Rendus Mathematique 346 (2008) 795 – 800. doi:http://dx.doi.org/10.1016/j.crma.2008.05.012. [10] M. Geveler, D. Ribbrock, S. Mallach, D. G¨oddeke, A simulation suite for Lattice-Boltzmann based real-time CFD applications exploiting multi-level parallelism on modern multiand many-core architectures, Journal of Computational Science 2 (2011) 113–123. [11] T.R. Hagen, J.M. Hjelmervik, K.A. Lie, J.R. Natvig, M.O. Henriksen, Visual simulations of shallow-water waves, Simulation Modelling Practice and Theory 13 (2005) 716–726. [12] Khronos OpenCL Working Group, The OpenCL Specification, Accessed November 2012. [13] M. Lastra, J.M. Mantas, C. Ure˜na, M.J. Castro, J.A. Garc´ıa, Simulation of shallow-water systems using graphics processing units, Mathematics and Computers in Simulation 80 (2009) 598–618. [14] NVIDIA, Cuda, http://www.nvidia.com/object/cuda_home_new.html, Accessed 2013. [15] NVIDIA, Nvidia, http://www.nvidia.com/, Accessed 2013. [16] NVIDIA, Kepler-The world’s fastest, most efficient HPC arquitecture, http://www.nvidia.com/object/nvidia-kepler. html, Accessed August 2012. [17] C. Par´es, Numerical methods for nonconservative hyperbolic systems: A theoretical framework, SIAM Journal on Numerical Analysis 44 (2006) pp. 300–321. [18] C. Par´es, M.J. Castro, On the well-balance property of Roe’s method for nonconservative hyperbolic systems. Applications to shallow-water systems, ESAIM: Mathematical Modelling and Numerical Analysis 38 (2004) 821–852. 19