DEVELOPMENT OF MASSIVELY PARALLEL NEAR PEAK PERFORMANCE SOLVERS FOR THREE-DIMENSIONAL GEODYNAMIC MODELLING

Academia and industry have not yet been able to adapt software and numerical algorithms to the rapid, fundamental changes in hardware évolution since the beginning of the 21th century. Many of today's applications and algorithms are as a resuit not well suited for the available hardware, what implies performances far below hardware's peak. We address in this thesis the current need to design new parallel algorithms and tools that ease the development of applications that are suited for today's and tomorrow's hardware. We present (1) the MATLAB HPC compiler HPC.m, which greatly simplifies the building of parallel high performance applications and (2) parallel algorithms for the 3D simulation of strongly nonlinear processes as mechanical and reactive porosity waves. To simulate mechanical porosity waves we employ a simple numerical algorithm that permits to resolve the déformation of fluid-filled viscoelastic porous média in 3D. The utilized mathematical model is based on Biot's poroelastic theory, extended to account for viscous déformation and plastic yielding. The algorithm is designed for massively parallel high performance Computing. It employs finite différence stencil calculations on a staggered grid to approximate spatial derivatives. Pseudo-translent itérations were utilized to formulate an explicit algorithm and Picard itérations to resolve the nonlinearities. The modelling results exhibit the impact of decompaction weakening on the formation of three-dimensional solitary-wave-like moving porosity channels. We evaluate the algorithm's suitability for the building of high performance massively parallel 3D solvers and compare the achievable performance and scalability to the ones of other state-of-the-art algorithms that were developed aiming at large-scale simulations. Our algorithm is found to be better suited for the building of high performance massively parallel 3D solvers than the other algorithms considered in this thesis. To simulate reactive porosity waves we use a solver for 3D déformation of fluid-filled reactive viscous porous média. The employed algorithm uses essentially the same numerical methods as the algorithm for the simulation of mechanical porosity waves. The Damkôhler number (Da) of the simulations is varied in order to estimate the respective rôles of viscous déformation (low Da) and reaction (high Da) on wave propagation. 3D waves are found to propagate independently of their source at constant speed by going through each other for ail the investigated Da. Soliton-like wave propagation as a resuit of metamorphic reaction provides an efficient mechanism for fluid flow in the Earth's crust. It is expected that this mechanism takes place at the meter-scale in the lower crust and at the kilometer-scale in the upper crust providing explanations for both metamorphic veins formation and fluid extraction in subduction zones. The here developed tool HPC.m transforms simple MATLAB scripts into parallel high performance applications for GPU-, CPU- and MIC-supercomputers, clusters or workstations. It is designed forstencil-based applications, in particular for iterative PDE solvers that use finite différences and a regular Cartesian grld. Its core is a source- to-source translator. We illustrate in this thesis the great performance and versatility of HPC.m by deploying it to generate solvers for a variety of physics across multiple Earth Science disciplines, Ail solvers run close to hardware's peak performance and scale linearly on the 80 GPUs of the Octopus cluster, hosted by the Institute of Earth Sciences at University of Lausanne (Lausanne, Switzerland). They achieve moreover a speedup over the fully vectorised MATLAB input script of about 250x to 500x on one GPU, of lOOOx to 2000x on one workstation with 4 GPUs and of 17 OOOx to 35 OOOx on 80 GPUs. We show in addition that our nonlinear poroviscoelastic two-phase flow solver scales also linearly on the 5000 GPUs of the Piz Daint supercomputer at the Swiss National Supercomputing Centre (CSCS, Lugano, Switzerland), achieving a speedup over the fully vectorised MATLAB input script of over 500 OOOx. We expect a similar scaling for ail the ten solvers. The source-to-source translator contained in HPC.m is the first known to the authors that can perform automatically ail tasks that are required for the génération of a massively parallel near peak performance supercomputing application from a code developed in classical prototyping environment as MATLAB. -- Académie et industrie n'ont pas encore été capables d'adapter les logiciels et les algorithmes numériques aux changements rapides et fondamentaux dans l'évolution du hardware depuis le début du 21ème siècle. Une grande partie des applications et algorithmes ne sont, par conséquent, pas bien adapté au hardware disponible ; cela implique des performances bien au-dessous du maximum atteignable par le hardware. Nous répondons dans cette thèse au besoin actuel de développer des nouveaux algorithmes parallèles et des outils qui simplifient le développement d'applications qui soient adaptées au hardware d'aujourd'hui et de demain. Nous présentons (1) le compilateur HPC MATLAB, HPCm, qui simplifie de manière importante la construction des applications parallèles de haute performance et (2) des algorithmes pour la simulation 3D de processus non-linéaires tel que des ondes de porosité mécaniques et réactives. Pour simuler des ondes de porosité mécaniques, nous utilisons un algorithme numérique simple qui permet de résoudre la déformation 3D de matière poreuse et viscoélastique, saturé de fluide. Le modèle mathématique utilisé est basé sur la théorie poroélastique de Biot, étendue pour tenir compte de déformation visqueuse et plastique. L'algorithme est conçu pour le calcul massivement parallèle de haute performance. Il utilise le calcul de différences finies pour approximer les dérivées spatiales. Des itérations pseudo-transitoires ont été utilisées pour formuler un algorithme explicit et des itérations de Picard pour résoudre des non-linéarités. Les résultats de modélisation montrent l'impact de l'affaiblissement en décompaction sur la formation de canaux de porosité, qui se déplacent de manière similaire à des solitons. Nous évaluons l'aptitude de l'algorithme pour la construction des solveurs massivement parallèles de haute performance et 3D ; de plus, nous comparons la performance et la scalabilité atteignable à ceux d'autres algorithmes qui ont été développés dans le but d'effectuer des simulations à grande échelle. Notre algorithme a été trouvé mieux adapté pour la construction des solveurs massivement parallèles de haute performance et 3D que les autres algorithmes considérés dans cette thèse. L'outil HPCm, développé dans cette thèse, transforme des scripts MATLAB en applications parallèles de haute performance pour des super-ordinateurs, des clusters et des ordinateurs desktops avec GPU, CPU ou MICs. Il est conçu pour des applications basées sur le calcul de stencil ; en particulier pour des solveurs PDE qui utilisent des différences finies et une grille Cartésien régulière. Le cœur de HPCm est un traducteur source-source. Dans cette thèse, nous illustrons la grande performance et versatilité de HPC m en en faisant usage pour la génération de solveurs pour une variété de physiques à travers de multiples disciplines dans les Sciences de la Terre. Tous les solveurs atteignent une performance qui est proche du maximum que le hardware peut fournir et montrent une scalabilité linéaire sur les 80 GPUs du cluster de calcul Octopus à l'Institut des Sciences de la Terre à l'Université de Lausanne (Lausanne, Suisse). De plus, ils atteignent un speedup sur le script MATLAB d'entrée entièrement vectorisé d'environ 250x à 500x sur un GPU, lOOOx à 2000x sur un ordinateur desktop avec 4 GPUs et 17 000x à 35 OOOx sur 80 GPUs. Par ailleurs, nous montrons que notre solveur non-linéaire de flux deux-phases poroviscoélastique atteigne une scalabilité linéaire sur 5000 GPUs du super-ordinateur Piz Daint au Centre Suisse du Calcul Scientifique (CSCS, Lugano, Switzerland) et arrive ainsi à un speedup sur le script MATLAB d'entrée entièrement vectorisé d'environ 500 OOOx. Nous nous attendons à une scalabilité similaire pour tous les dix solveurs. Le traducteur source- source qui est contenu dans HPCm est le premier connu aux auteurs qui peut effectuer tous les taches nécessaires pour la génération d'une application massivement parallèle de haute performance à partir d'un code développé dans un environnement de prototypage classique tel que MATLAB