Mehrgitterverfahren

Übersicht:

Schnellere Simulationen durch neue Lösertechnologien

Computersimulationen in der Industrie

Engpässe bei der numerischen Simulation

Schnellere Simulationen erfordern hierarchische Verfahren

Kombination hierarchischer Verfahren mit parallelen Rechnerarchitekturen

Entwicklung der Mehrgitterverfahren

Bett-Federrost

EMGM - Ein Mehrgitterverfahren zur Berechnung von Resonanzschwingungen am Beispiel des Bodensees

Schnellere Simulationen durch neue Lösertechnologien

In der Designphase von Automobil- und Flugzeugbau setzen die Hersteller immer feinere und genauere Computersimulationen ein. Gleichzeitig wächst auch die Zahl der erforderlichen Simulationsläufe stark. Das Problem bei aufwendigen Computersimulationen ist die Geschwindigkeit, mit der die Programme große Gleichungssysteme lösen können. Diese Gleichungssysteme mit vielen Millionen Unbekannten sind es, die komplexen Simulationen zu Grunde liegen und sie so aufwendig machen. Kombiniert man nun moderne Löser mit Parallelrechnertechnik, so kann in vielen Anwendungen die Zeit für eine Simulation auf einen Bruchteil sinken.

Computersimulationen in der Industrie

Strömungsmechanik und Strukturanalyse sind zwei typische Anwendungen, in denen die Industrie immer häufiger auf Computersimulationen setzt. Dadurch ist es möglich immense Summen einzusparen, die reale Tests verschlingen würden und durch den Zeitgewinn auch die langjährigen Designphasen entsprechend zu verkürzen. Darüber hinaus lassen sich durch reale Experimente längst nicht mehr alle wesentlichen physikalischen Parameter hinreichend genau bestimmen. So sind die Möglichkeiten von Windkanalsimulationen so weit ausgeschöpft, dass im Automobilsektor weitere Designverbesserungen allein durch diese Experimente nicht mehr praktikabel sind.

Bei einer Simulation im Computer repräsentieren komplexe Gitter die zu analysierenden Strukturen (ein Beispiel aus dem Automobilbereich zeigt Bild 1) . Die Komplexität der Struktur und die Anforderungen an die Genauigkeit bestimmen wesentlich die Auflösung dieser Gitter: Je höher die geforderte Simulationsgenauigkeit, desto feiner ist das zugrundeliegende Gitter, d.h. desto größer ist die Zahl der Gitterzellen. Heute führt dies oft zu Gittermodellen mit enorm feiner Auflösung. Bild 2 zeigt das Gittermodell eines Autos aus etwa 10 Millionen Gitterzellen. Dies ist erforderlich, damit eine numerische Windkanalsimulation die praktisch erforderliche Mindestgenauigkeit erzielen kann.

Engpässe bei der numerischen Simulation

Die Genauigkeitsanforderungen an Computersimulationen werden in Zukunft weiter steigen. Gleichzeitig wächst bei der Entwicklung von Prototypen die Zahl der erforderlichen Simulationen stark. Ein Grund hierfür liegt auch in den erweiterten Möglichkeiten, die Computersimulationen bieten. Physikalische Parameter, die in Experimenten nicht oder nur schwer bestimmt werden konnten, sind nun durch eine Berechnung am Computer zugänglich. Komplette Parameterstudien zur Designoptimierung sind überhaupt erst durch Computersimulationen möglich geworden.
Ob Computersimulationen in diesem Ausmaß aber praktikabel sind, hängt nun ganz wesentlich davon ab, wie lange die Berechnung einzelner Simulationen dauert. Dieser Aufwand wiederum ist direkt abhängig davon, wie fein das zugehörige Gittermodell ist. Dies bestimmt nämlich die Größe des Gleichungssystems, das der Computer lösen muss. Die Geschwindigkeit, mit der sich die oft riesigen Gleichungssysteme lösen lassen, begrenzt die Möglichkeiten einer numerischen Simulation ganz wesentlich.

Schnellere Simulationen erfordern hierarchische Verfahren

Numerische Standardverfahren stoßen bei Gleichungssystemen mit hunderttausend Unbekannten an die Grenzen der Effizienz. Bereits heute sind aber Unbekannte in der Größenordnung von vielen Millionen die Regel. Tendenz: stark steigend. Die Forschung der letzten Jahre hat klar gezeigt, dass eine effiziente Lösung derart großer Systeme nur mit hierarchisch operierenden Verfahren möglich ist. Der klassische Mehrgitteransatz bildet den historisch ersten hierarchischen Zugang: Anstatt nur auf dem gegebenen (extrem feinen) Gitter zu arbeiten, arbeiten Mehrgitterverfahren mit mehreren immer gröberen Gittern, um das gegebene Problem in optimaler Geschwindigkeit zu lösen. Optimal bedeutet hier, dass eine Verdopplung der Problemgröße auch die Rechenzeit lediglich verdoppelt. Die mathematisch viel einfacheren Standardverfahren sind nämlich meist weit davon entfernt.

Kombination hierarchischer Verfahren mit parallelen Rechnerarchitekturen

In der Vergangenheit haben sich viele Entwickler von Simulationssoftware sehr stark darauf verlassen, dass die Simulationsgeschwindigkeit durch schnellere, insbesondere parallele Computer wächst. In der Tat sind die erzielten Leistungssprünge in der Hardwareentwicklung sehr beeindruckend. Dem steht allerdings in vielen wichtigen Anwendungsbereichen eine dramatische Verfeinerung der Simulationsmodelle gegenüber. Der Einsatz numerisch nicht-optimaler Lösungsmethoden zur Simulation bedeutet, dass sich der Gewinn an reiner Hardwareleistung (relativ) immer weniger in einem wirklichen Gewinn für die Simulationsgeschwindigkeit ausdrückt.

Wirklich relevante Fortschritte bei der Computersimulation fordern daher eine effiziente Kombination von Fortschritten in der numerischen Algorithmik und der Computerentwicklung.

Entwicklung der Mehrgitterverfahren

Die wesentliche Idee der Mehrgitterverfahren stammt aus dem Wechselspiel zwischen der physikalischen Betrachtung des Problems und der mathematischen Abstraktion. Deshalb betrachten wir zuerst die physikalische Seite der Aufgabenstellung.

Der Anwendungsbereich der Mehrgitterverfahren bezieht sich vor allem auf Schwingungen komplexer Systeme, zB Stabilitätsuntersuchungen an großen Bauwerken wie Brücken, Türmen oder leichten Flächentragwerken, etwa dem Dach eines großen Stadions, aber auch auf die Bestimmung des Luftwiderstandes von Automobilkarosserien und Flugzeugen sowie auf die Bewegung großer Gewässer.

Die Gemeinsamkeit all dieser Systeme ist ihre Schwingfähigkeit. Folglich gibt es in jedem dieser Fälle:

• einen Gleichgewichtszustand:
  
  glatt ruht der See, die Brücke lagert bewegungslos, der Turm steht gerade, ohne durch Belastungen seitens des Winds oder durch Temperaturunterschiede gestört zu werden.
  

• eine rücktreibende Kraft:
  
  eine Kraft, die das System wieder auf den Gleichgewichtszustand hintreibt, wenn es angeregt worden ist. In Bauwerken sind das die elastischen Kräfte des Materials, im See die Schwerkraft und die Druckkräfte, welche den See wieder glätten.
  

• die Trägheit der im System enthaltenen Massen:
  
  Sie sind der Grund dafür, dass das System trotz rücktreibender Kraft nicht gleich wieder in seinen Gleichgewichtszustand zurückgeht, sondern über diesen hinausschießt und zu schwingen beginnt.
  

• äußere Kräfte:
  
  zum Beispiel Wind und Temperaturunterschiede, die dem System die Energie zuführen, wodurch es überhaupt erst in Bewegung versetzt wird.
  

• Reibungskräfte:
  
  Kräfte, die dafür verantwortlich sind, dass elastische und kinetische Energie in Wärme umgesetzt werden, wodurch das System letztendlich wieder den Zustand der geringsten Energie, seinen Gleichgewichtszustand, wieder findet.

Alle diese Eigenschaften finden sich aber auch in viel einfacheren Systemen wie einem Trommelfell, einem Trampolin oder einem Bett mit einem alten Federrost. Die mathematische Darstellung dieser Eigenschaften ist jedoch in allen Fällen von gleicher Art. Folglich sind im Allgemeinen die selben Mittel geeignet, das Durchhängen des Federrosts, die Deformation der Brücke oder den Tidenhub in einer Meeresbucht zu beschreiben.

Bett - Federrost

Im Folgenden soll deshalb das Prinzip des Verfahrens an dem einfachen Modell – Problem, dem Bett – Federrost, erläutert werden.

Wir stellen uns den Federrost als ein gleichmäßiges Rechtecksgitter aus Punkten vor, die mit ihren Nachbarpunkten (rechts, links, oben, unten) durch Federn gekoppelt sind. Dabei stellt die Funktion u(x, y), die Auslenkung des Federrostes gegen die Ebene des Bettrahmens dar. Nun gilt es die Kraft auf jeden einzelnen Punkt zu berechnen. Jede der an einem Punkt ansetzenden Federn trägt eine Kraft bei, die proportional der Differenz der Auslenkungen in ihren beiden Endpunkten ist. Es ergibt sich, wenn die Proportionalitätsfaktoren vernachlässigt werden:

F = (u_rechts - u) + (u_links - u) + (u_oben - u) + (u_unten - u) = u_rechts + u_links + u_oben + u_unten - 4u

Bei dieser Gleichung handelt es sich um eine Differenzengleichung. Für F = 0, also den Gleichgewichtszustand, ergibt sich:

u = (u_rechts + u_links + u_oben + u_unten) / 4

Man erkennt, dass die Auslenkung jedes Punktes das Mittel der Auslenkungen der vier Nachbarpunkte ist.

Betrachtet man nun einen derart ausgelenkten Federrost aus immer größerer Entfernung, so erscheint das Gitter immer feiner und gleicht erst einem Nylonstrumpf, dann einer biologischen Membran, wie zum Beispiel dem Trommelfell, in dem biologische Fasern die Federn des Federrosts darstellen und schließlich einer Seifenhaut, wo es die Kohäsionskräfte zwischen den Molekülen der Seifenlösung sind. In der Mathematik vollzieht man nun einen Grenzübergang, bei dem man den Federrost aus unendlicher Entfernung betrachtet, wodurch sich dieser als unendlich fein gesponnenes Federnetz darstellt. 

Dazu wird der obige Ausdruck für die Kraft ein wenig anders geschrieben und durch h² dividiert, wobei h der Abstand zwischen zwei benachbarten Gitterpunkten ist.

 (1 / h²)(u_rechts - 2u + u_links) + (1 / h²)(u_oben - 2u + u_unten)

  Bildet man nun den Grenzwert für h gegen 0, so führt das schließlich zum Laplace-Operator:

Dadurch ist aus unserer Differenzengleichung eine Differenzialgleichung geworden, also eine Gleichung, in der die unbekannte Größe eine Funktion ist und die eine Beziehung zwischen der unbekannten Funktion und ihren Ableitungen herstellt. Aus unserem diskreten Problem wurde dadurch ein kontinuierliches.

Der Laplace-Operator dient nicht nur der mathematischen Modellierung von Schwingungs-, sondern auch von Diffusions- und Wärmeleitungsprozessen und ist auch in der Quantenmechanik (Schrödinger-Gleichung) sowie in der Elektrodynamik (Maxwell'sche Gleichungen) zu finden.

Für die numerische Berechnung auf dem Computer muss man den umgekehrten Weg gehen und aus der Differenzialgleichung wieder eine Differenzengleichung machen (diskretisieren), da der Computer eine Funktion nur durch endlich viele Zahlen darstellen kann. Das entspricht also der Darstellung einer Membran durch einen Federrost. Dessen elastisches Verhalten ist ähnlich dem eines Nylonstrumpfes oder auch einer Seifenhaut, was mathematisch heißt, dass man das Verhalten eines groben Gitters als Näherung für das Verhalten eines feineren Gitters benutzen kann, welches seinerseits wieder als brauchbare Näherung für das Verhalten des mathematischen Systems selbst geeignet ist. Dieses mathematische System ist es ja eigentlich, was uns interessiert und welches durch den Laplace-Operator beschrieben wird. Dabei ist es günstig, wenn die Diskretisierung physikalisch interpretierbar ist wie in dem Beispiel mit dem Federrost.

Zur Lösung der Differenzengleichungen benötigt man Verfahren zur Lösung großer, linearer und dünnbesetzter Gleichungssysteme. Zwar gibt es dafür bereits das Gauß-Seidel-Verfahren sowie das Jacobi-Verfahren, allerdings verläuft die Konvergenz bei diesen Verfahren verhältnismäßig langsam; denn je größer n wird, desto mehr Schritte sind notwendig, um eine vorgegebene Genauigkeit zu erzielen.

Physikalisch lässt sich dieser Sachverhalt durch den Wärmeausgleichsprozess plausibel machen: Wird ein Stab an beiden Enden auf verschiedenen Temperaturen konstant gehalten, so stellt sich nach einer gewissen Zeit eine stabile Temperaturverteilung ein. Bei diesem Vergleich steht die lokale Temperatur für die Auslenkung u, die Wärme, die im Rahmen des Ausgleichsprozesses auf den Punkt zu oder von ihm wegfließt, repräsentiert die Korrektur durch ein numerisches Verfahren.

In jedem Zeitschritt findet ein Austausch von Wärme (oder auch Information) immer nur unter nächsten Nachbarn statt. Deswegen kann es sehr lange dauern, bis die Information über die Temperatur am Rande des Gebietes sich 100 Gitterpunkte weiter zu einem Punkt im Inneren des Gebietes herumgesprochen hat.

Die Konvergenzgeschwindigkeit eines Systems kann man zwar dadurch erhöhen, dass man die Korrektur jeweils mit einem geeigneten Faktor multipliziert (SOR - Verfahren), aber auch das bringt keine entscheidende Verbesserung; denn in erster Linie reduziert das Iterationsverfahren nicht den Fehler, sondern glättet ihn lediglich.

Hier setzt das Mehrgitterverfahren an, das durch gezieltes Wechseln von Gittern verschiedener Maschenweite eine enorme Konvergenzbeschleunigung erzielt. Man wendet das normale Iterationsverfahren zunächst im feinen Ausgangsgitter an und glättet in 2 bis 3 Schritten den Fehler. Die so erhaltenen Werte werden auf ein gröberes Gitter (Restriktion) übertragen und anschließend das Iterationsverfahren wiederholt. Das gröbere Gitter kann einfach die Teilmenge des feineren sein, bei der jede zweite Gitterzeile und -spalte gestrichen werden. Der Informationsaustausch findet jetzt über größere Entfernungen statt und ist dadurch effektiver.

Außerdem geht auf dieser Stufe das Rechnen schneller, da die Anzahl der Gitterpunkte gesunken ist. Allerdings darf man das Gitter nicht beliebig vergröbern, da zwar eine glatte, also langsam variierende Funktion, noch einigermaßen korrekt dargestellt wird, hingegen eine schnell variierende Funktion stark verfälscht werden kann.

Die Kunst besteht also darin, abzuschätzen, in welchen Bereichen man das Gitter vergrößern kann und in welchen Bereichen man auf dem feinsten Gitter rechnen muss.

Man restringiert wiederum und führt einige Iterationsschritte im gröbsten Gitter aus, wobei wieder einige Gitterpunkte unberücksichtigt gelassen werden. Da die Anzahl der Variablen inzwischen ziemlich klein geworden ist, kann man die zugehörige Gleichung direkt (nicht iterativ) lösen. Die hier errechneten Werte werden beim Übergang in das mittlere Gitter (Prolongation) beibehalten. Die Werte der Gitterpunkte, die im gröbsten Gitter nicht berechnet worden sind, werden jetzt als Mittel der Nachbarpunkte berechnet. Ebenso verfährt man nach Übergang ins feinste Gitter, so dass letztendlich für alle Gitterpunkte ein approximierter Wert berechnet worden ist. Alle bisher durchgeführten Vorgänge auf den drei Gitterstufen bilden zusammen genau einen Schritt des Mehrgitterverfahrens. Typisch sind drei bis sieben Gitterstufen.

Das Mehrgitterverfahren ist wesentlich effizienter als das Gauß-Seidel- oder das Jacobi-Verfahren; denn aufgrund der flexiblen Gitterwahl ist ein besserer Informationsaustausch und damit eine schnellere Konvergenz gewährleistet. Auch ist der Mehraufwand nicht so groß, wie man anfangs vermutet; er liegt lediglich bei Faktor 3, da die Anzahl der zu berücksichtigenden Gitterpunkte bei allen Grobgittern zusammen meist kleiner als die gesamte Anzahl aller Gitterpunkte des feinsten Gitters ist.

EMGM-Ein Mehrgitterverfahren zur Berechnung von Resonanzschwingungen am Beispiel des Bodensees

Physikalischer Hintergrund und mathematische Modellierung

In der Limnologie sind Resonanzschwingungen von Seen ein wohl bekanntes Phänomen. Diese so genannten „Seiches“ sind für die inneren Strömungen in Seen verantwortlich und haben daher eine entscheidende Bedeutung für die mikro-biologischen Verhältnisse in Seen. Auch für die Arbeit der Fischer ist die Kenntnis dieser inneren Strömungsverhältnisse wichtig.

Physiker und Biologen beschäftigen sich schon seit langem mit der Modellierung dieser Resonanzschwingungen. Wegen der viel zu groben Diskretisierung der zugrunde liegenden Gleichungen lieferten diese Berechnungen nur unzureichende Ergebnisse vor allem bei höherfrequenten Schwingungen.

Das Programm EMGM (Eigenvalue Problems with Multi-Grid Methods)

Das Programm EMGM löst Resonanzprobleme in zwei Raumdimensionen, d. h. im Detail, es berechnet die Eigenwerte und Eigenfunktionen elliptischer Differenzialoperatoren zweiter Ordnung.

Verwendung von Mehrgittermethoden

Die entwickelten Mehrgitterverfahren zur Lösung von Eigenwertproblemen sind vom Aufwand her optimal und herkömmlichen Verfahren an Rechenzeit und Speicheraufwand überlegen. Damit ist es möglich, Eigenwerte/vektoren von 10000x10000 Matrizen in wenigen Minuten auf Personal Computern zu lösen.

Konstruktion problemangepasster Gitter für Eigenwertprobleme

Um Resonanzprobleme auf einem Computer zu lösen, müssen die Eigenwerte/vektoren sehr großer Matrizen berechnet werden. Diese Matrizen bestehen fast nur aus Nullelementen. Nullelemente sind von entscheidender Bedeutung für die Geschwindigkeit der Mehrgitterverfahren.

(Beispiel einer derartigen Matrix für das oben genannte Bettrostbeispiel)

Durch spezielle Nummerierungstechniken entstehen Matrizen mit sehr regelmäßigem Besetzungsmuster. Auf Grund dieser regelmäßigen, vektorartigen Datenstruktur lässt sich EMGM sehr effizient auf Vektorrechnern implementieren.

Verwendung spezieller Techniken zur Rand/Uferanpassung der Gitter

Wegen der extrem komplizierten (fraktalen!) Struktur des Ufers benötigt man spezielle Techniken zur Anpassung des Gitters an den Rand.

Mit Hilfe von EMGM konnten die ersten zehn Oberschwingungen des Bodensees berechnet werden:

Erste, zweite, sechste und die zehnte Eigenschwingung des Bodensees (von links oben nach rechts unten).

Die zehnte Eigenschwingung ist historisch als Wasserwunder von Konstanz (1549) belegt.

Weiterführende Links:

http://hodgson.pi.tu-berlin.de/Vorlesungen/tfd_skript/node63.html

ausführliche mathematische Beschreibung

 

http://www.numerik.uni-kiel.de/~tpr/main.pdf

Dissertationsarbeit

 

Quelle:

Spektrum der Wissenschaft Digest - Wissenschaftliches Rechnen; Mehrgitterverfahren, Seite 42 ff.