Fast solvers for Tokamak fluid models with PETSC

Dit artikel presenteert een nieuwe semi-vergrotings geometrische multigrid-oplosser geïmplementeerd in PETSc voor de M3D-C1 tokamak-code, die de convergentiebeperkingen van de bestaande blok-Jacobi-voorwaardeoplosser aanpakt door gebruik te maken van de toroidale roosterstructuur om superieure robuustheid en prestaties te bereiken op complexe magnetohydrodynamische modellen.

De kern

Stel je voor dat je probeert te voorspellen hoe een wervelende, super hete soep van geladen deeltjes (plasma) zich gedraagt binnen een gigantische, donutvormige machine die een tokamak wordt genoemd. Deze machine is ontworpen om kernfusie-energie te creëren, net als de kracht van de zon. Deze "soep" is echter ongelooflijk chaotisch. Als je probeert de beweging stap voor stap met een computer te berekenen, wordt de wiskunde zo ingewikkeld dat de computer vastloopt, of het duurt zo lang dat het antwoord nutteloos is tegen de tijd dat het arriveert.

Dit artikel gaat over het bouwen van een slimmere, snellere rekenmachine voor dit specifieke type probleem.

Hier is de uiteenzetting van wat de auteurs hebben gedaan, met behulp van eenvoudige analogieën:

1. Het Probleem: De "Verkeersopstopping" van Wiskunde

De computercode die ze gebruiken (genaamd M3D-C1) probeert vergelijkingen op te lossen die beschrijven hoe het plasma beweegt. Om dit te doen, moet het een enorm puzzel miljoenen keren oplossen.

• De Oude Methode (Block Jacobi): Stel je voor dat je een enorme kaart van een stad met verkeersopstoppingen hebt. De oude methode was als het vragen aan een ander persoon om het verkeer op slechts één straat tegelijk te repareren, terwijl de rest van de stad wordt genegeerd. Als de stad klein is, werkt dit. Maar naarmate de stad groter wordt (meer "vlakken" of schijven van de donutvorm), kunnen de mensen die de straten repareren niet snel genoeg met elkaar praten. De verkeersopstoppingen worden erger en de oplossing vertraagt of stopt volledig.
• De Specifieke Uitdaging: Het plasma in deze machines is "anisotroop". Denk aan een stapel papier. Het is heel gemakkelijk om een vel papier over het oppervlak te laten glijden (makkelijke richting), maar heel moeilijk om het door de stapel te duwen (moeilijke richting). De oude wiskundige oplosser begreep deze "stapel papier"-structuur niet, dus probeerde het de moeilijke richting en de makkelijke richting op dezelfde onhandige manier op te lossen.

2. De Oplossing: De "Multigrid" Lift

De auteurs bouwden een nieuwe oplosser met behulp van een methode die Multigrid (MG) wordt genoemd.

• De Analogie: Stel je voor dat je probeert een verloren speelgoed te vinden in een enorm, meervoudig verdiept herenhuis.
  • De Oude Methode: Je controleert elke kamer, elke lade en elke hoek op de begane grond voordat je naar boven gaat. Het duurt eeuwen.
  • De Multigrid Methode: Je kijkt eerst naar een miniatuurmodel van het hele herenhuis vanuit een vogelperspectief. Je ziet snel het algemene gebied waar het speelgoed ontbreekt (het "grove" rooster). Vervolgens zoom je in op een middelgrote kaart om het te verkleinen. Tot slot ga je naar de eigenlijke kamer (het "fijne" rooster) om het speelgoed op te pakken.
  • Door het probleem eerst op de "grote lijn"-niveaus op te lossen, weet de oplosser precies waar hij moet zoeken als hij de kleine details bereikt. Dit maakt het ongelooflijk snel.

3. De "Geheime Ingrediënt": Semi-Coarsening

De auteurs realiseerden zich dat de tokamak lijkt op een stapel 2D-schijven (poloidale vlakken) die in een 3D-donut zijn gedraaid.

• Ze pasten hun "Multigrid Lift" specifiek toe op de stapelrichting (de toroidale richting).
• In plaats van te proberen de hele 3D-messigheid in één keer te vereenvoudigen, hielden ze de 2D-schijven gedetailleerd (omdat de vorm van de tokamakwand complex is), maar maakten ze de stapel schijven eenvoudiger naarmate ze hoger kwamen.
• Dit is als het nemen van een dik boek en alleen het aantal pagina's te verminderen terwijl de tekst op elke pagina duidelijk blijft. Het is een perfecte match voor de vorm van de machine.

4. De Resultaten: Snelheid en Betrouwbaarheid

Het team testte deze nieuwe oplosser op twee zeer moeilijke scenario's:

• Scenario A: De "Runaway Electron" (SPARC): Dit simuleert een gevaarlijk evenement waarbij deeltjes oncontroleerbaar versnellen.
  • Resultaat: De nieuwe oplosser was concurrerend met de oude op kleinere opstellingen en veel sneller op de grootste, meest complexe opstellingen. Het loste het probleem op in minder stappen, waardoor tijd werd bespaard.
• Scenario B: De "Stellarator" (Een andere, meer gedraaide machine): Deze geometrie is nog meer gedraaid en onregelmatig dan een standaard donut.
  • Resultaat: De oude oplosser faalde volledig en kon geen antwoord vinden. De nieuwe Multigrid-oplosser slaagde. Het was robuust genoeg om de gedraaide geometrie aan te kunnen die de oude tool brak.

5. De Hardware: Het Gebruik van Supercomputers

Ze draaiden deze tests op Perlmutter, een van 's werelds snelste supercomputers, die zowel krachtige CPU's als GPU's (grafische kaarten) gebruikt.

• Ze ontdekten dat hoewel de "opstelling" (het bouwen van de miniatuurmodellen) duur was, het daadwerkelijke oplossen ongelooflijk snel was op de GPU's.
• Ze ontdekten dat voor de moeilijkste problemen een "zware" smoother (een specifieke wiskundige truc) nodig was om te voorkomen dat de oplosser vastliep, wat iets meer rekenkracht vereiste maar zich terugbetaalde in snelheid.

Samenvatting

Het artikel beweert dat ze door het specifieke "stapel papier"-vorm van fusieplasma-machines te begrijpen, een nieuw wiskundig hulpmiddel (Multigrid) hebben gecreëerd dat:

1. Problemen sneller oplost dan de huidige standaardmethode op grote, complexe simulaties.
2. Niet crasht op gedraaide, complexe vormen waar de oude methode faalt.
3. Een cruciale eerste stap is naar het praktisch en snel genoeg maken van fusie-energiesimulaties om echte krachtcentrales te helpen ontwerpen.

Ze hebben niet beweerd dat dit kernfusie-energie zelf oplost, maar eerder dat het de snelle, betrouwbare rekenmachine biedt die nodig is om de fysica te simuleren die uiteindelijk zal leiden tot kernfusie-energie.

Technische samenvatting

Technische Samenvatting: Snelle Oplossers voor Tokamak-Fluidummodellen met PETSc

Probleemstelling
De M3D-C1-code is een primair hulpmiddel voor het simuleren van uitgebreide magnetohydrodynamische (MHD) plasma's in fusie-energieonderzoek, specifiek gericht op grootschalige instabiliteiten en disrupties in tokamaks en stellaratoren. De code maakt gebruik van een semi-impliciete tijdsintegrator die het oplossen van talrijke lineaire systemen per tijdstap vereist. Onder deze systemen wordt de impliciete voortgang van de impulsvergelijking geïdentificeerd als het meest uitdagende en minst robuuste onderdeel.

De huidige productieoplosser maakt gebruik van een Krylov-methode (FGMRES/GMRES) voorvoorwaardelijk gesteld door een Block Jacobi (BJ) methode, waarbij blokken vrijheidsgraden groeperen op vlakken met een constante toroïdale coördinaat. Hoewel effectief voor kleinere configuraties, verslechtert de convergentie van de BJ-oplosser aanzienlijk naarmate het aantal toroïdale vlakken toeneemt, als gevolg van de spectrale eigenschappen van de voorvoorwaardelijk gestelde matrix. Bovendien faalt de BJ-oplosser volledig om te convergeren op niet-axiale stellaratorconfiguraties, wat de toepasbaarheid van de code op complexe fusie-inrichtingen beperkt.

Methodologie
Dit werk ontwikkelt een geometrische multigrid (MG) oplossing die is toegesneden op de specifieke discretisatie en rasterstructuur van M3D-C1. De methodologie maakt gebruik van de volgende structurele kenmerken:

• Rasterstructuur: M3D-C1 gebruikt ongestructureerde 2D-rasters in het poloidale vlak (uitgezet rond de torus) gecombineerd met een regelmatig 1D-raster in de toroïdale richting.
• Semi-vergroving: De oplossing maakt gebruik van geometrische semi-vergroving specifiek in de toroïdale richting. Deze aanpak is gekozen omdat het toroïdale raster gestructureerd is, wat standaard geometrische MG-technieken (2:1-vergroving) mogelijk maakt, terwijl het complexe poloidale vlak onvergroefd blijft voor toekomstig werk.
• Algebraïsche interface: Hoewel de toepassing Jacobianen berekent op fijne rasters, worden de grof-rasteroperatoren algebraïsch geconstrueerd met behulp van een Galerkin-proces ($A_{i+1} = P^T A_i P$), waarbij $P$ een 3-punts stencil prolongatie-operator is. Dit creëert een hybride geometrisch/algebraïsche oplossing met een puur algebraïsche interface die compatibel is met PETSc.
• Gladders: Twee gladders worden onderzocht:
  1. Point-Block Jacobi (PBJ): Berekent expliciete dichte inversen van diagonaablokken (grootte 36 per vertex). Dit is snel en geheugenefficiënt, geschikt voor GPU-versnelling.
  2. Block Jacobi (BJ): De bestaande productieoplosser wordt hergebruikt als gladder. Dit omvat volledige LU-factoren op poloidale vlakken. Hoewel duur, wordt het gebruikt als post-gladder in een $V(0,1)$-cyclus voor de moeilijkste problemen waarbij PBJ stagneert.
• Implementatie: De oplossing is geïmplementeerd binnen de PETSc-bibliotheek, met gebruikmaking van diens Krylov-oplossers en abstracte interfaces voor lineaire algebra-back-ends (inclusief MUMPS en SuperLU voor directe oplossingen).

Belangrijkste Resultaten
De oplossing werd geëvalueerd op twee testgevallen met de Perlmutter-supercomputer (CPU- en GPU-nodes): een runaway-elektronmodel van een SPARC-tokamak-disruptie en een quasi-axiale stellaratormodel.

1. SPARC Tokamak (Schaalbaarheid en Snelheid):

   • Op configuraties met 4 tot 16 toroïdale vlakken convergerde de MG-oplosser (met PBJ-gladding) in één iteratie wanneer de voorvoorwaarde werd ververst, wat uitstekende algoritmische schaalbaarheid demonstreert.
   • Op de grootste configuratie (64 vlakken) nam de probleemcomplexiteit toe, wat een strakkere oplossertolerantie ($rtol = 10^{-7}$) en een agressievere $V(6,6)$-cyclus met frequente verversingen vereiste.
   • Prestaties: De MG-oplosser behaalde concurrerende prestaties in vergelijking met de BJ-oplosser. Op het 64-vlakken geval reduceerde de MG-oplosser het totale aantal Krylov-iteraties met een factor van ongeveer 16x ten opzichte van de BJ-oplosser. Hoewel de totale oplostijd op GPU's ongeveer vergelijkbaar was, toonde de MG-oplosser aanzienlijk snellere oplostijden op CPU's.
   • Complexiteit: Theoretische analyse suggereert dat de MG-oplosser een asymptotische werkcomplexiteit heeft van ongeveer 2x die van de BJ-oplosser als gevolg van de geometrische som van werk over rasterniveaus. De waargenomen prestaties geven aan dat grof-rasteroplossingen momenteel slecht presteren, wat ruimte voor optimalisatie suggereert.

2. Stellarator (Robuustheid):

   • Op een niet-axiale stellarator testgeval faalde de standaard Block Jacobi-oplosser volledig om te convergeren.
   • Daarentegen convergerde de MG-oplosser (met BJ als gladder) succesvol. Dit demonstreert een cruciaal robuustheidsvoordeel van de MG-aanpak voor niet-axiale geometrieën waar de spectrale eigenschappen van de matrix de BJ-voorvoorwaarde verslaan.

Betekenis en Claims
Het artikel claimt de eerste stap te presenteren in het toepassen van multigrid-methoden op complexe, wetenschappelijk relevante MHD-tokamakmodellen binnen de M3D-C1-code. De primaire bijdragen zijn:

• Demonstratie van Haalbaarheid: Bewijzen dat multigrid-methoden een pad kunnen bieden naar snellere, schaalbare en robuuste oplossers voor sterk anisotrope, slecht geconditioneerde gemagnetiseerde plasma-modellen.
• Toroïdale Semi-vergroving: Vaststellen dat semi-vergroving in de toroïdale richting een nuttige en effectieve eerste stap is voor gemeenschappelijke tokamak-codes met uitgezette rasterstructuren.
• Robuustheid: Aantonen dat MG-oplossers problemen kunnen oplossen (specifiek stellarator-disrupties) waar de huidige productie Block Jacobi-oplosser faalt.

Toekomstig Werk en Beperkingen
De auteurs blijven bescheiden over de huidige staat van de oplossing, met aandacht voor verschillende verbeteringsgebieden:

• Grof-rasterprestaties: De huidige prestaties worden beperkt door inefficiënte grof-rasteroplossingen, deels als gevolg van redundante parallelisme-eisen en MPI-fouten in directe oplossers (MUMPS/SuperLU) op grove rasters.
• Volledige 3D Multigrid: De meest impactvolle toekomstige optimalisatie is de toevoeging van 2D-multigrid in het poloidale vlak om een volledige 3D MG-oplosser te creëren. Dit zou het hanteren van ongestructureerde geometrische grove rasters vereisen.
• Gladders: Toekomstig werk omvat het optimaliseren van blokgroottes voor PBJ-gladders om beter te passen bij GPU-threadgroepen en het onderzoeken van FieldSplit-voorvoorwaarden die de drie ontkoppelde MHD-golven afzonderlijk behandelen.
• Verversingsstrategieën: Het huidige verversingsmodel (periodieke LU-factoren) is eenvoudig; slimmere algoritmen zijn nodig om quasi-stationaire toestanden efficiënter te hanteren.

Kortom, dit werk vestigt een levensvatbaar, robuust multigrid-kader voor M3D-C1 dat de bestaande oplossing overtreft in termen van iteratieaantallen en robuustheid op moeilijke geometrieën, terwijl het specifieke algoritmische en implementatiehobbels identificeert die moeten worden aangepakt om volledige textbook-multigrid-efficiëntie te bereiken.