A Tuned and Scalable Fast Multipole Method as a Preeminent Algorithm for Exascale Systems

Among the algorithms that are likely to play a major role in future exascale computing, the fast multipole method (FMM) appears as a rising star. Our previous recent work showed scaling of an FMM on GPU clusters, with problem sizes in the order of billions of unknowns. That work led to an extremely parallel FMM, scaling to thousands of GPUs or tens of thousands of CPUs. This paper reports on a a campaign of performance tuning and scalability studies using multi-core CPUs, on the Kraken supercomputer. All kernels in the FMM were parallelized using OpenMP, and a test using 10^7 particles randomly distributed in a cube showed 78% efficiency on 8 threads. Tuning of the particle-to-particle kernel using SIMD instructions resulted in 4x speed-up of the overall algorithm on single-core tests with 10^3 - 10^7 particles. Parallel scalability was studied in both strong and weak scaling. The strong scaling test used 10^8 particles and resulted in 93% parallel efficiency on 2048 processes for the non-SIMD code and 54% for the SIMD-optimized code (which was still 2x faster). The weak scaling test used 10^6 particles per process, and resulted in 72% efficiency on 32,768 processes, with the largest calculation taking about 40 seconds to evaluate more than 32 billion unknowns. This work builds up evidence for our view that FMM is poised to play a leading role in exascale computing, and we end the paper with a discussion of the features that make it a particularly favorable algorithm for the emerging heterogeneous and massively parallel architectural landscape

FMM-based vortex method for simulation of isotropic turbulence on GPUs, compared with a spectral method

The Lagrangian vortex method offers an alternative numerical approach for direct numerical simulation of turbulence. The fact that it uses the fast multipole method (FMM)--a hierarchical algorithm for N-body problems with highly scalable parallel implementations--as numerical engine makes it a potentially good candidate for exascale systems. However, there have been few validation studies of Lagrangian vortex simulations and the insufficient comparisons against standard DNS codes has left ample room for skepticism. This paper presents a comparison between a Lagrangian vortex method and a pseudo-spectral method for the simulation of decaying homogeneous isotropic turbulence. This flow field is chosen despite the fact that it is not the most favorable flow problem for particle methods (which shine in wake flows or where vorticity is compact), due to the fact that it is ideal for the quantitative validation of DNS codes. We use a 256^3 grid with Re_lambda=50 and 100 and look at the turbulence statistics, including high-order moments. The focus is on the effect of the various parameters in the vortex method, e.g., order of FMM series expansion, frequency of reinitialization, overlap ratio and time step. The vortex method uses an FMM code (exaFMM) that runs on GPU hardware using CUDA, while the spectral code (hit3d) runs on CPU only. Results indicate that, for this application (and with the current code implementations), the spectral method is an order of magnitude faster than the vortex method when using a single GPU for the FMM and six CPU cores for the FFT

Etude de l'adéquation des machines Exascale pour les algorithmes implémentant la méthode du Reverse Time Migation

As we are expecting Exascale systems for the 2018-2020 time frame, performance analysis and characterization of applications for new processor architectures and large scale systems are important tasks that permit to anticipate the required changes to efficiently exploit the future HPC systems. This thesis focuses on seismic imaging applications used for modeling complex physical phenomena, in particular the depth imaging application called Reverse Time Migration (RTM). My first contribution consists in characterizing and modeling the performance of the computational core of RTM which is based on finite-difference time-domain (FDTD) computations. I identify and explore the major tuning parameters influencing performance and the interaction between the architecture and the application. The second contribution is an analysis to identify the challenges for a hybrid and heterogeneous implementation of FDTD for manycore architectures. We target Intel’s first Xeon Phi co-processor, the Knights Corner. This architecture is an interesting proxy for our study since it contains some of the expected features of an Exascale system: concurrency and heterogeneity.My third contribution is an extension of the performance analysis and modeling to the full RTM. This adds communications and IOs to the computation part. RTM is a data intensive application and requires the storage of intermediate values of the computational field resulting in expensive IO accesses. My fourth contribution is the final measurement and model validation of my hybrid RTM implementation on a large system. This has been done on Stampede, a machine of the Texas Advanced Computing Center (TACC), which allows us to test the scalability up to 64 nodes each containing one 61-core Xeon Phi and two 8-core CPUs for a total close to 5000 heterogeneous coresLa caractérisation des applications en vue de les préparer pour les nouvelles architectures et les porter sur des systèmes très étendus est une étape importante pour pouvoir anticiper les modifications nécessaires. Comme les machines Exascale sont prévues pour la période 2018-2020, l'étude des applications et leur préparation pour ces machines s'avèrent donc essentielles. Nous nous intéressons aux applications d'imagerie sismique et en particulier à l'application Reverse Time Migration (RTM) car elle est très utilisée par les pétroliers dans le cadre de l'exploration sismique.La première partie de nos travaux a porté sur l'étude du cœur de calcul de l'application RTM qui consiste en un calcul de différences finies dans le domaine temporel (FDTD). Nous avons caractérisé cette partie de l'application en soulevant les aspects architecturaux des machines actuelles ayant un fort impact sur la performance, notamment les caches, les bandes passantes et le prefetching. Cette étude a abouti à l'élaboration d'un modèle de performance permettant de prédire le trafic DRAM des FDTD. La deuxième partie de la thèse se focalise sur l'impact de l'hétérogénéité et le parallélisme sur la FDTD et sur RTM. Nous avons choisi l'architecture manycore d’Intel, Xeon Phi, et nous avons étudié une implémentation "native" et une implémentation hétérogène et hybride, la version "symmetric". Enfin, nous avons porté l'application RTM sur un cluster hétérogène, Stampede du Texas Advanced Computing Center (TACC), où nous avons effectué des tests de scalabilité allant jusqu'à 64 nœuds contenant des coprocesseurs Xeon Phi et des processeurs Sandy Bridge ce qui correspond à presque 5000 cœur

Ein HPC-tauglicher Spektralelemente-Löser auf der Grundlage von statischer Kondensation und Mehrgittermethoden

Arbeitstitel: Erweiterte mathematische Methoden zur Simulation von turbulenten Strömungsvorgängen auf parallelen Rechnern:Inhaltsverzeichnis 1 Einleitung 3 2 Numerische Simulation physikalischer Prozesse 6 2.1 Königsdisziplin - Turbulente Strömungssimulation 6 2.2 Vom mathematischen Modell zur numerischen Lösung 9 2.2.1 Räumliche und zeitliche Diskretisierung 9 2.2.2 Allgemeine Reduktion auf Poisson- und Helmholtz-Gleichungen 11 2.3 Anforderungen an effiziente Lösungsverfahren 12 3 Basiskomponenten des entwickelten Verfahrens 16 3.1 Spektralelemente-Methode 16 3.1.1 Grundlagen 17 3.1.2 Gewählte Ansatzfunktionen und Stützstellen 20 3.1.3 Struktur des linearen Operators 24 3.2 Statische Kondensation 25 3.3 Geometrisches Mehrgitterverfahren 26 4 Das Spektralelemente basierte Mehrgitterverfahren auf kondensierten Gittern 31 4.1 Stand der Forschung 31 4.2 Mehrgitterverfahren auf kondensierten Gittern 32 4.2.1 Konzeption wirkungsvoller Glätter 34 4.3 Nachweis optimaler Eigenschaften 41 4.3.1 Lineare Komplexität 41 4.3.2 Ausgezeichnete Konvergenzgeschwindigkeit 43 4.3.3 Robustheit gegenüber Gitterverfeinerung 46 5 Konzeption des parallelen Mehrgitterlösers 49 5.1 Parallelrechner und Leistungsbewertungskriterien 49 5.2 Stand der Forschung 52 5.3 Grundlegende Struktur und Parallelisierung 54 5.3.1 Analyse des Speicherbedarfs 54 5.3.2 Zwei- und dreidimensionale Zerlegung 58 5.3.3 Parallelisierung und Kommunikation 62 6 Ergebnisse 65 6.1 Implementierung des Lösers 65 6.2 Hardwarespezifikation des Testsystems 66 6.3 Bewertung der Implementierung 68 6.3.1 Sequentieller Anwendungsfall 68 6.3.2 Nachweis der Skalierbarkeit im parallelen Anwendungsfall 76 6.3.3 Vergleich mit etablierten Lösungsansätzen bzw. Lösern 87 7 Zusammenfassung und Ausblick 89 Abbildungsverzeichnis 92 Tabellenverzeichnis 94 Abkürzungsverzeichnis 95 Symbolverzeichnis 96 Literaturverzeichnis 98 A Weiterführende Messergebnisse 106 A.1 Relative Mehrkosten der parallelen Implementierung 106 A.2 Sequentielle Lösungszeiten ohne Nachglättung im 2D-Fall 107 A.3 Sequentielle Lösungszeiten ohne Nachglättung im 3D-Fall 108Die rechnergestützte Simulation physikalischer Prozesse ist ein fester Bestandteil im Alltag von Wissenschaftlern aus den unterschiedlichsten Wissensbereichen. Unabhängig davon, ob das Ziel die Vorhersage des Wetters von morgen, die Konzentrationsbestimmung von Fluidteilchen in Mischprozessen oder die Erschaffung von Werkstoffen mit optimalen Materialeigenschaften ist, ohne den Einsatz von leistungsfähigen Rechnern ist dieses Ziel nicht zu erreichen. Aus dieser inhärenten Kopplung lässt sich eine grundlegende Aussage bzgl. der Laufzeit durchzuführender Simulationen ableiten. Schnellere Rechentechnik reduziert automatisch die Laufzeit einer bereits bestehenden Simulation und somit auch die Wartezeit auf die potentiell zu erwartenden Erkenntnisse. Zeitgleich ist die so erreichte Reduktion der Berechnungszeit auch ein Maß für die mögliche Erhöhung des Detailgrades einer bestehenden Simulation und somit auch ein Indikator für den zusätzlich zu erwartenden Erkenntnisgewinn. Ein Blick auf die seit 1993 herausgegebene Top500-Liste der schnellsten Supercomputer zeigt ein annähernd gleichbleibend exponentielles Wachstum der Rechenleistung. Dieser durch eine Interpretation von „Moores-Law“ beschriebene Sachverhalt wird laut aktuellen Prognosen auch in den nächsten Jahren bestehen bleiben. Für die im Bereich der Simulation tätigen Wissenschaftler gleicht dies einem Versprechen, dass ohne deren Zutun auch in Zukunft mit stetig kürzeren Simulationszeiten zu rechnen ist. Immer vorausgesetzt, es können genug finanzielle Mittel für die neue Hardware akquiriert werden. Doch dieser Schein trügt. Eine genauere Analyse der Entwicklung der Rechentechnik der letzten Jahre zeigt zwei maßgebliche Veränderungen. Zum einen stagniert die maximale Taktrate einer einzelnen CPU seit Erreichen der 4 GHz Grenze im Jahr 2004 und zum anderen wird, insbesondere seit der Einführung der ersten Dual Core CPU’s 2005, gesteigerte Rechenleistung fast gänzlich durch die Verwendung einer Vielzahl von Rechenkernen erreicht. Das aktuell mit mehr als zehn Millionen Rechenkernen an Position 1 der Top500-Liste geführte System TaihuLight (deu. Licht der Göttlichkeit) verdeutlicht die Dimensionen dieser Entwicklung. Die für dieses System in Aussicht gestellte maximale Rechenleistung von circa 125 Billiarden gleitkommaoperationen pro Sekunde, kann dabei nur von einer perfekt parallelisierten Simulationsrechnung erreicht werden. „Amdahls-Law“ zeigt jedoch, dass dieser perfekte Zustand, aufgrund von unvermeidlichen sequentiellen Abschnitten in den einzelnen im Programm verwendeten Algorithmen, nicht zu erreichen ist. Die genaue Abweichung vom vollparallelisierten Idealzustand wird dabei durch die sogenannte parallele Effizienz quantifiziert. Deren Wert muss hierbei per Definition zwischen Null und Eins liegen. Dem Paradigma „eine hohe parallele Effizienz ergibt eine hohe Rechenleistung und dies führt zur kürzestmöglichen Simulationslaufzeit“ folgend, wurden in den letzten Jahren die unterschiedlichsten Simulationsprogramme auf eben diese Effizienz getrimmt. In den meisten Fällen wurden hierfür Codes verwendet, die auf eine sehr lange Historie zurückgreifen, so dass alte bestehende Strukturen und Algorithmen unabhängig von deren wirklicher Eignung parallelisiert wurden. Diese Entwicklung führt jedoch mehr und mehr dazu, dass die Entwickler den Blick für die Vielseitigkeit der Faktoren, die zu einer akzeptablen Simulationslaufzeit führen, verlieren. Werden zum Beispiel Methoden niederer Ordnung, wie dies etwa bei den Standard Finite-Differenzen-Verfahren der Fall ist, zur Diskretisierung des Simulationsgebietes eingesetzt, steigt die Zahl der für kleine Lösungsfehler benötigten Gitterpunkte so schnell an, dass jedweder Arbeitsspeicher vor Erreichen der benötigten Genauigkeit aufgebraucht ist. Im Gegensatz dazu sind Methoden höherer Ordnung, wie dies etwa bei den Standard Finite-Elemente-Verfahren der Fall ist, aufgrund ihrer suboptimalen numerischen Komplexität kaum besser geeignet. Ein ähnliches Bild ergibt sich bei den Algorithmen, mit denen die Gleichungssysteme in den einzelnen Simulationsschritten gelöst werden. Stellvertretend sei hier das Jacobi-Verfahren genannt, welches sich zwar durch eine parallele Effizienz nahe Eins auszeichnet, jedoch zum einen eine nicht optimale quadratische numerische Komplexität und zum anderen eine von der Auflösung des Simulationsgitters abhängige maximale Iterationszahl besitzt. Sofern die Anwender der etablierten Simulationsprogramme keine Kosten für den Zugang zu Hochleistungsrechnern zu erwarten haben und diese Rechner immer wieder massiv ausgebaut werden, stellen die genannten Einschränkungen fürs Erste nur bedingt ein Problem dar. Denn, eine Simulation die nach Hinzunahme einer bestimmten Zahl von Rechenkernen um annähernd diesen Faktor beschleunigt wird ist etwas Ausgezeichnetes. Werden den Anwendern jedoch, wie bereits von immer mehr Universitätsrechenzentren diskutiert und in der Industrie bereits gängige Praxis, die Kosten für den Energieverbrauch in Rechnung gestellt, ergibt sich ein gänzlich anderes Bild. Ein Bild, in dem der Effizienz, die die angewandten Methoden bzw. die eingesetzten Algorithmen erreichen, die größte Bedeutung zufällt. Die Effizienz einer Methode wird hierbei ungenauerweise oft nur anhand deren Implementierung als Algorithmus bestimmt. Jedoch kann eine effizient implementierte Methode mit numerisch ungünstigen Eigenschaften einer nicht effizient implementierten Methode mit numerisch optimalen Eigenschaften deutlich unterlegen sein. Demnach ist es offensichtlich, dass nur für eine effizient implementierte Methode mit optimalen numerischen Eigenschaften die kürzestmögliche Simulationslaufzeit erreicht werden kann. Der Fokus der vorliegenden Arbeit liegt deshalb zu allererst auf dem Nachweis der optimalen numerisch/mathematischen Eigenschaften der entwickelten Methode. Diese Eigenschaften sind: lineare numerische Komplexität, Robustheit des Verfahrens gegenüber Gitterverfeinerungen im Simulationsgebiet und eine bisher unerreichte Konvergenzrate. Abschließend wird zusätzlich die Eignung der Methoden bzgl. deren Verwendung auf aktuellen Hochleistungsrechnern unter Verwendung von Zehntausenden von Rechenkernen belegt und auch deren effiziente Implementierung bzw. Umsetzung dargelegt. Ziel dieser Arbeit ist die Entwicklung effizienter mathematischer Methoden zur numerischen Simulation von physikalischen Prozessen und deren hochskalierende Implementierung auf Hochleistungsrechnern. Unter allen denkbaren Aufgabenstellungen zählen hierbei insbesondere diejenigen zu den herausforderndsten, die der Strömungsmechanik zugeordnet sind. Besonders die direkte numerische Simulation (DNS), welche zur Analyse von turbulenten Strömungsphänomenen eingesetzt wird, stellt hierbei höchste Ansprüche an die eingesetzten numerischen Verfahren. Die Entwicklung und Umsetzung der im Rahmen dieser Arbeit vorgestellten Methoden ist deshalb auf die Anwendung im Rahmen der turbulenten Strömungssimulation ausgerichtet. Diese Fokussierung dient jedoch allein dem Beleg der Leistungsfähigkeit und stellt keine prinzipielle Einschränkung der Methode dar.:Inhaltsverzeichnis 1 Einleitung 3 2 Numerische Simulation physikalischer Prozesse 6 2.1 Königsdisziplin - Turbulente Strömungssimulation 6 2.2 Vom mathematischen Modell zur numerischen Lösung 9 2.2.1 Räumliche und zeitliche Diskretisierung 9 2.2.2 Allgemeine Reduktion auf Poisson- und Helmholtz-Gleichungen 11 2.3 Anforderungen an effiziente Lösungsverfahren 12 3 Basiskomponenten des entwickelten Verfahrens 16 3.1 Spektralelemente-Methode 16 3.1.1 Grundlagen 17 3.1.2 Gewählte Ansatzfunktionen und Stützstellen 20 3.1.3 Struktur des linearen Operators 24 3.2 Statische Kondensation 25 3.3 Geometrisches Mehrgitterverfahren 26 4 Das Spektralelemente basierte Mehrgitterverfahren auf kondensierten Gittern 31 4.1 Stand der Forschung 31 4.2 Mehrgitterverfahren auf kondensierten Gittern 32 4.2.1 Konzeption wirkungsvoller Glätter 34 4.3 Nachweis optimaler Eigenschaften 41 4.3.1 Lineare Komplexität 41 4.3.2 Ausgezeichnete Konvergenzgeschwindigkeit 43 4.3.3 Robustheit gegenüber Gitterverfeinerung 46 5 Konzeption des parallelen Mehrgitterlösers 49 5.1 Parallelrechner und Leistungsbewertungskriterien 49 5.2 Stand der Forschung 52 5.3 Grundlegende Struktur und Parallelisierung 54 5.3.1 Analyse des Speicherbedarfs 54 5.3.2 Zwei- und dreidimensionale Zerlegung 58 5.3.3 Parallelisierung und Kommunikation 62 6 Ergebnisse 65 6.1 Implementierung des Lösers 65 6.2 Hardwarespezifikation des Testsystems 66 6.3 Bewertung der Implementierung 68 6.3.1 Sequentieller Anwendungsfall 68 6.3.2 Nachweis der Skalierbarkeit im parallelen Anwendungsfall 76 6.3.3 Vergleich mit etablierten Lösungsansätzen bzw. Lösern 87 7 Zusammenfassung und Ausblick 89 Abbildungsverzeichnis 92 Tabellenverzeichnis 94 Abkürzungsverzeichnis 95 Symbolverzeichnis 96 Literaturverzeichnis 98 A Weiterführende Messergebnisse 106 A.1 Relative Mehrkosten der parallelen Implementierung 106 A.2 Sequentielle Lösungszeiten ohne Nachglättung im 2D-Fall 107 A.3 Sequentielle Lösungszeiten ohne Nachglättung im 3D-Fall 10

High-order hybrid methods using Green's functions and finite elements

Green's functions provide an elegant mathematical framework for linear partial differential boundary value problems by directly defining solutions in integral form. For a Green's function that satisfies the desired boundary conditions, this form is a convolution of the Green's function with the right hand side. However, straightforward numerical computation of the convolution generally results in a prohibitively costly algorithm. This dissertation focuses on combining Green's function-based computations with finite element (FE) methods to alleviate their computational challenges while enhancing the representation capability of standard FE solutions on structured or semi-structured meshes. Two main areas of application are considered: a new finite-element-based method for N-body calculations and a new high-order method for elliptic interface problems on non-conforming volume meshes. First, the FE-based particle-particle--particle-mesh (FE-P3M) method for N-body problems is introduced. P3M methods are designed to compute quantities (e.g. potentials or forces) that are sums of discrete Green's functions and represent the effects of N interacting bodies. To avoid an N^2 operation when the quantity is desired at all body locations (as is often the case), P3M methods decompose the potential into short-range interactions between near neighbors and long-range interactions that are smooth and readily solved by mesh-based numerical methods. Instead of taking the traditional "Ewald-based" approach of using Gaussian screen functions to accomplish this decomposition, the FE-P3M method builds specially designed polynomial screens on an introduced finite element mesh. Due to this form of the screens, the long-range component of the potential is accurately resolved with a finite element method. When compared to classic P3M methods, which rely on the fast Fourier transform (FFT), the FE-P3M method allows for more flexible boundary conditions and possibly much lower communication costs in a parallel implementation. The method is described in detail, its cost and memory requirements are discussed, and its accuracy is demonstrated. The remaining chapters consider Green's functions in the context of layer potential solutions to homogeneous boundary value problems and the resulting integral equations (IEs) formulated on the boundary of the domain. A high-order coupling of finite elements and integral equations is presented for the solution of elliptic interface problems on embedded domains. This FE-IE method requires no special basis functions nor modifications to handle homogeneous or non-homogeneous jump conditions, and retains high-order convergence close to the embedded interface. The interface or embedded boundary conditions are enforced directly at the discretization points of the embedded surface mesh. High-order numerical convergence is demonstrated for both interior and exterior embedded domain problems, with a novel splitting defined and analyzed for the latter. Theoretical error convergence, existence and uniqueness of the decomposition, and numerical solution considerations are discussed. The new FE-IE interface method is then shown to be built from interior and exterior FE-IE subproblems. Finally, an extension to Stokes flow around embedded objects in two and three dimensions is presented