Unterabschnitte

• Datenstruktur für die Systemmatrix
• Gleichungslösung durch Gauß'sche Elimination
• Der CG-Gleichungslöser

4.6 Lösung des Gleichungssystems

Zur Berechnung des unbekannten Lösungsvektors (Potenzial, Temperatur) ist es notwendig, ein lineares algebraisches Gleichungssystem zu lösen, das allgemein in folgender Form angeschrieben werden kann:

$$[A]\cdot\{x\}=\{b\} \,.$$ (4.111)

Dabei ist $\{x\}$ der gesuchte Lösungsvektor, $\{b\}$ der sogenannte Rechte-Seiten-Vektor und $[A]$ die Systemmatrix.

Datenstruktur für die Systemmatrix

Der Rang der Matrix $[A]$ ist im Allgemeinen sehr hoch, er entspricht der Anzahl der Gitterknoten abzüglich jener, die auf einem Dirichlet'schen Rand liegen, sowie der bei schwebenden Randbedingungen durch das Zusammenfassen wegfallenden Knoten. Bei typischen dreidimensionalen Simulationen kann der Rang der Matrix ohne weiteres eine Größenordnung von $10^5\ldots 10^6$ erreichen. Andererseits ist die Besetzungsdichte (das Verhältnis zwischen der Anzahl der Matrixeinträge, die ungleich Null sind, zur Gesamtanzahl) sehr gering. In einer Zeile der Matrix sind nämlich genau soviele Einträge ungleich Null wie der Grad^4.7 des zugeordneten Gitterknotens. Da die Knotengrade im Mittel nur von der Art des Gitters, nicht aber von der Größe abhängen, wird mit zunehmenden Rang die Besetzungsdichte immer geringer. Typische Mittelwerte der Knotengrade liegen bei Tetraedergittern mit linearen Ansatzfunktionen etwa bei 11, bei quadratischem Ansatz ungefähr bei 25. So kann man beispielsweise bei einem Gitter mit 250000 Knoten und quadratischem Ansatz eine Besetzungsdichte von etwa 0.01% erwarten. Wenn man nun die gesamte Matrix aus obigem Beispiel in einem zweidimensionalen Feld abspeichern wollte, so bräuchte man dazu, unter der Annahme, dass für eine Gleitkommazahl 8 Byte benötigt werden, einen Speicherplatz von 466 Gigabyte. Angesichts dieses Beispiels erkennt man sofort die Notwendigkeit, die Systemmatrix in einer komprimierten Datenstruktur abzulegen. Würde man in optimaler Weise ausschließlich die Werte ungleich Null in einem eindimensionalen Feld abspeichern, so käme man mit 48 MB aus. Allerdings fehlt dann die Information in welcher Zeile bzw. Spalte der jeweilige Eintrag liegt. Zu diesem Zweck werden noch zwei zusätzliche Indexvektoren angelegt--der erste ordnet jedem Eintrag die Spaltennummer zu, die Elemente des zweiten Vektors zeigen auf die Zeilenanfänge innerhalb des Feldes. Dieses Datenkompressionsverfahren wird MCSR (Modified Compressed Sparse Row) genannt. Der zusätzliche Speicheraufwand beträgt 50%, sodass man in diesem Beispiel auf einen Gesamtspeicherbedarf von 72 MB kommt. Weiters besteht aufgrund der Symmetrie der Systemmatrix die Möglichkeit nur eine Hälfte (Dreiecksmatrix) abzuspeichern, was den Speicherbedarf nochmals auf die Hälfte reduziert (hier 36 MB).

Um gutes Laufzeitverhalten zu gewährleisten, muss die Datenstruktur so aufgebaut sein, dass der verwendete Gleichungslösungsalgorithmus effizient darauf zugreifen kann. Dabei ist es wichtig, dass beim Zugriff auf Matrixelemente keine Verzögerung durch langes Suchen entstehen und unnötige Operationen mit Null-Einträgen vermieden werden. Mit dem MCSR-Format kann eine Matrix$\times$Vektor-Multiplikation sehr effizient durchgeführt werden, wodurch es sich für die Anwendung beim Konjugierten Gradientenverfahren (siehe unten) besonders eignet.

Anmerkung: Theoretisch wäre eine solche Matrix$\times$Vektor-Multiplikation auch ohne vorherige Assemblierung der Systemmatrix, sondern direkt als Summierung über die Elementmatrizen möglich. Diese Methode kommt mit minimalem Speicherbedarf aus, hat aber einen ungleich höheren Rechenaufwand. Es können jedoch Teilgebiete des Gitters zu Teilmatrizen zusammengefasst werden, und die Matrix$\times$Vektor-Multiplikation getrennt über diese Teilmatrizen ausgeführt und anschließend aufsummiert werden. Dieser Algorithmus kann zur Parallelisierung der Finite Elemente Methode ausgenutzt werden, indem man jedem Teilbereich einen eigenen Prozess zuordnet. Die Aufteilung sollte so geschehen, dass die Anzahl der Kopplungen zwischen den Teilbereichen minimal ist, damit wird auch die Kommunikation zwischen den Prozessen minimiert.

Neben dem Konjugierten Gradientenverfahren wird auch das Gauß'sche Eliminationsverfahren zur Lösung des Gleichungssystems herangezogen. Dafür benötigt man ein Matrixformat, das in jeder Zeile zwischen dem ersten und dem letzten von Null verschiedenen Eintrag einen Speicherplatz reserviert, da diese Einträge später bei der Matrixfaktorisierung verwendet werden. Bei einer symmetrischen Matrix reicht es aus, wenn vom ersten von Null verschiedenen Element bis zur Diagonale Speicherplatz reserviert wird. Man spricht dabei von einer hüllenorientierten Struktur, deren Speicherbedarf von den Bandbreiten der einzelnen Zeilen abhängig ist und im Regelfall einem Vielfachen des MCSR Speicherbedarfs entspricht. Eine Reduktion der mittleren Bandbreite lässt sich durch Umordnen der Knotennummerierung erzielen, was in der Matrix ein Vertauschen von Zeilen und Spalten bewirkt (und natürlich auch in den $\{x\}$ und $\{b\}$ Vektoren). Dafür haben sich das Cuthill-McKee-Verfahren (CM) bzw. das Reverse-CuthillMcKee-Verfahren (RCM) etabliert [117]. Die Wirkungsweise des RCM Verfahrens ist in Abb. 4.5 anhand eines Beispiels dargestellt.

Abbildung 4.5: Einträge in der Systemmatrix einer einfachen zweidimensionalen Feldberechnung mit einem Dreiecksgitter und quadratischen Ansatzfunktionen unmittelbar nach der Assemblierung (a) und nach Anwendung des Reverse-Cuthill-McKee Verfahrens (b).

Wie weit die mittlere Bandbreite einer Matrix durch Umsortieren minimiert werden kann, hängt nicht nur vom verwendeten Verfahren ab, sondern in viel stärkerem Ausmaße von der Struktur der Matrix (bzw. des Gitters) selbst. Beispielsweise lassen sich bei der dreidimensionalen Widerstandsberechnung bei sehr langen und dünnen Leitungen sehr geringe mittlere Bandbreiten (60 und darunter) erreichen.

Gleichungslösung durch Gauß'sche Elimination

Als erster Schritt wird die Systemmatrix in ein Produkt aus einer unteren $[L]$ und einer oberen $[U]$ Dreiecksmatrix faktorisiert:

$$[A]=[L]\cdot[U]\,.$$ (4.112)

Danach kann die Lösung $\{x\}$ in zwei Teilschritten (Vorwärtseinsetzen und Rückwärtseinsetzen) berechnet werden:

$$[L]\cdot\{y\} =\{b\}$$ (4.113)

$$[U]\cdot\{x\} =\{y\}$$ (4.114)

Der rechenintensivste Teil ist die Faktorisierung, sie hat bei vollbesetzten Matrizen einen Aufwand von $O(n^3)$, Vorwärtseinsetzen und Rückwärtseinsetzen lediglich $O(n^2)$. Mit geringerer Bandbreite reduziert sich natürlich auch der Rechenaufwand entsprechend.

Für die Berechnung von Kapazitäten, Widerständen, oder transienten Vorgängen muss oft ein System mit gleicher Matrix $[A]$ aber verschiedener rechter Seite gelöst werden. In diesem Fall braucht die (langsame) Faktorisierung nur beim ersten Mal durchgeführt werden und die Lösung $\{x\}$ kann bei jedem weiteren Schritt sehr effizient mittels Vorwärtseinsetzen und Rückwärtseinsetzen berechnet werden.

Zur weiteren Reduktion des Speicherbedarfs kann die Assemblierung und Elimination kombiniert werden: In diesem Fall wird die Systemmatrix nicht in Teilmatrizen aufgespalten sondern es wird komplett mit allen Randknoten gearbeitet. Diese Methode bietet sich besonders an, wenn man nicht an einer Verteilung des Potenzials im Inneren des Simulationsgebietes, sondern lediglich an einer Kapazitäts- oder Widerstandsmatrix interessiert ist. Sobald sämtliche Elemente, die einen bestimmten Gitterknoten gemeinsam haben, assembliert worden sind, kann dieser Knoten aus dem Gleichungssystem eliminiert werden. Randknoten, die zu einer Elektrode gehören, können, da sie auf dem selben Potenzial liegen, ebenfalls zusammengefasst werden. In diesem Fall erspart man sich die Rücksubstitution und erhält direkt die gewünschte (Kapazitäts- bzw. Leitwert-) Matrix.

Der CG-Gleichungslöser

Das konjugierte Gradientenverfahren (CG) ist eine sehr häufig verwendete Methode lineare Gleichungssysteme iterativ zu lösen. Es zeichnet sich besonders dadurch aus, dass es nahezu keinen zusätzlichen Speicherplatz benötigt und im Regelfall mit deutlich weniger Rechenoperationen auskommt als die Gauß'sche Elimination. Allerdings handelt es sich bei der erhaltenen Lösung nur um eine Näherung, da die Iteration abgebrochen wird sobald ein Konvergenzkriterium erfüllt ist^4.8. Voraussetzung ist eine symmetrische und positiv definite Systemmatrix. Der Rechenaufwand des CG-Verfahrens ist proportional zur Wurzel der spektralen Konditionszahl der Matrix, das ist das Verhältnis zwischen größtem und kleinstem Eigenwert.

Näher soll im Rahmen dieser Arbeit nicht auf das CG-Verfahren eingegangen werden--es sei an dieser Stelle auf die ausführliche Darstellung in [118] verwiesen.

Zur Beschleunigung der Konvergenz besteht die Möglichkeit der Vorkonditionierung. Anstatt (4.112) wird das System

$$[M]^{-1}\cdot[A]\cdot\{x\}=[M]^{-1}\cdot\{b\}$$ (4.115)

verwendet. Die Matrix $[M]$ versucht man nun so zu wählen, dass die Konditionszahl des Produktes $[M]^{-1}\cdot [A]$ kleiner ist als die von $[A]$ alleine, etwa durch eine unvollständige Cholesky-Zerlegung [119,82,120].

Der CG-Algorithmus kann auch erweitert werden, sodass für mehrere rechte Seiten ohne großen zusätzlichen Aufwand gleichzeitig Lösungen berechnet werden können [121,122].

Fußnoten

... Grad^4.7

Der Grad eines Knotens gibt an, mit wievielen Nachbarknoten er durch Gitterelemente direkt verbunden ist.

... ist^4.8

Theoretisch liefert das CG-Verfahren die exakte Lösung in höchstens $N$ Schritten, praktisch erhält man aber schon viel früher eine hinreichend genaue Näherungslösung.