A GEMM-based direct solver for finite-difference Poisson problems in non-uniform grids

Deze paper introduceert een robuuste directe Poisson-oplosser voor niet-uniforme roosters die gebruikmaakt van een tensorformulering en GEMM-bewerkingen om de rekentijd te minimaliseren en hoge parallelle efficiëntie te bereiken op moderne CPU- en GPU-systemen.

De kern

Stel je voor dat je een enorme, driedimensionale puzzel moet oplossen. In de wereld van computerwetenschappen heet dit het oplossen van de Poisson-vergelijking. Dit is een wiskundige klus die nodig is om te begrijpen hoe druk zich verspreidt in een vloeistof, zoals lucht of water. Als je wilt simuleren hoe een storm ontstaat, hoe een vliegtuig vliegt of hoe bloed door aderen stroomt, moet je deze vergelijking duizenden keren per seconde oplossen.

Het probleem? De meeste computersimulaties werken het beste met een gelijkmatig raster (een rooster van gelijke blokjes), net als een perfect gestreken laken. Maar in de echte wereld is dat niet altijd handig.

Het probleem: De "Lekke" Netten

Stel je voor dat je een zwembad simuleert.

• Bij de randen (waar de muren zijn) moet je heel precies zijn, omdat daar de stroming verandert. Je hebt daar dus heel kleine, nauwkeurige blokjes nodig.
• In het midden van het zwembad stroomt het water rustig. Daar zijn grote, grove blokjes genoeg.

Als je een gelijkmatig raster gebruikt, moet je overal de kleine blokjes gebruiken, omdat je de randen niet kunt verwaarlozen. Dat betekent dat je een enorm aantal blokjes nodig hebt, zelfs waar ze niet nodig zijn. Dit maakt de berekening traag en duur.

De oude methoden om dit op te lossen (zoals "Multigrid") waren als een trage, ingewikkelde machine die op de nieuwe, snelle grafische kaarten (GPUs) van vandaag niet goed werkt. Ze werden vaak "verstikt" door de complexiteit.

De Oplossing: Een Slimme "GEMM"-Machine

De auteurs van dit paper (Pedro Costa en zijn team) hebben een nieuwe, directe methode bedacht. Ze noemen het een GEMM-gebaseerde solver.

Laten we dit uitleggen met een analogie:

1. De Oude Manier (FFT): De Perfecte Dans
Vroeger gebruikten wetenschappers een methode die leek op een perfecte dans (Fourier-transformaties). Als je rooster gelijkmatig was, konden ze de puzzel oplossen door de data in een specifieke volgorde te draaien en te spiegelen. Dit was razendsnel, maar het werkte alleen als het rooster perfect gelijkmatig was. Als je het rooster vervormde (omdat je kleine blokjes bij de wanden wilde), viel de dans uit elkaar en werkte het niet meer.

2. De Nieuwe Manier (GEMM): De Krachtige Vrachtwagen
Deze nieuwe methode is als het vervangen van die delicate dans door een krachtige, zware vrachtwagen (GEMM staat voor General Matrix-Matrix Multiplication).

• In plaats van te proberen een perfecte dans te doen met een vervormd rooster, pakken ze de data en gooien ze deze in een enorme, gestructureerde stapel.
• Vervolgens gebruiken ze de kracht van moderne computers (vooral de krachtige grafische kaarten van NVIDIA) om deze stapel in één keer te bewerken.
• Het is alsof je in plaats van één voor één blokken te verplaatsen, een hele pallet met blokken in één keer verplaatst.

Het Geniale Detail: De Hybrid
Het meest slimme aan deze nieuwe methode is dat je het hybride kunt maken.

• Je kunt de delen van je simulatie waar het rooster gelijkmatig is, laten doen door de snelle "dans" (FFT).
• Je kunt de delen waar het rooster vervormd is (bij de wanden), laten doen door de krachtige "vrachtwagen" (GEMM).

Ze kunnen deze twee methoden door elkaar gebruiken in dezelfde simulatie, zonder dat de computer in de war raakt.

Waarom is dit belangrijk?

Stel je voor dat je een auto wilt bouwen.

• De oude methode: Je bouwt de hele auto van ultra-dun, duur titanium, omdat je bang bent dat de wielen te dun zijn. Het is sterk, maar extreem duur en zwaar.
• Deze nieuwe methode: Je gebruikt titanium alleen waar het nodig is (bij de wielen) en goedkoop staal waar het kan (de carrosserie). Je gebruikt een slimme machine om deze verschillende materialen perfect aan elkaar te lassen.

De resultaten:

• Snelheid: Op één computercore is deze nieuwe methode tot 100 keer sneller dan de oude, trage methoden op vervormde roosters.
• Schalen: Als je duizenden computers (of GPU's) tegelijk gebruikt, werkt deze methode veel beter. De "vrachtwagen" (GEMM) is zo zwaar en efficiënt, dat de tijd die verloren gaat aan het organiseren van de data (communicatie) minder opvalt.
• Toekomst: Hierdoor kunnen wetenschappers veel gedetailleerdere simulaties doen van turbulente stromingen (zoals bij vliegtuigen of in de atmosfeer) zonder dat hun computers het laten afweten.

Kortom:
Ze hebben een manier gevonden om de wiskundige "puzzel" van vloeistoffen op te lossen, zelfs als het rooster ongelijkmatig is. Ze doen dit door de delicate dans (FFT) te combineren met de brute kracht van moderne rekenmachines (GEMM), waardoor simulaties sneller, goedkoper en realistischer worden.

Technische samenvatting

Hier is een gedetailleerde technische samenvatting van het paper in het Nederlands:

Titel: Een op GEMM gebaseerde directe solver voor eindige-differentie Poisson-problemen op niet-uniforme roosters

1. Het Probleem

Het numeriek oplossen van de driedimensionale Poisson-vergelijking is een fundamenteel en vaak kostbaar onderdeel van simulaties van onsamendrukbare stromingen (Navier-Stokes vergelijkingen), waarbij de druk moet worden bepaald om aan de divergentievrije voorwaarde te voldoen.

• Huidige uitdaging: Traditionele directe solvers die gebaseerd zijn op snelle Fourier-transformaties (FFT) vereisen uniforme roosters in de richtingen waar de transformatie wordt toegepast. Voor simulaties met sterk vervormde roosters (bijv. om wanden op te lossen in turbulente stromingen) zijn deze methoden niet direct toepasbaar.
• Alternatieven: Geometrische multigrid-methoden kunnen worden gebruikt voor niet-uniforme roosters, maar ze worden vaak inefficiënt op niet-uniforme netten (coarsening werkt minder goed) en lenen zich slecht voor moderne GPU-architecturen vanwege de inherente reductiestappen die het aantal onbekenden per niveau verkleinen, wat leidt tot onderbenutting van GPU-resources.
• Doel: Er is behoefte aan een solver die de flexibiliteit van niet-uniforme roosters combineert met de efficiëntie van directe methoden en goed schaalbaar is op moderne, heterogene hardware (CPU en GPU).

2. Methodologie

De auteurs presenteren een directe solver die een tensorgebaseerd raamwerk gebruikt, waarbij de Poisson-operator wordt gediagonaliseerd in twee richtingen en direct wordt opgelost in de derde richting.

• Eigendecompositie van 1D-operatoren:
  • Voor niet-uniforme roosters is de discretiseringsmatrix (tridiagonaal) niet-symmetrisch. Om efficiënte eigendecompositie mogelijk te maken, symmetriseren de auteurs de operator via een diagonale schaling (similariteitstransformatie).
  • De operator $T$ wordt getransformeerd naar een symmetrische matrix $\tilde{T} = D^{1/2} T D^{-1/2}$, waarbij $D$ een diagonale matrix is met de celbreedtes.
  • Omdat $\tilde{T}$ symmetrisch is, kan deze efficiënt worden ontbonden in $\tilde{T} = \tilde{Q}\Lambda\tilde{Q}^T$. De oorspronkelijke eigenvectoren worden vervolgens hersteld via $Q = D^{-1/2}\tilde{Q}$.
• 3D Solver Strategie:
  • De 3D Poisson-vergelijking wordt ontbonden in twee richtingen (bijv. $x$ en $y$) en direct opgelost in de derde ($z$) met een tridiagonale solver (TDMA).
  • De transformaties in de $x$- en $y$-richtingen worden uitgevoerd als GEMM (General Matrix-Matrix Multiplication) operaties. In plaats van batches van 1D FFT's, worden de onafhankelijke 1D transformaties samengevoegd tot één grote matrix-matrix vermenigvuldiging.
  • Hybride aanpak: Het algoritme ondersteunt een hybride configuratie: uniforme richtingen kunnen FFT's gebruiken, terwijl niet-uniforme richtingen GEMM's gebruiken. De domeindecompositie en communicatiepatronen blijven hierdoor identiek aan bestaande FFT-solvers.
• Implementatie:
  • Geïmplementeerd als een extensie van de CaNS solver (geschreven in Fortran met MPI en OpenACC/CUDA).
  • Gebruikt bibliotheekfuncties zoals LAPACK (voor eigendecompositie), cuBLAS/OpenBLAS (voor GEMM) en cuDecomp/2DECOMP&FFT (voor transposes).

3. Belangrijkste Bijdragen

• Generalisatie van FFT-methoden: De methode generaliseert de klassieke eigenfunctie-expansie naar niet-uniforme roosters zonder de parallelle structuur (pencil-decompositie) te veranderen.
• GEMM-gebaseerde transformaties: Door de transformaties te herschrijven als dense matrix-matrix vermenigvuldigingen, wordt de rekenintensiteit verhoogd. Dit is cruciaal voor GPU-efficiëntie, omdat GEMM-kernels zeer goed geoptimaliseerd zijn en een hoge rekenintensiteit hebben, waardoor communicatie- en transposetijden beter worden geamortiseerd.
• Hybride flexibiliteit: De solver kan per richting kiezen tussen FFT (voor uniforme roosters) en GEMM (voor niet-uniforme roosters), wat een optimale balans biedt tussen snelheid en flexibiliteit.

4. Resultaten

De prestaties zijn getest op CPU's (Snellius supercomputer) en GPU's (NVIDIA GB200), vergeleken met geometrische multigrid (MG) en block cyclic reduction (BCR).

• Single-Core CPU:
  • De directe methode (GEMM/FFT) is 1 tot 2 orders van grootte sneller dan geometrische multigrid op sterk gestreepte roosters.
  • Multigrid presteert slecht bij sterke roostergradiënten omdat de convergentie vertraagt.
• Strong Scaling (CPU):
  • De GEMM-varianten (GG) behouden een hogere parallelle efficiëntie dan de FFT-varianten (FF) bij hoge core-aantallen.
  • De berekeningsintensiteit van GEMM amortiseert de communicatiekosten (transposes) effectiever. Bij 8192 cores is de efficiëntie van de GEMM-variant aanzienlijk beter dan die van de FFT-variant.
• Weak Scaling (CPU):
  • Er is een duidelijke afweging: FFT-varianten vertonen een milde toename in tijd bij groeiende roosters ($O(N \log N)$), terwijl GEMM-varianten sneller groeien ($O(N^2)$).
  • Echter, door het gebruik van niet-uniforme roosters kan het totale aantal roosterpunten met een factor 2-3 worden verminderd, wat de extra kosten van GEMM compenseert.
• GPU Prestaties:
  • Op een enkele GPU is de volledig niet-uniforme configuratie (GG) ongeveer 2.8x langzamer voor de Poisson-solver dan de uniforme (FFT) versie, maar voor de volledige Navier-Stokes-stap is de vertraging slechts 1.8x.
  • Bij strong scaling op 64 GPUs (via NVLink) behouden de GEMM-varianten een parallelle efficiëntie van 45-66% en presteren ze beter dan de FFT-varianten door de hogere rekenintensiteit.
  • De methode schaal goed tot grote schalen, waarbij de GEMM-varianten minder gevoelig zijn voor communicatie-overhead.

5. Betekenis en Conclusie

De paper demonstreert dat het gebruik van dense lineaire algebra (GEMM) in plaats van FFT's een krachtige strategie is voor het oplossen van Poisson-vergelijkingen op niet-uniforme roosters in moderne HPC-omgevingen.

• Robuustheid: De methode biedt een robuust alternatief voor multigrid op sterk vervormde roosters, waar multigrid vaak faalt of inefficiënt wordt.
• Hardware-efficiëntie: Door de verschuiving naar GEMM wordt de solver beter afgestemd op de architectuur van moderne GPU's, die zijn ontworpen voor hoge rekenintensiteit.
• Praktische toepassing: Voor turbulente stromingen (DNS) waarbij wanden opgelost moeten worden, maakt deze methode het mogelijk om roosters te gebruiken die lokaal fijn zijn waar nodig en grof elders, zonder de enorme kosten van een volledig uniform fijn rooster. De "break-even" punt ligt vaak al bij een reductie van het totale aantal roosterpunten met een factor 2-3.

Kortom, deze solver combineert de nauwkeurigheid en snelheid van directe methoden met de flexibiliteit van niet-uniforme roosters, en is geoptimaliseerd voor de heterogene hardware van de toekomst.