20 Jahre SeisSol: Erdbebenwellen und mehr berechnen

Erdbeben entstehen in der Regel dort, wo unterirdisch tektonische Platten in Bewegung sind. Können diese nicht aneinander vorbeigleiten, baut sich an den Rändern Spannung auf. Gesteinsschichten geben nach und die gespeicherte Energie abrupt frei. In diesem Fall sprechen Geophysikerinnen von Transformbeben. Zerstörerischer wirken Subduktionsstörungen, dabei schiebt sich nach einer Kollision eine Platte unter die andere. Wer die Ursachen dieser Phänomene erforscht, erkundet neben seismischen Wellen auch Spannung und Bruchenergien, außerdem die Verformungen von Gestein, Temperatur, Untergrund und viele Voraussetzungen mehr. Und für all diese Phänomene gibt es Formeln: „SeisSol verwendet eine Discontinuous Galerkin-Diskretisierung, die bis 2006 zwar in der numerischen Strömungsmechanik, aber noch nicht für geophysikalische Anwendungen eingesetzt wurde“, erläutert die Geophysikerin Alice-Agnes Gabriel, heute Professorin am Scripps Institution of Oceanography an der University California in San Diego.

Sie stieß 2012 nach ihrer Promotion zum SeisSol-Team an der LMU, erarbeitete verschiedene Simulationen, integrierte und optimierte dabei grundlegende Formeln und habilitierte sich 2022 in diesem Themenfeld. Ursprünglich geht der Code jedoch auf den Mathematiker und Geophysiker Martin Käser zurück, der gemeinsam mit dem Numeriker Michael Dumbser an der LMU die Ausbreitung seismischer Wellen rechnerisch analysierte. Neben Differentialgleichungen, mit denen sie Fragen der Wellen- und Erdbebenbruchausbreitung lösten, setzten sie dabei auf die Diskontinuierliche Galerkin-Methode (DG) und kombinierten diese mit einem neuen Zeitschrittverfahren, den Arbitrary-High-Order Derivates. Hinter der ADER-DG-Methode stecken Gleichungen zur Darstellung von seismischen Wellen: Diese Gleichungen werden erstens für eine präzisere Berechnung von Wellen in viele kleine räumliche und zeitliche Abschnitte gestückelt; so lassen sich, zweitens, die vielen Rechenschritte effizient auf tausenden Knoten von High-Performance Computing-Systemen (HPC) verteilen. Den Unterschied in den Ergebnissen von der DG- sowie der ADER-DG-Methode können sich Nicht-Mathematiker vorstellen wie der zwischen Foto und Film: Liefern DG-Rechnungen viele Einzelbilder, schafft ADER-DG fließende Übergänge. Wellen werden realistischer dargestellt, es wird sichtbar, was Erschütterungen beschleunigt: „Es geht bei SeisSol darum, anhand von Sensordaten und Messungen die grundlegende Physik von Erdbeben und Tsunamis zu simulieren und die dabei ablaufenden Prozesse immer besser zu verstehen“, fasst Bader zusammen. „Zum Start meiner Professur 2011 ergab sich die Möglichkeit an SeisSol zu arbeiten. Mir gefiel von Anfang an die Kombination aus High Order Numerik, Höchstleistungsrechnen und Geophysik, das motivierte, mich intensiver mit SeisSol zu beschäftigen.“

In den ersten Jahren machte SeisSol als numerisches Verfahren von sich reden. Dieses wurde mit immer mehr Parametern komplexer, stellte Wellenausbreitungen auch dreidimensional dar: „2012 war SeisSol die akkurateste, aber teuerste, computergestützte Methode zur Lösung der Wellengleichung“, berichtet Gabriel. Der Code lieferte Präzision, brauchte aber viel Rechenzeit. Das sollte sich ändern, als sich das interdisziplinäre Team von LMU, TUM und LRZ – gefördert von der Deutschen Forschungsgesellschaft – für die Überarbeitung fand. Die Geophysikerinnen der LMU überlegten damals, wie sie die Entwicklung von Erdbrüchen modellieren und in SeisSol integrieren könnten. Die Informatik-Teams von TUM und LRZ planten den Code außerdem als Open Source-Programm, damit mehr Forschende SeisSol einsetzen konnten: „Der geschlossene, numerisch geprägte Code“, so Gabriel, „wurde zu einem anwendungsorientierten, mit dem wir erstmals ein echtes Erdbeben modellierten.“

Dafür erweiterte das Geophysik-Team das Set der Formeln und Algorithmen um Bruchprozesse. Währenddessen optimierten die Informatiker der TUM mit den Spezialistinnen des LRZ den Code: Mit Hilfe des adaptiven Zeitschrittverfahrens wurden Algorithmen verbessert, leistungskritische Rechnungen indes durch neue Programmierung ersetzt. Die Ein- und Ausgabe-Routinen von SeisSol wurden modernisiert, so dass schrittweise aus dem ehemaligen Fortran-Code bis 2023 endgültig ein C++-Paket mit Spuren von Python wurde. Als erste Simulation auf Petascale-Niveau wurden die Erschütterungen von Landers berechnet: SeisSol erwies sich als äußerst skalierbar und effizient, erledigte auf SuperMUC mehrere Billiarden Gleitkomma-Rechnungen pro Sekunde (PetaFLOPs). Mit dieser Arbeit erreichte das Team 2014 das Finale des renommierten Gordon-Bell-Preises der Association for Computing Machinery. Den Best Paper Award der Supercomputing-Konferenz (SC) in den USA gewann es drei Jahre später mit der bislang längsten und umfangreichsten Simulation des Sumatra-Andamanen-Erdbebens, die ebenfalls SuperMUC berechnete.