Repositorio Institucional de Documentos

EINA Universidad Zaragoza Máster Universitario en Meánia Apliada Trabajo Fin de Máster 2D Shallow Flow Simulation Using GPU Technologies Asier Lacasta Soto Grupo de Hidráulica Computacional - [email protected] Directora: Pilar García Navarro1 Co-Director: Javier Murillo Castarlenas1 (1) Area de Mecánica de Fluidos Escuela de Ingeniería y Arquitectura Universidad Zaragoza Curso 2011/2012 2 Aknowledgements I would like to express my appreiation to Dra. Pilar Garía Navarro for her valuable and onstrutive suggestions during the planning and development of this researh work. My grateful thanks are also extended to Dr. Javier Murillo for his help with the mo del solving my doubts. I would also like to thank the other memb ers of the group for their advies and p oints of view when neessary. I would like sp eially mention to Hetor Ratia. Finally, I wish to thank nVidia for their partial supp ort providing us with a Hardware part under the Nvidia Aademi Partnership program. This work has b een develop ed under pro jet CENIT-TECOAGUA CEN-20091028. i ii Resumen Los mo delos matemátios y méto dos numérios impliados en la simulaión de ujos on sup er- ie libre han sido estudiados durante tiemp o en el Grup o de Hidráulia Computaional de la Universidad de Zaragoza. Estos mo delos son la base de nuevos desarrollos omo el transp orte de sedimento, el mo delado de interaión on puentes o el aoplamiento hidrológio. A p esar de la alidad de estos méto dos, el oste omputaional es muy alto y en gran parte esto se deb e a la tenología numéria que requieren. Con la nalidad de sup erar esta limitaión, este traba jo estudia la implementaión de un ó digo de simulaión hidráulia orientada a ejeuión en GPU, p ermitiendo simular un amplio onjunto de situaiones transitorias en gran esala temp oral, on un tiemp o de simulaión razonable. El oste omputaional de éste tip o de herramientas ha sido reduido, tradiionalmente, utilizando ténias de paralelismo, impliando un alto número de pro esadores para reduir el tiemp o de álulo al máximo. En los últimos años, las freuenias de los pro esadores pareen hab er alanzado su límite (Figura 1 extraida de [9℄) p or lo que las ténias de paralelismo en pro esadores masivos son una nueva op ión. Figure 1: Evoluión de las freuenias de CPU desde 1985 hasta 2011 iii En este traba jo, se analiza el rendimiento del ó digo implementado en GPU, omparándolo on su equivalente en CPU. Este segudo, viene siendo desarrollado, en su totalidad, en Fortran mientras que el primero, ha sido desarrollado utilizando el lengua je de programaión C, ompartiendo el pro esamiento geométrio on la versión CPU. Las fuionalidades implementadas en la versión GPU, ubre una gran parte de situaiones de interés, tales omo el avane de una inundaión, los ambios de fondo y friión y algunas ondiiones de ontorno de entrada y de salidas. La implementaión del méto do en GPU no es trivial y requiere de un ono imiento en profundidad del funionamiento de esta tenología a ba jo nivel. Los b eneios de la versión GPU serán analizados a través de la aeleraión rep eto a la versión CPU en diferentes tip os de aso. EL rendimiento del ó digo GPU además, será medido teniendo en uenta el uso de mallas no estruturadas, las uales suelen ser neesarias en muhos o digos de CFD. Para su simulaión, se utilizará la GPU Tesla 2075 de nVidia. Además se utilizará el estándar CUDA, que hae la programaión más senilla que otros estándar en programión GPU, p ermietiendo al programador exprimir los b eneios de esta tenología. iv Abstrat The mathematial mo dels and numerial metho ds implied in the resolution of free surfae ows have b een studied for a long time within the Computational Hydrauli Group at the Universidad Zaragoza. They supp ort new developments suh as sediment transp ort, bridges mo deling or hydrologial oupling. Despite the quality that the numerial solvers prop osed by the group oer, the omputational ost of these metho ds is very high, due to the omplexity of the numerial to ols required. In order to avoid this limitation, the present work studies the implementation of a sienti hydrauli simulation to ol oriented to b e run on GPU, allowing to simulate a wide range of situations over large time sale problems, that otherwise an not b e omputed at an aordable ost. The omputational ost has b een traditionally redued by using parallel tehniques, involving a large numb er of pro essors in order to redue the simulation time as muh as p ossible. Sine CPU frequenies seem to b e reahing their maximum apaity (Figure 2 extrated from [9℄), nowadays Many-Core parallel tehniques app ear to b e an interesting option. Figure 2: CPU Frequeny evolution sine 1985 until 2011 The p erformane of the GPU version is analyzed omparing b oth CPU and GPU versions of the same o de. While the former was fully develop ed in Fortran language, the numerial v kernel of the new GPU version has b een written in C, sharing the geometrial prepro essing mo dule with the CPU version. The funtionalities implemented in the GPU version over a wide range of situations as they inlude all the harateristis that are desirable in the ontext of shallow ow simulation: o o ding advane, frition and b ed slop e soure-terms as weel as inlet and outlet b oundary onditions. The implementation of these requirements in the ontext of realisti simulations is not straightforward. This is explained when onsidering that, ontrary to other programming languages, the GPU version requires a go o d omprehension of the low level op erations, that do es not allow a diret onventional implementation. The b enets of the GPU version will b e analyzed in depth fo using on sp eed-up gain in omplex ases. The p erformane of the GPU o de is analyzed in depth to ensure not only the eieny but also the p ossibilities of GPU programming when using unstrutured meshes, that are often required in CFD o des. A Tesla 2075 nVidia GPU has b een used in the present study. Moreover, it has b een develop ed using nVidia-CUDA standard, whih makes friendly the programming for general purp ose appliations, allowing the programmer to exploit the many-ore paradigm. vi Contents 1 Intro dution 1 1.1 Context and assumptions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.2 Struture of the rep ort . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 2 Mathematial Mo del and Numerial Metho d 3 2.1 Approximate Riemann solution . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 2.2 Appliation to the 2D Shallow Water equations . . . . . . . . . . . . . . . . . . . 6 2.3 Numerial resolution . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7 3 CUDA Tehnology Overview 11 3.1 GPU Tehnology history . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11 3.2 nVidia CUDA tehnology . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12 3.3 CUDA development . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13 3.3.1 Example of implementation in a 1D ase . . . . . . . . . . . . . . . . . . . 14 3.4 Results Al l that glitters is not gold . . . . . . . . . . . . . . . . . . . . . . . . . . 16 4 Implementation 19 4.1 Mo del overview . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19 4.2 Memory oalesing . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21 4.3 Gathering data avoiding b ottlenek . . . . . . . . . . . . . . . . . . . . . . . . . . 23 4.4 Writing output les . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24 4.5 Compilation and other issues . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26 5 Results 29 5.1 Preision: A test-ase with analytial solution . . . . . . . . . . . . . . . . . . . . 29 5.2 Performane: A large-sale simulation at Júar River . . . . . . . . . . . . . . . . 33 5.3 Comparing with a distributed memory parallel implementation . . . . . . . . . . 39 6 Conlusions and future work 45 Bibliography 47 vii 2.1. APPROXIMATE RIEMANN SOLUTION - - 6 ?@ @ @ @ @ @ @ @           HHHHHHHH H      @ @ @ @ @ x′ lk 0 Un i Un j nk Figure 2.1: Riemann problem in 2D along the normal diretion to a ell side. where following previous work, [28℄ Z+X′ −X′ S(x′,0) dx′= (Tn)n k (2.13) where T is a suitable numerial soure matrix. This enables the following formulation of (2.8) ∂ ∂t ZΩ UdΩi+ NE X k=1 (δE−T)knklk= 0 (2.14) that is approximated by using the following linear problem ∂ ∂t RΩ ˆ UdΩi+PNE k=1 J∗ n,kδˆ Uklk= 0 ˆ U(x′,0)k=(Uiif x′<0 Ujif x′>0 (2.15) Integrating 2.15 over the same ontrol volume as b efore the following expression is obtained for eah k edge Z+X′ −X′ ˆ U(x′,1) dx′=X(Ui+Uj)−J∗(Uj−Ui) (2.16) and sine we want to satisfy (2.12), the onstraint that follows is: (δE−T)knk=˜ J∗(Uj−Ui) (2.17) Due to the non-linear harater of the ux matrix E , the denition of an approximated Jaobian matrix, e Jn,k , allows for a lo al linearization δ(En)k=e Jn,kδUk (2.18) and is exploited here [24℄. This approah provides a set of three real eigenvalues e λm k and eigenve- tors eem k . Then, it is p ossible to dene two matries e P= (ee1,ee2,ee3) and e P−1 with the following prop erty e Jn,k =e Pke Λke P−1 k (2.19) 5 CHAPTER 2. MATHEMATICAL MODEL AND NUMERICAL METHOD The dierene in vetor U aross the grid edge and the soure term are pro jeted onto the matrix eigenvetors basis: δUk=e PkAk(Tn)k=e PkBk (2.20) with Ak=α1α2α3T k and Bk=β1β2β3T k . Expressing all terms more ompatly: δ(E·n)k−(T·n)k= Nλ X m=1 e λ θαeem k (2.21) with θm k=1−β e λαm k (2.22) Finally, it is p ossible to dene the desired matrix in (2.17) e J∗ k= (e Pe Λ∗e P−1)k (2.23) with e Λ∗=e ΛΘ , where e Λk is a diagonal matrix with eigenvalues e λm,∗ k in the main diagonal and Θk is a diagonal matrix with θm k in the main diagonal: e Λk=  e λ10 0 0e λ20 0 0 e λ3  k Θk=   θ10 0 0θ20 0 0 θ3  k (2.24) 2.2 Appliation to the 2D Shallow Water equations The two-dimensional shallow water equations, whih represent depth averaged mass and momentum onservation, an b e obtained from the Navier-Stokes equations. Negleting diusion of momentum due to visosity and turbulene, wind eets and the Coriolis term, they form a system of equations [2℄ as in (2.1), where U= (h, qx, qy)T (2.25) are the onserved variables with h representing the water depth, qx=hu and qy=hv , with (u, v) the depth averaged omp onents of the velo ity vetor u along the (x, y) o ordinates resp etively. The uxes of these variables are given by: F=qx,q2 x h+1 2gh2,qxqy hT ,G= qy,qxqy h q2 y h+1 2gh2!T (2.26) where g is the aeleration of the gravity. The soure terms of the system are the b ed slop e and the frition terms: S=0,pb,x ρw−τb,x ρw ,pb,y ρw−τb,y ρwT (2.27) 6 2.3. NUMERICAL RESOLUTION where the b ed slop es of the b ottom level z are pb,x ρw =−gh∂z ∂x,pb,y ρw =−gh∂z ∂y (2.28) and the frition losses are written in terms of the Manning's roughness o eient n : τb,x ρw =ghSfx Sfx =n2u√u2+v2 h4/3,τb,y ρw =ghSfy Sfy =n2v√u2+v2 h4/3 (2.29) 2.3 Numerial resolution Following Go dunov's metho d, the solutions of the RP's are evolved for a time equal to the time step and the resulting solution is ell-averaged. The volume integral in the ell at time tn+1 leads to the up dating numerial sheme as: Un+1 iAi=Un iAi− NE X k=1 3 X m=1 (e λ−θαee)m klk∆t (2.30) with e λ±,m k=1 2(e λ±|e λ|)m k . When applied to the shallow water system presented in setion 2.2 the approximate Jaobian e Jn,k for the homogeneous part is onstruted with the following averaged variables [24℄ euk=ui√hi+ujphj √hi+phj ,evk=vi√hi+vjphj √hi+phj ,eck=rghi+hj 2 (2.31) leading to e λ1 k= (e un −ec)k,e λ2 k= (e un)k,e λ3 k= (e un +ec)k (2.32) and ee1 k=   1 eu−ecnx ev−ecny   k ,ee2 k=   0 −ecny ecnx   k ,ee3 k=   1 eu+ecnx ev+ecny   k (2.33) When ell averaging the solution in the 1D dimensional ase the time step ∆t is taken small enough so that there is no interation of waves from neighb ouring Riemann problems, attending to a distane ∆x/2 . In the 2D framework, onsidering unstrutured meshes, the equivalent distane to ∆x , that will b e referred to as χi in eah ell i must onsider the volume of the ell and the length of the shared k edges. χi=Ai maxk=1,NE lk (2.34) Considering that eah k RP is used to deliver information b etween eah pair of neighb ouring ells of dierent size, the asso iated distane min(Ai, Aj)/lk is relevant, so in ase that ˆ h(x′, t)≥0 in all k RP's the time step is limited by ∆t≤CFL ∆t e λ∆t e λ=min(χi, χj) maxm=1,2,3|e λm| (2.35) 7 CHAPTER 2. MATHEMATICAL MODEL AND NUMERICAL METHOD The previous stability ondition is insuient in presene of relatively imp ortant soure terms. The systemati ontrol of numerial stability in those ases has b een a matter of reent researh in the group as it is related with the appliability of the sheme to real situations. A simple generalization of the CFL ondition paying attention to the existene of the soure terms an lead to extremely small values of ∆t various orders of magnitude smaller than the value ditated by the homogeneous ondition, hene rendering the metho d impratial. This an b e avoided by means of a reonstrution of the approximate solution ˆ U(x′, t) that is not detailed here for the sake of oniseness. The strategy prop osed is based on enforing p ositive values of auxiliary quantities h∗ i h∗ i=hn i+α1 k−β e λ1 k≥0 (2.36) and h∗∗∗ j h∗∗∗ j=hn j−α3 k+β e λ3 k≥0 (2.37) so that, when they b eome negative, the numerial soure term is redued instead of reduing the time step size. For more details, see [21, 18℄. Furthermore, following the unied disretization in [6℄ the non-onservative term (Tn)k in (2.13) at a ell edge is written [20 ℄ as: (Tn)k=    0 pb ρw−τb ρwnx pb ρw−τb ρwny    k (2.38) where pb ρw and τb ρw attend to the pressure and frition exerted on the b ed resp etively. In this work the following expression for the thrust term pb ρw is prop osed: pb ρwk =     max pb ρwa,pb ρwbk if δd δz ≥0 and (e un)δz > 0 pb ρwb k otherwise (2.39) where d= (h+z) and pb ρwa k =−g(ehδz)kpb ρwb k =−ghr−|δz′| 2δz′ (2.40) with r=(i if δz ≥0 j if δz < 0δz′=     hi if δz ≥0 and di< zj hj if δz < 0 and dj< zi δz otherwise (2.41) The disretization of the frition term based on [21℄ is applied τb ρwk =g(ehSf)kdnSf,k =n2e un|e u| max(hi, hj)4/3k (2.42) 8 2.3. NUMERICAL RESOLUTION with dn the normal distane b etween neighb or ell enters. 9 10 3 CUDA Tehnology Overview Nowadays, GPU tehnologies start to onquer from ordinary business appliations to sieniti appliations. This general purpose orientation is denomined GPGPU 1 , allowing its develop ers to reah higher p erformane than in oventional arhitetures (Single Instrution Single Data) where the op erations are urrently p erformed sequentially. In the ase of sienti omputation, the GPGPU paradigm p erforms the numerial metho ds. nVidia has b een working in the improvement of the GPGPU paradigm, reating the CUDA to olkit. CUDA to olkit is a parallel arhiteture for graphi pro essing whih implements an intrution-set oriented to the GPU memory aess and op erations in C. Other more general implementations have b een p erformed through op en-soure platforms suh as Op enCL and others like PGI-Cuda as propietary-soure. Op enCL has the main advantage of b eing hardwareindep endent. It implies that the same o de ould b e exeuted on b oth nVidia and ATI GPUs. The main disadvantage is that the learning-urve is harder than for the CUDA to olkit. The other option is PGI-Cuda. It has the main advantage in the supp ort of CUDA primitives for Fortran but the disadvantage is the ost of it. So, as we are interested in simulating at nVidia GPUs, the implementation of the o de has b een develop ed using CUDA-To olkit. 3.1 GPU Tehnology history Sine the advent of Op enGL, GPUs added programmable shading to their apabilities. Eah pixel ould inop orate its pro essing as a program to b e shown on sreen after applying it. nVidia was the rst to pro due a hip apable of programmable shading. In 2002, ATI develop ed the rst Diret3D 9.0 aelerator, whih implemented lo oping and lengthy oating p oint math, b e- oming as exible as CPU and orders of magnitude faster for image-array op erations. Abstrating the graphial purp ose and taking a double-p oint array as if it were a vertexarray, the same op erations were able to b e applied, so with the nVidia CUDA To olkit, a new programming mo del for GPU omputing was stablished. After its app earane, Op enCL b eame broadly supp orted allowing develop ers o ding for AMD/ATI GPUs. 1 General Purp ose Graphi Pro essor Unit 11 CHAPTER 3. CUDA TECHNOLOGY OVERVIEW 3.2 nVidia CUDA tehnology The present work has b een develop ed using an nVidia Tesla GPU. The partiular organization and how it works is explained b elow and has followed [11 ℄. Most of the details are ommon with the previous GPU generations and it is previsible that will b e ommon with future generations to o. There are two main p oints of view when explaining how CUDA works. The rst is based on the hardware arhiteture. The minimum unit is the Streaming Pro essor (SP), where a single thread is exeuted. A group of SP's form the Streaming Multipro essor (SM), tipially with 32 SP's. Finally, a GPU is omp osed by b etween 2 and 16 SM's. The seond p oint of view is based on the way CUDA appliations are develop ed. The minimum unit is alled Thread. Threads are identied by lab els ranging b etween 0 and blokDim . The group of Threads is alled Blo k, and it ontains a (reommended) 32 multiple numb er of Threads. Finally any group of Blo ks is alled Grid. These elements are illustrated on Figure 3.1. Block 0 Block 1 Block 2 Thread Block Grid Figure 3.1: thread, blok, grid sheme omp osition Atual nVidia GPU's p erforms the threads sheduling inside the SM in groups of 32 alled Warps (we also reommend [15℄ for future onsiderations). Eah SM features two Warp shedulers and two instrution dispath units, allowing two Warps to b e issued and exeuted onurrently. Fermi's dual Warp sheduler selets two Warps, and issues one instrution from eah Warp to a group of sixteen ores, sixteen load/store units, or four SFU's. Beause Warps exeute indep endently, Fermi's sheduler do es not need to hek for dep endenies from within the instrution stream. Using this elegant mo del of dual-issue, Fermi ahieves near p eak hardware p erformane. Most instrutions an b e dual issued; two integer instrutions, two oating instrutions, or a mix of integer, oating p oint, load, store, and SFU instrutions an b e issued onurrently. Double preision instrutions do not supp ort dual dispath with any other op eration. Figure 3.2 shows how the SP are distributed inside the SM and how the multipro essors are distributed inside the GPU. Furthermore, Figure 3.3 shows the temp oral evolution inside the SM and how it works for a blo k with 256 elements ( warp = 256/32 = 8 elements). Any Thread an b e lab elled using blokDim , blokId and threadId . In an example with 14 Blo ks and 256 Threads/Blo k (3584 elements), we nd that for element 23 in Blo k 4, the lab els inside the o de are 12 3.3. CUDA DEVELOPMENT Streaming Multiprocessor Instruction Cache Register File WS/DU WS/DU Interconnected Network Shared Memory/L1 Cache Uniform Cache (a) GF100 Streaming Multipro essor (SM) Streaming Multiprocessor Streaming Multiprocessor Streaming Multiprocessor Streaming Multiprocessor Streaming Multiprocessor Streaming Multiprocessor Streaming Multiprocessor Streaming Multiprocessor Streaming MultiprocessorStreaming Multiprocessor Streaming MultiprocessorStreaming Multiprocessor Streaming MultiprocessorStreaming Multiprocessor L2 Cache Host Interface / GigaThread Engine Memory Controller Memory Controller (b) 14-SM based Fermi Arhiteure detail Figure 3.2: Desription of our Fermi 2075 GPU based on GF100/GF110 Arhiteture. Warp Scheduler Warp Scheduler Inst. Disp. Unit Warp 8 Instruction 5 Warp 2 Instruction 1 Warp 6 Instruction 17 Warp 8 Instruction 3 Warp 2 Instruction 2 ... ... Block 4 Warp 5 Instruction 14 Warp 7 Instruction 5 Warp 1 Instruction 2 Warp 5 Instruction 15 Warp 1 Instruction 3 Warp 7 Instruction 6 ... Inst. Disp. Unit Time Warp 6 Instruction 16 Figure 3.3: Exeution pip eline for a Strem Multipro essor (left) whih pro ess blo k number 4 (right) • blokDim=256 • blokId=4 • threadId=23 and then, the typial aess pattern, p oints to i=threadId+blokDim*blokId=23+256*4=1047 3.3 CUDA development The CUDA main funtions are related to the memory interation b etween CPU and GPU, in partiular, udaMempy with the dierent ags to stablish the way of the transfer. It is imp ortant to remark that these interations or data transfers b etween GPU and CPU are extremely slow and should b e minimized. Moreover, the allo ation and memory freeing op erations ould b e p erformed using their equivalenes in CUDA as shown in listing 3.1 13 CHAPTER 3. CUDA TECHNOLOGY OVERVIEW Listing 3.1: CUDA Most imp ortant funtions 1 // GPU Memory alloation 2 udaMallo(...,size); 3 // GPU Memory free 4 udaFree(..); 5 // Copy Host To devie 6 udaMempy(...,udaMempyHostToDevie); 7 // Copy Devie To Host 8 udaMempy(...,udaMempyDevieToHost); 9 // Copy Devie To devie 10 udaMempy(...,udaMempyDevieToDevie); The advantage of using GPU for programming numerial metho ds, omes from the High- Level Single Instrution Multiple Data (SIMD) or as nVidia alls, Single Intrutions Multiple Threads (SIMT) paradigm. Any op eration an b e exeuted in onurrene with many others allowing any CUDA Thread to aess to a partiular p osition while any other is aessing to another one. 3.3.1 Example of implementation in a 1D ase Consider, for example, the 1D transp ort equation: ∂u ∂t +c∂u ∂x = 0 (3.1) with c > 0 , and its initial and b oundary onditions u(x, 0) = f(x) u(0, t) = U0 applying the temp oral disretization with forward Euler and the upwind sheme: ∆ui ∆t=−ui−ui−1 δx (3.2) writing its as un+1 i=un i−un i−un i−1 δx ∆t·c (3.3) and the pro edure ould b e written in Standard C as follows Listing 3.2: Simple 1D transp ort equation in C 1 void upwindStepCPU(double *fn,double *fnmas1,double DELTAX){ 2 int i; 3 for (i=1; i<1/DELTAX; i++) { 4 fnmas1[i℄=fn[i℄+*DELTAT*(fn[i-1℄-fn[i℄)/(DELTAX); 5 } 6 } 14 4.2. MEMORY COALESCING 1 3 97 2 6 8 4 11 4 1 8 5 2 3 67 10 14 17 21 24 20 23 19 22 18 15 12 9 13 16 5 Figure 4.2: Strutured mesh with Cell Numb ering detail (Right) and Wall Numbering detail (Left) example 4.2 Memory oalesing Memory oalesing is the way the memory is ordered allowing half-Warp to aess global memory at the same time (using only 1 yle to p erform the load op eration). This means that, if a Thread (The rst one) in a Warp aesses to a partiular memory address and it aess pattern is suh that aess to the next address ( i, i+1, i+2.... ) the following 31 Threads do not need to read the memory again. Otherwise, two or more aesses are needed to allow eah Thread the aess to data. Memory oalesing is one of the most imp ortant things to take into aount when programming GPU's. Reent works [29℄ have demonstrated the eieny of oalesing tehniques, b eing this implementation b etter in some ases than shared memory strategies. Although there exist works dealing with the prots of using this strategy, the way to pro eed when using unstru- tured meshes is not lear. This topi will b e disussed in the next May 2012 GPU Tehnology Conferene [8℄ and some improvements are detailed in [25℄. In our ase, the p erfet memory oalesing tehnique ould b e implemented, [5 ℄ [7℄, if using strutured meshes. As it app ears in Figure 4.2, ell lab elling implies that the aess pattern for a Blo k of (in this ase 9) ells allows the programmer to make the p erfet math aess into a Warp. In other words, for any group of ells within a Warp, all the variables are aessible in only a oalesed reading. Being the present work oriented to a general implementation of the nite volume sheme on b oth strutured and unstrutured grids, the memory optimization is not as easy as desrib ed ab ove. Aording to the general up dating formula 2.30, this sheme works with the ell edge uxes or inter-ell elements through whih the Rienmann Problem is solved. In the ase of the strutured mesh, this ux takes plae into the left, right, upside and downside ell to a given ell, so all the op erations ould b e p erformed lo oping by ells. In unstrutured grids, this onept is dierent 21 CHAPTER 4. IMPLEMENTATION ... ... ... ... 32 Threads Warp Cell Data Array ... ... ... ... 32 Threads Warp Cell Data Array ... ... ... ... 32 Threads Warp Cell Data Array ... ... ... ... 32 Threads Warp Cell Data Array Time Figure 4.3: Misaligned and Coalesed aess pattern to ompute the ux variation for any group of elements following the sheme of Right, Left, Down, Up for W data (Stored by ell) in a mesh ordered as Figure 4.2. Light oloured orresp ond to the pro essed element 5, wih implies ells 2, 4, 6 and 8. 16 74 8 61 23 9 75 12 11 16 33 9 24 21 17 20 14 10 45 31 87 56 8 22 18 49 2359 27 30 Figure 4.4: Unstrutured mesh with Cell Numb ering detail (Right) and Wall Numb ering detail (Left) example and it is go o d idea to make the ux alulations by walls and then, to assign them to eah ell with the need to keep trae via a onnetivity matrix. For the general unstrutured ase it is imp ortant to deide how to stablish the order of the variables. It an b e p erformed through ells or through walls. Using as example Figure 4.2, the op erations of applying the variation to the ell (8) has no a oalesed pattern. There exists the need of searhing the neighb ouring ells (74,61,16) and alulating the ux through walls (33,9,16). Skething these op erations in an example for wall 33 (i=33, 1=8, 2=74) we have: Listing 4.2: Aess pattern for the main ux variation op eration. 1 alulateWallFluxes(...){ 2 // Loop by wall 3 int i = threadIdx.x+(blokIdx.x*blokDim.x); 4 if(i<nWall){ 22 4.3. GATHERING DATA AVOIDING BOTTLENECK ... ... ... ... ... ... ... ... ... ... ... ... 32 Threads Warp Wall Neighbouring Vector Cell Data Array Figure 4.5: Unoalesed aess pattern to get W data (Stored by ell). Pro essing wall 9 is light oloured when it aesses to ell 8 (i=9, 1=8) 5 1=wall[i℄; 6 2=wall[i+1℄; 7 // [COALESCED℄ Aess to the variables of the wall 8 // Normal Vetor 9 // Length of the side 10 // ... 11 ... 12 // [UNCOALESCED℄ Aess to the variables of 1 and 2 13 // Primitives variables 14 // Area of the ells 15 // ... 16 ... 17 // Store the value of the flux for the wall i 18 } 19 } Although this is the main funtion where the ux is alulated and it involves many unoalesed aesses to the variables, there are some op erations whose aess ould b e p erformed through the ells. 4.3 Gathering data avoiding b ottlenek One of the troubles when trying to make all the op erations inside the GPU is the identiation of global quantities suh as the minimum value of a vetor. As the Many-Core paradigm is not designed to share information b etween elements, redution op erations like min, max, sum... are p erformed at ublas library [23℄. ublas library has high-level funtions that work retrieving results to GPU or to CPU. When interested in using them without taking out the data from the GPU, that must b e sp eied. This ould b e done through ublasSetPointerMode_v2(handle, CUBLAS_POINTER_MODE_DEVICE) , stating that all results have to b e returned to the GPU memory. In our ase, it is essential that the algorithm alulates the minimum ∆t following the CFL ondition when running along all the ell edges. Then, following Figure 4.3 sheme, the minimum among all of them is seleted. Details are shown in Listing 4.3. 23 CHAPTER 4. IMPLEMENTATION 1 2 3 n-2 n-1 n 0.13 0.45 0.05 0.62 0.78 0.11 cublasIdamin() 3 dt[1..n] Δt=dt[3] Figure 4.6: Gathering minimum ∆t for all the domain Listing 4.3: Gathering ∆t op eration 1 __global__ void newDt(double *dt,double *vDt, int *id){ 2 // As ublasIdamin returns it value following 3 // 1-based indexing, we must to substrate 1 4 *dt=vDt[*id-1℄; 5 } 6 .... 7 ublasIdamin(handle,*npared,vDt,1,id); 8 newDt<<<1,1>>>(dt,vDt,d_id); 9 ... While alulation is ontrolled by host, it is neessary to transfer the up dated tn+1 . After δt is alulated, the up dating op eration an b e p erform as 4.4 and then, you an transfer the up dated value of tn+1 to CPU. Listing 4.4: Up dating ∆t 1 __global__ void updateT(double *dt,double *t){ 2 int i; 3 *t=*t+*dt; 4 } In order to alulate the global mass error, there is a sum of mass inside the mesh and the balane b etween the inlet and outlet b oundaries M =ρXhiAi (4.1) and then, it alulates the error as ǫ= M n+1 − M n+ M in − M out M n+1 (4.2) The sums are p erformed using ublasDasum where all elements are added within a vetor and the results stored in a variable, working similar to ublasIdamin . 4.4 Writing output les The feature of the newest CUDA mo dels allowing for simultaneous exeution and opy streams an b e used to hide delays aused by writing data to disk. 24 4.4. WRITING OUTPUT FILES Figure 4.7: Asynhronus dumping data diagram. Traditional udaMempy p erforms a synhronous opy, i.e., the all do es not return until the opy is omplete. However, alls to the new family of asynhronous funtions like udaMem- pyAsyn may return b efore the opy is omplete. Furthermore, the opy may b e assigned to a stream. In this way it is p ossible for the CPU host o de to all udaMempyAsyn and assign it to a opy stream, then launh kernels in an exeution stream. Both streams are pro essed simultaneously by the GPU. It is not p ossible to use udaMempyAsyn diretly to opy simulation results to Host memory in the ase of shallow ow simulation b eause the onurrent simulation would alter the values in the variables b eing opied. It is neessary to make a synhronous opy to a buer in GPU memory rst (Figure 4.7). One the opy of the results to the buer is omplete, a all to udaMempyAsyn is made whih opies the buer to host memory, and the simulation kernels are launhed simultaneously op erating on the usual variables. 25 CHAPTER 4. IMPLEMENTATION This sheme requires that the CPU launhes kernels after the all to the asynhronous opy. It is neessary to intro due a parallel CPU thread that waits for the opy to nish and then writes the results to disk. Thus, the main CPU thread will rst all udaMempyAsyn, then spawn a write thread and ontinue launhing kernels to advane the simulation. The rst task for the writing CPU thread will b e to wait for the opy stream to nish, then pro eed to write the results in host memory to disk. The limitation in this sheme is that the omputation time b etween dumps to disk has to b e greater than the writing time to disk itself. If that is not the ase, gains an still b e ahieved from using this sheme but further barriers are required. One of them is that the main CPU thread has to wait for the writing thread to nish b efore alling udaMempyAsyn. Dep ending on the problem, further gains an b e made e.g. using multibuering. 4.5 Compilation and other issues In the original Fortran version of the o de there are several funtions related to the prepro ess and p ostpro ess as skethed on gure 4.8. To b e more eient, the programming of that part of the o de in C has b een ommitted and the work has fo used on the eient programming of the numerial asp ets. So the prepro ess is p erformed through the Fortran version and the omputing kernel is p erformed using C/CUDA. To work with the two o des at the same time, they have b een ompiled together. The tehnique used is based on making a standard C interfae whih interop erates with CUDA and is alled from Fortran as shown in [1℄. The most ompliated and interesting detail of this op eration is the way of ompiling them. It is shown in Listing 4.5. Listing 4.5: Makele Sript 1 2 NVCC = nv 3 FORT = gfortran 4 5 FORTFLAGS = -w -O3 6 CUFLAGS = -g -w -O3 -m64 -arh sm_21 -Xptxas -dlm=a -I$(EXTRAE_HOME)/inlude 7 LDFLAGS = -L/opt/uda/4.0/lib64 -L$(EXTRAE_HOME)/lib -ludatrae -luda -ludart -lstd++ -lublas -lrt -lm -lpthread 8 OBJ = uda_bloks2mf.o SFS2Dv01_64.o 9 BIN = sfsGPU 10 11 $(BIN): $(OBJ) 12 $(FORT) $(FORTFLAGS) $(OBJ) $(LDFLAGS) -o $ 13 14 lean: 15 $(RM) $(OBJ) 26 4.5. COMPILATION AND OTHER ISSUES 16 17 leanEx: 18 $(RM) $(OBJ) $(BIN) 19 20 uda_atualiza.o: uda_atualiza.u 21 $(NVCC) $(CUFLAGS) $< - -o $ 22 23 uda_bloks2mf.o: uda_bloks2mf.u 24 $(NVCC) $(CUFLAGS) $< - -o $ 25 26 SFS2Dv01_64.o: SFS2Dv01_64.for 27 $(FORT) $(FORTFLAGS) $< - -o $ Bearing in mind that all the strutures are reated as Vetors in Fortran and Fortran indexing are 1-based (C uses 0-Based) an sp eial aess is required (Eq (4.5), (4.5) and (4.5)). Furthermore, Fortran stores the elements following Column-Ma jor Order while C storing is Row-Ma jor Order based. These two asp ets imply that: • The aess to the partiular p osition i of array V[M] is made, in C, as V(i) = V[i−1] (4.3) • The aess to the partiular p osition i, j of array V[MxN] is made in C as V(i, j) = V[(j−1) ·M+i−1] (4.4) • The aess to the partiular p osition i, j, k of array V[MxNxO] is made in C as V(i, j, k) = V[(k−1) ·M·N+ (j−1)M+i−1] (4.5) Load the mesh Load Initial Conditions Load BCs Stablish Sim. Length Stablish Sim. Length Compute Results Dump Data Free Resources t<tsim? Calc. dW Calc. dt Wet/Dry Correction Sync dt Update W No Yes Sync dW t=t+dt ⊗ Figure 4.8: Flux diagram for the appliation. Green-highlighted is the p orted slie of the o de 27 28 5 Results The ases hosen to show the results are fo used on how similar are the GPU numerial results to the ones obtained from the original CPU version (preision) and how eient this implementation an b e (p erformane). To ahieve this, two examples have b eed seleted. First, an aademi ase of unsteady ow with soure terms with analytial solution and seond a real life inundation ow of hydrauli interest. Furthermore, the GPU p erfomane has b een ompared with that of a distributed-parallel version of the CPU o de at [14℄ using a dam-break ow simulation with a large numb er of ells. 5.1 Preision: A test-ase with analytial solution This ase has b een used to minimize the dierenes b etween the results provided by the CPU and the GPU versions. The ase simulates the evolution of a mass of water ontained in a fritionless parab oloid. Test Case 1 orresp onds to zero initial velo ity and a urved initial free surfae shap e (Figure 5.1). As times go es on, the p otential energy transforms into kineti energy. It is a go o d ase b eause it has analytial solution [27℄ and there exists a hallenging wet/dry b oundary all the time. Figure 5.1: Left: Bed level and initial water depth state for test ase 1. As shown in Figure 5.2 and Figure 5.3 there are not visible dierenes b etween b oth simulations. In order to quantify the preision of the GPU implementation with resp et to the CPU, 29 CHAPTER 5. RESULTS Figure 5.2: Test ase 1. Left: GPU Simulated results for h and Right: CPU Simulated results for h at t= 42.03s. Figure 5.3: Test ase 1. Left: GPU Simulated result for |v| and Right: CPU Simulated results for |v| at t= 42.03s. the L1 , L2 and L∞ norm of the error in water depth at dierent times has b een alulated. Test Case 1 shows aeptable dierenes. This agrees with the error in the alulation reahing ma- hine preission ( O(−14) ) in b oth versions of the o de. The most sensitive region is the wet/dry b oundary where b oth the water depth and velo ity are very small. Test Case 2 orresp onds to the same fritionless ontainer but with dierent initial data orresp onding to a at surfae with velo ity. Although the visual omparison is also favorable, the detailed evaluation of the L1 , L2 and L∞ norm of the error in water depth at dierent times shows unaeptable dierenes whih ome from the preision of the double oating p oint data typ e, reahing O(L∞) = −4 . Studying the pro edene of the dierenes we nd the problem at the rst time step (See Figure 5.1). Following the numerial sheme, we found that: h∗∗∗ j=hn j−α1 k+β ˜ λ1 k≥0 (5.1) Attending to the new state for the seond time-step, we found the values for ell 65399 as app ears in Table 5.2 30 5.2. PERFORMANCE: A LARGE-SCALE SIMULATION AT JÚCAR RIVER if go o d preditions an b e obtained using redued domains of the study area or if, otherwise, it is preferable to dene large domains at the ost of less denition for the top ographi data if extremely long omputational times are to b e avoided. This hydrograph is syntheti sine no atual disharge reords exist [3℄. As the numerial domain D2 is lo ated 4 km downstrean of the Tous Dam it is p ossible to ompute a new dis- harge urve by reording the rate of ow disharge at an appropriate setion in D1 . Due to the huge magnitude of the o o ding the dierene b etween the two disharge urves is merely a lag time of a few hours. Both are displayed in Figure 5.14. Considering this, and the fat that no reords of the o o d wave arrival time exist, the same original disharge urve was set as inlet b oundary ondition in domains D1 and D2 when p erforming numerial simulations. At the oulet b oundary, downstream of the domains, the ow was let to exit freely without imp osing any onditions, as no information was provided. The initial depth of water in the river reah prior to the rain events is unknown. Taking into aount that the base ow of Júar River is roughly 50 m3s−1 whih is totally negligible in omparison with the sale of Tous outow hydrograph, the valley was assumed initially dry. Following [3℄ a Manning o eient of 0.030 sm−1/3 was used for the whole river b ed reah. Other zones of inreased Manning o eient are inluded. As the ground in the town area was fully paved with onrete, the o o d did not ero de it. Regarding reorded hydrauli data of the o o ding of the town of Sumaárel, a range for the maximum water elevation marks was olleted at 21 lo ations within or very lose to Sumaárel village. In b oth alulations a total time of 39 h was simulated with a omputational time of 5.5 h in the D1 domain and 22.3 h in the D2 domain. These gauging p oints are shown in Figure 5.11. Some gauges (numb ers 5, 9, 15, 17, 18 and 21) show no o o ding (zero or near zero maximum water depth) and orresp ond to lo ations just barely reahed by the o o ding so that they represent a sort of shore line of the o o d within the town. Table 1 ontains a summary of prob e lo ations, estimated maximum water depths and omputed maximum water depths on the two omputational domains. The values of the water depth at gauges 1 and 2, plaed in the lower part of the village indiate that the numerial solutions provided by b oth grids are a go o d predition of the maximum water level reahed by the o o ding at b oth stations. Both gauges register almost the same water level surfae evolution, as exp eted due to their proximity. Go o d agreement b etween maximum water elevation marks and predited data is also found for gauge lo ations 3 and 4, of similar b ed level elevation, and lo ated within the village. The results in table 1 show also a go o d agreement for gauge 5 that remains dry aording to the eld observations, despite it b eing lose to the river b ed. The elevation at gauge 6, within Sumaárel, is overestimated in approximately 1 m . The water depth at gauge 7 agrees well with the maximum water elevation mark, whilst water depth in gauge 8 is overestimated in approximately 1 m . 37 CHAPTER 5. RESULTS Gauge x(m)y(m) Est. max. h(m) Comp. max. h(m)D1 Comp. max. h(m)D2 1 2410 3290 17.5-19 18.613149 18.684626 2 2400 3335 8.0-9.0 10.181195 9.806911 3 2355 3315 7.0-8.0 7.270638 7.386148 4 2345 3380 7 6.775814 6.895801 5 2335 3175 0.2 0.000 0.00 6 2335 3420 5.0-6.0 7.464109 7.615280 7 2330 3365 6 6.101556 6.143140 8 2315 3450 5 6.561674 6.679546 9 2310 3590 0 0.304004 0.119698 10 2303 3255 4 3.887516 3.979779 11 2285 3425 2 3.039008 3.194761 12 2285 3500 5.0-6.0 4.772985 4.909878 13 2280 3280 2.5-3.0 4.186196 4.330580 14 2266 3550 2 3.549098 3.122085 15 2265 3400 0 1.928118 2.134662 16 2259 3530 3.0-4.0 3.698947 3.802850 17 2250 3440 0 0.661666 0.901334 18 2230 3525 0 1.041024 1.215631 19 2205 3445 2.0-3.0 2.026697 2.257170 20 2195 3440 2 1.857008 2.096829 21 2190 3485 0 0.000 0.00 Table 5.5: Gauges p osition, estimated maximum water depth and simulated water depth The results for gauges 9 and 10 show go o d agreement with eld observations. Gauge lo ation 9 remained dry along the o o ding and the simulation provides a maximum water depth in the sale of the entimeters. The numerial results for gauge 11 indiate an overestimation of the eld water depth estimation of approximately 1 m , whilst very go o d agreement is found for gauge 12. The simulations at the gauge lo ations 13 and 14 overestimate eld observations in approximately 1 m . Gauge lo ation 15 remained dry along the o o ding whereas the numerial simulation did not. On the other hand the results for gauge lo ation 16 are in go o d aordane with the observed eld data. Gauges 17 and 18 remained dry but the simulation estimates a maximum depth of nearly 1 m . The results for gauge lo ations 19 and 20 and 21 are in aordane with eld observations. The evolution of the omputed o o ding an b e seen in plan view in Figure 5.10 for times t= 5, 10, 15, 20, 25, and 30 hours. The omputed ow advanes and passes around the buildings but always moving inside the limit given by that line. Although mesh D1 has larger ells than D2 the numerial preditions from b oth grids are in general in agreement with observed data. It is remarkable that for this extreme event, despite the dierent lo ations of the inlet disharge setions and the dierent size of the ells in D1 and D2 , the water depth results for D1 are only slightly inferior than the ones obtained with D2 . 38 5.3. COMPARING WITH A DISTRIBUTED MEMORY PARALLEL IMPLEMENTATION It is very useful, when an exhaustive study is required, to rene the mesh in the area of interest. In this ase, the main trouble is to stablish the input hydrograph. Although water depth has no many signiative dierenes Figure 5.12 and 5.13, velo ity has not the same b ehaviour (Gauges 5 and 21 have b een ommited b eause of b oth have alulated the dry state). In Figure 5.15 is p ossible to appreiate the dierenes where the simulation p erformed with the oarse mesh makes a higher estimation of the velo ity. As displayed by the results of the water level time evolution at the gauges, the mesh renement in the zone of interest improves the quality of the preditions. The GPU simulation of the omputation on the rened mesh was 22 hours and 20 minutes (more than 28 days of simulation using CPU) and that for the oarse mesh was 5 hours and 30 minutes. The oarse mesh was a go o d aproximation of how the o o d advanes but not always an b e used to study the details in a partiular area. 5.3 Comparing with a distributed memory parallel implementation 28-Core ⋄ 1-Core GPU Cells 106648 CFL 0.9 tn 400.0 h0 5 - 0 Comp. Load (s.) 363.2 9383.83 250.79 Sup 25.84 37.41 Table 5.6: Computational load for a Dam-Break simulation (400 s.) with the mono-ore version, the MPI paralellized version and the new CUDA version. ⋄ Eah ore omes from an Intel i7 CPU 860  2.80 GHz This ase simulates the evolution of two onneted b oxes where one of them ontains 5 m. of water level and the other one is dry. The initial onditions and geometry are shown at 5.16. The reason to inlude this additional test ase is that it was run previously with a CPU version of the metho d paralellized through distrubuted mahines paradigm using Standard MPI. The simulation was run during 400 s. dumping data eah 200 time-steps. Furthermore, it has b een used CFL=0.9 and a manning o eient of m= 0.03 . The results show that the p ower of omputing of the GPU is omparable with the p ower of more than 30 omputers working at the same time using the Distrubuted Computing paradigm. Although CUDA programming is not as easy as MPI programming and it is imp ortant to note that not every implementations supp ort b oth kind of implementations, the p erformane of the rst tehnique is muh b etter. 39 CHAPTER 5. RESULTS Figure 5.12: Simulated and estimated water depth in 1-11 Gauges. 40 5.3. COMPARING WITH A DISTRIBUTED MEMORY PARALLEL IMPLEMENTATION Figure 5.13: Simulated and estimated water depth in 12-21 Gauges. 41 CHAPTER 5. RESULTS 0 2000 4000 6000 8000 10000 12000 14000 16000 0 20000 40000 60000 80000 100000 120000 140000 160000 Discharge (m3/s) t (s) 0 2000 4000 6000 8000 10000 12000 14000 16000 0 20000 40000 60000 80000 100000 120000 140000 160000 Discharge (m3/s) t (s) Figure 5.14: Tous syntheti hydrograph for D1 (Right) and D2 (Left) Figure 5.15: Comparison of Left: Coarse mesh velo ity mo dule and Righ: Rened mesh veloity module at t= 13h Figure 5.16: Initial onditions of water depth and mesh plot 42 5.3. COMPARING WITH A DISTRIBUTED MEMORY PARALLEL IMPLEMENTATION Figure 5.17: 5-0 Dam-Break simulation for (Right-Left, Top-Down) t=5, 10, 15, 20, 25, 30 seonds 43 44 6 Conlusions and future work A rst order nite volume sheme to disretize the Shallow Water equations on unstrutured meshes has b een implemented using GPUs. The asso iated sp eed-up has b een studied when solving dierent problems with nVidia Tesla Series 2070. The diulties generated by the use of unstrutured meshes have b een identied and partially overome so that our results show that it is p ossible to solve many dierent problems 30 times faster than a ommon CPU version on a single pro essor. Furthermore, only mahine preision dierenes are enountered b etween b oth implementations, so it is imp ortant to note that the sp eed of the simulation do es not aet the preision of the numerial metho d. Communiating data b etween CPU and GPU has a very exp ensive ost. An interesting strategy to redue the impat of the ommuniation has b een prop osed. The only neessity of ommuniation is the elapsed simulation time so that the CPU shedules the op erations. Previous work related to reduing the omputational ost by means of parallel CPU programming has b een ompared, showing that a GPU ould b e faster than 30 CPU ores involving less investment and less energy onsumption. The values of 50-100x sp eed-up announed in the related literature have not b een reahed in our implementation. Our interpretation is that it is not p ossible to b e more than 42 times faster than a CPU pro essor when working with double preision data and serious and areful sp eed-up omparisons are required in any ase. Although it is very ompliated to reah the theorial p erformane p eak, b oth implementations ould reah a reasonable p ower, so if b oth implementations are mostly optimized, sp eed ups like the related in this work are aeptable. As further work, it is interesting to explore the Multi-GPU paradigms, simulating with many GPUs and to study other implementations whih p erform the memory aess pattern under unstrutured meshes. 45 46