Development of linear solvers for large-scale CFD simulations on hybrid supercomputers

Abstract

(English) Divergence constraints are present in the governing equations of many physical phenomena, and they usually lead to a Poisson equation whose solution is one of the most challenging parts of scientific simulation codes. Indeed, it is the main bottleneck of incompressible computational fluid dynamics (CFD) simulations, and developing effective Poisson solvers is a critical task.

On the other hand, the fact that hardware's memory bandwidth tends to grow much slower than its peak performance has led to strongly memory-bound codes. Additionally, extreme-scale problems have comparable memory requirements, which can become problematic given the faster but smaller memory generally available in massively parallel accelerators. Using more computing resources is a natural way to avoid such limitations and tackle larger problems. However, this generally comes at the price of losing parallel efficiency and is exacerbated by the widening gap between memory and network bandwidths in fat nodes.

In this context, this thesis presents a strategy to mitigate several of the challenges above. It consists of exploiting s spatial symmetries to increase the arithmetic intensity of the simulations and reduce their memory requirements. Remarkably enough the proposed strategy is identical regardless of the mesh connectivity (structured or unstructured) and the geometric complexity of the problem, therefore applying to a wide range of academic and industrial configurations.

First of all, by proving the existence of an inexpensive block diagonalisation that transforms the original Poisson equation into a set of 2ŝ^s fully decoupled subsystems, we can accelerate significantly the iterative solvers' convergence. Additionally, imposing an adequate grid points' ordering allows reducing the operators' footprint and replacing the standard sparse matrix-vector product (SpMV) with the more compute-intensive sparse matrix-matrix product (SpMM). Numerical experiments making Algebraic Multigrid (AMG) and preconditioned Krylov subspace methods exploit a varying number of symmetries showed reductions of up to 23% and 72% in the iterations ' count, resulting in up to 1.7x and 5.6x overall speedups, respectively.

This approach is successfully extended to the preconditioning stage, giving rise to LRCFSAI(k), an enhanced variant of the Factored Sparse Approximate Inverse (FSAI) preconditioner. LRCFSAI(k) arises from introducing another level of approximation by taking advantage of the decoupled subsystems' close similarity to apply the same FSAI to all of them, therefore replacing SpMV with SpMM and leading to notable memory savings and increases in the arithmetic intensity. Of course, recycling the same preconditioner on all the subsystems worsens its convergence. Although being this effect much smaller than expected, it motivated the introduction of relatively cheap but very effective low-rank corrections. A key feature of these corrections is that thanks to being applied to each subsystem independently, the more symmetries being exploited, the more effective they become, leading to up to 5.7x faster convergences than the standard FSAI. Numerical experiments on up to 1.07 billion grids confirm the quality of LRCFSAI( k), which, despite being 2.6x lighter, outperforms the standard FSAI by a factor of up to 4.4x.

Finally, the advantages of exploiting symmetries and replacing SpMV with SpMM are extended throughout the simulations. As a result, matrix multiplications are accelerated, whereas their memory footprint is reduced, making massive simulations more affordable. As an example of practical application, we consider the numerical simulation of turbulent incompressible flows using a low-dissipation discretisation of unstructured collocated grids on both CPU and GPU. All the required matrices are classified into three sparsity patterns, showing up to 5.0x speed-ups and 8.0x memory savings.

(Català) Les equacions de conservació estan presents en molts fenòmens físics, i sovint condueixen a una equació de Poisson la solució de la qual és una de les parts més costoses dels codis de simulació numèrica. De fet, és el principal coll d'ampolla en les simulacions de dinàmica de fluids computacionals (CFD) incompressibles, i desenvolupar solucionadors de Poisson és una tasca crítica.

D'altra banda, el fet que l'amplada de banda de la memòria tendeixi a créixer molt més a poc a poc que el rendiment màxim ha portat a codis fortament memory-bound. A més, els problemes a exaescala tenen requisits de memòria comparables, la qual cosa pot ser problemàtica vist que la memòria en els acceleradors massivament paral·lels és més ràpida, però generalment més limitada. Utilitzar més nodes de càlcul és una manera natural d'evitar aquestes limitacions i abordar problemes més grans. Tanmateix, generalment comporta una pèrdua d'eficiència exacerbada per l'augment de la bretxa entre amplada de banda de memòria i xarxa als fat nodes.

En aquest context, aquesta tesi presenta una estratègia per mitigar diversos dels reptes anteriorment. Consisteix a explotar s simetries espacials per augmentar la intensitat aritmètica de les simulacions i reduir els seus requisits de memòria. L'estratègia proposada és idèntica independentment de la connectivitat de la malla (estructurada o no estructurada) i de la complexitat geomètrica del problema, i per tant és aplicable a una àmplia gamma de configuracions acadèmiques i industrials.

En primer lloc, demostrant l'existència d'una diagonalització per blocs que transforma l'equació de Poisson en un conjunt de 2^s subsistemes totalment desacoblats, podem accelerar significativament la convergència dels mètodes iteratius. A més, imposant un ordenament adequat de la malla podem reduir la memòria ocupada pels operadors i substituir el producte matriu-vector (SpMV) pel producte de matriu-matriu (SpMM), d'intensitat aritmètica major. Experiments numèrics basats en Algebraic Multigrid (AMG) i mètodes de Krylov explotant un nombre variable de simetries mostren reduccions de fins al 23% i 72% en el nombre d'iteracions, amb acceleracions generals de fins a 1,7x i 5,6x, respectivament.

Aquesta estratègia és estesa als precondicionadors donant lloc a LRCFSAI(k), una variant millorada del precondicionador Factored Sparse Approximate Inverse (FSAI). LRCFSAI(k) sorgeix d'introduir un nou nivell d'aproximació aprofitant la similitud dels 2^s subsistemes per aplicar el mateix FSAI a tots ells, substituint així SpMV per SpMM i garantint estalvis notables de memòria i increments en la intensitat aritmètica. És clar que reciclar el mateix precondicionador a tots els subsistemes empitjora la convergència. Tot i ser aquest efecte molt menor de l'esperat, va motivar la introducció de correccions de baix rang molt efectives. Una característica clau d'aquestes correccions és que, gràcies a ser aplicades a cada subsistema per separat, com més simetries s'exploten, més efectives esdevenen, donant lloc a convergències fins a 5,7x més ràpides que el FSAI estàndard. Experiments numèrics en malles de fins a 1.070 milions confirmen la qualitat de LRCFSAI(k), el qual, tot i ser 2,6x més lleuger, accelera el FSAI estàndard fins a 4,4x vegades.

Finalment, els avantatges d'aprofitar simetries i substituir SpMV per SpMM s'estenen a tota la simulació. Com a resultat, les multiplicacions per matrius s'acceleren, mentre que el consum de memòria es redueix, fent més assequibles les simulacions massives. Com a exemple d'aplicació pràctica, considerem la simulació numèrica de fluxos incompressibles turbulents mitjançant una discretització poc dissipativa de malles col·locades no estructurades tant en CPU com GPU. Finalment, classificant tots els operadors discrets en tres tipus podem generalitzar els guanys obtinguts, amb acceleracions de fins a 5,0x i estalvis de memòria de fins a 8,0x.

(Español) Las ecuaciones de conservación están presentes en muchos fenómenos físicos, y a menudo conducen a una ecuación de Poisson la solución de la cual es una de las partes más costosas de los códigos de simulación numérica. De hecho, es el principal cuello de botella en las simulaciones de dinámica de fluidos computacional (CFD) incompresibles, y desarrollar solvers de Poisson es una tarea crítica. Por otro lado, el hecho que el ancho de banda de la memoria tienda a crecer mucho más despacio que el rendimiento máximo ha llevado a códigos fuertemente memory-bound. Además, los problemas a exaescala tienen requisitos de memoria comparables, lo cual puede ser problemático ya que la memoria en los aceleradores masivamente paralelos es más rápida, pero generalmente más limitada. Utilizar más nodos de cálculo es una forma natural de evitar tales limitaciones para abordar problemas más grandes, pero generalmente comporta una pérdida de eficiencia. En este contexto, esta tesis presenta una estrategia para mitigar varios de los retos anteriores. Consiste en explotar s simetrías espaciales para aumentar la intensidad aritmética de las simulaciones y reducir sus requisitos de memoria. La estrategia propuesta es idéntica independientemente de la conectividad de la malla (estructurada o no estructurada) y de la complejidad geométrica del problema, siendo así aplicable a una amplia gama de configuraciones académicas e industriales. Demostrando la existencia de una diagonalización por bloques que transforma la ecuación de Poisson en 2^s subsistemas totalmente desacoplados podemos acelerar significativamente la convergencia de los métodos iterativos. Además, imponiendo un ordenamiento adecuado de la malla podemos reducir la memoria ocupada por los operadores y sustituir el producto matriz-vector (SpMV) por el producto matriz-matriz (SpMM), de intensidad aritmética mayor. Experimentos numéricos explotando un número variable de simetrías con Algebraic Multigrid (AMG) y con métodos de Krylov muestran reducciones de hasta el 23% y 72% en el número de iteraciones, con aceleraciones de hasta 1,7x y 5,6x, respectivamente. Esta estrategia es extendida a los precondicionadores dando lugar a LRCFSAI(k), una variante mejorada del preacondicionador Factored Sparse Approximate Inverse (FSAI). LRCFSAI(k) surge de introducir otro nivel de aproximación aprovechando la similitud de los 2^s subsistemas para aplicar el mismo FSAI en todos ellos, sustituyendo así SpMV por SpMM y garantizando ahorros de memoria e incrementos de intensidad aritmética. Está claro que reciclar el mismo preacondicionador empeora la convergencia. Aun ser este efecto muy menor de lo esperado, motivó la introducción de correcciones de bajo rango muy efectivas. Una característica de estas correcciones es que, gracias a ser aplicadas a cada subsistema por separado, se tornan más efectivas al explotar más simetrías, dando lugar a convergencias hasta 5,7x más rápidas que con el FSAI estándar. Experimentos numéricos en mallas de hasta 1.070 millones confirman la calidad de LRCFSAI(k), el cual, a pesar de ser 2,6x más ligero, acelera el FSAI estándar hasta 4,4x. Finalmente, las ventajas de aprovechar simetrías y sustituir SpMV por SpMM se extienden a toda la simulación. Como resultado, las multiplicaciones por matrices se aceleran, mientras que el consumo de memoria se reduce, haciendo más asequibles las simulaciones masivas. Como ejemplo de aplicación práctica, consideramos la simulación numérica de flujos incompresibles turbulentos tanto en CPU como GPU. Así clasificando los operadores discretos en tres tipos podemos generalizar las ganancias obtenidas, con aceleraciones de hasta 5,0x y ahorros de memoria de hasta 8,0x.