Sie sind hier: AMDiS > Die Poissongleichung

Die Poissongleichung

Sonntag 16. Januar 2011 von
Simon Praetorius
Die Poissongleichung gilt als eines der einfachsten Beispiele partieller Differentialgleichungen, die häufig auch als Testfall für verschiedene Algorithmen Verwendung findet. Für diese Gleichung kann man (unter ein paar Voraussetzungen) auch analytische Lösungen angeben, so dass eine Verifizierung der numerischen Ergebnisse möglich ist. Die Poissongleichung soll als erstes Problem in AMDiS implementiert werden.
Dieses einfachste Beispiel einer elliptischen partiellen Differentialgleichung ist benannt nach Simeon Denis Poisson, und lässt sich aufschreiben in der Form
-Δu=f
dabei bezeichnet Δ den Laplace-Operator, u die unbekannte Funktion und f:ℝ→ℝ einen Quellterm (auch bezeichnet als rechte Seite) des Problems. Detailliertere Beschreibung zu Lösung und Anwendung der Gleichung kann auf Wikipedia nachgelesen werden. Zur Diskretisierung mittels FEM muss zuerst ein Gebiet festgelegt werden, in dem die Gleichung definiert wird. Dann sind Bedingungen an die Lösung am Rand des Gebietes zu stellen, um Eindeutigkeit und Existenz einer Lösung zu gewährleisten. Danach muss man die Gleichung in eine schwache Formulierung überführen.

Sei Ω das Gebiet, in dem die PDGL gelten soll, dann bezeichne ∂Ω den Rand des Gebietes. Dort verlangen wird, dass eine Dirichlet-Bedingung gilt, d.h. u|∂Ω=g für eine Funktion g:∂Ω→ℝ. Der Einfachheit halber sei Ω das Einheitsquadrat [0,1]2. Allgemein muss man noch fordern, dass der Rand des Gebietes Ω möglichst (stückweise) glatt ist und auch die Funktionen f und g sollten einigen Glattheitsbedingungen genügen. Dies ist notwendig, um zeigen zu können, dass eine Lösung der Differenzialgleichung existiert und es für die praktische Umsetzung sehr entscheidend. Näheres zu den Bedingungen an die Daten kann in entsprechender Literatur nachgelesen werden (z.B. Christian Großmann, Hans-Görg Roos: "Numerik partieller Differentialgleichungen", Vieweg+Teubner, Auflage 3, 2005)

Zur Beschreibung des Gebietes Ω kann die Macro-Mesh Datei macro.stand.2d verwendet, die sich im Order demos/macro der Standard AMDiS-Installation befindet. Die Poisson-Gleichung in schwacher Form erhält man durch Multiplikation der Gleichung mit einer Testfunktion, Integration über das Gebiet Ω und partieller Integration der zweiten Ableitungen. Wegen der Dirichlet-Randbedingung ist die Testfunktion v so gewählt, dass sie am Rand von Ω verschwindet, also v|∂Ω=0. Es ergibt sich
(∇u, ∇v)Ω = (f,v)Ω, ∀v mit v|∂Ω=0
wobei das Skalarprodukt definiert ist durch (a,b)Ω:=∫Ωa⋅b dx. Die linke Seite der Gleichung wird oft als Operatormatrix und die rechte Seite als Lastvektor bezeichnet.

Die Implementierung der Poisson-Gleichung in AMDiS wird mit folgendem Hauptprogramm realisiert (poisson.cc)
c++ Code
  • #include "AMDiS.h"
  • #include "main.h"
  •  
  • using namespace AMDiS;
  • using namespace std;
  •  
  • int main(int argc, char* argv[])
  • {
  • FUNCNAME("main");
  •  
  • // ===== Name einer Parameter-Datei als erstes Argument an das Programm erwarten =====
  • TEST_EXIT(argc == 2)("usage: ellipt initfilen");
  •  
  • // ===== Parameter des Problems aus Datei einlesen =====
  • Parameters::init(true, argv[1]);
  •  
  • // ===== Ein skalares Problem (nur eine unbekannte, eine Gleichung) definieren =====
  • ProblemScal ellipt("ellipt");
  • ellipt.initialize(INIT_ALL);
  •  
  • // === Für automatische Gitteradaption verwendete Objekte ====
  • AdaptInfo adaptInfo("ellipt->adapt");
  • AdaptStationary adapt("ellipt->adapt", ellipt, adaptInfo);
  •  
  • // ===== Operatormatrix aufbauen =====
  • Operator matrixOperator(ellipt.getFeSpace());
  • matrixOperator.addSecondOrderTerm(new Laplace_SOT);
  • ellipt.addMatrixOperator(matrixOperator);
  •  
  • // ===== Lastvektor aufbauen =====
  • Operator rhsOperator(ellipt.getFeSpace());
  • rhsOperator.addZeroOrderTerm(new CoordsAtQP_ZOT(new F));
  • ellipt.addVectorOperator(rhsOperator);
  •  
  • // ===== Randbedingungen definieren =====
  • ellipt.addDirichletBC(1, new G);
  •  
  • // ===== Berechnung starten =====
  • adapt.adapt();
  •  
  • // ===== Ausgabe des Ergebnisses =====
  • ellipt.writeFiles(adaptInfo, true);
  • }
Zu allererst muss die AMDiS-Bibliothek eingebunden werden mittels
#include "AMDiS.h"
Dabei muss beim Kompilieren der entsprechende include-Pfad zu AMDiS.h gesetzt sein. Durch AMDiS werden einige C-Macros bereit gestellt, die hier und da Anwendung finden. Die gängigsten sind FUNCNAME("..."), wodurch der Funktionsname der aktuellen Methode gesetzt wird und bei Meldungen/Ausgaben durch andere Macros entsprechend mit ausgegeben wird. Ein Macro zur Ausgabe ist TEST_EXIT(bedingung)(text), das 'text' ausgibt, wenn 'bedingung' NICHT erfüllt ist. Bei Nichterfüllung der Bedingung wird außerdem das Programm abgebrochen. Im obigen Fall wird überprüft, ob genau zwei Komandozeilenargumente (Dateiname und erstes Argument = Name der Parameterdatei) angegeben sind. Weitere Macros in AMDiS sind im Artikel 'C-Macros zur Vereinfachung der Programmierung' beschrieben.

Nach der Initialisierung der Parameter (Parameterdatei wird eingelesen. Darin stehen Angaben zu Polynomgrad, Gitter, Verfeinerung, Gleichungssystemlöser, ...), wird ein Problem definiert. Das erste Argument ist der Name des Problems, der dann auch als Schlüsselwort in der Parameterdatei verwendet wird. Der Aufruf von initialize(INIT_ALL) erzeugt ein Gitter, den Finite-Elemente-Raum, den Lösungsvektor und die Systemmatrix/Rechte-Seite-Vektor, einen Fehlerschätzer und Markierer für Gitteradaption, auf Grundlage der Werte in der Parameterdatei für das definierte Problem. Bei mehreren (gekoppelten) Problemen, kann für jedes Problem individuelle Objekte erzeugt, oder auch mit dem ersten Problem Daten geteilt werden, z.B. durch
problem1.initialize(INIT_ALL);
problem2.initialize(init_flag, problem1, adopt_flag)
Wobei in 'init_flag' die Flags gesetzt werden, welche definieren was nur für das problem2 erzeugt werden sollen und durch 'adopt_flag' angegeben wird, was von problem1 übernommen wird. Damit kann man z.B. erzwingen, dass beide Probleme auf dem selben Gitter und mit den selben Finite-Elemente-Räumen rechnen, aber z.B. die Löser unterschiedlich sind. Die möglichen Flags findet man in der AMDiS-Datei 'ProblemStatBase.h' im AMDiS-include Ordner.

Die Objekte AdaptInfo und AdaptStationary speichern alle Informationen über den Lösungsprozess (auch aus der Parameterdatei eingelesen), z.B. Toleranzen, maximale Iterationszahlen für Gitteradaption, im Fall von instationären Problemen auch noch Zeitschrittweiten und die aktuelle 'Zeit' des Problems. Diese Objekte werden benötigt, auch wenn keine Gitteradaption oder nur ein stationäres Problem gelöst wird.

Nach der Initialisierung des Programms wird nun das Poisson-Problem implementiert. AMDiS macht dies in Form von Operatoren, die in die Matrix oder die rechte Seite der Gleichung eingefügt werden. Es gibt eine ganze Liste von implementierten Operatoren, die in der Latex-Datei operatorTerms.tex im Ordern doc dokumentiert sind. mit den Befehl
# pfclatex operatorTerms.tex
kann unter Linux eine PDF-Dokumentation der Operatoren erzeugt werden. Beispiele für Operatoren sind der Laplace-Operator, ein Transport-Term oder eine lineare Funktion der Lösungsvariablen. Bei der Implementierung geht man, wie für FEM-Formulierungen üblich, von der schwachen Form der Gleichungen aus. Wird ein Laplace-Operator hinzugefügt, dann heißt das, dass das Skalarprodukt (∇u, ∇v)Ω implementiert wird. Man unterscheidet zwischen zero-order, first-order und second-order Termen. Bei zero-order Termen steht die Unbekannte u und die Testfunktion ohne Ableitung. Es kann dann die Multiplikation der Unbekannten mit irgendeiner Funktion, z.B. der Term (h(x)*u, v)Ω, implementiert werden. Wird ein zero-order Term auf die rechte Seite (den Lastvektor) eingefügt, dann entfällt die Unbekannte u im Skalarprodukt, im Beispiel also (h(x), v)Ω, außer man setzt mit dem Befehl
rhsOperator.setUhOld(dofvector);
einen (Pointer auf einen) DOFVector, der anstelle der Unbekannten u eingesetzt wird, z.B. in einem instationäre Problem die Lösung aus dem letzten Zeitschritt.
First-order Terme enthalten eine Ableitung auf der Unbekannten u ODER auf der Testfunktion v, z.B. (n*grad(u), v)Ω oder (n*u, grad(v))Ω, wobei ein zusätzlicher Flag entscheidet, was von beiden implementiert werden soll: GRD_PHI entspricht der Ableitung auf der Unbekannten u, GRD_PSI auf der Testfunktion v. (Näheres dazu in einem Artikel, wenn ein solcher Term in einem Problem auftaucht)
Second-order Terme enthalten eine Ableitung auf der Unbekannten und eine auf der Testfunktion, also eine schwache zweite Ableitung. Als Beispiel sei hier der Laplace_SOT, wie oben beschrieben, genannt.

Um eine Randbedingung für das Problem zu definieren, hier Dirichlet-Randbedingung mit Funktion g, wird die Funktion addDirichletBC(RandNummer, Funktor) aufgerufen. Die RandNummer entspricht der Nummer, die in der Macro-Datei für den Rand gesetzt wurde. Damit sind verschiedene Randbedingungen an verschiedenen Rändern möglich, wenn unterschiedliche Nummern gesetzt sind. Der Funktor entspricht einer AbstractFunction mit Argumenttype WorldVector und Rückgabetyp double, also einer Funktion g:ℝn→ℝ, die dann aber nur am Rand ausgewertet wird. Wie andere Randbedingungen implementiert werden können, wird im Artikel 'Einbindung von Randbedingungen'.

Mit dem Aufruf adapt.adapt() des Objekts AdaptStationary wird der stationäre Lösungsprozess gestartet.

Die Funktionen (bzw. Funktoren) F und G werden ausgelagert in eine zweite Datei (main.h)
c++ Code
  • // Dirichlet Randfunktion g
  • class G : public AbstractFunction<double, WorldVector<double> >
  • {
  • public:
  • double operator()(const WorldVector<double>& x) const
  • {
  • return exp(-10.0 * (x * x));
  • }
  • };
  •  
  • // Function f auf der rechten Seite der Poisson-Gleichung
  • class F : public AbstractFunction<double, WorldVector<double> >
  • {
  • public:
  •  
  • F : AbstractFunction<double, WorldVector<double> >(4) {}
  •  
  • double operator()(const WorldVector<double>& x) const
  • {
  • int dow = x.getSize();
  • double r2 = (x * x);
  • double ux = exp(-10.0 * r2);
  • return -(400.0 * r2 - 20.0 * dow) * ux;
  • }
  • };
Eine AbstractFunction ist hier ein Funktion mit Argumenttyp WorldVector (das entspricht einem Element des n) und Rückgabetyp double. Die Argumente müssen als const deklariert werden und auch der überladene Operator 'operator()' darf die Klasse nicht verändern (Schlüsselwort const). Der Konstruktor der Klasse AbstractFunction kann noch ein Argument aufnehmen. Wie in der Klasse F wird ein int-Wert übergeben (hier AbstractFunction<...>(4)), der angibt, welche Quadratur-Formel bei der Auswertung verwendet werden soll, genauer: bis zu welchem Polynomgrad die Quadraturformel Polynome exakt integrieren soll. Der default-Wert ist 0.

Je nachdem, wie AMDiS kompiliert wurde, muss man nun einfach ein MakeFile aus den Demos anpassen, oder die CMakeFiles.txt in den Ordner von poisson.cc kopieren und dann "ccmake ." ausführen, um die Kompiliervorgang zu konfigurieren.

Die Lösung des oben beschriebenen und implementierten Problems sieht nach Ausführung von "./poisson init/ellipd.dat.2d" folgendermaßen aus (Parameterdatei aus den Beispielen verwendet):
Lösung der Poissongleichung
Lösung der Poissongleichung mit Adaptivem Gitter
Darstellung der Lösung der Poissongleichung als Warp
Der Befehl "ellipt.writeFiles(adaptInfo, true);" erzeugt eine ".vtu"-Datei für die Lösung der PDGL. Diese kann dann z.B. mit dem freien Tool Paraview geöffnet werden.
Besucher: 6089 | Permalink | Kategorie: AMDiS
Tags: , ,
Kommentar hinzufügen