Diplomarbeit von Lars Mentrup - M3/Allgemeines · Technische Universit at Munchen Fakult at f ur...

Technische Universitat MunchenFakultat fur Mathematik

Entwicklung einermassenerhaltenden

Semi-Lagrange-Methode zurSimulation von

Spurenstofftransport in derAtmosphare auf einem adaptiven

dreidimensionalen Gitter

Diplomarbeit von Lars Mentrup

Aufgabensteller: Prof. Dr. F. Bornemann

Betreuer: Dr. J. Behrens

Abgabetermin: 5. September 2003

Ich erklare hiermit, dass ich die Diplomarbeit selbstandig und nur mit denangegebenen Hilfsmitteln angefertigt habe.

Munchen, den 5. September 2003

Danksagung:Ich mochte an dieser Stelle Herrn Dr. Jorn Behrens fur seine eingehendeBetreuung und die interessante Themenstellung herzlich danken.

Inhaltsverzeichnis

1 Einleitung 1

2 Numerische Grundlagen 32.1 Raumdiskretisierung mit der Finite-Elemente-Methode . . . . 32.2 Adaptivitat . . . . . . . . . . . . . . . . . . . . . . . . . . . . 92.3 Die Advektionsgleichung . . . . . . . . . . . . . . . . . . . . . 112.4 Diskretisierung mit der Semi-Lagrange-Methode . . . . . . . 132.5 Massenerhaltung . . . . . . . . . . . . . . . . . . . . . . . . . 17

3 Massepakete-Semi-Lagrange-Methode 193.1 Definitionen und Begriffe . . . . . . . . . . . . . . . . . . . . 193.2 Beschreibung der Massepakete-Semi-Lagrange-Methode . . . 20

3.2.1 Die Massepakete-Semi-Lagrange-Methode in 1D . . . 213.2.2 Initialisierung der Masse . . . . . . . . . . . . . . . . . 243.2.3 Unterteilung eines Tetraeders in Massepakete . . . . . 253.2.4 Verfolgung der Triangulierung des nachsten Zeitschritts

stromaufwarts . . . . . . . . . . . . . . . . . . . . . . 313.2.5 Massenzuweisung an die Eckpunkte der Tetraeder des

neuen Gitters . . . . . . . . . . . . . . . . . . . . . . . 313.2.6 Konzentrationsberechnung . . . . . . . . . . . . . . . . 33

3.3 Implementierung: Datenstrukturen und Algorithmen . . . . . 343.3.1 Datenstrukturen . . . . . . . . . . . . . . . . . . . . . 363.3.2 Algorithmus: Unterteilung eines Tetraeders in Masse-

pakete . . . . . . . . . . . . . . . . . . . . . . . . . . . 393.3.3 Algorithmus: Horizontale Suche . . . . . . . . . . . . . 48

4 Ergebnisse: Tests anhand von Modellproblemen 514.1 Darstellung der Testumgebung . . . . . . . . . . . . . . . . . 514.2 Modellproblem I: Lineares Windfeld . . . . . . . . . . . . . . 544.3 Modellproblem II: Zirkulares Windfeld . . . . . . . . . . . . . 594.4 Modellproblem III: Konvergentes Windfeld . . . . . . . . . . 634.5 Fazit . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67

5 Zusammenfassung und Ausblick 68

Notation 70

Literatur 71

1 EINLEITUNG 1

1 Einleitung

Viele physikalische Prozesse lassen sich durch partielle Differentialgleichun-gen modellieren. Hierzu zahlen exemplarisch die Potential-, die Wellen- oderdie Warmeleitungsgleichung. Um den Spurenstofftransport in der Atmo-sphare zu modellieren, verwenden wir hier die Advektionsgleichung, die furlaminare Stromungen eine einfache und gute Modellierung darstellt. Turbu-lente Stromungen sind in diesem Modell nicht darstellbar.

Die Modellgleichungen fur physikalische Prozesse lassen sich meist nichtanalytisch losen. Hier kommt die Numerik ins Spiel. Durch Diskretisierungder Gleichungen uberfuhrt man die unendlich-dimensionalen Probleme inendlich-dimensionale und mithin durch computergestutzte Verfahren losba-re Probleme uber. Bei der Diskretisierung ist darauf zu achten, moglichstviele physikalische Eigenschaften zu ”retten“. Im Fall des Spurenstofftrans-ports ist es z.B. erwunscht, das Prinzip der Massenerhaltung zu beachten.

Um die Diskretisierung durchzufuhren, wird die Semi-Lagrange-Methode(SLM) eingesetzt. Sie stellt eine Kombination aus Zeitdiskretisierung undFinite-Elemente-Methode (FEM) dar. Die SLM ist in dem Programmpaketamatos3d/flash3d von Behrens realisiert und in [Beh96] dokumentiert. Indiesem Programmpaket werden Stromungen der Atmosphare im Raum simu-liert. amatos3d ist ein Finite-Elemente-Gittergenerator fur das dreidimensio-nale Gitter und flash3d die Simulationssoftware, die auf amatos3d aufsetzt.flash3d setzt die Semi-Lagrange-Methode ein, um die Advektionsgleichungzu losen. Im allgemeinen ist die SLM jedoch nicht massenerhaltend. Aufgabedieser Diplomarbeit ist es, ein massenerhaltendes Semi-Lagrange-Verfahrenzu entwickeln und in das Programmpaket einzubinden.Grundlegende Idee des neuen Verfahrens ist, die gesuchten Konzentrati-onswerte nicht mehr direkt wie bei der SLM zu berechnen. Vielmehr wirddas Simulationsgebiet in Masse-Volumen-Einheiten, sog. Massepakete, un-terteilt und deren Bewegung verfolgt, um dann am stromabwarts gelegenenGitterpunkt die Konzentration anhand der Massen und Volumen, die umden Punkt liegen, zu berechnen. Dieses Verfahren bezeichnet der Autor alsMassepakete-Semi-Lagrange-Methode (MPSLM).

Uberblick uber die Arbeit Zunachst wird in Kapitel 2 auf die Diskre-tisierung mit der Finite-Elemente-Methode eingegangen, um darauf hin dieAdvektionsgleichung und ihre Diskretisierung mittels der Semi-Lagrange-Methode vorzustellen. Außerdem wird das Prinzip der Massenerhaltung er-lautert. Im Kapitel 3 wird zunachst ausfuhrlich die Diskretisierung mit derMassepakete-Semi-Lagrange-Methode beschrieben und danach auf Aspekteder Implementierung eingegangen. Anhand von Modellproblemen werden in

1 EINLEITUNG 2

Kapitel 4 Testsimulationen fur die SLM bzw. MPSLM durchgefuhrt und dieErgebnisse diskutiert. In Kapitel 5 schließlich werden die Ergebnisse zusam-mengefasst und ein Ausblick auf mogliche Verbesserungen des neuen Verfah-rens gegeben. Auf der beiliegenden CD-ROM finden sich die Quellcodes derSemi-Lagrange-Methode und der Massepakete-Semi-Lagrange-Methode furdie durchgefuhrten Testsimulationen. Außerdem sind die entsprechenden Er-gebnisdateien enthalten, die unter anderem den Fehler in der L2-Norm unddie Eigenschaft der Massenerhaltung protokollieren. Zusatzlich finden sichdarauf Animationen einiger Tests.

2 NUMERISCHE GRUNDLAGEN 3

2 Numerische Grundlagen

Dieses Kapitel bietet einen kurzen Uberblick uber die in dieser Arbeit ver-wendeten numerischen Methoden und Konzepte. Anfangs wird in Abschnitt2.1 die Diskretisierung durch die Finite-Elemente-Methode (FEM) exempla-risch anhand einer elliptischen Differentialgleichung dargestellt. In 2.2 wirdein kurzer Uberblick zur adaptiven Steuerung des in der FEM verwendetenGitters gegeben. Anschließend wird in 2.3 auf die fur die Simulation vonStromungen in der Atmosphare verwendete, zeitabhangige Advektionsglei-chung eingegangen. Die Diskretisierung der Advektionsgleichung mit Hilfeder Semi-Lagrange-Methode wird in Abschnitt 2.4 behandelt. Die Forderungnach Massenerhaltung wird in 2.5 naher erortert.

2.1 Raumdiskretisierung mit der Finite-Elemente-Methode

Im Folgenden seien nur die wesentlichen Schritte erlautert, um die Finite-Elemente-Methode einzufuhren: Von der klassischen kontinuierlichen Dar-stellung einer Differentialgleichung geht man uber zur Variationsformulie-rung, die daraufhin nach Galerkin diskretisiert wird. Aus Effizienzgrundenwerden bei der Diskretisierung Finite Elemente verwendet. Fur genauereDarstellungen sei hier auf die Literatur verwiesen (z.B. [QV94], [GR94],[KA00], [Bra97]).Als Beispiel diene uns die Poisson-Gleichung, eine elliptische Differential-gleichung zweiter Ordnung, mit homogenen Randbedingungen (Dirichletbe-dingung):

−∆u = f auf Ω (2.1)u = 0 auf dem Rand Γ = ∂Ω (2.2)

Hierbei stellt ∆ den Differentialoperator zweiter Ordnung ∆ = ( ∂2

∂x21, . . . , ∂

2

∂x2d)T

mit x = (x1, . . . , xd) ∈ Ω ⊂ Rd dar. Ω sei ein offenes, beschranktes Gebietund so beschaffen, dass der Gauß’sche Divergenzsatz anwendbar ist.Man betrachtet nun das Skalarprodukt von f mit beliebiger Testfunktionv ∈ C∞0 (Ω):

〈f, v〉 :=∫Ω

f(x)v(x) dx(2.1)= −

∫Ω

∇ · (∇u)(x)v(x) dx (2.3)

=∫Ω

∇u(x) · ∇v(x) dx−∫∂Ω

∇u(x) · ν(x)v(x) dx (2.4)

=∫Ω

∇u(x) · ∇v(x) dx (2.5)

Der Ubergang von (2.3) zu (2.4) ergibt sich aus dem Divergenzsatz. Wegenv|∂Ω = 0 verschwindet das Randintegral in (2.4).


Wir definieren nun die Bilinearform a : V × V → R

a(u, v) := 〈f, v〉 (2.6)

mit u, v ∈ V . Dabei ist

V := u : Ω→ R | u ∈ C(Ω), ∂iu existiert und ist stuckweisestetig fur alle i = 1, . . . , d, u = 0 auf ∂Ω

Definition 2.1. a(u, v) = 〈f, v〉 ∀v ∈ V heißt Variationsgleichung.

Definition 2.2. u ∈ V heißt schwache Losung der Differentialgleichung(2.1), wenn gilt:

a(u, v) = 〈f, v〉0 ∀v ∈ V

Bemerkung 2.3. Fur V wird im Fall der Poisson-Gleichung der Sobolev-Raum H1(Ω) gewahlt:

H1(Ω) := u : Ω→ R | u ∈ L2(Ω), u hat schwache Ableitungen∂iu ∈ L2(Ω) fur alle i = 1, . . . , d

Falls zusatzlich u|∂Ω = 0 gilt, ergibt sich:

H10 (Ω) := u ∈ H1(Ω) | u = 0 auf ∂Ω

Es sei erwahnt, dass auf H1(Ω) ein Skalarprodukt definiert ist mit

〈u, v〉1 :=∫Ω

u(x)v(x) dx+∫Ω

∇u(x) · ∇v(x) dx

sowie der davon erzeugten Norm

‖u‖1 :=√< u, u >1 =

∫Ω

|u(x)|2 dx+∫Ω

|∇u(x)|2 dx1/2

2

Sei nun a wie in (2.6) definiert und b : V → R eine Linearform, so kannman folgende zu (2.1) aquivalente Variationsgleichung aufstellen und derenLosung suchen:

Gesucht u ∈ V mit a(u, v) = b(v) ∀v ∈ V

Nach (2.5) wird a(u, v) =∫

Ω∇u · ∇v dx und b(v) =∫

Ω fv dx gewahlt. Diesstellt ein unendlich-dimensionales Problem dar, da dimV = ∞ gilt. DieLosung u der Differentialgleichung ist in den meisten Fallen jedoch nicht


in geschlossener Form darstellbar. Es gilt eine numerische Methode zu ent-wickeln, um eine diskretisierte Losung uh zu finden, fur die der Diskretisie-rungsfehler gegen 0 strebt:

limh→0‖u− uh‖V = 0 (2.7)

mit einem Funktionenraum V mit Norm ‖·‖V und mit dimV = ∞. Gilt(2.7), spricht man von Konvergenz der Methode. Dabei ist zu beachten, dassdie diskretisierte Losung uh in einem endlich-dimensionalen Unterraum Vhliegt, der insbesondere selbst ein Funktionenraum mit der Norm ‖·‖V ist.Fur Vh fordern wir in diesem Zusammenhang:

uh ∈ Vh ⊂ V (2.8)

Gilt (2.8), spricht man von einer konformen Methode. Eine weitere gefor-derte Eigenschaft einer Methode ist die Stabilitat. Das bedeutet, dass dieLosung uh stetig von den Anfangswerten abhangt und kleine Storungen derAnfangswerte nur kleine, beschrankte Storungen in der Losung uh zur Folgehaben:

‖unh‖ ≤M∥∥u0

h

∥∥ mit festem M > 0

Um eine numerische Losung zu gewinnen, wird der Diskretisierungsansatznach Galerkin verfolgt.Man ersetze V durch einen Teilraum Vh ⊂ V mit dimVh = N <∞ und losedie Aufgabe:

Gesucht uh ∈ Vh mit a(uh, v) = b(v) ∀v ∈ Vh (2.9)

Fur den endlich-dimensionalen Raum Vh, den sog. Ansatzraum, lasst sich ei-ne Basis aus endlich vielen linear unabhangigen Funktionen ϕi ∈ Vh wahlen,die Vh aufspannen:

Vh := linϕ1, . . . , ϕN (2.10)

Jede Funktion w ∈ Vh laßt sich also durch Linearkombination der Basis vonVh darstellen:

w =N∑i=1

ciϕi

mit geeigneten Koeffizienten ci ∈ R. Insbesondere gilt also: uh =∑N

i=1 ξiϕi.Man kann nun (2.9) wie folgt formulieren:

Gesucht uh =N∑i=1

ξiϕi ∈ Vh, mit a(uh, v) = b(v) ∀v ∈ Vh (2.11)

Da v ∈ Vh ebenfalls durch die gegebene Basis darstellbar ist, kann man(2.11) auch durch

a(uh, ϕi) = b(ϕi) ∀i = 1, . . . , N


und schließlichN∑i=1

a(ϕj , ϕi)ξj = b(ϕi) ∀j = 1, . . . , N

ausdrucken. Fur das Modellproblem der Poisson-Gleichung ist

a(ϕj , ϕi) =∫Ω

∇ϕj∇ϕi dx ,

b(ϕi) =∫Ω

fϕi dx .

Daraus ergibt sich mit

Ah = (a(ϕj , ϕi))ij ∈ RN×N

ξ = (ξ1, . . . , ξN )T ∈ RN

qh = (b(ϕi))i ∈ RN

das zum Galerkin-Verfahren zugehorige Gleichungssystem

Ahξ = qh . (2.12)

Wir fuhren nun beispielhaft den Teilraum Vh ein, der in dieser Arbeit zurAnwendung kommt und die geforderten Anspruche an Konvergenz und Sta-bilitat erfullt.Sei Ω ein polygonal berandetes Gebiet und der Rand Γ aus endlich vielenGeradenstucken zusammengesetzt. Dann wahlt man eine Triangulierung Thwie folgt:

Definition 2.4. Eine Triangulierung Th einer Menge Ω ⊂ Rd besteht ausendlich vielen Teilmengen T von Ω, mit folgenden Eigenschaften:(T1) Jedes T ∈ Th ist abgeschlossen.(T2) Fur jedes T ∈ Th ist sein nichtleeres Inneres int(T ) ein Lipschitz-Gebiet.(T3) Ω = ∪T∈ThT .(T4) Fur verschiedene T1 und T2 aus Th ist der Schnitt von int(T1) undint(T2) leer.

In der vorliegenden Arbeit wird mit Tetraedern als Elemente T1, . . . , Tn ge-arbeitet. Die in einem Punkt zusammenfallenden Eckpunkte mehrerer Tetra-eder werden als Knoten K1, . . . ,Km mit Koordinaten k1, . . . , km bezeichnet.Als Approximation der Randwertaufgabe (2.1), (2.2) mit linearen finitenElementen auf einer gegebenen Triangulierung Th von Ω versteht man dieWahl:

Vh := u ∈ C(Ω)| u|T ∈ P1(T ) ∀ T ∈ Th, u = 0 auf ∂Ω


Dabei stellt P1(T ) die Menge der Polynome 1. Grades auf T dar. Wahltman nun die nodale Basis

ϕi(Kj) = δij ∀j = 1, . . . ,m mit

δij = 1, wenn i = j,

δij = 0, sonst,(2.13)

kann man uh ∈ Vh in eindeutiger Weise durch die Knotenwerte u(Ki),i = 1, . . . ,m bestimmen. Es gilt also:

uh(x) =m∑i=1

u(Ki)ϕi(x) x ∈ Ω

Die Funktionswerte u(Ki) an den Knoten Ki werden Freiheitsgrade vonu ∈ Vh genannt.Um eine Basis von Vh = ϕ1, . . . , ϕN fur lineare finite Elemente anzugeben,fuhren wir baryzentrische Koordinaten ein.Sei T ein Tetraeder mit Eckpunkten K1, K2, K3, K4, dann laßt sich T wiefolgt beschreiben:

T :=

x =

4∑i=1

λiki | 0 ≤ λi ≤ 1,4∑i=1

λi = 1

=

x = k1 +

4∑i=2

λi(ki − k1) | λi ≥ 0,4∑i=2

λi ≤ 1

Die Parameter λi werden baryzentrische Koordinaten auf T genannt. Zuihrer Berechnung stellt man folgendes Gleichungssystem auf:

Bλ =(x1

)⇔

4∑j=1

kijλj = xi

4∑j=1

λj = 1

mit

B =

k11 k12 k13 k14

k21 k22 k23 k24

k31 k32 k33 k34

1 1 1 1

Man erhalt ein λi, indem man die i-te Spalte von B durch ( x1 ) ersetzt undmit der Cramer’schen Regel folgendes berechnet:

λi(x) =1

det(B)det

k11 . . . x1 . . . k14

k21 . . . x2 . . . k24

k31 . . . x3 . . . k34

1 . . . 1 . . . 1

(2.14)

Die lineare Finite-Elemente-Basisfunktion kann also folgendermaßen defi-niert werden:


Definition 2.5 (Lineare FEM-Basisfunktion auf T ∈ Th). Die Basis-funktion ϕi(x) zum Knoten Ki ist auf dem Element T mit den EckpunktenKi,Kj ,Kk,Kl durch

ϕi(x) = λi(x)

gegeben. λi(x) wird mit Hilfe von (2.14) berechnet, wobei die Koordinatender Eckpunkte Ki,Kj ,Kk,Kl eingesetzt werden.

Nach Definition der baryzentrischen Koordinaten wird (2.13) erfullt. Wegender Form, die die lineare Finite-Elemente-Basisfunktion um einen Knotender Triangulierung erzeugt, wird sie auch Pyramidenfunktion genannt.

Die Bedingung der nodalen Basis aus (2.13) wird allgemein an alle Finite-Elemente-Ansatzfunktionen gestellt. Damit gilt fur die Matrixeintrage vonAh aus (2.12) zu einer beliebigen Finite-Elemente-Ansatzfunktion:

(aij) =∫Ω

∇ϕi · ∇ϕj 6= 0 ⇔ Ki,Kj ∈ T fur ein T ∈ Th

Dies ist genau dann der Fall, wenn Ki und Kj Nachbarn sind. Da dies furdie wenigsten Eintrage von Ah zutrifft, ist Ah uberwiegend durch Nullele-mente besetzt und heißt deshalb auch dunn besetzt. Die Eigenschaft derFinite-Elemente-Methode, dass die erzeugten Matrizen dunn besetzt sind,ist ein Grund dafur, dass die FEM in der Anwendung sehr beliebt ist. Mit ge-eigneten Algorithmen lassen sich Gleichungssysteme, die zu dunn besetztenMatrizen fuhren, wesentlich effizienter losen als voll besetzte Matrizen.

Betrachtungen zur Konvergenz und Stabilitat der FEMIm folgenden wird nur knapp auf Ergebnisse der Analyse der FEM eingegan-gen. Tiefergehendere Darstellungen finden sich z.B. in [KA00] oder [QV94].Die Bilinearform a heißt stetig bzgl. ‖·‖, falls

|a(u, v)| ≤M ‖u‖ ‖v‖

fur alle u, v ∈ V mit einem M > 0 und elliptisch auf V , wenn

|a(u, u)| ≥ α ‖u‖2

fur u ∈ V mit α > 0. Es gilt:

Satz 2.6 (Stabilitat des Galerkin-Verfahren). Die Losung uh nach(2.9) ist stabil im folgenden Sinn

‖uh‖ ≤1α‖b‖ unabhangig von h,

wobei‖b‖ := sup

v∈V,v 6=0

|b(v)|‖v‖

.


Den Beweis findet man z.B. in [KA00].

Definition 2.7. Sei u ∈ V die Losung des betrachteten kontinuierlichenProblems und sei uh ∈ Vh die durch die Galerkin-Diskretisierung gewonneneLosung, dann bezeichnet

e = ‖u− uh‖Vden Diskretisierungsfehler des Verfahrens.

Satz 2.8 (Lemma von Cea). Unter der Voraussetzung, dass die Biline-arform a bzgl. einer Norm ‖·‖ auf V stetig und elliptisch ist, gilt fur denDiskretisierungsfehler folgende Abschatzung:

‖u− uh‖V ≤ C infv∈Vh‖u− v‖V , mit C > 0, fest.

Den Beweis hierzu findet man z.B. in [KA00].

Unter der Voraussetzung, dass die kontinuierliche Losung u existiert, be-trachtet man nun eine Folge von Teilraumen Vl ⊂ V mit dimVl = Nl <∞,l ∈ N und Vl ⊂ Vm fur l < m ∈ N so dass

liml→∞

Vl → V

Dann gilt nach [Beh96]:

Satz 2.9 (Approximation durch diskrete Unterraume).Ist uh ∈ Vl ⊂ V eine Losung im Sinne der Galerkin-Methode, so konvergiertmit l→∞ uh gegen u:

liml→∞‖u− uh‖ → 0

Das bedeutet also: Erhoht man die Zahl der Freiheitsgrade, verringert sichder Diskretisierungsfehler. Damit stellt sich auch die Frage: Gibt es eine Stra-tegie, gezielt die Freiheitsgrade zu erhohen, um den Diskretisierungsfehlerzu verringern? Damit waren wir bei Fragen der adaptiven Gittersteuerung,die im folgenden Abschnitt naher erlautert werden.

2.2 Adaptivitat

Ist eine Triangulierung Th von Ω im gesamten Gebiet gleichmaßig fein, sospricht man von einem uniformen Gitter. Es stellt sich heraus, dass mandurch gezielte Erhohung der Freiheitsgrade in Teilgebieten von Ω den Dis-kretisierungsfehler e = ‖u− uh‖V erheblich verringern kann, ohne auf demgesamten Gebiet mit derselben Verfeinerungsstufe des Gitters rechnen zumussen. Verwendet man im Bezug auf die Verfeinerungsstufe heterogeneGitter, so spricht man von einem adaptiven Gitter. Adaptive Gitter eroffnendie Moglichkeit, die zur Verfugung stehenden Rechnerkapazitaten wesentlich


effizienter zu nutzen.Einerseits kann man bei gleicher Fehlertoleranz Rechenzeiten und Spei-cherbedarf stark senken, indem nur noch in Teilgebieten mit voller Feinheitgerechnet und das Gitter in anderen Gebieten ausgedunnt wird. Ist manalso mit der bisher auf dem uniformen Gitter erreichten Genauigkeit zufrie-den, erreicht man eine Beschleunigung des Verfahrens. Ist man andererseitsmit der Rechenzeit und dem Speicherbedarf des Verfahrens zufrieden, kannman mit niedrigerer Knotendichte in den ”uninteressanten“ bzw. mit hoher-er Knotendichte in den ”interessanten“ Gebieten rechnen. So erreicht manmit adaptiven Gittern eine hohere Genauigkeit bei gleichem Ressourcenbe-darf im Vergleich zu uniformen Gittern. Welches sind also die interessantenGebiete?

Abbildung 2.1: Friedrichs-Keller-Triangulierung des Einheitswurfels

Das Gebiet Ω wird durch finite Elemente gemaß (2.4) trianguliert. Als Bei-spiel fur ein Gebiet Ω dient uns der Einheitswurfel Ω = E = [0, 1]3. Hier ver-wenden wir zur Anfangstriangulierung die Friedrichs-Keller-Triangulierungdes Einheitswurfel (siehe Abbildung 2.1). Die Triangulierung wird anhandder gegebenen Anfangsbedingungen verfeinert und im weiteren Verlauf dernumerischen Losung einer Differentialgleichung in jedem Zeitschritt den Ge-gebenheiten angepasst. Es wird also eine adaptive Gittersteuerung einge-setzt. Damit soll erreicht werden, dass nur an Knoten K Werte uh berechnetwerden, an denen es sich lohnt.Es wird eine Fehlerschranke TOL vorgegeben. Ergibt der Fehlerschatzer,dass der Fehler lokal auf dem Tetraeder T unterhalb der Fehlerschranke liegt,wird das Gitter nicht verfeinert. Wird die Schranke uberschritten, wird lokal


verfeinert. Wird sie stark unterschritten, wird das Gitter lokal vergrobert.

Algorithmus 2.10 (Adaptive Steuerung des Gitters).

‖u− uh‖V |T

TOL Vergrobere> TOL Verfeinere

Letztendlich entscheidend ist die Gute des Fehlerschatzers auf T . Er orien-tiert sich im verwendeten Programmpaket amatos3D/flash3D von Behrensam Gradienten am Knoten Ki. Als Verfeinerungs- bzw. Vergroberungsalgo-rithmus kommt der in [Ban91] beschriebene Algorithmus durch Bisektionder Elemente T ∈ Th zur Anwendung.So entstehen Gebiete mit unterschiedlich großer Gitterweite, also lokal hoher-er Auflosung (siehe Abbildung 2.2), mit der man erreicht, dass der Diskre-tisierungsfehler lokal geringer ausfallt. Die adaptive Gittersteuerung stellteinen eigenen Forschungszweig innerhalb der Numerik partieller Differen-tialgleichungen dar. Zur Vertiefung dieser Thematik sei der Leser auf dieLiteratur verwiesen (z.B. [FG00]).

Abbildung 2.2: Beispiel fur eine verfeinerte Anfangstriangulierung des Ein-heitswurfels mit Schnitt bei z=0,5. Gut erkennbar ist die in einem Bereich auf-tretende hohere Knotendichte.

2.3 Die Advektionsgleichung

Die passive Advektionsgleichung ist eine hyperbolische Differentialgleichungerster Ordnung. Sie ist gut geeignet, um einfache Stromungen zu modellie-ren. Ein weiterer Vorteil ist, dass interessante Testfalle auch analytisch losbarsind. Das bietet die Moglichkeit, die numerisch berechneten Ergebnisse mit


der analytischen Losung zu vergleichen und somit den Diskretisierungsfehlerdirekt zu bestimmen. Das nutzen wir im Kapitel 4 aus, um die Diskretisie-rungsfehler der Semi-Lagrange-Methode mit der in Kapitel 3 entwickeltenMassepakete-Semi-Lagrange-Methode zu vergleichen.Die passive Advektionsgleichung in Euler-Koordinaten mit divergenzfreierStromung ∇ · a = 0 lautet:

∂u

∂t+ a · ∇u = 0 (2.15)

Dabei stellt u : Ω× [0, T ]→ R die gesuchte Losung, a : Ω× [0, T ]→ Rd ein

orts- und zeitabhangiges Windfeld und ∇ = ( ∂∂x1

, . . . , ∂∂xd

)T den Differen-tialoperator erster Ordnung bzgl. der Ortskoordinaten dar.Wir werden nun einen Koordinatenwechsel durchfuhren und leiten die Ad-vektionsgleichung in Lagrange-Koordinaten her (s. [QV94], [KA00]). DieIdee dahinter ist, dass die Konzentration entlang der Charakteristik kon-stant ist. Den Koordinatenwechsel kann man so verstehen, dass in Euler-Koordinaten der Standpunkt des Beobachters fest gewahlt ist und die Par-tikel der Stromung vorbeiziehen. Die Lagrange-Koordinaten werden dage-gen so gewahlt, dass sich die Position des Beobachters mit der Stromungmitbewegt. Um die Koordinatentransformation von Euler- zu Lagrange-Koordinaten durchzufuhren, wird folgendes Hilfsproblem betrachtet:

Finde zu (x, s) ∈ Ω× [0, T ] ein Vektorfeld X : Ω× [0, T ]2 → Rd mit

dXdt (x, s, t) = a(X(x, s, t), t), t ∈ [0, T ],X(x, s, s) = x.

(2.16)

Die Trajektorien X(x, s, ·) sind die eben genannten Charakteristiken. Ist astetig auf Ω×[0, T ] und Lipschitz-stetig auf Ω fur festes t ∈ [0, T ], so existierteine eindeutige Losung X = X(x, s, t).Ist nun u eine Losung der Advektionsgleichung (2.15), so folgt mit

u(x, t) := u(X(x, s, t), t) bei festem s ∈ [0, T ] (2.17)

und der Anwendung der Kettenregel die Beziehung

du

dt(x, t) =

(∂u

∂t+ a · ∇u

)(X(x, s, t), t)

Die Differentialgleichung fur passive Advektion lautet also in Lagrange-Koordinaten:

du

dt(x, t) = 0

und damit schließlich fur divergenzfreie Stromung

du

dt(x, t) = 0 (2.18)


Dem Zusammenhang (2.17) kommt zentrale Bedeutung zu, denn insbeson-dere gilt, dass der Wert von u durch den Wert u an den stromaufwartsgelegenen Punkten entlang der Charakteristik X(x, s, t) festgelegt ist. EinVerfahren, dass (2.18) numerisch lost, muss nun erstens die gewohnlicheDifferentialgleichung (2.16) approximieren und zweitens per Interpolationden Wert u an der stromaufwarts durch (2.16) gegebenen Stelle bestimmen.Diese Aufgabe leistet die Semi-Lagrange-Methode (SLM).

2.4 Diskretisierung mit der Semi-Lagrange-Methode

Die Diskretisierung mit Hilfe der Semi-Lagrange-Methode hat sich im Be-reich der Simulation von Ozean- und Atmospharenmodellen etabliert. Derentscheidende Vorteil des Verfahrens ist, dass man nicht an die Courant-Friedrichs-Lewy-Bedingung gebunden ist, die durch die Modellgleichung inEuler-Koordinaten impliziert ist. Diese Stabilitatsbedingung laßt nur kleineZeitschrittweiten zu.

∆

t

−

x

t

xm xm− 2α mm m

t

− α

∆

n

A’

B

C

A

+t

t n

n

Abbildung 2.3: Schematische Darstellung der Drei-Schritt-SLM fur eindimensio-nale Advektion

Um die Diskretisierung mit der SLM zu veranschaulichen, zeigen wir an-hand von Abbildung 2.3 die Vorgehensweise mittels ruckwarts genommenerDifferenzenquotienten im eindimensionalen Fall (siehe auch [Beh96]).Die Werte von u seien zu den Zeitpunkten tn und tn−∆t an allen Gitterpunk-ten xi bekannt. Exemplarisch werde nun der Wert von u bei (xm, tn + ∆t),in Abbildung 2.3 mit C bezeichnet, gesucht. Dazu wird die exakte Trajek-torie (AC) (durchgehende Linie) verfolgt. Der Wert bei A ware der exakteWert fur u(xm, tn + ∆t). Die Trajektorie (AC) ist jedoch im allgemeinen


unbekannt, sie wird durch die lineare Trajektorie (A′C) (gestrichelt) appro-ximiert. Das Schema, um die Advektionsgleichung (2.18) zu diskretisieren,lautet damit:

u(xm, tn + ∆t)− u(xm − 2αm, tn −∆t)2∆t

= 0 (2.19)

Dabei wird zur Vereinfachung der Notation αm = xm −X(xm, tn + ∆t, tn)verwendet.Die Position des stromaufwarts gelegenen Punktes A′ wird durch numeri-sche Losung der gewohnlichen Differentialgleichung (2.16) berechnet. DieGenauigkeit der Approximation der Trajektorie ist demnach von zentra-ler Bedeutung. Bezogen auf (2.16) wird zur Diskretisierung von dX ∆x =xm − (xm − 2αm) = 2αm und von dt tn + ∆t − tn − ∆t = 2∆t gewahlt.Daraus ergibt sich eine Naherung mit O(∆t2):

αm = a(xm − αm, tn)∆t (2.20)

Die implizite Formel wird iterativ mit einem geeigneten Startwert α(0)m gelost.

Verfahren hoherer Genauigkeit findet man in [Beh96], [KA00] und [QV94].Im Anschluss an die nun folgende allgemeine Formulierung derSemi-Lagrange-Methode als Algorithmus wird auf die in dieser Arbeit ge-wahlte Form der Interpolation eingegangen.

Das Drei-Schritt-Semi-Lagrange-Verfahren fuhrt man wie folgt durch:

Algorithmus 2.11. (Drei-Schritt-Semi-Lagrange-Verfahren)

1. Berechne fur alle xi zum Zeitpunkt tn + ∆t die Trajektorien-Approxi-mation αi nach (2.20)

2. Berechne u zur Zeit tn − ∆t mittels geeigneter Interpolation an denstromaufwarts gelegenen Punkten xi − 2αi fur alle Knoten Ki ∈ Thdes Zeitschritts tn + ∆t

3. Berechne die Werte u an den Gitterpunkten xi zur Zeit tn + ∆t mitHilfe von (2.19)

Es fallt auf, dass in das Drei-Schritt-Verfahren nur Werte u in den Zeit-schichten tn − ∆t bzw. tn + ∆t verwendet werden. Lediglich das Windfelda zum Zeitpunkt tn geht in die Berechnung der Trajektorien ein. Die Aus-wertung von u zum Zeitpunkt tn kann deshalb weggelassen werden. Manerreicht so eine Verdopplung der Zeitschrittlange.Es liegt also nahe, durch Notationsanderung von tn − ∆t → tn und tn →tn + ∆t

2 das folgende Zwei-Schritt-Semi-Lagrange-Verfahren zu formulieren:


Algorithmus 2.12. (Zwei-Schritt-Semi-Lagrange-Verfahren)

1. Berechne fur alle xi zum Zeitpunkt tn + ∆t2 mit Hilfe der Auswertung

von a die Trajektorien-Approximation αm nach (2.20)

2. Berechne u zur Zeit tn mittels geeigneter Interpolation an den strom-aufwarts gelegenen Punkten xi − αi fur alle Knoten Ki ∈ Th des Zeit-schritts tn + ∆t

3. Berechne u an den Gitterpunkten xi zur Zeit tn + ∆t mit Hilfe von(2.19)

Fur die Approximation der Trajektorien αm wurde mit (2.20) bereits eineMoglichkeit angegeben. Sei noch die Form der Interpolation auf dem Finite-Elemente-Gitter mit linearen finiten Elementen nachgetragen, die u.a. imweiteren Verlauf dieser Arbeit benutzten Programmpaket amatos3d/flash3dverwendet wird.Fur die Interpolation der Werte u an den Knoten Ki der TriangulierungT +h zum Zeitpunkt tn + ∆t seien die Trajektorien αi gegeben. Außerdem sei

das Element T ∈ T −h bekannt, in dem der jeweils stromaufwarts gelegenePunkt zu Ki und αi liegt. Die mit − bzw. + versehene Triangulierung Thist das adaptiv verfeinerte Gitter zum Zeitschritt tn respektive tn + ∆t.Dann wird bei der benutzten Konfiguration von amatos3d/flash3d folgendeInterpolation verwendet:

Algorithmus 2.13. (Lineare Interpolation mit Clipping)

1. Berechne die lineare Interpolation ulin(xi−αi) =∑

i,Ki∈Tu(Ki)λi(x) von

u an der Position xi−αi in T ∈ T −h . Dabei sind λi die baryzentrischenKoordinaten von xi − αi in T .

2. Setze

u(x) =

umax , falls umax < ulin(x)umin , falls umin > ulin(x)ulin(x) , sonst

Dabei ist umax = maxKi∈T

u(Ki) und umin analog das Minimum uber alle

Eckpunkte des Tetraeders T.

Weitere Moglichkeiten der Interpolation finden sich in [Beh96].

Bemerkung 2.14. Die Suche, in welchem Element T des Finite-Elemente-Gitters T −h der stromaufwarts gelegene Punkt zu einem xi liegt, ist nichttrivial und daruberhinaus ein gewichtiger Part bei der benotigten Rechenzeitzur Losung der diskretisierten Advektionsgleichung.


Betrachtungen zur Konvergenz und Stabilitat der SLMEs seien hier nur die zentralen Satze aus dem Paper von Falcone und Fer-retti [FF98] genannt. Sie behandeln die Konvergenz und Stabilitat der Semi-Lagrange-Methode angewandt auf lineare hyperbolische Differentialgleichun-gen erster Ordnung. Angepasst auf unser Problem der Advektionsgleichung∂u∂t + a · ∇u = 0 lautet das SLM-Schema in der Notation von Falcone undFerretti:u

nj =

∑m∈Q

un−1m ϕm (xj + hΦa(xj)) fur alle j ∈ Q

u0j = u0(xj)

(2.21)

Hierbei ist unj der Wert am Knoten xj zum Zeitschritt tn undQ = 1, . . . ,Mdie Menge der Knotenindizes. ϕ1, . . . , ϕM sei die nodale Basis auf denKnoten xj . a ∈ Cp(Rd) sei ein gegebenes Windfeld und Φa ein konsistentesVerfahren p-ter Ordnung zur Approximation der gewohnlichen Differential-gleichung, die die Trajektorien bzgl. a beschreibt. h ist die Zeitschrittweiteund k die Schrittweite bzgl. des Ortes.Sei nun Pk ein Projektionsoperator auf der nodalen Basis ϕ1, . . . , ϕM furden mit der Definition

Pkw(x, t) :=∑m∈Q

w(xm, t)ϕ(x) (2.22)

geltelimk→0‖w − Pkw‖∞ = 0 fur ein w ∈W 1,∞(Ω)

Dabei ist W 1,∞(Ω) die Menge der in Ω einmal schwach differenzierbarenFunktionen aus L∞(Ω) mit Ableitungen in L∞(Ω).Φa stellt also die Approximation bzgl. t und Pk die Approximation bzgl. desOrtes dar. Sei nun noch u0 ∈ W 1,∞(Ω) und ‖u(t)− Pku(t)‖∞ ≤ E(k) furein beliebiges t ∈ [0, T ]. Dann gilt nach [FF98] fur ein festes C > 0:

Satz 2.15 (Konvergenz der Semi-Lagrange-Methode).Wenn

∑j∈Q |ϕj(x)| ≤ 1, dann gilt fur eine beliebige Zeitschrittweite h < h0∥∥∥ukh(T )− u(T )

∥∥∥∞≤ C

(hp +

1hE(k)

).

Falcone und Ferretti zeigen unter weiteren Bedingungen, dass dieselbe Ab-schatzung auch fur die L2-Norm gilt. Speziell fur lineare finite Elemente undder Pyramidenfunktion als Basisfunktion gilt, dass

∑j∈Q |ϕj(x)| = 1 und

damit der obige Satz anwendbar ist.

Die Autoren kommen zu folgendem Ergebnis zur L∞-Stabilitat im Bezugauf die Zeitschrittweite h des Semi-Lagrange-Verfahrens:


Satz 2.16 (L∞-Stabilitat der Semi-Lagrange-Methode). Sei die Orts-schrittweite k fest und es gelte (2.22) fur den Projektionsoperator Pk auf dernodalen Basis, dann ist das Schema (2.21) fur ein beliebiges h > 0 stabil.

Beweise sowohl fur die Stabilitats- als auch die Konvergenzaussage findensich in [FF98].

Bemerkung 2.17. Die Abschatzung aus Satz 2.15 zeigt, dass es nur Sinnmacht, eine geeignete Kombination der Genauigkeit der Interpolation aufden Knotenwerten Pk und der Genauigkeit des Verfahrens Φ zur Approxi-mation der Trajektorien zu wahlen.

Unter den gestellten Bedingungen lasst die SLM also prinzipiell beliebigweite Zeitschritte zu, ließe man den Diskretisierungsfehler außer Acht. EinNachteil des Verfahrens ist bisher jedoch, dass die Forderung nach Massen-erhaltung im allgemeinen nicht erfullt wird, wie wir im nachsten Abschnittsehen.

2.5 Massenerhaltung

Die mit u bezeichnete Funktion stellt in den spateren Abschnitten dieser Ar-beit eine Massendichtegroße (Konzentration) dar. Das mit a(x, t) bezeich-nete Windfeld gibt im Gebiet Ω eine Stromung vor. Fur Stromungen inGasen und Flussigkeiten gilt ein Grundprinzip der Physik, die so genannteKontinuitatsgleichung :

∂u

∂t+ div(ua) = 0

Herleitung und genauere Darstellung finden sich z.B. in [LL91]. Die Konti-nuitatsgleichung besagt, dass sich eine Anderung der Große u uber die Zeitt nur durch ein- oder ausstromendes Gas ergeben kann. Integriert man dieDichte u uber das Gebiet Ω, erhalt man die im Gebiet vorhandene Masse

m =∫Ω

u dx

Es folgt sofort das Prinzip der Massenerhaltung

dm

dt= 0

falls kein Fluss in das bzw. aus dem Gebiet Ω stattfindet. Erwunscht ist nunauch, dass das Prinzip der Massenerhaltung bei der Semi-Lagrange-Finite-Elemente-Diskretisierung der Advektionsgleichung beachtet wird. Laßt manfur einen Moment Stromungen in das bzw. aus dem Gebiet Ω außer acht, so


ergibt sich fur triangulierte Gebiete Ω ⊂ R3 mit linearen finiten Elementenfolgende Forderung:∑

i

uk(Kki )

14V (P ki ) =

∑i

u0(K0i )

14V (P 0

i ) fur k = 1, . . . , N (2.23)

Dabei ist uk(Kki ) die Konzentration am Knoten Kk

i zum Zeitschritt k undV (P ki ) das Volumen des Finite-Elemente-Patches um Kk

i :

Pi =⋃

Ki∈Tj

Tj

Da die Lage der Ki und Pi durch die adaptive Veranderung des Gittersnicht fest sind, wird auch hier mit k indiziert. k = 0 stellt die Startwertedar. Da die Werte uk(Kk

i ) durch Interpolation bestimmt werden, sind siemit einem Fehler behaftet, der sich aus der Interpolationsordnung ergibt.Damit folgt aus (2.23) sofort, dass die Semi-Lagrange-Methode das Prin-zip der Massenerhaltung verletzt. Mit Hilfe der Massepakete-SLM wird inKapitel 3 eine Semi-Lagrange-Methode entwickelt, die durch Konstruktionexakt massenerhaltend gemaß (2.23) bis auf Massenverlust uber den Randdes Simulationsgebiets Ω ist.

3 MASSEPAKETE-SEMI-LAGRANGE-METHODE 19

3 Massepakete-Semi-Lagrange-Methode

Wie in Abschnitt 2.5 beschrieben ist die Semi-Lagrange-Methode mit dervorliegenden Interpolation nicht massenerhaltend. In diesem Kapitel wirdein Verfahren eingefuhrt, dass massenerhaltend ist. Dabei wird die Idee bei-behalten, die Trajektorien zuruck zu verfolgen. Jedoch wird sie insofern mo-difiziert, dass nicht einzelne Knoten, sondern ganze Elemente T des GittersT n+1h zum Zeitschritt tn+1 stromaufwarts verfolgt werden. Die Interpolation

an den stromaufwarts gelegenen Punkten hingegen wird durch einen neuenAlgorithmus ersetzt. Entscheidender Gedanke hierbei ist, nicht mehr kon-zentrationsbasiert zu rechnen, sondern die Konzentrationsverteilung in denTetraedern T des alten Gitters T nh zum Zeitschritt tn durch Masse-Volumen-Einheiten, sog. Massepakete, zu approximieren. Das Gitter T nh wird so durchdie Massepakete ersetzt. Dann wird die Position der Massepakete in denstromaufwarts advehierten Elementen des Gitters T n+1

h festgestellt und dieMasse den entsprechenden Elementen des Gitters T n+1

h zugeteilt. Mit derso bestimmten Masse und mit Hilfe des Volumens der T ∈ T n+1

h konnen dieneuen Konzentrationswerte u(Ki), Ki ∈ T n+1

h bestimmt werden.Durch die Berechnung der Konzentration an den Knoten von T n+1

h uber den

”Umweg“ der Massepakete ist sichergestellt, dass die in T nh vorhandene Mas-se exakt auf das Gitter des nachsten Zeitschritts T n+1

h ubertragen wird. Dadie Massepakete die zentrale Idee des neuen Verfahrens darstellen, nennt derAutor das Verfahren die Massepakete-Semi-Lagrange-Methode (MPSLM).Zunachst werden nun in Abschnitt 3.1 einige Definitionen und Begriffe erlau-tert. In Abschnitt 3.2 wird die Massepakete-Semi-Lagrange-Methode sche-matisch an einem Beispiel im Eindimensionalen erklart, um daraufhin denAlgorithmus und dessen Bestandteile im Detail zu erlautern. Abschließendwird in 3.3 auf Aspekte der Implementierung eingegangen. Darunter sinddie verwendeten Datenstrukturen und zwei zentrale Algorithmen der imple-mentierten MPSLM.

3.1 Definitionen und Begriffe

Die Begriffe Knoten Ki bzw. Tetraeder Tj einer Triangulierung Th wurdenbereits in Kapitel 2 eingefuhrt. In den weiteren Abschnitten dieses Kapitelswerden zum Verstandnis zusatzlich die folgenden Begriffe benotigt.Mit Ekj , k = 1, . . . , 4 bezeichnen wir die vier Eckpunkte des Tetraeders Tjeiner Triangulierung Th. Das Volumen an dem Eckpunkt Ekj ist:

V (Ekj ) =14V (Tj)

wobei V (Tj) das Volumen des Tetraeders Tj ist.Mit dem Knotenvolumen V (Ki) wird das Volumen aller Eckpunkte in Ki


bezeichnet. Es gilt:V (Ki) =

∑Ki=Ekj

V (Ekj )

Mit der Masse m(Ekj ) am Eckpunkt Ekj ist die mit dem Eckpunkt assoziierteTeilmasse des Tetraeders Tj gemeint. Es gilt:

m(Tj) =4∑

k=1

m(Ekj )

Die Knotenmasse m(Ki) am Knoten Ki wird wie folgt definiert:

m(Ki) =∑

Ki=Ekj

m(Ekj )

Es wird also streng zwischen den Eckpunkten Ekj eines Tetraeders Tj undden Knoten Ki einer Triangulierung Th unterschieden. Wie zuvor bereitsangedeutet, werden die Tetraeder Tj in Teilmengen, sog. Massepakete, un-terteilt. Mit Massepaket M l

j wird ein Teilgebiet des Tetraeders Tj mit Indexl ∈ 1, . . . , L bezeichnet. L ist hierbei die Anzahl der Teilmengen in Tj .Die Unterteilung in Massepakete wird so vorgenommen, dass das Innere derMassepakete int(M l

j), l = 1, . . . , L paarweise disjunkt ist und es gilt:

Tj =⋃l

M lj

Daraus ergibt sichV (Tj) =

∑l

V (M lj)

Ein Tetraeder wird also vollstandig in Massepakete unterteilt.Den Massepaketen M l

j wird eine Masse und ein Volumen entsprechend demGrad der gewunschten Unterteilung, der Massenverteilung in Tj und demVolumen von Tj zugeordnet. Die Zuordnungsvorschriften werden im Laufedes Abschnitts 3.2.3 eingefuhrt. Die Motivation fur die gewahlten Begriffeund deren Volumen bzw. Masse wird in den folgenden Abschnitten klar.Der Grad der Unterteilung eines Tetraeder Tj in Massepakete wird im fol-genden mit dem Begriff Verfeinerungsstufe des Tetraeders Tj bezeichnet.Um klar von dem Begriff der Verfeinerungsstufe bei der adaptiven Unter-teilung der Triangulierung Th zu trennen, wird im Zusammenhang mit deradaptiven Gitterverfeinerung der Begriff Gitterstufe verwendet.

3.2 Beschreibung der Massepakete-Semi-Lagrange-Methode

Im kommenden Teilabschnitt wird zunachst anhand der Massepakete-Semi-Lagrange-Methode (MPSLM) im Eindimensionalen das Prinzip der MPSLM


veranschaulicht. In den weiteren Teilabschnitten werden die Einzelheiten derMPSLM erlautert. Hierbei sei betont, dass die MPSLM unter der Vorausset-zung linearer finiter Elemente entwickelt wird und nicht ohne weiteres auffinite Elemente hoherer Ordnung ubertragbar ist.

3.2.1 Die Massepakete-Semi-Lagrange-Methode in 1D

Analog der Beschreibung der SLM in Abschnitt 2.4 wird nun zunachst dieVorgehensweise bei der MPSLM an einem Beispiel mit eindimensionalerOrtskoordinate gezeigt.Sei die Masse m(Ekj ) an den Eckpunkten Ekj aller Elemente Tj ∈ T nh zumZeitschritt tn bekannt. Durch T nh = (x0, x1), . . . , (xM−1, xM ) ist im Ein-dimensionalen eine Triangulierung durch Intervalle auf Ω ⊂ R gegeben. Mist die Anzahl der Intervalle. Jedes Element Tj besitzt also zwei EckpunkteE1j , E

2j . In jedem Knoten mit Koordinaten Ki = (xi, tn) mit i = 0, . . . ,M

sind also zwei Eckpunkte des linken bzw. rechten Intervalls vorhanden. Nurdie Randknoten K0,KM fallen nur mit einem Eckpunkt zusammen.Die Elemente Tj werden nun in einzelne Teilintervalle geteilt. Dies sind diesog. Massepakete. Je nach Verfeinerungsgrad ergeben sich mehr oder wenigerMassepakete M l

j ∈ Tj . Die genaue Einteilung der Elemente Tj in Massepa-kete wird in 3.2.3 erlautert. Innerhalb eines Massepaketes wird homogeneMassenverteilung angenommen. Der geometrische Mittelpunkt jedes Masse-paketes M l

j ist deshalb zugleich dessen Massenschwerpunkt. Es werden diebaryzentrischen Koordinaten λl1, λ

l2 des Massenschwerpunkts jedes Masse-

paketes M lj bezogen auf die beiden Eckpunkte Ekj , k = 1, 2 des Elements

Tj berechnet (siehe Abschnitt 2.4). Die Massen der Eckpunkte m(Ekj ) wer-den gemaß den baryzentrischen Koordinaten λl1, λ

l2 und dem Volumen des

Massepaketes V (M lj) auf die Massepakete verteilt:

m(M lj) =

2∑i=1

λlim(Eij)V (M lj)

Dieses Verfahren wird in 3.2.3 erlautert.Damit haben wir die Elemente Tj ∈ (x1, x2), . . . , (xM−1, xM ) des GittersT nh in Massepakete M l

j unterteilt, sowie die Masse und den Massenschwer-punkt jedes Massepaketes bestimmt.

In Abbildung 3.1 sind exemplarisch einige Elemente Tj im Zeitschritt tn farb-lich unterschiedlich markiert. Hiermit wird illustriert, dass nicht die Knotenwie bei der SLM, sondern die Elemente von Bedeutung sind. Pro ElementTj wird beispielhaft eine Unterteilung in vier Massepakete M l

j , l = 1, . . . , 4vorgenommen. Die Unterteilung der Elemente in Massepakete ist ebenfallsdurch farbliche Trennstriche dargestellt. Die Massenschwerpunkte der Mas-sepakete M l

j sind eingetragen.


U m−1U U

m−1D D Dm m+1

m+1m

t n

− αm+1

xm+1− α

mx

m− αm−1

xm−1

x x xm−1 m m+1

t n+1

Abbildung 3.1: Schematische Darstellung der Massepakete-SLM fur eindimensio-nale Advektion

Wir berechnen nun exemplarisch die Konzentration u(Dm) am Punkt Dm

fur den Zeitschritt tn+1 = tn + ∆t.Dazu werden die in Abbildung 3.1 mit

Dm−1 = (xm−1, tn+1), Dm = (xm, tn+1), Dm+1 = (xm+1, tn+1)

bezeichneten stromabwarts (downstream) gelegenen Knoten des Gitters T n+1h

betrachtet. Sie bilden die Intervalle (Dm−1, Dm) bzw. (Dm, Dm+1), die zuden Intervallen (Um−1, Um) bzw. (Um, Um+1) entlang der approximiertenTrajektorien (gestrichelt) stromaufwarts verfolgt werden. Dabei bezeichnenwir mit

Um−1 = (xm−1−αm−1, tn), Um = (xm−αm, tn), Um+1 = (xm+1−αm+1, tn)

die stromaufwarts (upstream) gelegenen Punkte. Die Werte αm−1, αm, αm+1

sind durch die Trajektorienapproximation aus Abschnitt 2.4 gegeben.Es wird nun festgestellt, welche Massepakete in das Intervall(Um−1, Um) fallen. Als Kriterium dient hier, ob der Massenschwerpunkt ei-nes Massepaketes in (Um−1, Um) liegt. Daraufhin werden die baryzentrischenKoordinaten der Massenschwerpunkte in (Um−1, Um) bezogen auf die Eck-punkte Um−1 und Um berechnet. Die Masse der Massepakete in (Um−1, Um)wird auf die Eckpunkte Um−1 bzw. Um entsprechend den baryzentrischenKoordinaten der Massenschwerpunkte aufgeteilt. Analog geht man fur dasIntervall (Dm, Dm+1) vor. (Eine genauere Darstellung findet sich spater in


Teilabschnitt 3.2.5.) Wir haben somit die Masse an den Eckpunkten der In-tervalle (Um−1, Um) bzw. (Um, Um+1) gegeben.Man nimmt nun an, dass die Masse, die Um von den Intervallen (Um−1, Um)bzw. (Um, Um+1) zugewiesen ist, durch Advektion entlang der Trajektoriestromabwarts in den Punkt Dm transportiert werden. Dem Knoten Dm wer-den die bereits berechneten Massen zugeordnet, die in Um liegen.Zusammen mit dem Knotenvolumen V (Dm), das sich wie in Abschnitt 3.1berechnet, erhalt man die Konzentration u(Dm) mit

u(Dm) =m(Dm)V (Dm)

Diese Berechnung fuhrt man fur alle Knoten aus T n+1h durch und erhalt

so fur alle Knoten Ki ∈ T n+1h die Werte u(Ki). Die Einzelmassen in den

Eckpunkten der Intervalle bleiben fur den nachsten Zeitschritt gespeichert.

Vor dem ersten Zeitschritt muss noch die Masse an den Eckpunkten Ekjbestimmt werden. Dies wird in Abschnitt 3.2.2 naher erlautert. Die MPSLMlaßt sich damit durch folgenden Algorithmus unabhangig von der Dimensionder Ortskoordinaten beschreiben:

Algorithmus 3.1. (Massepakete-Semi-Lagrange-Methode)

1. Initialisierung: Bestimme die Masse m(Eki ) an den Eckpunkten Ekialler Elemente Ti der Triangulierung T 0

h zum Zeitpunkt t0 = 0.

2. Verteile die Masse der Eckpunkte m(Eki ) auf die Massepakete M li der

Elemente Ti. Der Massenanteil je Massepaket errechnet sich aus denbaryzentrischen Koordinaten des Massenschwerpunkts λlk und dem Vo-lumen des Massepaketes V (M l

i ).

3. Verfolge die Elemente des Gitters Tn+1j ∈ T n+1

h zum Zeitschritt tn+1

anhand der berechneten Trajektorien stromaufwarts. Die stromaufwartsverfolgten Elemente seien mit Tnj bezeichnet.

4. Liegt der Massenschwerpunkt von M li in einem Element Tnj des strom-

aufwarts verfolgten Gitters, dann weise den Eckpunkten des ElementsTnj die Masse des Massepaketes M l

i entsprechend den baryzentrischenKoordinaten des Massenschwerpunkts µlk von M l

i im neuen ElementTnj zu. Die Masse der Eckpunkte m(Ekj ) von Tn+1

j entspricht der Mas-se an den Eckpunkten von Tnj . Speichere die Masse an den Eckpunktenvon Tn+1

j fur den nachsten Zeitschritt.

5. Fur alle Knoten Ki aus T n+1h : Berechne die Konzentration u(Ki) an

Ki mit Hilfe der Knotenmasse m(Ki) und dem Knotenvolumen V (Ki).

6. Im nachsten Zeitschritt: Gehe zu 2.


Die nachsten Abschnitte behandeln jeden Einzelschritt der MPSLM genau-er. Dabei wird auf eine Triangulierung durch Tetraeder in 3D eingegangen,jedoch zur Anschauung auch Dreiecke in 2D betrachtet.

3.2.2 Initialisierung der Masse

Zunachst sind durch die Anfangsbedingungen nur Konzentrationswerte u(Ki)an den Knoten Ki der Anfangstriangulierung T 0

h gegeben. Die Berechnun-gen erfolgen jedoch massenbasiert, um Massenerhaltung zu erreichen. Da-zu wird fur jeden Eckpunkt Ekj eines Elements Tj die zugeordnete Massem(Ekj ) berechnet. Um das zu Ekj zugehorige Volumen V (Ekj ) zu berechnen,teilt man den Tetraeder in vier Teilmengen: Durch jede Kante des Tetra-eders legt man eine Schnittebene, die gemeinsam mit dem Mittelpunkt dergegenuberliegenden Kante eindeutig bestimmt ist. Das Tetraeder wird alsodurch jede Schnittebene halbiert. Es ergeben sich entsprechend den sechsKanten sechs Schnittebenen. Die Teilmengen, die einen Eckpunkt gemein-sam haben, werden vereinigt. Diese Teilmengen stellen gerade die grobsteStufe der Unterteilung von Tj in Massepakete M l

j mit l = 1, 2, 3, 4 dar. DieseMassepakete werden Massepakete der 1. Verfeinerungsstufe genannt.Die Massepakete der 1. Verfeinerungsstufe lassen sich durch die baryzentri-schen Koordinaten in Tj beschreiben:

M lj = x ∈ Tj | λj(x) < λk(x), j 6= k

Das Volumen der Eckpunkte Ekj fallt mit dem Volumen der Massepaketeder 1. Verfeinerungsstufe zusammen und es gilt:

V (M lj) = V (Ekj ) =

14V (Tj)

Gut zu erkennen sind die Volumina in Abbildung 3.2 im Abschnitt 3.2.3. Ge-meinsam mit der gegebenen Konzentration u(Ki) an den Knoten Ki ergibtsich die den Eckpunkten zugeordnete Masse mit

m(Ekj ) =14V (Tj)u(Ki) fur Ekj = Ki

Bemerkung 3.2. Um den Speicheraufwand zu verringern, ware es auchmoglich die Masse statt eckpunktbasiert nur knotenbasiert zu speichern, al-so nur die Gesamtmasse, die dem Knoten Ki vom Patch aller Tetraederum Ki zugewiesen wird. Dies fuhrt jedoch zu hoherer numerischer Diffu-sion. Das Gebiet, das durch ein dorthin advehiertes Massepaket betroffenist, ware der gesamte Patch von Tetraedern um den Knoten Ki. An einemBeispiel kann man das gut veranschaulichen: Wird ein Massepaket in einTetraeder advehiert, wird die Masse auf die Eckpunkte des Tetraeders ver-teilt. Die Gesamtmasse an den Knoten berechnet sich aus der Gesamtmassean den Eckpunkten in Ki. Daraus bestimmt man mit dem Knotenvolumen


die Konzentration. Im nachsten Zeitschritt werden die Massepakete neu er-zeugt. Bei knotenbasierter Speicherung der Masse ist nicht mehr feststellbar,aus welchem Tetraeder des Patches um Ki die Masse kam. Die Masse mussdann im nachsten Zeitschritt auf alle Massepakete aller Tetraeder des Pat-ches um Ki verteilt werden. Durch die eckenbasierte Speicherung der Massekann man die Masse einem Tetraeder zuordnen. Die Masse wird nur auf dieMassepakete dieses Tetraeders verteilt. Dadurch wird also die numerischeDiffusion pro Zeitschritt auf ein Tetraeder beschrankt.

3.2.3 Unterteilung eines Tetraeders in Massepakete

Zur Unterteilung eines Tetraeders in Massepakete der 1. Verfeinerungsstufewird die bereits in 3.2.2 geschilderte Art der Unterteilung durch Schnittebe-nen verwendet. Die Abbildung 3.2 zeigt das Ergebnis dieser Unterteilung.

Bemerkung 3.3 (Zusammenhang mit Finite-Volumen-Methoden).Betrachtet sei die Vereinigung aller Massepakete der 1. Verfeinerungsstufeum jeden Knoten Ki. Die Familie aller so entstandenen Vereinigungsmen-gen wird Donald-Diagramm genannt. Es wird zur Unterteilung von bereitsvorhandenen Triangulierungen verwendet, um sog. Kontrollvolumina ein-zufuhren. Die Kontrollvolumina dienen dazu, den Fluss uber die Randereines Kontrollvolumens zu bilanzieren. Das Donald-Diagramm ist verwandtmit dem bekannteren Voronoi-Diagramm, das eine Unterteilung eines Ge-biets anhand einer gegebenen Punktmenge vornimmt. Naheres hierzu findetsich in [KA00].

Gerade bei adaptiven Gittern reicht die grobe Unterteilung durch Masse-pakete der 1. Verfeinerungsstufe nicht aus: Landet ein Massenschwerpunkteines Massepakets, das aus einem Gebiet niedriger Gitterstufe stammt, ineinem Tetraeder eines Gebiets hoher Gitterstufe, ist das Volumen des Mas-sepaketes viel großer als das des Tetraeders im hoch verfeinerten Gebiet. DieMasse wird nach dem hier verwendeten Verfahren jedoch nur diesem einenTetraeder zugewiesen. Dadurch entsteht ein untragbar hoher Fehler bei derApproximation der Konzentrationsverteilung im Gebiet um den Tetraederder hohen Gitterstufe. Daruberhinaus erhalten die benachbarten Tetraedergar keine Masse zugewiesen, da das Massepaket sie zwar ebenso uberdeckt,jedoch nur die Position des Massenschwerpunktes fur die Zuweisung vonMasse ausschlaggebend ist. Wir werden deshalb nun eine beliebig feine Un-terteilung eines Tetraeders einfuhren.

Zunachst erlautern wir das Prinzip der Unterteilung an einem zweidimensio-nalen Beispiel mit Dreiecken. Seien mit K sog. Kantenteilungspunkte undmit F Flachenteilungspunkte benannt. Kurz: Kantenteiler K bzw. Flachen-teiler F . Da die Unterscheidung verschiedener Teilungspunkte einer Ka-tegorie nicht weiter relevant ist, wird auf die Indizierung verzichtet. In


Abbildung 3.2: Unterteilung eines Tetraeders in Massepakete, hier: Verfeinerungs-stufe 1. Im Bild sind vier Teilmengen des Tetraeders zu erkennen. Gemeinsam mitder Masse und dem Volumen der Teilmengen wird so ein Massepaket definiert.Außerdem sind im Bild die Massenschwerpunkte der Massepakete zu sehen.

Abbildung 3.3 erkennt man zunachst die zu Abbildung 3.2 aquivalente Ver-sion der Unterteilung erster Stufe fur Dreiecke (links). Sie entsteht, indemman das Dreieck durch die Seitenhalbierenden in sechs Teilmengen teilt.Jedem Eckpunkt El, l = 1, 2, 3 des Dreiecks wird die Vereinigungsmengeder an ihm anliegenden Teilmengen zugeordnet. Dies sind die Massepake-te M l, l = 1, 2, 3 erster Verfeinerungsstufe in Dreiecken. In baryzentrischenKoordinaten beschreibt man die Menge um einen Eckpunkt mit

M l = x | λl(x) < λk(x), l 6= k

Wir konstruieren nun die 2. Verfeinerungsstufe: Zu den Eckpunkten trittnun ein Kantenteiler K fur jede Kante im Dreieck. Die Kantenteiler werdenauf den Kanten aquidistant verteilt. Die Kantenteiler der 2. Verfeinerungs-stufe sind also in den Kantenmittelpunkten platziert. Durch die Kantenteilerwerden nun parallel zu allen Seitenhalbierenden weitere Schnittgeraden ge-legt. Diese unterteilen mit den bereits vorhandenen Seitenhalbierenden dasDreieck in kleinere Teilmengen. Die Massepakete der 2. Verfeinerungsstu-fe entstehen nun durch die Vereinigung der Teilmengen, die einen Eckpunktbzw. einen Kantenteiler gemeinsam haben (siehe Abbildung 3.3: zweite Dar-stellung von links).


In der 3. Verfeinerungsstufe werden zwei Kantenteiler aquidistant auf dieKanten verteilt. Es tritt nun zusatzlich ein Flachenteiler hinzu. Auch hierwird wiederum mittels parallel zu allen Seitenhalbierenden durch die Kan-tenteiler gelegten Schnittgeraden das Dreieck in Teilmengen unterteilt. DieVereinigungsmenge der Teilmengen, die einen Eckpunkt, einen Kantenteileroder einen Flachenteiler gemeinsam haben, ergeben die Massepakete der 3.Verfeinerungsstufe (siehe Abbildung 3.3: zweite Darstellung von rechts). Furdie Massepakete der 4. Verfeinerungsstufe ergibt sich die in Abbildung 3.3ganz rechts dargestellte Aufteilung.Folgt man diesem Schema, dann kann man Schritt fur Schritt eine belie-big feine Unterteilung des Dreiecks erzeugen. Dabei wachst die Anzahl derKantenteiler und deren assoziierte Massepakete linear und die Anzahl derFlachenteiler und deren assoziierte Massepakete quadratisch.

Massepaket an einem Eckpunkt Massepaket an einem Kantenteiler Massepaket an einem Flächenteiler

1. VS 2. VS 3.VS 4.VS

Abbildung 3.3: Darstellung der Unterteilung in Massepakete in 2D mit Verfeine-rungsstufen 1,2,3,4 von links nach rechts.

Bei Tetraedern im Dreidimensionalen werden statt der Seitenhalbierendenim Dreieck die Flachenhalbierenden und deren parallele Ebenen durch dieKantenteiler und Flachenteiler gelegt, um den Tetraeder in Teilmengen zuunterteilen. Eine Flachenhalbierende verlauft durch eine Kante des Tetra-eders. Gemeinsam mit dem Kantenmittelpunkt der gegenuberliegenden Kan-te ist die Flachenhalbierende eindeutig bestimmt. Entsprechend der Anzahlvon sechs Kanten existieren im Tetraeder sechs Flachenhalbierende. Auchim Tetraeder werden die Massepakete durch die Vereinigungsmenge der anden Eckpunkten, Kantenteilern bzw. Flachenteilern anliegenden Teilmengendefiniert. Zusatzlich treten ab der Verfeinerungsstufe 4 noch Raumteiler Rhinzu. Die Anzahl der Raumteiler wachst kubisch.

Um die Massepakete-Semi-Lagrange-Methode durchfuhren zu konnen, istdie Berechnung des Volumens eines Massepaketes und der baryzentrischenKoordinaten des Massenschwerpunktes notwendig. Mittels der Bestimmungder Eckpunkte der Vereinigungsmengen um die Eckpunkte, Kanten-, Flachen-und Raumteiler kann man die Massenschwerpunkte und die Volumina derMassepakete berechnen. Die Massenschwerpunkte ergeben sich als geome-


Massepakettyp i Position Vi Anzahl ni1-Quader Eckpunkt 1

4k3V (T ) 4 ∀k ≥ 1

2-Teilkorper eines rh. D. Kantenteiler 1412k3V (T ) 6(k − 1) ∀k ≥ 2

3-Halbiertes rh. D. Flachenteiler 3k3V (T ) 4

∑kj=2(j − 2) ∀k ≥ 3

4-Rhombischer Dodekaeder Raumteiler 6k3V (T )

∑kj=3

(j−3)(j−2)2 ∀k ≥ 4

Tabelle 3.1: Ubersicht der Massepakettypen i eines Tetraeders fur beliebige Ver-feinerungsstufe k. Zusatzlich sind einige Eigenschaften der verschiedenen Massepa-kettypen aufgefuhrt. Mit rh. D. ist rhombischer Dodekaeder abgekurzt. V(T) istdas Volumen des Tetraeders, der von den Massepaketen unterteilt wird.

trischer Mittelpunkt der Massenpakete gemaß der Annahme, dass die Masseinnerhalb eines Massepaketes homogen verteilt ist. Die Bestimmung der Eck-punkte der Vereinigungsmengen ist rein technischer Natur und wird ausge-spart. In Tabelle 3.1 sind die Ergebnisse dieser Berechnungen bzgl. der Volu-mina der unterschiedlichen Massepakete fur beliebige Verfeinerungsstufen kaufgefuhrt. Zusatzlich ist noch angegeben, welche Anzahl ni an Massepake-ten in einem Tetraeder einer Verfeinerungsstufe k vorhanden ist. Außerdemwird angegeben, welcher geometrische Korper sich durch die Vereinigungs-mengen bzgl. der Position eines Massepaketes im Tetraeder ergibt. Bei derWahl eines regelmaßigen Tetraeder als Referenzelement sind alle Massepake-te um einen Raumteiler rhombische Dodekaeder. Fur Flachenteiler ergebensich halbierte rhombische Dodekaeder. In Abbildung 3.4 wird ein rhombi-sches Dodekaeder dargestellt, sowie dessen Halbierung. Fur die Kantenteilerergibt sich ein Teilkorper des rhombischen Dodekaeder. Fur die Massepaketean den Eckpunkten ergeben sich Quader.

Die gewahlte Art der Unterteilung durch Halbierende ist vor allem deshalbvorteilhaft, da sich die Teilungsverhaltnisse und Volumina der Massepaketedes hier gewahlten Referenzelements auf alle Elemente des Gitters Th durchaffin-lineare Abbildungen ubertragen lassen. Dies hat zur Folge, dass die ge-suchten Schwerpunkte und Volumina der Massepakete nur einmalig fur dasReferenzelement zu berechnen sind und daraufhin auf jedes (nicht degene-rierte) Element in der Triangulierung angewendet werden konnen (siehe auch[KA00]). Durch die regelmaßige Unterteilung des Tetraeders mit Flachen-halbierenden bzw. der zu den Flachenhalbierenden parallelen Schnittebenenergibt sich, dass alle Massepakete, die einem Eckpunkt zugeordnet sind, dasgleiche Volumen besitzen. Ebenso gilt fur Massepakete, die mit Kantentei-lern assoziiert sind, dass sie das gleiche Volumen besitzen. Analog gilt dies


a2a3

a3

32a3a

Abbildung 3.4: Darstellung eines rhombischen Dodekaeders und dessen Halbie-rung, die sich durch Teilung an der Seitenflache eines Tetraeders ergibt.

auch fur die mit Flachen- bzw. Raumteilern assoziierten Massepakete. InTabelle 3.1 finden sich die Volumina der unterschiedlichen Massepakettypenwieder.

Sei das Volumen des Referenzelement auf 1 normiert und seien mit M l

die Massepakete beliebiger Verfeinerungsstufe des Tetraeders T bezeichnet.Durch die Art der beschriebenen Unterteilung ergibt sich, dass das Innereder Massepakete int(M l) paarweise disjunkt ist. Es gilt fur jede Verfeine-rungsstufe des Tetraeders T :∑

l

V (M l) = V (T ) = 1

Wir kehren nun zu der Notation aus Abschnitt 3.2.1 zuruck. Wir konnennun das Volumen eines beliebigen Massepaketes V (M l

j) in einem TetraederTj ∈ Th angeben. Zusatzlich sind die Koordinaten des Massenschwerpunktesbekannt.Einem Massepaket M l

j wird proportional zu den baryzentrischen Koordina-ten seines Massenschwerpunktes λlk, k = 1, 2, 3, 4 bezogen auf die EckpunkteEkj und zu dem Volumen des Massepaketes V (M l

j) eine Teilmasse des Te-traeders Tj wie folgt zugeordnet:

m(M lj) =

∑k

λlkm(Ekj )V (M lj) (3.1)

m(Ekj ) ist die an den Eckpunkten Ekj vorliegende Masse des Tetraeders Tj .Aus der disjunkten vollstandigen Unterteilung von Tj in Massepakete M l

j


folgt∑

l V (M lj) = 1 zusammen mit der Eigenschaft der baryzentrischen

Koordinaten∑

k λlk = 1 ergibt sich, dass∑

l

m(M lj) =

∑k

m(Ekj )

gilt. Die Masse der Eckpunkte m(Ekj ) eines Tetraeders Tj wird also exaktauf die Massepakete M l

j von Tj abgebildet.

Bemerkung 3.4 (Approximation exakter Schnitte zwischen Tetra-edern). Die Unterteilung der Tetraeder in Massepakete und die Zuweisungder Massen auf die Eckpunkte des stromaufwarts advehierten Gitters stellteine Approximation der Berechnung der exakten Schnitte zwischen der Tri-angulierung zum alten Zeitschritt und dem stromaufwarts advehierten Git-ter zum neuen Zeitschritt dar. Laßt man die Verfeinerungsstufe l aller Te-traeder des alten Gitters gegen unendlich gehen und ubertragt mittels derMassepakete die Masse auf die Tetraeder des neuen Gitters ist die Berech-nung der Bestimmung der exakten Schnitte zwischen dem alten und demneuen stromaufwarts advehierten Gitter aquivalent. Hier besteht der Zu-sammenhang zu den cell-integrated Semi-Lagrange-Methoden (CISL), diedie Massenerhaltung uber die Berechnung der exakten Schnitte erreichen.

Bei der Wahl der Verfeinerungsstufe im Tetraeder ist abzuwagen zwischender Genauigkeit der Approximation der Konzentrationsverteilung in Tj undder Durchfuhrbarkeit des Verfahrens hinsichtlich der vorhandenen Ressour-cen. Zu niedrige Verfeinerung laßt den Approximationsfehler zu groß werden.Mit zu hoher Verfeinerung erreicht man wegen des kubischen Anwachsensder Anzahl der Massepakete pro Tetraeder schnell die Grenze des Speichers.Als in der Anwendung praktikabel hat sich die folgende a-priori-Wahl her-ausgestellt: Es wird einerseits die Gitterstufe des zu unterteilenden Tetra-eders im Vergleich zur grobsten bzw. feinsten Gitterstufe in dem adaptivenGitter als Kriterium gewahlt und andererseits die Konzentration im Te-traeder im Vergleich zur niedrigsten bzw. hochsten Konzentration in derTriangulierung. Die beiden Kriterien werden multiplikativ verknupft. Dasbedeutet, dass ein Tetraeder der grobsten Gitterstufe und der hochsten Kon-zentration mit der maximal vorgegebenen Verfeinerungsstufe in Massepaketeunterteilt wird. Tetraeder der feinsten Gitterstufe und niedriger Konzentra-tion werden dagegen mit der minimalen vorgegebenen Verfeinerungsstufeunterteilt. Fur die verwendeten Testfalle aus Kapitel 4 mit 7 GitterstufenUnterschied zwischen grobstem und feinstem Gitter haben sich Verfeine-rungsstufe 4 bis Verfeinerungsstufe 14 als gute Werte herausgestellt. Beidieser Gelegenheit sei darauf hingewiesen, dass die adaptive Verfeinerung derTriangulierung pro Gitterstufe eine sukzessive Halbierung des Volumens zurFolge hat, also mit exponentiellem Abfall des Volumens einhergeht, wahrenddas Volumen der Massepakete mit jeder Verfeinerungsstufe nur kubisch ge-ringer wird.


Optional kann man zusatzlich noch eine Erhohung der Verfeinerungsstufeder Tetraeder in Massepakete zuschalten. Das verwendete Kriterium, umzu entscheiden ob ein Tetraeder durch mehr Massepakete unterteilt wird,ist das Gesamtvolumen, das die Massepakete besitzen, die einem Tetraederdes stromaufwarts gelegenen Gitters zugewiesen werden. Ist das Volumenwesentlich zu hoch, werden die Tetraeder des alten Gitters, aus denen dieMassepakete stammen, feiner unterteilt. Dies kann auch mit mehreren Ite-rationen durchgefuhrt werden.

Bemerkung 3.5 (MPSLM und Adaptive Gitter). Ein Teil des Ge-schwindigkeitsgewinns, der durch Adaptivitat erzielt wird, geht durch diestarkere Unterteilung von Tetraedern grober Gitterstufen in Massepaketewieder verloren.

3.2.4 Verfolgung der Triangulierung des nachsten Zeitschrittsstromaufwarts

Dieser Teil der Semi-Lagrange-Methode wird beibehalten, allerdings wer-den nun ganze Elemente T ∈ T n+1

h stromaufwarts verfolgt. Dies geschiehtmit den Trajektorien α, die fur jeden Knoten K der Triangulierung T n+1

h

berechnet werden. Ein Massepaket M , das aus einem Tetraeder der Tri-angulierung T nh hervorgegangen ist, ist genau dann in einem Element Tder stromaufwarts verfolgten Triangulierung T n+1

h , wenn der Massenschwer-punkt von M in T liegt. Der Suchalgorithmus, der bestimmt, in welchem Tder Massenschwerpunkt liegt, wird in Abschnitt 3.3.3 beschrieben.

3.2.5 Massenzuweisung an die Eckpunkte der Tetraeder des neu-en Gitters

Ist ein Massepaket M li , das aus Ti ∈ T nh hervorgegangen ist, in einem strom-

aufwarts gelegenen Tetraeder Tnj gefunden, werden die baryzentrischen Ko-ordinaten µlk, k = 1, . . . , 4 des Massenschwerpunkts im Tetraeder Tnj be-rechnet. Zu beachten ist, dass die baryzentrischen Koordinaten des Massen-schwerpunkts λ auf Eckpunkte eines Tetraeders Ti im alten Gitter und die µauf die Eckpunkte eines stromaufwarts advehierten Tetraeders Tnj aus demneuen Gitter bezogen sind. Die Masse wird dann gemaß diesen Koordinatenµlk den Eckpunkten Ekj des Tetraeders Tnj zugewiesen. Die Gesamtmassem(Ekj ) an den Eckpunkten des Tetraeders Tnj ergibt sich wie folgt:

m(Ekj ) =∑

M li∈Tnj

µlkm(M li ) fur alle k = 1, . . . , 4 (3.2)

Durch die Eigenschaft der baryzentrischen Koordinaten mit4∑

k=1

µlk = 1,∀l ist

sicher gestellt, dass die gesamte Masse des Massepaketes auf die Eckpunkte


Ekj verteilt wird. Die den Eckpunkten zugewiesene Gesamtmasse wird ge-speichert, um im nachsten Zeitschritt wieder von neuem auf Massepaketeverteilt zu werden.

Satz 3.6 (Globale Massenerhaltung der MPSLM). Wird die in Ab-schnitt 3.2.3 vorgestellte Aufteilung der Masse auf Massepakete und die indiesem Abschnitt dargestellte Zuweisung der Masse der Massepakete auf dieTriangulierung des nachstes Zeitschritts durchgefuhrt, so ist die Massepakete-Semi-Lagrange-Methode bis auf Massenverlust uber den Gebietsrand Γ mas-senerhaltend fur alle Zeitschritte tn, n = 1, . . . , N .

Beweis:Die Gesamtmasse m(n)(Ω) in Ω zum Zeitschritt tn ist die Summe aller Mas-sen der Eckpunkte Eki der Tetraeder Ti in der Triangulierung T nh von Ω:

m(n)(Ω) =∑i,k

m(Eki ) (3.3)

Die Verteilung der Massen m(Eki ) auf Massepakete M li gewichtet nach bary-

zentrischen Koordinaten λlk, k = 1, . . . , 4 und den Volumina der MassepaketeV (M l

j) erfolgt nach (3.1) mit

m(M li ) =

∑k

λlkm(Eki )V (M li ) fur alle l

Diese Formel garantiert wie in Abschnitt 3.2.3 gezeigt, dass fur die Gesamt-masse aller Massepakete in Ω gilt:∑

i,l

m(M li ) =

∑i,k

m(Eki )(3.3)= m(n)(Ω) (3.4)

Die Massen der Massepakete m(M li ) werden nach der Vorschrift (3.2) auf

die Eckpunkte Ekj der stromaufwarts gelegenen Tetraeder Tnj der Triangu-lierung T n+1

h verteilt.Unter Berucksichtigung der Eigenschaft der baryzentrischen Koordinaten∑

k µk = 1 ergibt sich, dass die gesamte Masse aus T nh auf die Triangulierungdes nachsten Zeitschritts T n+1

h ubertragen wird. Dies gilt bis auf Massepa-kete M l

i , die außerhalb der stromaufwarts verfolgten Triangulierung T n+1h

liegen: ∑j,k

m(Ekj ) =∑i,l

m(M li )−

∑M li∩T

n+1h =∅

m(M li )

Setzt man die bisher erzielten Ergebnisse ineinander ein, ergibt sich dieGesamtmasse m(n+1)(Ω) in Ω zur Zeit tn+1 wie folgt:

m(n+1)(Ω) =∑j,k

m(Ekj )(3.4)= m(n)(Ω)−

∑M li∩T

n+1h =∅

m(M li ) (3.5)


Damit ist gezeigt, dass bis auf Massenverlust uber den Gebietsrand dieMasse in Ω der Triangulierung T nh exakt auf die Triangulierung T n+1

h desnachsten Zeitschritts ubertragen wird. Sei nun mit m(n)

Verlust der Massenver-lust

∑M li∩T

n+1h =∅m(M l

i ) bei der Ausfuhrung des n-ten Zeitschritts bezeich-net. Dann gilt fur die Masse in Ω zum Endzeitschritt N unter Berucksichti-gung von (3.5):

m(N)(Ω) = m(0)(Ω)−∑n

m(n)Verlust mit n = 1, . . . , N

Dabei ist m(0)(Ω) die Startmasse.Damit ist gezeigt, dass die Massepakete-Semi-Lagrange-Methode massener-haltend bis auf Massenverluste uber den Rand Γ des Gebiets Ω ist.

2

Ohne es hier zu zeigen ist anschaulich klar, dass die MPSLM sogar lokalmassenerhaltend ist: Die Masse eines Tetraeders zum alten Zeitschritt wirdexakt auf die Massepakete verteilt, um dann wiederum auf die stromaufwartsadvehierten Tetraeder des neuen Gitters exakt verteilt zu werden. Je nachGroße der Massepakete (also der Verfeinerungsstufe und Gitterstufe einesTetraeders) konnen Fehler auftreten, falls ein Massepaket nur teilweise ineinem Tetraeder liegt. Jedoch ist der auftretende Fehler lokal begrenzt. BeiWindstille, also a(x, t) = 0 fur alle x ∈ Ω, t > 0, gilt sogar exakte lokaleMassenerhaltung im Tetraeder.

3.2.6 Konzentrationsberechnung

Um die Konzentration am KnotenKi der Triangulierung T n+1h zu berechnen,

summiert man die Masse aller Eckpunkte in Ki:

m(Ki) =∑

Ekj =Ki

m(Ekj )

Zusammen mit dem Knotenvolumen um Ki

V (Ki) =∑

Ki∈Tn+1j

14V (Tn+1

j )

erhalt man schließlich die Konzentrationswert fur den Knoten Ki

u(Ki) =m(Ki)V (Ki)

Dabei wird davon ausgegangen, dass unter den Tetraedern Tn+1j der Trian-

gulierung T n+1h keine degenerierten Tetraeder mit V (Tn+1

j ) = 0 vorkommen.


3.3 Implementierung: Datenstrukturen und Algorithmen

Die hier beschriebene Massepakete-Semi-Lagrange-Methode ist eingebettetin das Programmpaket amatos3d/flash3d von Jorn Behrens. Fur die 2D-Version dieses Pakets findet sich eine Beschreibung in [Beh96]. amatos3dstellt den adaptiven Finite-Elemente-Gittergenerator zur Verfugung undflash3d setzt mit dem Semi-Lagrange-Verfahren und der Steuerung der Ad-aptivitat mittels eines Fehlerschatzers auf amatos3d auf.

Im Flussdiagramm (siehe Abbildung 3.5) ist die Einbettung der MPSLMin das Programmpaket amatos3d/flash3d von Behrens dargestellt. Die vor-genommenen Anderungen und Erganzungen, um diese Einbettung zu reali-sieren, sind grau unterlegt. Wie bereits in den vorhergehenden Abschnittenerlautert, wird die Semi-Lagrange-Methode durch die vom Autor entwickel-te massenerhaltende Massepakete-Semi-Lagrange-Methode ersetzt. Hierbeiist vor allem der Programmabschnitt ”Neue Werte berechnen“ betroffen.Der genaue Ablauf der MPSLM ist in der Darstellung 3.5 rechts. Erkennbarist die optionale starkere Verfeinerung der Tetraeder durch Massepakete,falls das Gesamtvolumen der Massepakete

∑V (M), die in einem Tetra-

eder gefunden wurden, das Volumen des jeweiligen Tetraeders V (T ) starkubersteigt. Diese Option ist bei den in Kapitel 4 durchgefuhrten Tests ausge-schaltet, da sich nur zu vernachlassigende Verbesserung des Fehlers ergeben.Zu den bereits in flash3d vorhandenen Datenstrukturen kommen noch dreiDatenstrukturen hinzu, die im folgenden Abschnitt 3.3.1 vorgestellt werden.Von zentraler Bedeutung ist die Berechnung der Parameter des Referenz-elements fur die verschiedenen Massepakete-Verfeinerungsstufen eines Te-traeders. Dieser Vorgang wird in Abschnitt 3.3.2 ausfuhrlich beschrieben.Die Suche der Massenschwerpunkte der Massepakete in den Tetraedern desstromaufwarts verfolgten Gitters wird abschließend in Abschnitt 3.3.3 dar-gestellt. Alle anderen Verfahrensschritte ergeben sich aus der im vorigen Ab-schnitt dargestellten Vorgehensweise und den dort angegebenen Formeln.Die hier vorgestellten Algorithmen sind als Fortran 90-Code realisiert undbefinden sich in den Dateien SLM interpolation.f90, SLM generateMP.f90.Die Datenstrukturen sind in SLM MPtypes.f90 definiert. Diese Dateien be-finden sich auf der beiliegenden CD zum Beispiel im Unterverzeichnis/MPSLM/lin/flash3d/src/flash.


MPSLM-Programm - Start

Initialisierung Preprocessing

Speichere altes Gitter

Trajektorien berechnen

Fehler schätzen Gitter anpassen

Gitter geändert?

Alle Zeitschritte durchgeführt?

Diagnose durchführen Grafik ausgeben

Postprocessing

MPSLM-Programm Stop

Neue Werte berechnen mit masseerhaltendem Verfahren

Informationen aus altem und neuem Gitter auslesen

Upstream-Koordinaten und -Volumen der Tetraeder des neuen Gitters berechnen

A priori-Schätzer zur Berechnung der Verfeinerungsstufe der Tetraeder im alten Gitter

Massepakete-Tabelle entsprechend den Tetraeder-Verfeinerungsstufen mittels MPK-Tabelle erzeugen

Horizontale Suche der Masse-pakete in den stromaufwärts adve-hierten Tetraedern im neuen Gitter

)()( TVMVTM

>>∑∈

Für alle Tetraeder: Masse der Massepakete auf die Eckpunkte der Tetraeder verteilen; Massen an Ecken in der Tetraedertabelle für nächsten Zeitschritt speichern

Volumen der Tetraeder im neuen Gitter berechnen

Neue Konzentrationswerte an den Knoten des neuen Gitters berechnen

Anlegen der Massepakete-koeffizienten-Tabelle zur späteren Berechnung der Massepakete

Anlegen der Tetraeder-tabelle zur Speicherung der Massen an den Eckpunkten über alle Zeitschritte hinweg

Exakter Fehler bezogen auf das Modellproblem berechnen

Erzeugte Tabellen deallokieren

ja

nein

Erhöhe Verfei-nerungsstufe der Ursprungstetraeder der Massepakete M

nein

ja

ja

nein

Legende

Übernommene Elemente aus amatos3d/flash3d

Modifizierte Elemente aus amatos3d/flash3d

Neue Elemente

Abbildung 3.5: Flussdiagramm zur Darstellung des Ablaufs von ama-tos3d/flash3d mit MPSLM. Grau unterlegt sind die vorgenommenen Anpassungen,um die MPSLM in den bisherigen Programmablauf einzubetten. Der bisherige SLM-Algorithmus zur Berechnung der Konzentrationswerte zum neuen Zeitschritt wirdkomplett ausgetauscht und durch den dargestellten MPSLM-Algorithmus ersetzt.


.

.

.

.

Tetraeder T1

Ti

.

.

.

.

Tn

.

.

.

.

Tetraeder T1

T j

.

.

.

.

.

.

Tm

Tetraedertabelle für Zeitschritt tn

Tetraedertabelle für Zeitschritt tn+1

i_MPperTetra: Anzahl der MP in Tj

i_MPindex[i]: Indizes der durch das Tetraeder erzeugten Massepakete

i_refinelvl: Gitterstufe des Tetraeders

l_valueassigned: Ist Masse an den Eckpunkten vorhanden?

l_searchtetra: Sollen die Massenpakete des Tetraeder gesucht werden?

l_refine: Soll der Tj stärker durch Massepakete verfeinert werden?

r_volume: Volumen von Tj

r_value[4]: Massen an den Eckpunkten von Tj

Abbildung 3.6: Aufbau der Tetraeder-Tabelle. Mit einem Beispiel fur ein Daten-element aus der Tabelle. Die Namen der Variablen im Datenelement stammen ausdem Programmcode.

3.3.1 Datenstrukturen

Zu den in amatos3d/flash3d bereits vorhandenen werden zusatzlich drei we-sentliche Datenstrukturen erzeugt:

1. Tetraedertabelle: Liste aller Tetraeder zu den Zeitschritten tn, tn+1

2. Massepakete-Koeffiziententabelle: Koeffizienten zur Erzeugung der Mas-sepakete vorgegebener Verfeinerungsstufe

3. Massepaketetabelle: Liste aller erzeugter Massepakete

In der Tetraedertabelle sind die Tetraeder Ti, i = 1, . . . ,M der Triangu-lierung T nh , sowie die Tetraeder Tj , j = 1, . . . , N der Triangulierung T n+1

h

gespeichert. Alle verfahrensrelevanten Daten, die durch die Tetraeder gege-ben sind, werden hier gespeichert. In Abbildung 3.6 ist die Datenstrukturdargestellt. Beispielhaft ist ein Datenelement aus der Tabelle angegeben.Darin werden die Namen der Variablen im Programmcode und ihre Funk-tion bezuglich des Programmablaufs dargestellt. Die Koordinaten der Eck-punkte eines Tetraeders sind hier nicht gespeichert. Sie werden bei Bedarfaus dem Gittergenerator amatos3d durch einen entsprechenden Funktions-aufruf ausgelesen.Die Variablen i MPperTetra, i MPindex[i] dienen zur Verknupfung der


i_origin: Ursprungsecke Projektion wird in dieser Richtung vorgenommen

Verfeinerungsstufe 1

.

.

.

.

VS i

.

.

.

.

VS n

Anzahl der Masse-pakete m in der VS i

Pointer auf Massepakete-Koeffizienten-Array

Pointer auf Massepakete-Koeffizienten M1

.

.

.

.

Pointer auf Mj

.

.

.

.

Pointer auf Mm

Pointerarray Pointerarray auf die Masse-pakete-Koeffizienten

Parameter des Massepaketes Mj in VS i

i_prvdir[3]: Vorhergehende Richtungen im Baum

i_celltype: Massepaket an Ecken-, Kanten-, Flächen- oder Raumteiler

r_mass[4]: Koeffizienten zur Zuweisung des Massenanteils von Mj im Tetraeder

r_volume: Volumen des Massepaketes Mj im Referenzelement

r_coeff[4]: Baryzentrische Koordinaten des Massen-schwerpunkts im Referenzelement

l_end: Technische Variable zur Erweiterung des Baumes bei Berechnung der VS i+1; Ist das Massepaket am Ende des Baumes?

.

.

.

.

Pointer auf Massepaket M1

Pointer auf Mj

.

.

.

.

Pointer auf Mn

i_MPinTetra: Index des Tetraeders, in den das Massepaket advehiert wird

i_MPbyTetra: Index des Tetraeders, von dem das Massepaket stammt

r_mass: Masse des Massepaketes

l_searchMP: Soll das Massepaket noch gesucht werden?

l_assignvalue: Soll die Masse auf die Eckpunkte des Tetraeder mit dem Index i_MPinTetra verteilt werden?

r_volume: Volumen des Massepaketes

r_coord[3]: Koordinaten des Massenschwerpunkts in globalen Koordinaten

Massepakete-Array Einträge des Massepakets Mj

Abbildung 3.7: Aufbau der Massepakete-Koeffiziententabelle. Die einmal berech-neten Koeffizienten dienen spater der Generierung der Massepakete. Der relativkomplizierte Aufbau uber mehrere Pointer ist durch die Forderung der dynami-schen Erzeugung fur beliebige Verfeinerungsstufen n begrundet.

Tetraedertabelle mit der spater erlauterten Massepaketetabelle. Sie sind nurbelegt, falls l valueassigned wahr ist. Massepakete werden fur einen Te-traeder also nur erzeugt, falls tatsachlich Masse an den Eckpunkten vor-liegen. Dies kann mittels der Variable r value[4] festgestellt werden, inder die Masse eines Tetraeders Ti ∈ T nh gespeichert ist. Eckige Klammernbedeuten, dass ein Array vorliegt. Ist keine konstante Zahl angegeben han-delt es sich um einen dynamisch allokierten Array, der je nach Bedurfnisangepasst wird. In der Variable i MPindex[i] z.B. werden die Indizes derMassepakete des Tetraeders gespeichert. Die Lange des Arrays variiert jenach Verfeinerungsstufe des Tetraeders. Im Gegensatz dazu bleibt die An-zahl der Eckpunkte eines Tetraeders konstant und somit auch die Großedes vorzuhaltenden Speichers fur die Massen an den Eckpunkten r value.Die Variable l refinelvl dient der Speicherung der Gitterstufe des Tetra-eders im adaptiven Gitter. Sie wird zur Berechnung der Verfeinerungsstufedes Tetraeders verwendet. l searchtetra, l refine dienen der optionalenNachverfeinerung des Tetraeders mit mehr Massepaketen. l refine zeigtan, ob der Tetraeder nachverfeinert werden soll. Bei der erneuten Sucheder Massepakete zeigt l searchtetra an, welche Tetraeder in der Sucheubergangen werden konnen, da sie nicht nachverfeinert wurden. r volumeschließlich speichert das Volumen des Tetraeders.


i_origin: Ursprungsecke Projektion wird in dieser Richtung vorgenommen

Verfeinerungsstufe 1

.

.

.

.

VS i

.

.

.

.

VS n

Anzahl der Masse-pakete m in der VS i

Pointer auf Massepakete-Koeffizienten-Array

Pointer auf Massepakete-Koeffizienten M1

.

.

.

.

Pointer auf Mj

.

.

.

.

Pointer auf Mm

Pointerarray Pointerarray auf die Masse-pakete-Koeffizienten

Parameter des Massepaketes Mj in VS i

i_prvdir[3]: Vorhergehende Richtungen im Baum

i_celltype: Massepaket an Ecken-, Kanten-, Flächen- oder Raumteiler

r_mass[4]: Koeffizienten zur Zuweisung des Massenanteils von Mj im Tetraeder

r_volume: Volumen des Massepaketes Mj im Referenzelement

r_coeff[4]: Baryzentrische Koordinaten des Massen-schwerpunkts im Referenzelement

l_end: Technische Variable zur Erweiterung des Baumes bei Berechnung der VS i+1; Ist das Massepaket am Ende des Baumes?

.

.

.

.

Pointer auf Massepaket M1

Pointer auf Mj

.

.

.

.

Pointer auf Mn

i_MPinTetra: Index des Tetraeders, in den das Massepaket advehiert wird

i_MPbyTetra: Index des Tetraeders, von dem das Massepaket stammt

r_mass: Masse des Massepaketes

l_searchMP: Soll das Massepaket noch gesucht werden?

l_assignvalue: Soll die Masse auf die Eckpunkte des Tetraeder mit dem Index i_MPinTetra verteilt werden?

r_volume: Volumen des Massepaketes

r_coord[3]: Koordinaten des Massenschwerpunkts in globalen Koordinaten

Massepakete-Array Einträge des Massepakets Mj

Abbildung 3.8: Aufbau der Massepakete-Tabelle. Mit einem Beispiel fur ein Da-tenelement in der Tabelle. Die Namen der Variablen stammen aus dem Programm-code.

Nun kommen wir zu der in Abbildung 3.7 dargestellten, etwas kompliziertaufgebauten Massepakete-Koeffiziententabelle (MPK-Tabelle). Sie speichertdie Koeffizienten zur Erzeugung der Massepakete. Wie bereits in Teilab-schnitt 3.2.3 dargestellt, ist es nur notig diese Koeffizienten einmalig fur dasReferenzelement zu berechnen und sie dann geeignet mit dem zu unterteilen-den Tetraeder zu verknupfen, um schließlich die Massepakete zu erzeugen.Der komplizierte Aufbau uber mehrere Pointer-Ebenen hat seinen Grundin der dynamischen, sukzessiven Generierung der Tabelle. Aus der Verfei-nerungsstufe 1 wird also die Verfeinerungsstufe 2 abgeleitet, 3 aus 2, biszur gewunschten hochsten Verfeinerungsstufe n. Der Generierungsprozessder Koeffizienten wird in Teilabschnitt 3.3.2 dargestellt. Dort wird auch dieVerwendung der Variablen naher erlautert.

In der Massepakete-Tabelle werden die notigen Variablen bezuglich der Mas-sepakete gespeichert (siehe Abbildung 3.8). Die Werte der Variablen erge-ben sich einerseits aus der Verknupfung der Tetraedertabelle und der MPK-Tabelle. Hierzu zahlt i MPbyTetra, die den Index wiedergibt, aus welchemTetraeder das Massepaket erzeugt wurde. r mass und r volume speichern diedem Massepaket zugeordnete Masse sowie Volumen. Die Variable r coordgibt die globalen Koordinaten des Massenschwerpunkts des Massepaketeswieder. Sie ergibt sich aus den globalen Koordinaten der Eckpunkte des Te-traeders und den baryzentrischen Koordinaten des Massepaketes, die in der


MPK-Tabelle zu finden sind. Andererseits sind in der Massepakete-TabelleVariablen zur Durchfuhrung der Suche, in welchem stromaufwarts adve-hierten Tetraeder das Massepaket zu finden ist, gespeichert. Dazu zahlendie bereits erwahnten globalen Koordinaten r coord. Hinzu kommt nochi MPinTetra. Diese Variable speichert den Index, in welchem stromaufwartsadvehierten Tetraeder das Massepaket liegt. Falls l searchMP wahr ist, mussdas Massepaket noch gesucht werden. l assignvalue gibt an, ob die Massedes Massepaketes einem Tetraeder zugewiesen werden kann. Dies ist dannnicht der Fall, wenn das Massepaket außerhalb des stromaufwarts advehier-ten Gebiets liegt.

3.3.2 Algorithmus: Unterteilung eines Tetraeders in Massepakete

Der Algorithmus dient zur Generierung der Massepakete-Koeffiziententabelle(MPK-Tabelle). Er berechnet die im Abschnitt 3.2.3 zur Durchfuhrung desMPSLM-Algorithmus benotigten Werte r coeff, r mass, r volume. Da-bei gibt r coeff die baryzentrischen Koordinaten eines Massenschwerpunk-tes wieder, r mass beinhaltet die Koeffizienten zur Berechnung des Massen-anteils eines Massepaketes und r volume wird zur Berechnung des Volumensdes Massepaketes verwendet. Die genannten Koeffizienten werden einmaligfur das Referenzelement berechnet und in der MPK-Tabelle gespeichert.Dieser Algorithmus wird wahrend der Initialisierung von flash3d durch-gefuhrt. Die Daten stehen dann fur alle Zeitschritte zur Verfugung.

1. VS 2. VS 3. VS 4. VS

Abbildung 3.9: Veranschaulichung der Erzeugung der Massepakete-Koeffiziententabelle. Die nachsthohere Verfeinerungsstufe entsteht durch hin-zukommende Massepakete und lasst sich durch geeignete Projektion auf dasReferenzelement zuruckfuhren. Hier in 2D mit Verfeinerungsstufen 1,2,3,4 vonlinks nach rechts.

Das grundlegende Prinzip, das der Generierung der Koeffizienten zu Grun-de liegt, soll zunachst anhand Abbildung 3.9 in 2D dargestellt werden. Wirbetrachten die Massepakete an den Eckpunkten der Dreiecke. Man erkennt,dass die Flache dieser Massepakete in jedem Dreieck der Abbildung identisch


ist. Die Anzahl der Massenpakete an Eckpunkten bleibt konstant. Auch dieFlachen der Massepakete an Kantenteilern ist identisch. Es tritt jedoch proVerfeinerungsstufe an jeder Kante jeweils ein Massepaket hinzu. Die Flacheder Massepakete an Flachenteilern ist ebenfalls identisch. In Verfeinerungs-stufe 3 tritt erst ein Flachenteiler-Massepaket hinzu und in der nachstenVerfeinerungsstufe dann zwei. In Verfeinerungsstufe 5 kommen drei Masse-pakete hinzu.Betrachtet man im dreidimensionalen Fall den Tetraeder ergeben sich nochRaumteiler-Massepakete. Das erste Massepaket dieser Art tritt in Verfeine-rungsstufe 4 auf. In Verfeinerungsstufe 5 kommen drei weitere hinzu und inVerfeinerungsstufe 6 dann bereits 6 Massepakete. Allgemein ergibt sich furTetraeder mit bekanntlich 4 Eckpunkten, 6 Kanten und 4 Seitenflachen fol-gende Massepakete-Anzahl n(i) an EckpunktenE, KantenteilernK, Flachen-teilern F und Raumteilern R fur die Verfeinerungsstufe i:

n(i)E = 4 ∀i ≥ 1n

(i)K = n

(i−1)K + 6 ∀i ≥ 2

n(i)F = n

(i−1)F + 4(i− 2) ∀i ≥ 3

n(i)R = n

(i−1)R + (i−3)(i−2)

2 ∀i ≥ 4

In jeder weiteren Verfeinerungsstufe kommen also zu den bereits bestehen-den Massepaketen neue hinzu.Letztendlich mochte man jedoch nur die Position der Massenschwerpunkteund das Volumen aller Massepakete im Referenzelement. Mittels Zentralpro-jektion des Tetraeders der Verfeinerungsstufe i bzw. der Massenschwerpunk-te innerhalb dieses Tetraeders auf seinen Vorganger der Verfeinerungsstufei−1 und der Festlegung des Tetraeders der Verfeinerungsstufe 1 als Referen-zelement kann man dies erreichen. Als Projektionspunkt dient ein Eckpunktdes Tetraeders. Es ergibt sich als Projektionsfaktor p(i):

p(i) =i− 1i

∀i > 1

Mittels des Strahlensatz ergibt sich bei Projektion in Richtung des Eckpunk-tes Ek fur baryzentrische Koordinaten im Tetraeder folgende Formel:

λ(i)j =

(1− p(i)) + p(i)λ

(i−1)j fur j = k,

p(i)λ(i−1)j sonst

(3.6)

Betrachtet man die Aneinanderreihung von l Projektionen ergibt sich alsProjektionsfaktor:

p(l) =l∏1

p(i) =1l


Damit erkennt man nun leicht, dass in (3.6)

liml→∞

λ(i)j =

1, fur j = k,

0, sonst

gilt. Ein Massepaket ruckt also mit jeder Projektion dem Eckpunkt Eknaher. Es sei daran erinnert, dass die Zentralprojektion winkeltreu, paralle-lentreu und streckenteilungstreu ist und damit der geometrische Mittelpunktdes projizierten Korpers durch Projektion des geometrischen Mittelpunktsdes Ausgangskorpers berechnet werden kann. Insbesondere gilt dies also furdie hier gesuchten Massenschwerpunkte der Massepakete. Der Volumenanteiljedes einzelnen Massepakete am Gesamtvolumen des Tetraeders vermindertsich mit jeder Projektion wie folgt:

V(i)M = (

1p(i)

)3V(i−1)M (3.7)

Nach diesen Voruberlegungen gehen wir im Folgenden auf die Generierungder MPK-Tabelle ein. Dabei wird die Position der Massenschwerpunkte ausder jeweils vorhergehenden Verfeinerungsstufe abgeleitet. Wir betrachtendas Vorgehen nach den einzelnen Massepaket-Typen getrennt.

Eckpunkt-Massepakete Zunachst betrachten wir die Massepakete an denEckpunkten E1, E2, E3, E4. Sie bilden die Verfeinerungsstufe 1. Fur die ba-ryzentrischen Koordinaten r coeff der Massenschwerpunkte ergibt sich:

λ(1)j =

4596 , falls j = k,1796 , sonst.

Der Volumenanteil der Massepakete am Gesamtvolumen des Referenzele-ments betragt in der 1. Verfeinerungsstufe VE = 1

4 . Er wird in der Va-riable r volume gespeichert. Dieser Anteil vermindert sich mit jeder weite-ren Verfeinerungsstufe gemaß (3.7). Die Position der Massenschwerpunkteder hoheren Verfeinerungsstufen ergeben sich durch die sukzessive Projekti-on der Massenschwerpunkte in Richtung ihrer Eckpunkte nach der Formel(3.6). Zu diesem Zweck wird in der Variable i origin der Eckpunkt, demdas Massepaket zugeordnet ist, gespeichert.

Kantenteiler-Massenpakete Diese Massepakete treten ab der 2. Verfei-nerungsstufe auf. In dieser haben sie folgende baryzentrischen Koordinatenr coeff:

λ(2)j =

2666 , falls j ∈ k,m766 , sonst.

wobei Ek, Em die Eckpunkte der Kante sind, an dem das Kantenteiler-Massepaket liegt. Der Volumenanteil r volume der Kantenteiler-Massepakete


4

3

2

1

0

K

K

K

E

K

K

K

K

K

F F

SS

F

SS

K

SS

E

K

K

K

K

K

F F

SS

F

SS

K

K

K

F F

SS

F

F F R

ZZZ

F

HHHHH

K

(((((((((((((

PPPPPP

E

((((((((((((((((((((

(((((((((((((

XXXXXXXX

E

Abbildung 3.10: Baumstruktur fur die Datenstruktur-Generierung derMassepakete-Koeffizienten. E, K, F, R stehen dabei fur die geometrische Positionder Massepakete im Tetraeder an einer Ecke E, an einem Kantenteiler K, an einerFlachenteiler F bzw. an einem Raumteiler R. Die Ziffern links des Baums stehenfur die Verfeinerungsstufe in der die Massepakete generiert werden. Sie werden ausden Massepaketen der vorherigen Verfeinerungsstufe abgeleitet. Die Massepaketeeiner Verfeinerungsstufe ist die Gesamtheit aller Massepakete der in dieser Stufeneu generierten und der bereits vorher generierten Massepakete.

am Gesamtvolumen des Referenzelements betragt anfangs VK = 1496 . Auch

hier wird jedem Massepaket in der Variable i origin ein Eckpunkt zuge-ordnet in dessen Richtung das Massepaket mit jeder Erhohung der Ver-feinerungsstufe projiziert wird. In jeder neuen Verfeinerungsstufe tritt proKante ein weiteres Kantenteiler-Massepaket hinzu. Die baryzentrischen Ko-ordinaten dieses Massepaketes ergibt sich aus der Uberlegung, dass die Pro-jektion der Kantenteiler-Massepakete entweder in Richtung des einen oderdes anderen Eckpunktes durchgefuhrt werden kann. Daraus ergibt sich eineSymmetrieforderung. Das hinzukommende Massepaket kann also aus denbaryzentrischen Koordinaten des in der vorhergehenden Verfeinerungsstufeneu generierten Massepakets erzeugt werden, indem man die Projektion inRichtung des zweiten, nicht in i origin stehenden Eckpunktes der Kan-te vornimmt. Die Information uber den zweiten Eckpunkt findet sich ini prvdir[1]. So wird in jeder Verfeinerungsstufe an jeder Kante jeweilsein Kantenteiler-Massepaket zu den bereits bestehenden hinzugefugt. Wel-ches Massepaket in der vorherigen Verfeinerungsstufe neu erzeugt wurde,ist uber die Variable l end auswertbar. Der Variablenname ruhrt daher, dasman sich die Erzeugung der neuen Massepakete wie die Erweiterung einesBaumgraphen vorstellen kann. Die zuletzt erstellten Massepakete sind amEnde eines Zweiges. Die Baumstruktur, die die Generierung aller Massepa-kete wiedergibt, ist in Abbildung 3.10 dargestellt.


Flachenteiler-Massepakete Sie treten erstmals ab der 3. Verfeinerungs-stufe an den vier Seitenflachen des Tetraeders auf. Ihre baryzentrischen Ko-ordinaten r coeff ergeben sich zu

λ(3)j =

413 , falls j ∈ k,m, n113 , sonst.

wobei Ek, Em, En die Eckpunkte der Flache sind, an dem das Flachenteiler-Massepaket liegt. Der Volumenanteil der vier Massepakete ist anfangs je-weils VF = 3

27 (r volume). Analog zu den bisherigen Massepaketen wird denFlachenteiler-Massepaketen ein Eckpunkt in i origin zugeordnet, in dessenRichtung sie mit jeder Erhohung der Verfeinerungsstufe projiziert werden.Wie bei den Kantenteiler-Massepaketen werden bei der Erhohung der Ver-feinerungsstufe neue Massepakete erzeugt, die von den in der vorherigenVerfeinerungsstufe neu erzeugten Massepaketen abgeleitet werden. Hier ge-staltet sich diese Generierung jedoch etwas komplizierter. Zur Erinnerung: injeder Verfeinerungsstufe kommen pro Seitenflache zu den zuletzt erzeugtena Massepaketen a + 1 weitere hinzu. Es werden also alle in der vorherigenVerfeinerungsstufe erzeugten Massepakete um ein weiteres erganzt. Bis aufein Massepaket M (i−1)

? , aus dem zwei neue Massepakete hervorgehen. Seimit Ek der Eckpunkt festgelegt, in dessen Richtung alle Massepakete einerSeitenflache bei der Erhohung der Verfeinerungsstufe von i−1 auf i projiziertwerden. Um die Koordinaten der neu hinzukommenden Schwerpunkte zu be-rechnen, wird die Zentralprojektion auf die Schwerpunkte der in der vorhe-rigen Verfeinerungsstufe i−1 erzeugten Massepakete in Richtung eines zwei-ten Eckpunktes Em der Seitenflache angewandt. (Welche Massepakete in derletzten Verfeinerungsstufe erzeugt wurden ist uber l end auswertbar.) Da-mit hatten wir alle Massenschwerpunkte der Verfeinerungsstufe i bis auf denvon M (i)

? generiert. Der Massenschwerpunkt fur M (i)? ergibt sich aus dem im

vorherigen Schritt erzeugten Massenschwerpunkt des Massenpaketes M (i−1)?

durch Projektion in Richtung des dritten Eckpunktes En der Seitenflache.Um die Information vorzuhalten, welches Massepaket zuletzt in Richtung Enprojiziert wurde und damit M (i−1)

? ist, wird in i prvdir[1], i prvdir[2]die Projektionsrichtungen der vorherigen beiden Verfeinerungsschritte ge-speichert. Fur die als erstes auftretenden Flachenteiler-Massepakete der 3.Verfeinerungsstufe werden die Variablen i prvdir[1], i prvdir[2] dabeiso vorbelegt, dass sie als M (3)

? ihrer Seitenflache fur den nachsten Verfeine-rungsschritt dienen.

Raumteiler-Massepakete Sie treten ab der 4. Verfeinerungsstufe auf. Da-bei ist es zunachst nur ein Massepaket, dessen Massenschwerpunkt

λ(4)j =

14, fur alle j


als baryzentrische Koordinaten r coeff besitzt. Der Volumenanteil des Mas-sepakets ist VR = 6

64 . Auch hier wird gemaß (3.7) der Volumenanteil r volumemit jeder Projektion, also Erhohung der Verfeinerungsstufe, geringer.Erhoht man die Verfeinerungsstufe werden alle bereits vorhandenen Masse-pakete Richtung ihres Eckpunktes i origin projiziert. Noch einen Schrittkomplizierter ist nun die Generierung neuer Massepakete bei der Erhohungder Verfeinerungsstufe. Dabei gibt es drei verschiedene Massepakete, wasderen ”Nachkommenschaft“ angeht: M (i−1),M

(i−1)? ,M

(i−1)?? .

Aus Massepaketen M (i−1) geht nur ein neues Massepaket hervor, das durchProjektion in Richtung eines zweiten Eckpunktes Ek des Tetraeders ent-steht. Aus Massepaketen M

(i−1)? gehen zwei Massepakete hervor, davon ist

eines vom Typ M (i) mit Projektion Richtung Ek, aus dem zweiten geht je-doch wieder eines vom Typ M

(i)? hervor und es wird Richtung des dritten

Eckpunktes Em des Tetraeders projiziert. Insofern haben wir noch Analo-gie zu den Flachenteiler-Massepaketen. Aus dem dritten Typ M (i−1)

?? jedochgehen drei Massepakete hervor, davon eines vom Typ M (i), das durch Pro-jektion Richtung Ek entsteht, eines vom Typ M

(i)? , das durch Projektion

Richtung des dritten Eckpunktes Em generiert wird und schließlich wirdpro Verfeinerungsstufe aus M (i−1)

?? ein neues Massepaket vom Typ M(i)?? er-

zeugt, das durch Projektion Richtung des vierten Eckpunktes En des Te-traeders bestimmt ist. Die Information uber den Typus der ”Nachkommen-schaft“ kann aus den Variablen i prvdir[1], i prvdir[2], i prvdir[3]fur jedes Raumteiler-Massepaket gewonnen werden. Das Raumteiler-Masse-paket, das in der Verfeinerungsstufe 4 als erstes erzeugt wird, ist dabei sozu wahlen, dass es vom M

(4)?? -Typ ist.

Zusammenfassend wird also so vorgegangen, um eine neue Verfeinerungs-stufe i zu generieren:

• Projiziere die baryzentrischen Koordinaten r coeff aller Massenschwer-punkte der Verfeinerungsstufe i− 1 mittels (3.6) in Richtung von demin i origin gespeicherten Eckpunkt und belege r volume gemaß (3.7)

• Die durch l end ausgezeichneten, also die in der Verfeinerungsstufei− 1 erzeugten Massepakete, werden je nach ihrem Typus, der durchi prvdir[1], i prvdir[2], i prvdir[3] bestimmt ist, auch in an-dere Richtungen projiziert. Damit werden die neuen Massepakete er-zeugt. Ihr Volumen r volume ergibt sich nach (3.7). Setze l end furdie neuen Massepakete auf wahr und fur alle anderen Massepakete auffalsch.

Dabei muss beachtet werden, dass Eckpunkt-Massepakete ab der 1. Verfei-nerungsstufe, Kantenteiler-Massepakete erst ab der 2. Verfeinerungsstufe,Flachenteiler-Massepakete ab der 3. und Raumteiler-Massepakete ab der 4.


Verfeinerungsstufe hinzukommen.

Jedes Massepaket ließe sich auch als Raumteiler-Massepaket interpretie-ren, wobei ein Eckpunkt-Massepaket das initiale M

(0)?? -Massepaket ware.

Dementsprechend ist der Baum in Abbildung 3.10 zusammenhangend dar-gestellt. Da jedoch die Massepakete an den Ecken, Kanten bzw. Flachenandere Korper darstellen und ihre Massenschwerpunkte deshalb nicht di-rekt voneinander abgeleitet werden konnen, ist die Klassifizierung nach denverschiedenen Massepaket-Typen sinnvoll. Eine Folgerung dessen ist, dassbeispielsweise die Fortsetzung des ganz rechten Zweiges in Abbildung 3.10ausgehend von dem Raumteiler-Massepaket R, wiederum der bereits darge-stellte Baum mit Ausgangsknoten E aus der initialen Verfeinerungsstufe 0selbst ist.Abschließend sei noch die Bedeutung der Variablen r mass[4] erlautert.In ihnen wird der Massenanteil eines Massepaketes gespeichert, den diesesMassepaket je Eckpunkt erhalt: r mass[i] = r coeff[i] ∗ r volume. DieWerte fur die Massenanteile pro Eckpunkt mussen wie die baryzentrischenKoordinaten des Massenschwerpunktes eines Massepaketes nur einmalig be-rechnet werden.

Man kann nun die Massepakete-Koeffiziententabelle mit dem so beschriebe-nen Vorgehen Verfeinerungsstufe um Verfeinerungsstufe sukzessive bis zumgewunschten Verfeinerungsgrad generieren. In Abbildung 3.11 ist der Ablaufdes oben dargestellten Algorithmus nochmals sehr komprimiert gezeigt. Inder darauf folgenden Abbildung 3.12 kann man das Ergebnis der Berechnungder baryzentrischen Koordinaten der Massenschwerpunkte im Referenzele-ment fur verschiedene Verfeinerungsstufen betrachten.


MPK-Tabelle dynamisch allokieren

Massepakete-Koeffizienten der Verfeinerungsstufe 1 erzeugen

Neue Verfeinerungsstufe anlegen: Zentralprojektion der alten Schwerpunkte der Massepakete in Richtung ihrer Herkunftsecke

Neue Massepakete-Koeffizienten an die Enden des Baumgraphen anhängen

Maximale Verfeinerungsstufe erzeugt?

Für alle Verfeinerungsstufen: Volumenkoeffizienten zuordnen

Für alle Verfeinerungsstufen: Massenkoeffizienten zuordnen

Initialisierung von flash3d Start: Initialisierung der Masse-pakete-Koeffiziententabelle

Optional: Für alle Verfeinerungs-stufen: Grafikausgabe der Massenschwerpunkte im normalisierten Tetraeder

Stop: Initialisierung der Masse-pakete-Koeffiziententabelle Weitere Initialisierung von flash3d

Verfeinerungsstufe 2, 3 oder 4 ?

V.stufe = 2: Erstes Auftreten von Kantenteilern; Kantenmassepakete erzeugen V.stufe = 3: Erstes Auftreten von Flächenteilern; Flächenmassepakete erzeugen V.stufe = 4: Erstes Auftreten von Raumteilern; Raummassepakete erzeugen

ja

nein

nein

ja

Abbildung 3.11: Schematische Darstellung des Algorithmus zur Erzeugung derMassepakete-Koeffiziententabelle als Flussdiagramm.


Abbildung 3.12: Verschiedene Verfeinerungsstufen i des Referenzelements mitden zugehorigen Massenschwerpunkten der Massepakete. Von links nach rechts:i = 1, 2, 3, 4, 5, 6, 11, 16, 21, 26, 31, 36


3.3.3 Algorithmus: Horizontale Suche

Aufgabe ist es, die mit Hilfe der Tetraedertabelle und der Massepakete-Koeffiziententabelle erzeugten Massepakete M1, . . . ,MN einem Tetraeder Tjim Gitter T n+1

h des nachsten Zeitschritts tn+1 zuzuordnen oder festzustellen,dass ein Massepaket Mi außerhalb von T n+1

h liegt, also von dem vorgege-benen Windfeld aus dem Simulationsgebiet advehiert wird. Dabei stellt Ndie Gesamtzahl der Massepakete, die aus den Tetraedern Ti ∈ T nh erzeugtwerden.Seien mit a1, a2, a3, a4 ∈ R3 die Koordinaten der Eckpunkte eines TetraedersT bezeichnet. Um zu prufen, ob ein Punkt p in einem Tetraeder T liegt, wirdfur jede durch drei Eckpunkte aufgespannte Seitenflache von T festgestellt,ob der Punkt p auf derselben Seite wie der vierte Eckpunkt von T liegt. Istdas der Fall liegt p in T , sonst nicht.Dies soll exemplarisch fur einen Fall dargestellt werden. Sei die Seitenflache,die durch a1, a2, a3 gebildet wird, betrachtet. Durch den Normalenvektor aufdie Seitenflache

n = (a2 − a1)× (a3 − a1)

wird eine Ebene beschrieben, die den R3 in zwei Halbraume teilt. Nun be-trachtet man das euklidische Skalarprodukt der Vektoren p−a1 bzw. a4−a1

mit n:

sp = (p− a1)Tnsa4 = (a4 − a1)Tn

Gilt nunsgn(sp) = sgn(sa4) (3.8)

so liegen p und a4 in demselben Halbraum bzw. fur sp = 0 liegt p in derEbene, die durch n beschrieben wird. Durch zyklisches Vertauschen der In-dizes der Eckpunkte fuhrt man obige Prufung fur alle Seitenflachen durch.Die Schnittmenge der durch den Normalenvektor und den jeweils viertenEckpunkt des Tetraeders festgelegten Halbraume ergibt das Innere von T .Ist also (3.8) wahr fur alle Seitenflachen von T (mit entsprechender An-passung des Index in (3.8)), liegt p in int(T ). Gilt fur eine oder mehrereSeitenflachen, dass sp = 0 ist, liegt p auf der Seitenflache, Kante bzw. einemEckpunkt von T und damit kann p ebenfalls dem Tetraeder T zugeordnetwerden.Wir konnen nun feststellen, ob ein Punkt p in einem Tetraeder liegt odernicht. Aber wie stellt man nun fest, in welchem Tetraeder der Punkt p liegt?Der bisher in amatos3d/flash3d verwendete Algorithmus beruht auf hier-archischer Suche. Hierbei wird von der grobsten Gitterstufe ausgehend, furjeden Tetraeder dieser Gitterstufe getestet, ob p innerhalb oder außerhalbeines Tetraeders ist. Sollte der Punkt in keinem Tetraeder liegen, ist manfertig. Er liegt außerhalb des Simulationsgebiets. Stellt man fest, er liegt


in einem Tetraeder, so muss man nun auf der nachsthoheren Gitterstufedieselbe Prozedur der Suche vollziehen. Nach und nach steigt man so Git-terstufe um Gitterstufe hoher, um schließlich bei einem Tetraeder zu enden,der nicht mehr verfeinert ist. Fur den Fall, dass die Verfeinerung der Tetra-eder mittels Bisektion (siehe Abschnitt 2.2) durchgefuhrt wird, ergibt sicheine minimale Zahl von n Prufungen nach obigen Schema und eine maxi-male Zahl von 6 + 2(n − 1) Prufungen, n sei die Zahl der Gitterstufe desTetraeders indem der Punkt p liegt. Auf der grobsten Gitterstufe sind beider hier verwendeten Friedrichs-Keller-Triangulierung des Einheitswurfels 6Tetraeder vorgegeben, auf den weiteren Gitterstufen existieren je 2 Tetra-eder.Die horizontale Suche verfahrt nach einem anderen Prinzip. Sei ein nichtmehr hoher verfeinertes Starttetraeder T im Gitter vorgegeben. Es wird nungepruft, ob p in T liegt. Falls ja, ist man bereits nach einer Prufung (3.8)fertig. Falls nein, untersuche den Nachbartetraeder von T , der in Richtungder Seitenflache liegt, die bei obiger Prufung ergab, dass p nicht in T liegt.Man wiederholt diese Prozedur solange, bis der gesuchte Tetraeder gefundenist bzw. kein Nachbartetraeder existiert, also p außerhalb des Simulationsge-biets liegt. Dieses Vorgehen erhalt den Namen ”horizontale Suche“, da mansich ausschließlich innerhalb der Tetraeder bewegt, die nicht mehr hoherverfeinert sind.Fur die Suche in welches Tetraeder der Massenschwerpunkt eines Massepa-kets advehiert wird, wahlt man als Starttetraeder das Tetraeder, aus dem dasMassepaket stammt. Nimmt man nun ein gutartiges Windfeld an, besteht dieAussicht nach weniger als n Schritten den gesuchten Tetraeder zu finden.Da die Suche bei der Massenpakete-Semi-Lagrange-Methode den großtenAnteil an der Rechenzeit einnimmt, ist also eine erhebliche Beschleunigungdes Verfahrens gegenuber der hierarchischen Suche zu erwarten. Sollte manjedoch bei der horizontalen Suche ein hoch verfeinertes Gebiet durchqueren,so schmilzt der Vorteil schnell dahin, denn mit jedem Sprung zum Nach-bartetraeder legt man nur eine kurze Strecke in Richtung des gesuchtenTetraeders zuruck. Das kann sich in einer erheblich hoheren Suchzeit nie-derschlagen. Eine scharfe Obergrenze von Suchoperationen lasst sich nichtangeben.

Die Simulationsergebnisse aus Kapitel 4 zeigen, dass zumindest fur dieseTestfalle die gewunschte Beschleunigung eintritt. Der implementierte Algo-rithmus ist schematisch in Abbildung 3.13 dargestellt.


Neue Werte an den Knoten berechnen: Start: Horizontale Suche der Massepakete im neuen Gitter

Für alle Massepakete

Massepaket bereits gefunden?

Für alle Seitenflächen des Tetraeders im neuen Gitter: Schwerpunkt des Massepakets auf derselben Seite wie der vierte Eckpunkt des Tetraeders?

Herkunftstetraeder als Starttetraeder der Suche verwenden

Nehme benachbartes Tetraeder in Richtung der Seitenfläche, die ergab, dass der Schwerpunkt nicht im Tetraeder ist

nein Ist das Nachbartetraeder außerhalb des Gitters?

Massepaket außerhalb des Gitters

Massenschwerpunkt im Tetraeder. Index des Tetraeders speichern

Alle Schwerpunkte der Massepakete gefunden?

nein

nein

Stop: Horizontale Suche der Massepakete im neuen Gitter Neue Werte an den Knoten berechnen

ja

ja

ja

nein

ja

Abbildung 3.13: Flussdiagramm zur Suche der Massepakete in den Tetraederndes stromaufwarts verfolgten Gitters T n+1

h

4 ERGEBNISSE: TESTS ANHAND VON MODELLPROBLEMEN 51

4 Ergebnisse: Tests anhand von Modellproblemen

Die hier entwickelte Massepakete-Semi-Lagrange-Methode (MPSLM) wirdin diesem Kapitel mit der bisherigen Semi-Lagrange-Methode (SLM) ver-glichen. In Abschnitt 4.1 wird die fur alle Modellprobleme gleiche Anfangs-bedingung dargestellt und einige Vorbemerkungen gegeben. Darauf folgendwerden die Modellprobleme

• Lineares Windfeld (4.2),

• Zirkulares Windfeld (4.3),

• Konvergentes Windfeld (4.4)

behandelt. Es werden die Faktoren

• Fehler in der L2-Norm,

• Massenerhaltung,

• Konvergenzverhalten

untersucht. Abschließend werden in Abschnitt 4.5 einige Schlussfolgerungenund zusammenfassende Bemerkungen zu den Testergebnissen dargestellt.

Hinweis: Fur die Durchfuhrung der Testsimulationen wird ein Rechner derFirma SGI mit vier RISC-Prozessoren MIPS R10000 CPU, R10010 FPU, 225MHz, 1024 MB RAM verwendet. Der Quellcode ist in Fortran 90 geschrie-ben. Als Compiler wird MIPSpro 7 Fortran 90 eingesetzt. Auf der beiliegen-den CD-ROM befinden sich in den Unterverzeichnissen /SLM bzw. /MPSLM dieQuellcodes von amatos3d/flash3d je nach eingesetztem Verfahren. Jeweils inden Unterverzeichnissen /lin, /circ, /conv ist das entsprechende Wind-feld vorgegeben. Außerdem enthalten die Verzeichnisse auch kompilierte Pro-grammdateien, die auf dem SGI-Rechner ausfuhrbar sind. In /SLM/resultsbzw. /MPSLM/results sind die Ergebnisdateien aller Testlaufe abgelegt, dieunter anderem den L2-Fehler und die relative Masse als Daten enthalten. In/SLM/animations bzw. /MPSLM/animations schließlich findet man animier-te GIF-Dateien von ausgewahlten Testlaufen jeweils auf einem uniformenund adaptiven Gitter fur alle drei Testbeispiele.

4.1 Darstellung der Testumgebung

Als Gebiet Ω dient der Einheitswurfel E = [0, 1]3 (vgl. 2.1). Es werdenAnfangswerte fur u wie folgt vorgegeben:

u(x, 0) =

1 x ∈ B0,125

((0,25 0,5 0,5)T

),

0 sonst.


Die durch die Anfangswerte definierte Kugel wird Tracer genannt.Das Gebiet wird durch finite Elemente diskretisiert. Je nach verwendetemGitter ergeben sich die minimale Gitterweite hmin bzw. die maximale Git-terweite hmax, die in Tabelle 4.1 aufgefuhrt sind.Alle Modellprobleme werden auf drei uniformen Gittern und funf adapti-ven Gittern getestet. Alle Knoten Ki, die sich in B0,125

((0,25 0,5 0,5)T

)befinden, erhalten den Anfangswert u(Ki, 0) = 1, die restlichen den An-fangswert u(Kj , 0) = 0.

Abbildung 4.1: Darstellung der Anfangswerte auf einem adaptiven Gitter. Schnittbei z = 0,5. Im Bild ist der Tracer erkennbar. Im Ubergangsbereich zwischen Tracer(u = 1) und dem sonstigen Gebiet (u = 0) ist die erhohte Knotendichte sichtbar.

Ein Beispiel fur eine Anfangstriangulierung mit Knotenwerten ist in Abbil-dung 4.1 zu finden: Bei adaptiven Gittern wird die Anfangstriangulierungder gegebenen Anfangsbedingung angepasst, indem in einer Umgebung desTracerrandes vor dem ersten Zeitschritt nachverfeinert wird.Bei der Wahl des Windfeldes wird bei den Modellproblemen auf analytischeLosbarkeit geachtet, um den Fehler in der L2-Norm ‖u− uh‖ explizit ange-ben zu konnen.Um die Verfahren auf Massenerhaltung zu testen, wird nach der Initialisie-rung der Gitterwerte die vorhandene Anfangsmasse in der Triangulierung


VS 13 13 14 14 15 15 8 15 9 16 10 17 11 18 12 19hmin 0,0625 0,0541 0,0313 0,0313 0.0313 0,0271 0,0156 0,0156hmax 0,1083 0,0884 0,0625 0,3536 0,2500 0,2165 0,1768 0,1250

Tabelle 4.1: Ubersicht uber minimale und maximale Kantenlange hmin bzw. hmax

fur verschieden gewahlte Verfeinerungsstufen VS des Gitters. Mit der Abkurzung8 15 ist ein Gitter charakterisiert mit der grobsten Verfeinerungsstufe 8 und der fein-sten Verfeinerungsstufe 15. 13 13, 14 14, 15 15 stellen uniforme Gitter dar, wahrend8 15, . . . , 12 19 adaptive Gitter sind.

Th von Ω bestimmt und zu Vergleichszwecken gespeichert.

m(k) =∑Ki∈Th

14u(Ki)V (Ki) k = 0, . . . , N (4.1)

Hierbei ist N die Anzahl der Zeitschritte. In den weiteren Zeitschritten k =1, . . . , N wird mit der Formel (4.1) die vorhandene Masse berechnet und dierelative Masse

m(k)rel =

m(k)

m(0)fur alle k = 1, . . . , N

gebildet. Fur alle k sollte wegen der Forderung des Massenerhalts m(k)rel = 1

gelten.

Zu jedem Modellproblem werden außerdem Konvergenzbetrachtungen durch-gefuhrt. Hierbei wird zwischen den Fallen SLM bzw. MPSLM und uniformenbzw. adaptiven Gitter unterschieden. Fur uniforme Gitter werden die Fehlerin der L2-Norm im 13 13- mit dem 15 15-Gitter und fur adaptive Gitter das9 16- mit dem 12 19-Gitter in Bezug gesetzt, um die Konvergenzrate r zuschatzen:

r(k) =log E

(k)g

E(k)f

log hgmin

hfmin

fur alle Zeitschritte k = 0, . . . , N (4.2)

Dabei stellt E(k)g bzw. E(k)

f den auftretenden Fehler im groberen bzw. fei-

neren Gitter zum Zeitschritt k dar. hgmin bzw. hfmin bezieht sich auf dieKantenlange der feinsten Gitterstufe im jeweils groberen bzw. feineren Git-ter. Die Werte r(k) aus (4.2) werden uber (fast) alle berechneten Zeitschrittesummiert und gemittelt:

r =1k

∑k

r(k)


So erhalt man eine empirisch ermittelte Konvergenzrate des Verfahrens furdas jeweilige Modellproblem. (Bei der Konvergenzbetrachtung werden in al-len Fallen die ersten zwei Zeitschritte nicht mit einbezogen, da sich hier zumTeil unbrauchbar hohe negative oder positive Konvergenzraten ergeben ha-ben. Sie hatten die Werte der gemittelten Konvergenzrate verfalscht.)

4.2 Modellproblem I: Lineares Windfeld

Das erste Modellproblem ist ein konstantes, lineares Windfeld parallel zurx-Achse:

a(x, t) =

c00

Eine Darstellung des Windfelds mit Geschwindigkeitsvektoren sieht man inAbbildung 4.2. Dieses Windfeld gewahrleistet, dass der Tracer in seiner Formund Große erhalten bleibt. Das erleichtert die Bestimmung des Fehlers in derL2-Norm, da die analytische Losung u an den Knoten einfach zu berechnenist.

Zunachst wird nun auf die numerischen Ergebnisse in Bezug auf den Fehlerin der L2-Norm eingegangen. Dabei stellt die Abbildung 4.3 die fur die SLM(oben) bzw. MPSLM (unten) erhaltenen Werte dar. Man erkennt, dass dieFehlerentwicklung ahnlich verlauft. Allerdings ist der Fehler im MPSLM-Verfahren leicht hoher als der des SLM-Verfahrens.Die Sprunge im Graphen der L2-Norm finden ihren Ursprung in der Berech-nung der L2-Norm. Die analytische Losung u ist mit der Wahl des Windfeldsa bekannt:

u(x, t) =

1 x ∈ B0,125

((0,25 + ct 0,5 0,5)T

),

0 sonst.

Fur die Berechnung des L2-Fehlers werden die Werte der analytischen undder numerischen Losung an den Gitterknoten herangezogen. Es hangt vomZeitschritt, dem Windfeld und dem Verfeinerungsgrad der Triangulierung Thab, ob einige Gitterknoten innerhalb oder außerhalb des Radius des Tracersder analytischen Losung liegen. Das bedeutet, dass der Fehler ‖u− uh‖2 jenach Lage des Tracers schwankt. Die Abhangigkeit vom Gitter erkennt mangut in den Graphen zum L2-Fehler (Abbildung 4.3). Mit zunehmend hohererVerfeinerung, also hoherer Knotendichte in einer Umgebung um den Randdes Tracers, fallen die Schwankungen geringer aus.Bemerkenswert ist, dass die Fehlerkurve des 8 15- bzw. 15 15-Gitters nahezuzusammenfallt und damit die Starke der adaptiven Verfeinerung nochmalsinnfallig macht.


Außerdem erkennt man, dass mit hoherer Verfeinerung des Gitters auchder Fehler in der L2-Norm geringer ausfallt. Im Zusammenhang mit derKonvergenzeigenschaft des SLM-Verfahrens ist dies auch zu erwarten.

Im oberen Diagramm der Abbildung 4.4 wird das Verhalten der SLM bzgl.Massenerhaltung gezeigt. Man erkennt sofort die Abweichung der Massevom geforderten Wert mrel = 1. Bezogen auf die Verfeinerungsstufen ist kei-ne Systematik zu erkennen. Die Werte fur die Massenerhaltung der MPSLMin Abbildung 4.4 zeigen dagegen, dass das Ziel der Massenerhaltung ein-gehalten wird. Betrachtet man den Graphen genau, sieht man ab dem 60.Zeitschritt, dass fur das 13 13-Gitter sehr geringe Massenverluste auftreten.Dies ist dem Phanomen der numerischen Diffusion zuzuschreiben. Die nume-rische Diffusion verursacht, dass nach einigen Zeitschritten geringe Spurenvon Masse durch numerische Fehler bis an den Rand des Gebiets Ω trans-portiert wurden. Das ist nicht gewollt, aber unvermeidbar. Die Massepaketenahe dem Rand ∂Ω werden bei entsprechendem Wind uber den Rand desGebiets hinaus getragen. Fur feinere Gitterstufen ergibt sich dieser Verlustin spateren Zeitschritten, da geringere numerische Diffusion auftritt. Dienumerische Diffusion laßt sich sehr gut in den Animationen auf der beilie-genden CD-ROM beobachten. Sie außert sich im Auseinandertreiben desTracers. Den Verlust uber den Rand des Gebiet kann man nicht erkennen,da die Konzentration der Masse am Rand nicht hinreichend groß ist.

Lineares Windfeld Uniformes Gitter Adaptives GitterSLM 0,98 1,24MPSLM 0,84 1,21

Tabelle 4.2: Konvergenzrate im Vergleich der Verfahren und der Gittersteuerungfur das lineare Windfeld

Die Konvergenzrate der beiden Verfahren, die in Tabelle 4.2 aufgefuhrt sind,sind in etwa vergleichbar und entsprechen jenen Erwartungen, die sich ausder Verwendung von linearen finiten Elementen ergeben.


Abbildung 4.2: Darstellung des linearen Windfeldes. Schnitt bei z = 0,5. Im Bildist das Tracerfeld mit seinen Konzentrationswerten erkennbar.


0,00E+00

1,00E-03

2,00E-03

3,00E-03

4,00E-03

5,00E-03

6,00E-03

7,00E-03

0 2 4 6 8 10 12 14 16 18 20 22 24 26 28 30 32 34 36 38 40 42 44 46 48 50 52 54 56 58 60 62 64 66 68 70

Zeitschritt

L2-N

orm

Feh

ler

13_1314_1415_158_159_1610_1711_1812_19

0,00E+00

1,00E-03

2,00E-03

3,00E-03

4,00E-03

5,00E-03

6,00E-03

7,00E-03

0 2 4 6 8 10 12 14 16 18 20 22 24 26 28 30 32 34 36 38 40 42 44 46 48 50 52 54 56 58 60 62 64 66 68 70

Zeitschritt

L2-N

orm

Feh

ler

13_1314_1415_158_159_1610_1711_1812_19

Abbildung 4.3: Fehler ‖u− uh‖2 fur das SLM- (oben) und das MPSLM-Verfahren(unten) angewandt auf das Modellproblem mit linearem Windfeld. In der Legendesind jeweils die grobste und feinste Gitterstufe angegeben.


0,0000E+00

2,0000E-01

4,0000E-01

6,0000E-01

8,0000E-01

1,0000E+00

1,2000E+00

1,4000E+00

0 2 4 6 8 10 12 14 16 18 20 22 24 26 28 30 32 34 36 38 40 42 44 46 48 50 52 54 56 58 60 62 64 66 68 70

Zeitschritt

Rel

ativ

e M

asse

13_1314_1415_158_159_1610_1711_1812_19

8,0000E-01

9,0000E-01

1,0000E+00

1,1000E+00

1,2000E+00

0 2 4 6 8 10 12 14 16 18 20 22 24 26 28 30 32 34 36 38 40 42 44 46 48 50 52 54 56 58 60 62 64 66 68 70

Zeitschritt

Rel

ativ

e M

asse

13_1314_1415_158_159_1610_1711_1812_19

Abbildung 4.4: Betrachtung der Massenerhaltung fur das SLM- (oben) und dasMPSLM-Verfahren (unten) angewandt auf das Modellproblem mit linearem Wind-feld. In der Legende sind jeweils die grobste und feinste Gitterstufe angegeben.


4.3 Modellproblem II: Zirkulares Windfeld

Das zirkulare Windfeld wird durch

a(x, t) =

−(x2 − 0,5)(x1 − 0,5)

0

beschrieben (siehe Abbildung 4.5). Es verursacht ein Kreisen um die durch(0,5 0,5 0,5)T verlaufende Achse parallel zur z-Achse. Das Windfeld ist sogewahlt, dass die Winkelgeschwindigkeit um das Zentrum des Windfeldskonstant ist. Dadurch ergibt sich ein Starkerwerden des Windfelds in x di-rekt proportional zum Abstand von x zum Zentrum. Der Tracer in seinerForm als Ball bleibt in der initialisierten Große erhalten. Folglich ist dieKonzentration im Ball konstant. Bei der verwendeten Programmkonfigura-tion kreist der Ball in 96 Zeitschritten einmal um die Achse und kehrt zumAusgangspunkt zuruck.Wie auch beim linearen Windfeld zeigt sich in Abbildung 4.6 das wiederkeh-rende Muster in der Schwankung des L2-Fehlers. Das Muster wiederholt sichjeweils nach einer Viertelumdrehung des Tracers, also nach 24 Zeitschritten.Analog zum linearen Windfeld zeigt sich bezogen auf den Graphenverlauf,dass die MPSLM mit der SLM vergleichbar ist. Allerdings sind die Fehlerdes MPSLM-Verfahrens bei vergleichbaren Gittern z.T. um das 1,75-fachehoher als die des SLM-Verfahrens. Das liegt in dem großeren Fehler, der beider Approximation der Werte uh durch die Massepakete entsteht.Wie beim linearen Windfeld ergibt sich keine Massenerhaltung fur das

SLM-Verfahren. Das MPSLM-Verfahren dagegen ist grundsatzlich massen-erhaltend, verliert jedoch durch numerische Diffusion nach ca. 48 Zeitschrit-ten beim grobsten Gitter signifikant an Masse. Feinere Gitter verlieren erstspater Masse uber den Rand von Ω.

Zirkulares Windfeld Uniformes Gitter Adaptives GitterSLM 1,26 1,17MPSLM 0,45 1,11

Tabelle 4.3: Konvergenzrate im Vergleich der Verfahren und der Gittersteuerungfur das zirkulare Windfeld

Auch in Tabelle 4.3 zeigt sich, dass die Konvergenzraten fur das Modell-problem auf adaptivem Gitter gleichauf liegen. Fur das uniforme Gitter istjedoch ein großer Unterschied von mehr als einer halben Ordnung festzu-stellen. Der Grund hierfur ist bisher nicht gefunden und bedurfte noch einernaheren Untersuchung.


Abbildung 4.5: Darstellung des zirkularen Windfeldes. Schnitt bei z = 0,5. ImBild ist das Tracerfeld mit seinen Konzentrationswerten erkennbar.


0,00E+00

1,00E-03

2,00E-03

3,00E-03

4,00E-03

5,00E-03

6,00E-03

7,00E-03

8,00E-03

9,00E-03

1,00E-02

0 8 16 24 32 40 48 56 64 72 80 88 96 104 112 120 128 136 144 152 160 168 176 184 192

Zeitschritt

L2-N

orm

Feh

ler

13_13

14_14

15_15

8_15

9_16

10_17

11_18

12_19

0,00E+00

1,00E-03

2,00E-03

3,00E-03

4,00E-03

5,00E-03

6,00E-03

7,00E-03

8,00E-03

9,00E-03

1,00E-02

0 8 16 24 32 40 48 56 64 72 80 88 96 104 112 120 128 136 144 152 160 168 176 184 192

Zeitschritt

L2-N

orm

Feh

ler

13_1314_1415_158_159_1610_1711_1812_19

Abbildung 4.6: Fehler ‖u− uh‖2 fur das SLM- (oben) und das MPSLM-Verfahren(unten) angewandt auf das Modellproblem mit zirkularem Windfeld. In der Legendesind jeweils die grobste und feinste Gitterstufe angegeben.


8,0000E-01

9,0000E-01

1,0000E+00

1,1000E+00

1,2000E+00

1,3000E+00

1,4000E+00

1,5000E+00

1,6000E+00

1,7000E+00

0 8 16 24 32 40 48 56 64 72 80 88 96 104 112 120 128 136 144 152 160 168 176 184 192

Zeitschritt

Rel

ativ

e M

asse

13_1314_1415_158_159_1610_1711_1812_19

8,0000E-01

9,0000E-01

1,0000E+00

1,1000E+00

1,2000E+00

0 8 16 24 32 40 48 56 64 72 80 88 96 104 112 120 128 136 144 152 160 168 176 184 192

Zeitschritt

Rel

ativ

e M

asse

13_1314_1415_158_159_1610_1711_1812_19

Abbildung 4.7: Betrachtung der Massenerhaltung fur das SLM- (oben) unddas MPSLM-Verfahren (unten) angewandt auf das Modellproblem mit zirkularemWindfeld. In der Legende sind jeweils die grobste und feinste Gitterstufe angegeben.


4.4 Modellproblem III: Konvergentes Windfeld

Das konvergente Windfeld

a(x, t) =

(x1 − 1,25)(x2 − 0,5)(x3 − 0,5)

treibt den Tracer dem Punkt (1,25 0,5 0,5)T zu. Dabei wird der Tracer zu-nehmend verdichtet. Auch hier wurde das Windfeld so gewahlt, dass derBall in seiner Form als Kugel erhalten bleibt, um damit die Berechnungder analytischen Losung zu erleichtern. Das kann man an einem Beispielim R

1 verdeutlichen: Sei K der Konvergenzpunkt auf den die Stomung zu-fließt. Seien A1, M1 und I1 drei Punkte mit der Beziehung A1 = M1 + r1

bzw. I1 = M1 − r1. r1 stellt den Radius dar. K sei so platziert, dass I1 amnachsten liegt. M1 ist der Mittelpunkt der betrachteten Kugel in R1. A1

liegt also außerhalb des Intervalls (M1,K) und I1 innerhalb. Analog zumoberen Windfeld ergibt sich nach einem Zeitschritt:

I2 = M1 − r1 + (M1 − r1 −K)∆t = M1(1 + ∆t)− r1(1 + ∆t)−K∆tM2 = M1 + (M1 −K)∆t = M1(1 + ∆t)−K∆tA2 = M1 + r1 + (M1 + r1 −K)∆t = M1(1 + ∆t) + r1(1 + ∆t)−K∆t

Sollte die Form der Kugel erhalten sein, musste r2 = A2−M2 und r2 = M2−I2 ubereinstimmen. Durch Einsetzen erhalt man: r2 = r2 = r2 = r1(1 + ∆t).Dieses Ergebnis ubertragt sich auch auf das verwendete Windfeld in R3.Wie bei den vorhergehenden Modellproblemen ist der L2-Fehler der bei-

den Verfahren vom Verlauf her vergleichbar (siehe Abbildung 4.9). Ebensoergibt sich ein hoherer Fehler fur die MPSLM. Auch das Kriterium derMassenerhaltung ergibt ein ahnliches Bild wie bei den vorhergehenden Test-problemen (siehe Abbildung 4.10). Die Semi-Lagrange-Methode weicht vonder Ausgangsmasse ab. Die Massepakete-Semi-Lagrange-Methode halt dasPrinzip der Massenerhaltung ein.

Konvergentes Windfeld Uniformes Gitter Adaptives GitterSLM 0,69 0,93MPSLM 0,70 0,90

Tabelle 4.4: Konvergenzrate im Vergleich der Verfahren und der Gittersteuerungfur das konvergente Windfeld

Die Konvergenzraten der beiden Verfahren sind insgesamt niedriger als beiden vorher betrachteten Modellproblemen. Sie stimmen fur SLM und MPSLMnahezu uberein.


Abbildung 4.8: Darstellung des konvergenten Windfeldes. Schnitt bei z = 0,5.ImBild ist das Tracerfeld mit seinen Konzentrationswerten erkennbar.


0,00E+00

5,00E-03

1,00E-02

1,50E-02

2,00E-02

2,50E-02

3,00E-02

3,50E-02

0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18

Zeitschritt

L2-N

orm

Feh

ler

13_1314_1415_158_159_1610_1711_1812_19

0,00E+00

5,00E-03

1,00E-02

1,50E-02

2,00E-02

2,50E-02

3,00E-02

3,50E-02

0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18

Zeitschritt

L2-N

orm

Feh

ler

13_1314_1415_158_159_1610_1711_1812_19

Abbildung 4.9: Fehler ‖u− uh‖2 fur das SLM- (oben) und das MPSLM-Verfahren(unten) angewandt auf das Modellproblem mit konvergentem Windfeld. In der Le-gende sind jeweils die grobste und feinste Gitterstufe angegeben.


8,0000E-01

9,0000E-01

1,0000E+00

1,1000E+00

1,2000E+00

1,3000E+00

1,4000E+00

1,5000E+00

1,6000E+00

1,7000E+00

0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18

Zeitschritt

Rel

ativ

e M

asse

13_1314_1415_158_159_1610_1711_1812_19

8,0000E-01

9,0000E-01

1,0000E+00

1,1000E+00

1,2000E+00

0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18

Zeitschritt

Rel

ativ

e M

asse

13_1314_1415_158_159_1610_1711_1812_19

Abbildung 4.10: Betrachtung der Massenerhaltung fur das SLM- (oben) und dasMPSLM-Verfahren (unten) angewandt auf das Modellproblem mit konvergentemWindfeld. In der Legende sind jeweils die grobste und feinste Gitterstufe angege-ben. Es wurden drei Tests auf uniformem Gitter und funf auf adaptiven Gitterdurchgefuhrt.


4.5 Fazit

Experimentell wurde nachgewiesen, dass das MPSLM-Verfahren bis auf nu-merische Diffusion Massenerhaltung garantiert. Wie man auch an den Ani-mationen in den Unterverzeichnissen /MPSLM/animations bzw./SLM/animations auf der beiliegenden CD erkennen kann, ist die nume-rische Diffusion bei der MPSLM geringfugig hoher als die der SLM.Der Fehler im Sinne der L2-Norm ist durchgehend hoher als der des bis-herigen Semi-Lagrange-Verfahrens. Der Fehler liegt je nach betrachtetemModellproblem um das 1,25- bis 1,75-fache hoher. Der Verlauf des L2-Fehler-graphen beider Methoden ist ahnlich (siehe Abbildungen 4.3, 4.6 und 4.9).Bis auf das zirkulare Modellproblem auf uniformen Gitter sind die Konver-genzraten fur das MPSLM-Verfahren und das SLM-Verfahren gleichwertig.Dies lasst aus experimenteller Sicht den Schluss zu, dass die Massepakete-Semi-Lagrange-Methode die wesentlichen Eigenschaften der Semi-Lagrange-Methode beibehalten konnte und zusatzlich die Massenerhaltung erreichtwird. Je nach der Anzahl der verwendeten Massepakete zur Diskretisierungder Konzentration, ist dies jedoch mit einer hoheren Laufzeit im Vergleichzur SLM verbunden.

5 ZUSAMMENFASSUNG UND AUSBLICK 68

5 Zusammenfassung und Ausblick

Aufgabe dieser Diplomarbeit war, eine massenerhaltende Semi-Lagrange-Methode zu entwickeln und diese in das Programmpaket amatos3d/flash3dvon Jorn Behrens einzubinden.Zunachst wurde in Kapitel 2 begrundet, warum die Semi-Lagrange-Methodenicht massenerhaltend sein kann. Als Ursache wurde die verwendete Formder Interpolation, die auf Konzentrationswerten basiert, festgestellt. In Ka-pitel 3 wurde die Massepakete-Semi-Lagrange-Methode eingefuhrt, die dasPrinzip der Trajekorienverfolgung beibehalt, die konzentrationsbasierte In-terpolation der SLM aber ersetzt. Grundlegender Gedanke der MPSLMist, die Konzentrationsverteilung durch Massepakete zu approximieren. DieMPSLM wurde eingehend in Kapitel 3 behandelt und ihre Implementie-rung und Einbindung in amatos3d/flash3d dargestellt. Es wurde bewiesen,dass die Massepakete-Semi-Lagrange-Methode durch Konstruktion bis aufMassenverluste uber den Rand des Simulationsgebiets massenerhaltend ist.Diese Eigenschaft wurde durch die Ergebnisse der Testsimulationen in Ka-pitel 4 bestatigt. Die weiteren Ergebnisse des Kapitels 4 zeigten, dass derFehler in der L2-Norm der MPSLM zwar hoher als der Fehler der SLMist, jedoch einen ahnlichen Verlauf wie die SLM hat. Die SLM und dieMPSLM erscheinen deshalb vergleichbar. Die Aufgabe, eine massenerhal-tende Semi-Lagrange-Methode zu entwickeln und in das Programmpaketamatos3d/flash3d von Behrens einzubinden, wurde also erfullt.

Eine weitergehende Aufgabe ware nun die theoretische Fundierung bzgl.Konvergenz und Stabilitat der MPSLM. Die Ergebnisse aus den Testsimu-lationen in Kapitel 4 erscheinen diesbezuglich ermutigend.

Um Moglichkeiten zu finden, die Effizienz der MPSLM zu erhohen, wareaußerdem interessant zu untersuchen, inwiefern man durch Anderungen amVerfahren die Anzahl der Massepakete reduzieren kann, ohne an Genauig-keit zu verlieren. Hierbei ware eine nahere Betrachtung der Wahl der Ver-feinerungsstufe mittels des in Abschnitt 3.2.3 erorterten a-priori-Schatzershilfreich. Ebenso erscheint dem Autor noch die nahere Betrachtung der be-reits implementierten adaptiven Nachverfeinerung der Tetraeder in Mas-sepakete attraktiv. Eine andere Moglichkeit konnte die Modifizierung derUnterteilung der Tetraeder in Massepakete sein, die bisher durch die nahezuaquidistante Platzierung der Massepakete gekennzeichnet ist. Eine konzen-trationsbezogene Platzierung der Massepakete eroffnet unter Umstanden dieMoglichkeit, weniger Massepakete erzeugen zu mussen als bei der herkomm-lichen MPSLM.

Die hier vorgestellte MPSLM ist fur Simulationen auf einem dreidimensio-nales Gitter entwickelt. Eine Ubertragung auf zweidimensionale Gitter ist

5 ZUSAMMENFASSUNG UND AUSBLICK 69

leicht moglich. Dies wurde auch die Beschrankung des Verfahrens hinsicht-lich des Speicherbedarfs wesentlich entscharfen, da die Anzahl der Massepa-kete pro Element der Triangulierung nicht mehr kubisch, sondern nur nochquadratisch anwachst. Es ergibt sich somit auch eine schnelle Ausfuhrungder Simulation in 2D. Bei starkerer Verfeinerung steht eine hohere Genau-igkeit in Aussicht.

Interessant erscheint dem Autor auch, die Genauigkeit der MPSLM auf drei-dimensionalen Gittern erhohen. Eine Moglichkeit, dies zu erreichen, ware dieAnwendung finiter Elemente hoherer Ordnung, gerade auch im Hinblick derErgebnisse von Falcone und Ferretti in [FF98], die in Abschnitt 2.4 erlautertwurden. Die MPSLM ist jedoch unter der Voraussetzung linearer finiter Ele-mente entwickelt worden. Die Erweiterung auf finite Elemente hoherer Ord-nung ware also noch zu untersuchen.

NOTATION 70

Notation

u Losung der partiellen Differentialgleichunguh Losung der diskretisierten partiellen Differentialgleichung

Ω Grundgebiet bezuglich der raumlichen Variablen x

Γ = ∂Ω Rand von ΩV BanachraumVh endlichdimensionaler Finite-Elemente-Raum

L2(Ω) Raum der integrierbaren Funktionen uber Ω, fur die gilt:∫Ω

|f |2 dx <∞

H1(Ω) Sobolev-Raum aller Funktionen in L2(Ω),deren erste Ableitung ebenfalls in L2(Ω) liegt.

a(x, t) Vektorfeld mit Windgeschwindigkeiten in x ∈ Ω ⊂ Rd, d = 2, 3und t ∈ [0, T ]

x Raumliche Variable x ∈ Ω ⊂ Rd, d = 2, 3xi Komponente i der raumlichen Variable xt Zeitliche Variable t ∈ [t0, T ]

tn, tn+1 Alter bzw. neuer Zeitschritt

m(·) MasseV (·) Volumen

Th Triangulierung des Gebiets ΩT nh , T n+1

h Triangulierung des Gebiets Ω zum alten bzw. neuen Zeitschritt.Im allgemeinen gilt: T nh 6= T n+1

h

Ti Finites Element i der Triangulierung ThEki Ecke k des Finiten Elements TiKi Knoten i ∈ ThM li Massepaket l, das bei der Unterteilung von Ti entsteht

λlk Baryzentrische Koordinaten eines Massepaketes lbezogen auf die Eckpunkte Eki eines Tetraeders Ti ∈ T nh

µlk Baryzentrische Koordinaten eines Massepaketes lbezogen auf die Eckpunkte Ekj eines Tetraeders Tj ∈ T n+1

h

LITERATUR 71

Literatur

[Ban91] Bansch, Eberhard: Local Mesh Refinement in 2 and 3 Dimen-sions. IMPACT Comput. Sci. Engrg., 3:181–191, 1991.

[Beh96] Behrens, Jorn: Adaptive Semi-Lagrange-Finite-Elemente-Methode zur Losung der Flachwassergleichungen: Implementie-rung und Parallelisierung. Nummer 217 in Berichte zur Polar-forschung. Alfred-Wegener-Institut, Bremerhaven, 1996.

[BFL99] Bern, Marshall W., Joseph E. Flaherty und Mitchell

Luskin: Grid Generation and Adaptive Algorithms. Springer-Verlag, New York, 2. Auflage, 1999.

[BR87] Bank, Randolph E. und Donald J. Rose: Some Error Esti-mates for the Box Method. SIAM J. Numer. Anal., 24(4):777–787,1987.

[Bra97] Braess, Dietrich: Finite Elemente. Springer-Verlag, Berlin,Heidelberg, New York, 2. Auflage, 1997.

[Ede01] Edelsbrunner, Herbert: Geometry and Topology for MeshGeneration. Cambridge University Press, Cambridge, 2001.

[FF98] Falcone, Maurizio und Roberto Ferretti: ConvergenceAnalysis for a Class of high-order semi-Lagrangian AdvectionSchemes. SIAM J. Numer. Anal., 35(3):909–940, 1998.

[FG00] Frey, Pascal Jean und Paul-Louis George: Mesh Genera-tion. Hermes Science, Oxford, 2. Auflage, 2000.

[GR94] Großmann, Chrisitan und Hans-Gorg Roos: Numerik par-tieller Differentialgleichungen. Teubner-Verlag, Stuttgart, 2. Auf-lage, 1994.

[Hac88] Hackbusch, W.: On First and Second Order Box Schemes.Computing, 41:277–296, 1988.

[HMS56] Hammer, P.C., O.P. Marlowe und A.H. Stroud: Numericalintegration over simplexes and cones. Math. Tables Aids Comp.,10(55):130–137, Juli 1956.

[KA00] Knabner, Peter und Lutz Angermann: Numerik partiellerDifferentialgleichungen. Springer-Verlag, Berlin Heidelberg NewYork, 2000.

[LeV94] LeVeque, Randall J.: Numerical Methods for ConservationLaws. Birkhauser-Verlag, Basel, 1994.

LITERATUR 72

[LL91] Landau, Lev D. und E. M. Lifschitz: Lehrbuch der theorti-schen Physik – Band VI Hydrodynamik. Akademie Verlag, Berlin,5. Auflage, 1991.

[NM02] Nair, Ramachandran N. und Bennert Machenhauer: TheMass-Conservative Cell-Integrated Semi-Lagrangian AdvectionScheme on the Sphere. Mon. Wea. Rev., 130:649–667, Marz 2002.

[OBSC00] Okabe, Atsuyuki, Barry Boots, Kokichi Sugihara undSung Nok Chiu: Spatial Tessellations: Concepts and Applicati-ons of Voronoi Diagrams. Wiley Series in Probability and Stati-stics. John Wiley & Sons Ltd., Chichester, UK, 2. Auflage, 2000.

[QV94] Quarteroni, Alfio und Alberto Valli: Numerical Approxi-mation of Partial Differential Equations. Springer-Verlag, Berlin,Heidelberg, New York, 1994.

[Ric02] Richter, Gerard R.: Stable Finite Element Box Methods forHyperbolic Equations. SIAM J. Numer. Anal., 39(5):1667–1683,2002.

[SC91] Staniforth, Andrew und Jean Cote: Semi-Lagrangian In-tegration Schemes for Atmospheric Models — A Review. Mon.Wea. Rev., 119:2206–2223, 1991.

[Zie84] Zienkiewicz, O. C.: Methode der finiten Elemente. Carl HanserVerlag, Munchen, Wien, 2. Auflage, 1984.

[Zie96] Zienkiewicz, O. C.: Origins, Milestones and Directions of theFinite Element Method – A Personal View, Band 4 der ReiheHandbook of Numerical Analysis, Kapitel 1, Seiten 3–67. ElsevierScience B.V., 1996.

Date post:	12-Aug-2019
Category:	Documents
Upload:	duongquynh
View:	213 times
Download:	0 times

Diplomarbeit von Lars Mentrup - M3/Allgemeines · Technische Universit at Munchen Fakult at f ur...

Documents