Adaptive Discontinuous-Galerkin-Verfahren zur L¨osung der ... · 1 Hiermit erkl¨are ich, dass ich...

Technische Universitat Munchen

Fakultat fur Mathematik

Adaptive

Discontinuous-Galerkin-Verfahren zur

Losung der Flachwassergleichungen

mit verschiedenen Flussfunktionen

Diplomarbeit von Robert Hochstetter

Aufgabensteller: Prof. Dr. Hans-Joachim Bungartz

Betreuer: Dr. Michael Bader

Abgabetermin: 1. Juni 2009

1

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

Munchen, den 2. Juni 2009

Inhaltsverzeichnis

1 Einleitung 8

2 Mathematische Grundlagen 10

2.1 Die Flachwasser-Gleichungen . . . . . . . . . . . . . . . . . . 102.2 Das Discontinuous-Galerkin-Verfahren . . . . . . . . . . . . . 14

2.2.1 Schwache Form . . . . . . . . . . . . . . . . . . . . . . 142.2.2 Diskrete Formulierung . . . . . . . . . . . . . . . . . . 152.2.3 Ansatz- und Testfunktionen . . . . . . . . . . . . . . . 152.2.4 Der numerische Fluss . . . . . . . . . . . . . . . . . . 18

3 Gitterstruktur und Stapelsystem 20

3.1 Gitterstruktur . . . . . . . . . . . . . . . . . . . . . . . . . . . 203.2 Traversierung mittels Sierpinski-Kurve . . . . . . . . . . . . . 233.3 Das Stapelsystem . . . . . . . . . . . . . . . . . . . . . . . . . 25

3.3.1 Die Stapel . . . . . . . . . . . . . . . . . . . . . . . . . 253.3.2 Die Datenstruktur . . . . . . . . . . . . . . . . . . . . 263.3.3 Typeinteilung und Kantenparameter . . . . . . . . . . 283.3.4 Behandlung der Randkanten . . . . . . . . . . . . . . 32

3.4 Adaptivitat . . . . . . . . . . . . . . . . . . . . . . . . . . . . 333.4.1 Verfeinerung . . . . . . . . . . . . . . . . . . . . . . . 343.4.2 Vergroberung . . . . . . . . . . . . . . . . . . . . . . . 353.4.3 Vorgehen bei der Adaption . . . . . . . . . . . . . . . 36

4 Numerische Flusse 38

4.1 Allgemeine Herleitung . . . . . . . . . . . . . . . . . . . . . . 384.2 Schockwellen . . . . . . . . . . . . . . . . . . . . . . . . . . . 404.3 Der Lax-Friedrichs-Fluss . . . . . . . . . . . . . . . . . . . . . 414.4 Der Fluss nach Roe . . . . . . . . . . . . . . . . . . . . . . . . 44

2

INHALTSVERZEICHNIS 3

4.4.1 Herleitung . . . . . . . . . . . . . . . . . . . . . . . . . 454.4.2 Berechnung der Matrix A . . . . . . . . . . . . . . . . 46

5 Details zur Implementierung 51

5.1 Grundlegende Annahmen . . . . . . . . . . . . . . . . . . . . 515.2 Anpassung des Algorithmus und der Datenstruktur . . . . . . 52

5.2.1 Der modifizierte Lax-Friedrichs-Fluss . . . . . . . . . . 525.2.2 Der lokale Lax-Friedrichs-Fluss . . . . . . . . . . . . . 535.2.3 Der globale Lax-Friedrichs-Fluss . . . . . . . . . . . . 55

5.3 Details zur Berechnung des Lax-Friedrichs -Parameters λ . . 575.3.1 Die Transportgleichung . . . . . . . . . . . . . . . . . 575.3.2 Die Flachwasser-Gleichungen mit konstanten Ansatz-

funktionen . . . . . . . . . . . . . . . . . . . . . . . . 575.3.3 Die Flachwasser-Gleichungen mit linearen Ansatzfunk-

tionen . . . . . . . . . . . . . . . . . . . . . . . . . . . 585.4 Behandlung der Randkanten . . . . . . . . . . . . . . . . . . . 625.5 Details zur Berechnung des Roe-Flusses . . . . . . . . . . . . 63

6 Numerische Ergebnisse 65

6.1 Rechenzeit- und Speicheranalyse . . . . . . . . . . . . . . . . 656.1.1 Speicheranalyse . . . . . . . . . . . . . . . . . . . . . . 666.1.2 Rechenzeitanalyse . . . . . . . . . . . . . . . . . . . . 66

6.2 Testrechnungen anhand des radialsymmetrischen Dammbruch-Problems . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 676.2.1 Konstante Ansatzfunktionen . . . . . . . . . . . . . . 686.2.2 Lineare Ansatzfunktionen . . . . . . . . . . . . . . . . 716.2.3 Der Parameter λ . . . . . . . . . . . . . . . . . . . . . 72

7 Abschließende Bemerkung 80

A Verwendete mathematische Verfahren 82

A.1 Explizites Euler-Verfahren . . . . . . . . . . . . . . . . . . . . 82A.2 Runge-Kutta-Verfahren zweiter Ordnung . . . . . . . . . . . . 82

B matlab-Funktionen 85

B.1 matlab-Funktion fur Daten in Matrixform . . . . . . . . . . . 85B.2 matlab-Funktion fur Daten in Vektorform . . . . . . . . . . . 87B.3 matlab-Funktion fur Bild-Plots . . . . . . . . . . . . . . . . . 89

INHALTSVERZEICHNIS 4

B.4 matlab-Funktion fur zweidimensionalen Schnitt . . . . . . . . 91

Abbildungsverzeichnis

2.1 Definition von ξ und hb der Flachwassergleichungen nach Ai-zinger er al. [5] (Bock [3]) . . . . . . . . . . . . . . . . . . . . 11

2.2 Definition von ξ und H in der konservativen Form (2.4) derFlachwassergleichungen (Bock [3]) . . . . . . . . . . . . . . . 13

2.3 Uberdeckung von Ω mit Dreiecken Tk (Bock [3]) . . . . . . . 152.4 Darstellung durch unstetige Polynome in 1D (Bock [3]) . . . . 162.5 Das Crouzeix-Raviart-Element . . . . . . . . . . . . . . . . . 17

3.1 Hangender Knoten . . . . . . . . . . . . . . . . . . . . . . . . 213.2 Initiale Triangulierung des Berechnungsgebietes Ω . . . . . . 223.3 Verfeinerungsbaum eines rekursiv verfeinerten Dreiecks (Schrauf-

stetter [12]) . . . . . . . . . . . . . . . . . . . . . . . . . . . . 223.4 Die ersten sechs Iterationen der Sierpinski-Kurve (Schraufs-

tetter [12]) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 243.5 Dreiecksgitter mit Sierpinski-Kurve und zugehoriger Verfei-

nerungsbaum (Schraufstetter [12]) . . . . . . . . . . . . . . . 243.6 Die Sierpinski-Kurve teilt die Knoten in grune und rote Kno-

ten (Schraufstetter [12]) . . . . . . . . . . . . . . . . . . . . . 263.7 Mogliche Falle bei der Traversierung eines Dreiecks (Schrauf-

stetter [12]) . . . . . . . . . . . . . . . . . . . . . . . . . . . . 283.8 Nummerierung der Ecken eines Dreiecks (Schraufstetter [12])

bisher (links) und beim Crouzeix-Raviart-Element . . . . . . 293.9 Rekursive Definition des Dreieckstyps mit Kantenparametern

alt/neu (Schraufstetter [12]) . . . . . . . . . . . . . . . . . . . 313.10 Dreieckstypen mit Farbkanten (Radzieowski [8]) . . . . . . . . 313.11 Rekursive Farbung der Knoten: Die Referenzfarbe wechselt

zwischen den Leveln (Schraufstetter [12]). . . . . . . . . . . . 32

5

ABBILDUNGSVERZEICHNIS 6

3.12 Unklarheit bei hangenden Knoten: Wie werden die Kanten-parameter aufgeteilt? . . . . . . . . . . . . . . . . . . . . . . . 34

3.13 Verfeinerung des Gitters (Radzieowski [8]) . . . . . . . . . . . 343.14 Die vier Falle der Verfeinerung (Radzieowski [8]) . . . . . . . 35

4.1 Das Riemann-Problem mit den zwei Initialzustanden q− undq+ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39

4.2 Ausbreitung der Schockfronten fur die eindimensionalen Flachwasser-Gleichungen nach Leveque [7] . . . . . . . . . . . . . . . . . . 40

4.3 Die Zustande q−, q∗ und q+ getrennt durch Unstetigkeitenuber die Charakteristiken λ1 und λ2. . . . . . . . . . . . . . . 41

4.4 Die zur Berechnung des Lax-Friedrichs-Flusses notwendigenKnoten q− und q+ auf der Kante e. . . . . . . . . . . . . . . . 42

4.5 Verschmieren der Lax-Friedrichs-Losung (blaue, gestrichelteLinie) im Vergleich zur exakten Losung (rote, durchgezogeneLinie). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43

4.6 Nichtlineare Charakteristiken λi fur nichtlineare hyperboli-sche Erhaltungsgleichungen . . . . . . . . . . . . . . . . . . . 44

5.1 Flusse in Richtung der außeren Normalen nei . . . . . . . . . . 515.2 Abschatzung von v2

x und v2y durch die Sekanten v∗x und v∗y . . . 61

6.1 Lokaler Lax-Friedrichs-Fluss, konstante Ansatzfunktionen . . 696.2 Schnitt durch den Losungsverlauf bei konstanten Ansatzfunk-

tionen . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 706.3 Nicht konvergentes Verfahren bei lokalem Lax-Friedrich-Fluss

mit linearen Ansatzfunktionen . . . . . . . . . . . . . . . . . 716.4 Lokaler Lax-Friedrichs-Fluss, lineare Ansatzfunktionen, λ mit

Maximumsnorm abgeschatzt. . . . . . . . . . . . . . . . . . . 736.5 Schnitt durch den Losungsverlauf bei linearen Ansatzfunktionen 746.6 Verlauf der maximalen λ-Werte pro Zeitschritt . . . . . . . . 756.7 Lokaler Lax-Friedrichs-Fluss und Lax-Friedrichs-Parameter λ 776.8 Verlauf der Losung und des Lax-Friedrichs-Parameters λ . . . 78

A.1 Berechnung eines Zeitschrittes beim expliziten Euler-Verfahren. 83A.2 Berechnung eines Zeitschrittes beim Runge-Kutta-Verfahren

zweiter Ordnung . . . . . . . . . . . . . . . . . . . . . . . . . 83

ABBILDUNGSVERZEICHNIS 7

Kapitel 1

Einleitung

In der Technik ist es von großer Wichtigkeit, Stromungsverhalten von Flussig-keiten und deren Auswirkungen auf ein System, wie zum Beispiel die Stro-mung in Turbinen (im Maschinenbauwesen) oder die Stromung durch einenKanal beziehungsweise um Hindernisse wie Bruckenpfeiler (im Bauinge-nieurwesen) simulieren zu konnen.

Um solche Stromungsverlaufe von Wasser, oder allgemeiner um die Stromungvon Fluiden zu modellieren, gibt es verschiedene Ansatze. Haufig angewen-det wird dabei die Beschreibung der Stromung durch ein System partiellerDifferentialgleichungen zweiter Ordnung, die Navier-Stokes-Gleichungen.

Aufgrund ihrer recht allgemeinen Natur lassen sich diese Navier-Stokes-Gleichungen fur bestimmte Anwendungen sehr vereinfachen. So kann zumBeispiel zur Simulation von Wasserstromungen die Kompressibilitat außerAcht gelassen werden, da Wasser als nahezu inkompressibel gilt. Die Glei-chungen werden folglich weniger umfangreich und somit weniger aufwandigzu losen.

Eine solche spezielle Vereinfachung der Navier-Stokes-Gleichungen sind dieFlachwasser-Gleichungen, mit denen Wasserstromungen in einer oder zweiRaumrichtungen simuliert werden konnen. Dabei wird angenommen, dassdas Verhaltnis der Wassertiefe zu den anderen beiden Raumdimensionen(der Breite und der Lange) sehr klein ist. Genaueres zu diesen Vereinfa-chungen findet man zum Beispiel in der Arbeit von Bock [3]. Diese Flach-wassergleichungen werden beispielsweise in der Ozean-Simulation eingesetzt

8

KAPITEL 1. EINLEITUNG 9

um die Wellenausbreitung nach einem Dammbruch oder nach einem Seebe-ben vorherzusagen.

Da die Berechnungsgebiete bei solchen Aufgabenstellungen naturgemaß sehrgroß sind, aber im Gegenzug die Berechungen an den entstehenden Schock-fronten recht genau sein mussen, werden bei der numerischen Losung an dieverwendeten Rechner bezuglich Rechenzeit und insbesondere Speicherplatzgroße Anforderungen gestellt. Es ist daher von großem Interesse, effizien-te Verfahren zur Losung dieses Problems zu konstruieren beziehungsweiseweiterzuentwickeln, welche besonderes Augenmerk auf diese beiden begren-zenden Faktoren richten.

In dieser Arbeit sollen die Flachwasser-Gleichungen mit einem Discontinuous-Galerkin-Ansatz gelost werden, welcher auf einem besonders speichereffizi-enten Ansatz mit raumfullenden Kurven als Datenstruktur arbeitet. DasKernstuck dieser Verfahren sind die Flussfunktionen, die den Transport derzu berechnenden Großen durch die Raumdimensionen beschreiben. In dieserArbeit sollen diese Flussfunktionen naher betrachtet werden und ihr Einflußauf das Verfahren untersucht werden.

Dazu wird in Kapitel 2 an die mathematischen Grundlagen wie das Dis-continuous-Galerkin-Verfahren und eine kurze Herleitung der Flachwasser-Gleichungen herangefuhrt. In Kapitel 3 wird dann naher auf den informa-tischen Teil, namlich die bisherige Implementierung, eingegangen, um denAblauf einer Traversierung mit den raumfullenden Kurven und das Stapel-system zum Datenaustausch zu erklaren. Auf den folgenden Kapiteln liegtdas Hauptaugenmerk dieser Arbeit. In Kapitel 4 werden die mathemati-schen Aspekte der Flussfunktionen und deren Unterschiede diskutiert. InKapitel 5 werden die sich durch deren Neuimplementierung ergebendenalgorithmischen Veranderungen erklart und in Kapitel 6 einige Ergebnissenumerischer Simulationen anhand einfacher Testfalle dargestellt. In Kapitel

7 wird schließlich ein kurzer Ausblick auf mogliche Erweiterungen gegeben.

Kapitel 2

Mathematische Grundlagen

In diesem Kapitel werden die mathematischen Grundlagen fur diese Arbeitzusammengefasst. Dabei werden zuerst die Flachwassergleichungen kurz vor-gestellt, anschließend wird die verwendete numerische Klasse der Disconti-nuous-Galerkin-Verfahren und deren mathematische Grundlagen, also schwa-che Form und die spezielle Wahl der Ansatz- und Testfunktionen beschrie-ben.

2.1 Die Flachwasser-Gleichungen

Die Flachwassergleichungen sind ein System partieller Differentialgleichun-gen in zwei Raumdimensionen in den Variablen

• ξ ∈ R, der Wellenhohe (also der Hohenabweichung vom Meeresspiegel),und

• v =

(vx

vy

)∈ R2, der Geschwindigkeit in x- und y-Richtung.

Nach Aizinger und Dawson [5] lautet eine Darstellung:

∂ξ

∂t+∇ · (vH) = 0, (2.1)

∂v∂t

+ v · ∇v + τbfv + fck× v + g∇ξ − νT

H∆(Hv) =

1H

F . (2.2)

10

KAPITEL 2. MATHEMATISCHE GRUNDLAGEN 11

hb

ξ

Abbildung 2.1: Definition von ξ und hb der Flachwassergleichungen nachAizinger er al. [5] (Bock [3])

Die ubrigen verwendeten Großen haben folgende Bedeutung:

• Die Zeit t,

• die Gesamtwassertiefe H = ξ + hb (siehe Abbildung (2.1)), dabei ist ξ

die Wellenhohe zum mittleren Wasserspiegel und hb die Meerestiefe,

• der Reibungskoeffizient fur die Reibung am Boden τbf ,

• der Parameter fur die Coriolis-Kraft fc und der zugehorige lokale Nor-malenvektor k,

• der Ortsfaktor g,

• die tiefengemittelte Zahigkeit der Flussigkeit νT und

• die Variable F, in der weitere Gegebenheiten wie der variable Druckder Atmosphare vorgenommen werden konnen.

(2.1) heißt Transportgleichung und modelliert die physikalische Eigenschaftder Massenerhaltung, denn da Wasser als inkompressibel gilt, kann die Was-serhohe mit der Masse identifiziert werden. Gleichung (2.2) beschreibt dieImpulserhaltung des Systems. Die Transportgleichung wird als erste ein-fache Testgleichung mit konstanter Geschwindigkeit v in dieser Form ver-wendet, bevor nach Neuerungen im Verfahren mit dem ganzen Flachwasser-Gleichungssystem gerechnet wurde.


Um den Umgang und die Berechnung zu erleichtern werden in dieser Ar-beit die meisten Einflussgroßen des Modells vernachlassigt, wodurch einekompaktere Darstellung der Flachwassergleichungen erreicht wird:

∂ξ

∂t+∇ · (vH) = 0,

∂v∂t

+ v · ∇v + g∇ξ = 0 . (2.3)

Um eine vollstandigere Herleitung der Flachwasser-Gleichungen einzusehen,sei auf Schwanenberg [14] verwiesen.

Zum Losen des kompletten Systems der Flachwasser-Gleichungen wird eineandere sehr kompakte Formulierung verwendet, die sich an der Darstellungvon Remacle et al. in [9] orientiert:

∂

∂t

ξ

ξvx

ξvy

+ div

ξv

ξvxv + 12gξ2ex

ξvyv + 12gξ2ey

=

0−gξ(∂xH + τbf )−gξ(∂yH + τbf )

. (2.4)

Oder vereinfacht:

∂

∂tq + div F(q) =

∂

∂t

q1

q2

q3

+ div

F1(q)F2(q)F3(q)

= r. (2.5)

Dabei sind ex und ey die Einheitsvektoren in x- und y-Richtung, q der drei-dimensionale Vektor der Zustandsvariablen und F ein sechsdimensionalerFlussvektor (da jedes v zweidmensional ist). Jedes Fi(q) ist also zweidimen-sional. Der Divergenz-Operator div wird dabei als Summe der partiellenAbleitungen in die beiden Raumrichtungen komponentenweise verstanden:

div F =

divF1(q)divF2(q)divF3(q)

=

∇ · F1(q)∇ · F2(q)∇ · F3(q)

(2.6)

Im Vergleich zur Darstellung (2.2) ist hier die Wasserhohe ξ unterschied-lich definiert. War in der ersten Darstellung noch die Abweichung der Was-serhohe vom mittleren Wasserspiegel hb als ξ gegeben, ist ξ nun die totale


ξ

H

z

x

Abbildung 2.2: Definition von ξ und H in der konservativen Form (2.4) derFlachwassergleichungen (Bock [3])

Wassertiefe (siehe Abbildung (2.2)) und H der Abstand des Wassergrundes,die Tiefe des Wasserbetts, zum Koordinatenlevel ”z = 0“. Man nennt dieseForm der Flachwassergleichungen auch die konservative Form, da das Pro-blem in den erhaltenen (konservativen) Variablen formuliert ist. Die weiterenVariablen der Gleichung (2.4) sind namlich nicht mehr die Geschwindigkei-ten vx und vy, sondern die Impulsdichten in beiden Raumrichtungen ξvx undξvy, also die Variablen, welche den Gesetzen der Impulserhaltung genugen.

Nimmt man vereinfachend einen ”flachen”Meeresgrund ohne jegliche Stei-gung an und vernachlassigt auch noch die Reibung des Wassers am Boden, soerhalt man die in dieser Arbeit geloste Form der Flachwasser-Gleichungen:

∂

∂tq + div F(q) = 0.


2.2 Das Discontinuous-Galerkin-Verfahren

2.2.1 Schwache Form

Sucht man bei der Behandlung partieller Differetialgleichungen nach ei-ner Losung, die die Gleichungen und ihre Anfangs- und Randbedingungenerfullt, stellt sich als erstes die Frage nach deren Existenz. Fur viele Proble-me ergeben sich dabei Schwierigkeiten: es existiert oftmals keine sogenannteklassische Losung, welche die Anforderungen (zum Beispiel Glattheitsanfor-derungen) erfullt.

Aus diesem Grunde hat man den Begriff einer schwachen Losung (im Ge-gensatz zur klassischen Losung) eingefuhrt. Diese schwache Losung muss derDifferentialgleichung nicht punktweise genugen, sondern stattdessen eine In-tegralbeziehung erfullen.

Man wahlt eine Funktion φ aus einem Funktionenraum V , genannt Test-raum, multipliziert die Gleichung mit dieser Funktion. Anschließend wird dieGleichung integriert und gefordert, dass die Gleichung fur alle φ ∈ V erfulltsei. Dies verringert die Glattheitsanforderungen an die Losung, ermoglichtdamit also die Losung großerer Gleichungsklassen. In der Vektorform (2.4)ergibt dies:Suche ein q fur das die Gleichung∫

Ω

∂q∂t

φ dxdy +∫

ΩdivF(q)φ dxdy = 0 (2.7)

fur alle φ ∈ V erfullt ist. Ω ist dabei das Berechnungsgebiet.Das Integral ist hier naturlich komponentenweise fur jede Gleichung zu ver-stehen:∫

Ω

∂qi

∂tφ dxdy +

∫Ω∇ · Fi(q)φ dxdy = 0 ∀φ ∈ V, i = 1 : 3. (2.8)

Nun wendet man auf das zweite Glied der linken Seite den Gauß‘schen In-tegralsatz an, und transformiert dieses Integral in ein Oberflachenintegral:∫

Ω

∂qi

∂tφdxdy−

∫Ω

Fi(q)∇φdxdy+∫

∂ΩFi(q)·nφds = 0 ∀φ ∈ V, i = 1 : 3.

(2.9)


Abbildung 2.3: Uberdeckung von Ω mit Dreiecken Tk (Bock [3])

Dabei ist n der außere Einheitsnormalenvektor auf dem Rand ∂Ω.

2.2.2 Diskrete Formulierung

Der Raum V der Testfunktionen φ ist ein unendlichdimensionaler Raum,denn die Funktionen sind durch ihre Werte an allen Punkten des Berech-nungsgebiets definiert. Um den Raum der moglichen Losungen auf eine end-liche Dimension zu beschranken, zerlegt man das Berechnungsgebiet Ω nunin disjunkte Teilgebiete Tk, in unserem Fall in gleichschenklig-rechtwinkligeDreiecke (siehe Abbildung (2.3)). Einer der Grunde dieser Wahl wird in Ka-pitel 3 erlautert.

Beim Discontinuous-Galerkin-Verfahren wird nun die Gleichung (2.9) nichtmehr geschlossen auf ganz Ω gelost, sondern auf jedem Tk.

∫Tk

∂qi

∂tφdxdy−

∫Tk

Fi(q)∇φdxdy+∫

∂Tk

Fi(q)·nφds = 0 ∀φ ∈ V, i = 1 : 3.

(2.10)Dabei braucht die Gesamtlosung nicht stetig zu sein, was sich im Namen

des Discontinuous-Galerkin-Verfahrens vererbt hat.

2.2.3 Ansatz- und Testfunktionen

Nun wird der Funktionenraum V durch einen endlichdimensionalen RaumVh ersetzt. Im Allgemeinen wahlt man Vh als Spann von Funktionen, die gut


p1 p

2p3

I1

I2 I

3

Abbildung 2.4: Darstellung durch unstetige Polynome in 1D (Bock [3])

beherrschbar sind. Dies sind speziell Polynome von begrenztem Grad.

V mh :=

p|p : R2 → R; p ist Polynom vom Grad ≤ m auf Tk

(2.11)

Die Aufgabenstellung lautet nun also wie folgt:Finde eine Approximation (qh)i ∈ V m

h , so dass folgende Gleichung erfulltist:∫

Tk

∂(qh)i

∂tφdxdy −

∫Tk

Fi(qh)∇φdxdy +∫

∂Tk

Fi(qh) · nφds = 0, (2.12)

∀φ ∈ V mh , i = 1 : 3, ∀Tk ∈ Ω.

Fur m = 0 sind die gesuchten Losungen (qh)i also stuckweise konstanteFunktionen, fur m = 1 eine stuckweise lineare Funktionen und so weiter,welche die verschiedenen Losungskomponenten des Differentialleichungssy-tems im Raum V m

h approximieren.

In dieser Arbeit werden fur V mh die Raume der konstanten und der linearen

Polynome V 0h und V 1

h verwendet, die Losung setzt sich also aus stuckweisekonstanten, beziehungsweise linearen Funktionen zusammen.Als Basis fur die Raume V m

h , m = 1, 2, auf jedem Tk werden die Lagrange-Polynome li(x, y) vom Grad m verwendet. Also kann beispielweise jedesElement (qh)i ∈ V 1

h , i = 1 : 3, als Linearkombination der Basisfunktionen


h

hx

y

23

1

Abbildung 2.5: Das Crouzeix-Raviart-Element mit den Freiheitsgraden aufden Kantenmittelpunkten (Bock [3])

dargestellt werden:

(qh)i(t, x, y) =d∑

j=1

qji (t)l

j(x, y). (2.13)

Es werden somit in jedem Zeitschritt konstante Koeffizienten qji (t) fur die

Basisfunktionen gesucht, um die Approximationslosung darzustellen.

Da das Gitter aus rechtwinklig-gleichschenkligen Dreiecken besteht (sieheAbbildung (2.3)), reicht es die Ansatzfunktionen fur ein Referenzelement,dem Crouzeix-Raviart-Element (siehe Abbildung (2.5)), zu definieren. Umdie Ansatzfunktionen fur diejenigen Elemente zu erhalten, die gedreht dazuin der Ebene liegen ist dann nur eine einfache Transformation notig. Dielinearen Ansatzfunktionen lj(x, y), j = 1 : 3 fur das Referenzelement lautenfolgendermaßen:

l1 = 1− 2h

y

l2 =2h

x +2h

y − 1, (2.14)

l3 = 1− 2h

x

Diese Wahl der Lagrange-Polynome als Ansatzfunktionen hat den Vorteil,dass der Wert der gesuchten Funktion auf den Kantenmittelpunkten jedesElements durch lediglich einen Koeffizienten bestimmt wird, da die beiden


anderen Ansatzfunktionen dort gleich Null sind. Fur eine genauere Beschrei-bung der linearen Ansatzfunktionen sei auf Schwaiger [13] verwiesen.

Die Forderung, dass Gleichung (2.10) fur alle φ ∈ V gilt, wird nun durch eineschwachere Forderung ersetzt. Es genugt nun, diese Gleichung fur alle φ ausdem Ansatzraum V m

h zu erfullen. Damit erhalt man die Aufgabenstellung:Finde (qh)i ∈ V m

h , so dass die Gleichungen∫Tk

∂(qh)i

∂tlj dxdy −

∫Tk

Fi(qh)∇lj dxdy +∫

∂Tk

Fi(qh) · nlj ds = 0 (2.15)

fur alle Basisfunktionen lj von V mh gelten.

2.2.4 Der numerische Fluss

Schließlich wird der dritte Summand, das Oberflachenintegral, in Gleichung(2.15) betrachtet. Da die Gesamtlosung aus den auf den Elementen kon-stanten beziehungsweise linearen Funktionen nicht notwendig stetig ist, istnicht klar, ob die Werte auf dem Rand des Elements Tk oder diejenigen desjeweiligen Nachbarelements benutzt werden sollen. Die Auswertung diesesRandintegrals ist damit nicht eindeutig.

Man behilft sich dabei mit einem Kunstgriff. Der Flussterm Fj(qh) · n imRandintegral wird durch einen numerischen Fluss (F · n)∗e uber die jeweili-gen Kanten e des Elements ersetzt. Dieser Fluss stellt die einzige Interaktionzwischen zwei benachbarten Elementen dar, muss also physikalisch gesehenden Massen- und Impulsaustausch moglichst gut beschreiben. Anschaulicherzeugt der numerische Fluss uber eine Kante e bei der Flachwassergleichungzum Beispiel das Hin- und Herschwappen von Wasser zwischen den zwei zurKante e benachbarten Elementen.

Die Aufgabenstellung hat sich damit folgendermaßen verandert:Finde Funktionen (qh)i ∈ V m

h , die die folgenden Gleichungen erfullt:∫Tk

∂(qh)i

∂tljdxdy −

∫Tk

Fi(qh)∇ljdxdy +∑

e

∫e(F · n)∗el

jds = 0 (2.16)

fur alle Basisfunktionen lj von V mh fur i = 1 : 3.


Damit bildet der numerische Fluss neben der Gitterzerlegung des Gebiets Ωeinen zentralen Teil des Discontinuous-Galerkin-Verfahrens. Die Flussfunk-tionen werden in Kapitel 4 mathematisch analysiert. In Kapitel 6 werdeneinige Ergebnisse numerischer Tests angegeben und bewertet.

Kapitel 3

Gitterstruktur und

Stapelsystem

Dieses Kapitel befasst sich mit der informatischen Seite des bisher implemen-tierten Verfahrens. Wie bereits angedeutet, wird dies durch eine speicher-effiziente Abarbeitungsreihenfolge basierend auf einer raumfullenden Kurve,der Sierpinski-Kurve, erreicht. Es werden hier die Anforderungen und Be-reitstellung des Gitters und die Abarbeitungsreihenfolge sowie insbesondereder Datenaustausch auf Basis eines Stapelsystems geklart. Am Ende desKapitels wir noch kurz auf die Durchfuhrung einer Adaption eingegangen.

3.1 Gitterstruktur

Wie bereits kurz erwahnt wird das Berechnungsgebiet Ω in dieser Arbeit ingleichschenklig-rechtwinklige Dreiecke zerlegt. Dreiecke werden oft als Ele-mente des Gitters gewahlt, weil sich hier polynomielle Ansatzfunktionen sodarstellen lassen, dass sich speziell die Werte auf den Ecken und Kantenmit-telpunkten sehr effizient berechnen lassen (siehe Crouzeix-Raviart-Element,Kapitel (2.2.3)).

An das Dreiecksgitter werden nun folgende Anforderungen gestellt:

• mathematische Anforderungen

Um sinnvolle Ergebnisse zu erhalten wird die Zulassigkeit des Gittersgefordert. Diese Forderung umschließt, dass Ω von der Triangulierungkomplett uberdeckt werden muss und dass alle Dreieckselemente bis

20

KAPITEL 3. GITTERSTRUKTUR UND STAPELSYSTEM 21

1

2

3

4

Abbildung 3.1: Hangender Knoten. Der (doppelt gezeichnete) mittlere Kno-ten ist ein Knoten der Dreiecke 2, 3 und 4, aber kein Knoten von Dreieck1.

auf ihre Kanten und Ecken disjunkt sein mussen. Außerdem darf keinKnoten als hangender Knoten existieren. Dies bedeutet, dass jederEckpunkt eines Dreiecks auch Eckpunkt des benachbarten Dreieckssein muss (siehe Abbildung (3.1)). Zur Erlauterung der Grunde dieserForderungen siehe Radzieowski [8] und Schraufstetter [12].

• informatische Anforderungen

Um die Rechenzeit- und insbesondere die Speichereffizienz zu gewahr-leisten soll das Gitter adaptiv verfeinerbar sein. Es sollen wahrend derBerechnung also nur diejenigen Dreiecke verfeinert werden, auf deneneine bestimmte mathematische Große, namlich der lokale Fehler, einenGrenzwert uberschreitet. Die restlichen Dreiecke sollen (bis auf die Ver-feinerung zur Vermeidung hangender Knoten) nicht verfeinert werden,sondern wenn moglich, sogar vergrobert werden, um speicherschonendzu arbeiten.

Unstrukturierte Dreiecksgitter erfordern einen wesentlichen Mehraufwandbeim Abspeichern der Daten, denn im Gegensatz zu regularen Dreiecks-gittern, wie dem hier verwendeten, ware es notig, die Nachbarschaftsbezie-hungen aller Dreieckselemente abzuspeichern. Dies ist ein weiterer Grund,ein regelmaßiges Gitter in Form von gleichschenklig-rechtwinkligen Dreie-cken zu verwenden. Die Speichereffizienz genießt namlich sehr hohe Prioritat(langere Rechenzeiten konnen gegebenfalls verkraftet werden, wahrend derSpeicherplatz hingegen begrenzt ist).


1

1

0

Abbildung 3.2: Initiale Triangulierung des Berechnungsgebietes Ω

1

2 3

24

5

1

2 3

4 5

Abbildung 3.3: Verfeinerungsbaum eines rekursiv verfeinerten Dreiecks(Schraufstetter [12])

Das Berechnungsgebiet Ω = [0; 1]× [0; 1] wird zur Gittergenerierung unter-teilt in 4 Initial-Dreiecke, die durch die Diagonalen des Einheitsquadrats be-grenzt werden (siehe Abbildung (3.2)). Soll das Gitter verfeinert werden, sowird rekursiv im zu verfeinernden Dreieck auf der Mitte der Hypotenuse einneuer Knoten und von diesem zum gegenuberliegenden Knoten eine neueKante eingefugt, bis der gewunschte Verfeinerungsgrad erreicht ist. Dabeimuss naturlich beachtet werden, dass keine hangenden Knoten entstehen,angrenzende Dreiecke mussen also gegebenenfalls auch verfeinert werden.


Auf diese Weise wird ein binarer Verfeinerungsbaum erzeugt, der als Wur-zel ein Initial-Dreieck besitzt. Jedes der beiden Kinder eines Knotens die-ses Baumes ist also ein verfeinertes Dreieck mit halbem Flacheninhalt, dieBlatter des Baumes sind die feinen Dreiecke, auf denen tatsachlich die Be-rechnung durchgefuhrt wird (Abbildung (3.3)).

Dieser Baum wird effizienterweise aber nicht als Baum gespeichert, sondernals linearisierter Bitstream. Dabei wird der Baum mittels Tiefensuche durch-laufen, die ”0“ fur jedes Blatt speichert und fur jeden anderen Knoten ”1“.Die Große des Dreiecks beziehungsweise des Blattes auf dem die Berechnungstattfindet, kann dabei durch dessen Level im Baum oder durch die Positionim Bitstream bestimmt werden.

3.2 Traversierung mittels Sierpinski-Kurve

Zur Berechnung der Losung mit dem hier angewandten Verfahren wird dieDifferentialgleichung zu diskreten Zeitpunkten auf jedem Dreieck Tk desGitters gelost. Es muss also fur jedes Blatt des Verfeinerungsbaumes, sprichfur jede ”0“ des Bitstreams, eine Losung berechnet werden. Dazu wird dieGleichung (2.16) mit einem numerischen Verfahren, namlich dem explizitenEuler- oder dem (ebenfalls expliziten) Runge-Kutta-Verfahren zweiter Ord-nung gelost (siehe Anhang). Es wird also nacheinander fur jedes Tk ∈ Ωein kleines Gleichungssystem gelost. Man spricht von einem ”zellorientier-ten Ansatz“ im Gegensatz zum ”globalen Ansatz“, bei dem ein großes Glei-chungssystem fur das ganze Berechnungsgebiet gelost wird. Der Vorteil liegtdabei in einer geringeren Rechenzeit.

Es ist nun notwendig, sicherzustellen, dass jedes der Dreiecke genau einmaldurchlaufen wird, am besten immer in der gleichen Reihenfolge. Dies wirdvom Verlauf einer raumfullenden Kurve vorgegeben. Eine raumfullende Kur-ve ist eine Abbildung von R → Rn (hier von [0; 1] → [0; 1]2), die surjektivund stetig ist. Fur allgemeinere Definitionen siehe Bader [1] und Schrauf-stetter [12].

In dieser Arbeit wird dazu die Sierpinski-Kurve (Abbildung (3.4)) verwen-


Abbildung 3.4: Die ersten sechs Iterationen der Sierpinski-Kurve (Schraufs-tetter [12])

Abbildung 3.5: Dreiecksgitter mit Sierpinski-Kurve und zugehoriger Verfei-nerungsbaum (Schraufstetter [12])

det, die vom namensgebenden polnischen Mathematiker Waclaw Sierpin-ski (1882-1969) im Jahre 1912 definiert wurde. Im Gegensatz zu anderengelaufigen raumfullenden Kurven wie der Peano-Kurve arbeitet sie auf gleichschenklig-rechtwinkligen Dreiecken und nicht auf Quadraten. Sie ist rekursiv definiertund entsteht, wenn man ein Ausgangsdreieck durch Bisektion halbiert. Manverbindet Mittelpunkte der neu entstandenen Dreiecke und verfeinert an-schließend auch diese Dreiecke weiter nach der oben genannten Vorschrift.Fur eine formale Definition und eine Beschreibung der erzeugenden Gram-matik siehe Schraufstetter [12].

Die Sierpinski-Kurve und die oben beschriebene Diskretisierung des Berech-nungsgebietes Ω werden also nach der selben Vorschrift erzeugt. Die Tra-versierung des Gitters auf Basis der Sierpinski-Kurve ist aquivalent zumDurchlauf des 0-1-Bitstreams und damit des Verfeinerungsbaumes mittelsTiefensuche (siehe Abbildung (3.5)). Als theoretisches Konstrukt hinter derDurchfuhrung der Traversierung garantiert die Sierpinski-Kurve also dieKorrektheit des Durchlaufs. Um schließlich das ganze Berechnungsgebietbestehend aus dem vier Initial-Dreiecken zu durchlaufen, werden die vier


Sierpinski-Kurven hintereinander ausgefuhrt. Die Frage der Abarbeitungs-reihenfolge ist damit nun geklart.

3.3 Das Stapelsystem

Beim zellorientierten Ansatz wird auf jedem Dreieck Tk ein kleines linea-res Gleichungssystem gelost. Bisher wurde dabei außer Acht gelassen, dassdie einzelnen Elemente Informationen austauschen mussen, was aber bishernoch nicht gewahrleistet ist. Diese auszutauschenden Großen werden im nu-merischen Verfahren (2.15) durch das Oberflachenintegral, beziehungsweisein (2.16) durch den numerischen Fluss (F ·n)∗, beschrieben. Der Austauschwird durch ein System von vier Speichern vom Typus Stapel (engl. stack)(siehe Bader et al. [2]) durchgefuhrt.

3.3.1 Die Stapel

Ein Stapel ist ein Speicher, auf dem Daten nach dem last-in-first-out-Prinzip(LIFO) behandelt werden. Es sind auf einem Stapel also nur zwei Operatio-nen moglich: das Ablegen von Daten auf dem Stapel, genannt push und dasHolen des obersten Stapelelements vom Stapel pop.

Zwei der vier benotigten Stapel werden Eingabestapel und Ausgabesta-pel genannt. Bei der Traversierung des Gitters werden in der durch dieSierpinski-Kurve festgelegten Reihenfolge die Stapelelemente vom Eingabe-stapel geholt (pop). Wenn ein Stapelelement schließlich fertig bearbeitet ist,wird es auf den Ausgabestapel abgelegt (push). Die Daten liegen also jetztauf dem Ausgabestapel in umgekehrter Reihenfolge. Nun werden Ein- undAusgabestapel vertauscht, das heißt der Eingabestapel wird zum Ausgabe-stapel und umgekehrt. Die nachste Traversierung wird folglich ruckwarts,das heißt auch entlang der Sierpinski-Kurve, aber in der Gegenrichtungdurchgefuhrt.

Wie in Abbildung (3.6) zu sehen, wird von der Sierpinski-Kurve jeder Knotenzweimal besucht. Knoten 6 wird beispielsweise das erste Mal nach Knoten4 und Knoten 5 besucht, das zweite Mal dann nach Knoten 10 und Knoten7. Damit hier die korrekten Daten verwendet werden, sind die beiden ande-


1 2

3

4 5

6

78

9

10

11

13

12 14

Abbildung 3.6: Die Sierpinski-Kurve teilt die Knoten in grune und rote Kno-ten (Schraufstetter [12])

ren Stapel, die sogenannten Farbstapel, notig. Wie man im Bild (3.6) sieht,werden die Knoten der Dreieckselemente in zwei Mengen, namlich in gruneund rote Knoten geteilt. Dabei liegen die roten Knoten alle auf der rechtenSeite, die grunen Knoten alle auf der linken Seite der Kurve.

Auf diesen Farbstapeln werden die Informationen eines erst ein Mal besuch-ten Knotens abgelegt (und zwar auf dem Stapel der jeweiligen Knotenfarbe).Wenn der Knoten das zweite Mal besucht wird, wird die Knoteninformati-on von diesem Farbstapel ausgelesen. Die Korrektheit dieses Stapelsystemswurde von Schraufstetter [12] bewiesen. Diese Farbstapel werden in Kapitel5 noch von entscheidender Bedeutung bei der Implementierung der neuenFlussterme sein.

3.3.2 Die Datenstruktur

Bei dem Durchlauf entlang der Sierpinski-Kurve werden von jedem Dreieckdie Werte der Zustandsvariablen (q)i und die Flusswerte (Fi · n)∗ auf ei-ner Kante zur Berechnung des Wertes im neuen Zeitschritt benotigt. DieInformationen befinden sich namlich, wie in Abschnitt (2.2.3) beschrieben,beim Crouzeix-Raviart-Element nicht auf den Ecken jedes Elements, son-dern auf den Kantenmittelpunkten. Es mussen also die Kanten nach Farbeneingeteilt werden. Wie man spater sehen wird ist es allerdings nur notig,


die Farbe derjenigen Kante zu kennen, durch welche die Sierpinski-Kurvedas Dreieck weder betritt noch verlasst. Diese Kante, genannt Farbkante,bekommt die Farbe der beiden Knoten zugordnet, die diese Kante beruhren(siehe Abbildung (3.10).

Auf dem Ein- beziehungsweise Ausgabestapel liegen die Dreiecke mit Daten,namlich den Werten der korrespondierenden Zustandsvariablen (q)i und dennumerischen Flusswerten (F · n)∗ aller Kanten des Dreiecks. Auf den Farb-stapeln konnen aber nur die Informationen von Kanten der jeweiligen Farbegespeichert werden, da sonst die Korrektheit des Stapelsystems nicht mehrgewahrleistet ware. Deshalb arbeiten Ein- und Ausgabestapel auf einer an-deren Datenstruktur als die Farbstapel.

Auf den Farbstapeln wird eine Datenstruktur abgelegt, die Color-Stack-Element genannt wird. Diese beinhaltet als benotigte Berechnungsgroßenlediglich die Flusswerte einer einzigen Kante, namlich der oben beschriebe-nen Farbkante. Auf die Zustandwerte (q)i einer Kante kann dabei verzich-tet werden, denn die Flusswerte sind die einzigen mit dem Nachbarelementinteragierenden Werte, wie bereits in Abschnitt (2.2.4) festgestellt wurde.Hiermit sind diese Flusswerte die einzigen Werte, die uber die Farbstapelweitergegeben werden mussen.Dem gegenuber liegen auf dem Ein- und Ausgabestapel Elemente, in denensowohl die Werte der Zustandsvariablen als auch Flusswerte aller Kanteneines Dreiecks gespeichert sind, hier Stack-Elemente genannt. Kurz gesagtwerden also auf dem Ein- und Ausgabestapel die Dreiecke samt Informatio-nen gespeichert, wahrend auf den Farbstapeln lediglich Kanteninformatio-nen abgelegt werden.

Zusammenfassend wird in diesem Verfahren folgendermaßen vorgegangen:

• Benotigt man wahrend der Traversierung Daten eines (oBdA grunen)Knoten, so werden die Knotendaten

– aus dem Stack-Element ausgelesen, das frisch vom Eingabestapelgeholt wurde. Dies geschieht, wenn der Knoten wahrend dieserTraversierung noch nicht besucht wurde.

– vom grunen Farbstapel geholt, falls das Element bereits einmalbesucht wurde.


H KV

Abbildung 3.7: Mogliche Falle bei der Traversierung eines Dreiecks (Schrauf-stetter [12])

• Die Operation wird durchgefuhrt. Hierbei handelt es sich entweder umdie Berechnung eines neuen Zeitschritts, um die Durchfuhrung einerAdaption oder um die Ausgabe des Wertes.

• Nach der Durchfuhrung der Operation werden die Daten der Kante

– im Stack-Element gespeichert, falls der Knoten bereits zum zwei-ten Mal besucht wurde. Wenn die Berechnungen fur das ganzeStack-Element abgeschlossen asind, wird dies auf dem Ausgabe-stapel abgelegt.

– auf dem grunen Farbstapel gespeichert, falls der Knoten erst zumersten Mal in dieser Traversierung benutzt wurde.

• Anschließend kann mit dem nachsten Stack-Element fortgefahren wer-den.

• Sind alle Elemente abgearbeitet, so werden Ein- und Ausgabestapelsowie der rote und grune Stapel vertauscht.

In jedem Color-Stack-Element werden fur die jeweilige Kante zwei Fluss-werte fur jedes (q)i gespeichert. Der alte Flusswert, der dazu benutzt wird,die neuen Zustandswerte zu berechnen, und der neue Flusswert, der ausden neu berechneten Zustandwerten berechnet wird. Dieser wird dann inder nachsten Traversierung als alter Flusswert zur Berechnung des neuenZweitschritts benutzt.

3.3.3 Typeinteilung und Kantenparameter

Nun mussen die Dreiecke in verschiedene Typen eingeteilt werden, die imProgrammcode alle gesondert behandelt werden. Es wird nun vorausgesetzt,dass die Sierpinski-Kurve alle Dreiecke entlang der Hypotenuse durchlauft.


1

2

3

1

2

3

Abbildung 3.8: Nummerierung der Ecken eines Dreiecks (Schraufstetter [12])bisher (links) und beim Crouzeix-Raviart-Element

Wenn man Abbildung (3.7) betrachtet, dann kann man erkennen, dass be-zuglich der Kante, uber die die Kurve in ein Dreieck eintritt oder austritt,drei verschiedene Moglichkeiten bestehen. Diese drei Moglichkeiten sind:

• Eintritt uber die Hypotenuse, Austritt uber die Kathete (Dreiecksbe-zeichnung H),

• Eintritt uber die Kathete, Austritt uber die andere Kathete (Bezeich-nung V ) und

• Eintritt uber die Kathete, Austritt uber die Hypotenuse (BezeichnungK).

Die Knoten seien dabei folgendermaßen bezeichnet: Der Knoten, der derEcke gegenuberliegt, an dem die Kurve in das Dreieck eintritt, wird immerals Knoten 1, der Knoten gegenuber der Ecke an dem das Dreieck verlassenwird immer als Knoten 3 und der Knoten gegenuber der Hypotenuse alsKnoten 2 definiert.

Nun mussen noch weitere Kantenparameter definiert werden, um festzule-gen auf welchen Stapel zugegriffen werden muss, um die richtigen Daten zuerhalten. So muss fur jeden Knoten gespeichert werden, ob er das erste oderdas zweite Mal besucht wird. Wird er das erste Mal besucht, so muss er auseinem Dreieckselement vom Eingabestapel ausgelesen werden. Nachdem dieOperation auf diesem Knoten durchgefuhrt wurde, mussen die Daten diesesKnotens auf dem Farbstapel mit der entsprechenden Farbe abgelegt werden.Wird ein Knoten hingegen das zweite Mal besucht, mussen seine Daten vomentprechenden Farbstapel geholt und anschließend in ein Dreieckselementauf dem Ausgabestapel abgelegt werden.


Zu diesem Zweck wurde eine Klassifizierung der Knoten in alte und neueKnoten vorgenommen. Ein Knoten ist vom Parameter alt, wenn er bereitseinmal besucht wurde und nun das zweite Mal besucht wird. Entsprechendist er ein neuer Knoten, wenn er noch nicht besucht wurde und in dieserTraversierung das erste Mal verwendet wird. Allgemein kann man Knoten,die frisch vom Eingabestapel ausgelesen werden, als neue, Knoten die vonFarbstapeln geholt werden als alte Knoten klassifizieren.

Es ergibt sich dabei folgende Vereinfachung:Knoten, die auf der Kante liegen, durch welche die Sierpinski-Kurve in dasElement eintritt, sind immer alte Knoten, weil dieser Knoten schon beimVorgangerelement benutzt wurde. Ebenso sind Knoten auf Austrittskantenimmer neue Knoten, denn sie werden das zweite Mal direkt im Anschluß beider Behandlung des Nachbarelements benutzt. Somit verbleibt nur noch eineinziger Knoten, namlich auf derjenigen Kante, durch den die Kurve wedereintritt noch austritt, um ein Dreieck naher zu spezifizieren. Der Parameteralt/neu dieser Kante definiert den Dreieckstyp. Als Typen sind moglich Ho,Hn, Vo, Vn, Ko und Kn. Der tiefgestellte Index o steht dabei fur ”alt“ (eng-lisch ”old“), der Index n fur ”neu“ (englisch ”new“).

Fur die Kinder eines Dreiecks im Verfeinerungsbaum sind damit folgendeAbleitungsregeln gegeben (siehe Abbildung (3.9)):

Ho −→ (Vo,Ko),

Hn −→ (Vn,Ko),

Vo −→ (Ho,Ko), (3.1)

Vn −→ (Hn,Kn),

Ko −→ (Hn, Vo),

Kn −→ (Hn, Vn).

Fur jede durchgefuhrte Berechung sind diese Parameter zusammen mit demFarbparameter rot/grun wichtig, um zu wissen, von welchem Stapel dieDaten geholt werden und auf welchem Stapel sie wieder abgelegt werden


alt

alt/neu

neu

alt alt

alt/neu alt/neu

neuneu

alt/neu

alt

alt

alt/neu

alt

alt

neu

neu

alt/neu

neu

neu

alt/neu

neu

neu

neu

alt

alt

alt

Abbildung 3.9: Rekursive Definition des Dreieckstyps mit Kantenparame-tern alt/neu (Schraufstetter [12])

H V K

Abbildung 3.10: Dreieckstypen mit Farbkanten (Radzieowski [8])

mussen.

Auf vergleichbare Weise kann bei der Knotenfarbe vorgegangen werden.Auch hier ist es nicht notig, fur jeden Knoten eines Dreiecks die zugehorigeKnotenfarbe zu speichern. Kennt man namlich die Farbe des Knotens aufder Hypotenuse, so kann man bei Kenntnis des Dreieckstypus die beiden an-deren Knotenfarben ableiten (siehe Abbildung (3.10)). An dieser Stelle seinoch einmal daran erinnert, dass beim Crouzeix-Raviart-Element die Kan-tenfarbe die Knotenfarbe des daraufliegenden Knoten definiert. Dies ist nurfur die Farbkanten, also die Kanten durch die das Dreieck weder betretennoch verlassen wird, relevant. Denn die Eintritts- beziehungsweise Austritts-kanten werden nie auf einem Farbstapel abgelegt, da sie entweder gerade im


Abbildung 3.11: Rekursive Farbung der Knoten: Die Referenzfarbe wechseltzwischen den Leveln (Schraufstetter [12]).

letzten Element verwendet wurden oder im nachsten Element verwendetwerden.

Mit diesen Regeln lasst sich ebenfalls wieder die Farbe der Kinder einesDreiecks ableiten (siehe Abbildung (3.11)). Dabei bezeichne der hochgestell-te Index g die Farbe ”grun“ des Dreiecks und entsprechend ”r“ die Fraberot. Es ergeben sich folgende Regeln:

Hg −→ (V r,Kr),

Hr −→ (V g,Kg),

V g −→ (Hr,Kr), (3.2)

V r −→ (Hg,Kg),

Kg −→ (Hr, V r),

Kr −→ (Hr, V r).

Die Farbe der Dreiecke andert sich also mit jedem Level, die Kinder von

”roten“ Dreiecken sind immer ”grun“, wahrend die Kinder von ”grunen“Dreiecken immer ”rot“ sind. Die neu eingefugte Ecke liegt namlich immerauf der anderen Seite der Kurve als die bisher farbgebende Ecke.

3.3.4 Behandlung der Randkanten

Ein Problem, das bisher noch nicht betrachte wurde, ist die Behandlungvon Randkanten beziehungsweise synonym dazu Randknoten. Randknotenwerden im Gegensatz zu inneren Knoten wahrend einer Traversierung nureinmal durchlaufen. Dies liegt daran, dass Randknoten nur ein einziges Drei-


eck beruhren. Uber Randkanten findet also keine Interaktion mit anderenElementen statt, die Berechnung der Werte auf Randkanten kann in jederTraversierung in einem einzigen Schritt erfolgen. Eine Beschreibung der Re-chenvorschriften bezuglich Auslaufrandbedingungen und Wandrandbedin-gungen findet man bei Schwaiger [13].

Wie alle anderen Knoten werden auch Randknoten in alt und neu unterschie-den, allerdings nach etwas modifizierten Kriterien. So ist wahrend der ers-ten Traversierung ein Randknoten, der auf der Eintrittskante liegt, ein alterKnoten, alle anderen Randknoten sind vom Parameter neu. Andert sich dieTraversierungsrichtung werden aus alten neue Randknoten und umgekehrt.Die Randknoten mit dem Parameter neu werden wie alle neuen Knotennach ihrer Verwendung auf den jeweiligen Farbstapel abgelegt. Dort bleibensie allerding im Gegensatz zu internen Knoten bis zum Ende der Traversie-rung, da sie nur ein Mal benutzt werden. Nach jeder Traversierung werden,wie bereits angedeutet, nicht nur Ein- und Ausgabestapel vertauscht, son-dern auch der rote und grune Farbstapel. Bei der nachsten Traversierung,die in umgekehrter Richtung die Sierpinski-Kurve entlanglauft, werden die-se Randknoten als alte Randknoten behandelt und nach ihrer Verwendungwieder in die Dreiecke der Datenstruktur ”Stack-Element“ gespeichert. Nachjeder zweiten Traversierung sind also alle Daten auf dem Ein- beziehungs-weise Ausgabestapel gespeichert.

3.4 Adaptivitat

Wie in Kapitel (3.1) gefordert, soll das Gitter adaptiv verfeinerbar sein, ei-nerseits um den lokalen Fehler durch Berechnung auf kleineren Elementenbegrenzen zu konnen und andererseits um dafur nicht mit ubermaßigemSpeicher- und Zeitaufwand bezahlen zu mussen. Allerdings kann das Gitteran manchen Stellen nicht beliebig verfeinert und an anderen Stellen ver-grobert werden. Denn dadurch konnten hangende Knoten (siehe Abbildung(3.1)) entstehen, mit welchen man die bisher vorhandene Code-Strukturnicht mehr nutzen konnte. Aus diesem Grund wird bei dem vorliegenden Ver-fahren der Forderung aus Kapitel 3.1 entsprochen, keine hangenden Knotenzuzulassen. Bei den Discontinuous-Galerkin-Verfahren sind hangende Kno-ten zwar grundsatzlich erlaubt, allerdings ist nicht automatisch klar, wie die


1

2

34

Abbildung 3.12: Wird die blaue Kante durch Verfeinerung hinzugefugt, sowird Knoten 1 geloscht und die Knoten 2, 3 und 4 hinzugefugt (beziehungs-weise bei Vergroberung Knoten 2, 3 und 4 geloscht und Knoten 1 hinzu-gefugt). Der Flussterm von Knoten 1 muss aus den Flusstermen von Knoten2 und 3 berechnet werden (beziehungsweise in Flussterme fur Knoten 2 und3 aufgespalten werden).

Abbildung 3.13: Verfeinerung des Gitters (Radzieowski [8])

Interaktion zwischen einer Kante, die mehrere Nachbarn besitzt (siehe Ab-bildung (3.12)), und eben diesen Nachbarkanten mathematisch auszufuhrenware. Moglicherweise wurde dies auch betrachtliche numerische Fehler nachsich ziehen. Des weiteren muss die Verfeinerung immer Vorrang vor der Ver-groberung haben, da die Verfeinerung eine Notwendigkeit ist, um den lokalenFehler akzeptabel klein zu halten, wahrend die Vergroberung als positiverEffekt an Stellen, an denen ”nicht viel passiert“, interpretiert werden kann,da hier Rechenzeit und Speicherkapazitat geschont werden konnen.

3.4.1 Verfeinerung

Aufgrund der Forderung, keine hangenden Knoten zuzulassen, zieht dieVerfeinerung eines Dreiecks moglicherweise auch die Verfeinerung von an-deren Dreiecken nach sich. Es kann sich dadurch (bezuglich des Flachen-inhaltes) ein Dreieck nur um den Faktor 2 beziehungsweise um den Faktor


Abbildung 3.14: Die vier Falle der Verfeinerung (Radzieowski [8])

√2 bezuglich der Seitenlange von einem Nachbardreieck unterscheiden. Ist

die Hypotenuse des Dreiecks gleich der Kathete eines Nachbardreiecks, hates den halben Flacheninhalt, oder ist die Kathete gleich der Hypotenusedes Nachbardreiecks besitzt den doppelten Flacheninhalt. Damit unterliegtdie Adaptivitat also einer naturlichen Begrenzung (siehe Abbildung (3.13)).

Wie in (Abbildung (3.14)) dargestellt, gibt es vier mogliche Falle der Verfei-nerung. Je nach dem, welche Kante fur eine Verfeinerung markiert wurde,muss einer der skizzierten vier Falle angewandt werden. So beschreibt dasBild links oben eine Verfeinerung der Hypotenuse, die Bilder rechts obenund links unten die Verfeinerung jeweils einer Kathete und das Bild rechtsunten ein Verfeinerung beider Katheten. Dazu ist neben dem jeweiligen Bildder entstehende Verfeinerungsbaum skizziert. Im Bitstream bedeutet zumBeispiel die Verfeinerung einer Kathete (in Abbildung (3.14) im Bild linksunten) dass “0“, die fur das alte Blatt steht, durch die Abfolge ”10100“ersetzt wird.

3.4.2 Vergroberung

In der Vergroberung werden zwei Dreieckselemente durch das Loschen einerKante zu einem Dreieck zusammengefasst. Dies ist dann moglich, wenn dielokalen Fehler auf den beiden Dreiecken klein genug sind, so dass es mathe-matisch unbedenklich ist, sie zu einem großeren Element zusammenzufassen.Dazu sind weitere Anforderungen notig:

• Beide Dreiecke mussen im Verfeinerungsbaum das gleiche Vaterdreieck


besitzen; sonst ware die Abarbeitungsreihenfolge gemaß der Sierpinski-Kurve nicht mehr gegeben.

• Verfeinern hat Vorrang vor Vergrobern; deshalb durfen in keinem derbeiden Dreiecke zur Verfeinerung markierte Kanten vorkommen.

• Von jedem der beiden Dreiecke musssen jeweils beide Katheten zurVergroberung vorgesehen sein.

• Keine der Katheten darf Hypotenuse eines angrenzenden Dreiecks sein.Es entstunde sonst durch Loschen der gemeinsamen Kante ein hangenderKnoten an der Hypotenuse des neuen Dreiecks.

Im Bitstream wird durch die Vergroberung zweier Dreiecke zu einem Dreieckdie Zeichenfolge ”100“ durch ”0“ ersetzt.

3.4.3 Vorgehen bei der Adaption

Eine vollstandige Adaption wird durch mehrere Traversierungen vollzogen.Dabei werden die eben festgelegten Regeln fur Verfeinerung und Vergroberungbeachtet. Sie lauten in chronologischer Reihenfolge:

• Markierung der Kanten

In einer ersten Traversierung wird auf Basis eines Adaptionskriteriumsfur samtliche Elemente gepruft, welche Kanten verfeinert und welchevergrobert werden durfen. Diese Kanten werden dabei fur die jeweiligeOperation markiert. Als Adaptionskriterium wird dabei die Anderungpro Zeitschritt als Große betrachtet. Dabei existiert ein oberer und einunterer Schwellenwert. Ist der berechnete Wert großer als der obereWert, so wird die Kante zur Verfeinerung markiert, ist der berechneteWert kleiner als der untere Wert, zur Vergroberung. Anschaulich be-deutet dies, dass auf Elementen, auf denen eine starke Anderung derZustandsgroßen beobachtet wird, eine hohere Auflosung des Gittersangestrebt wird.

• Uberprufung der Konformitat

In der nachsten Traversierung wird uberpruft, welche Elemente nochmodifiziert werden mussen, um die Konformitatsbedingung (siehe (3.1))zu erfullen. Existieren nach dieser Traversierung immer noch hangende


Knoten, so werden weitere Traversierungen durchgefuhrt, bis ein kon-formes Gitter erreicht ist.

• Eigentliche Adaption

Zur eigentlichen Adaption wird nun eine letzte Traversierung durch-gefuhrt. Dabei mussen fur alle verfeinerten oder vergroberten DreieckeWerte fur die Knoten auf den neuen Kanten berechnet werden. Beivergroberten Elementen muss ein Wert geloscht und aus zwei altenWerten ein neuer Knotenwert berechnet werden. Wird ein Elementverfeinert, so mussen aus einem bestehenden Wert Knotenwerte furdie neuentstandenen Knoten generiert werden und in der richtigenReihenfolge auf die richtigen Stapel abgelegt werden.

Diese Adaptionsoperationen werden der Arbeit von Radzieowski [8] geauerbeschrieben.

Im Gesamtablauf einer Traversierung wird immer zuerst ein neuer Zeitschrittberechnet. Auf Basis der so erhaltenen neuen Werte wird dann das Adap-tionsverfahren durchgefuhrt. Der durchgefuhrte Zeitschritt wird dabei nieverworfen, selbst wenn sich herausstellt, dass auf einem zu groben Gittergerechnet wurde. Deshalb hinkt die Adaption dem verursachenden Effekt im-mer etwas hinterher. Dies wird allerdings durch das Konformitatskriteriumetwas abgemildert, da auch in der Umgebung eines sehr feinen Gitters immernoch eine relativ hohe Auflosung vorhanden ist.

Kapitel 4

Numerische Flusse

Dieses Kapitel befasst sich mit dem mathematischen Hauptteil der vorliegen-den Arbeit. Dabei werden die physikalischen und numerischen Grundlagenfur die Verwendung eines numerischen Flussterms geschaffen. Des weiterenwerden zwei Vertreter der Flussfunktionen, namlich der Lax-Friedrichs-Flussund der Roe-Fluss eingefuhrt und auf das Verfahren angewendet.

4.1 Allgemeine Herleitung

Bei der mathematischen Herleitung des Discontinuous-Galerkin-Verfahrenswurde in Kapitel 2.2.4 bereits angesprochen, aus welchem Grund das Ober-flachenintegral in Gleichung (2.15) durch einen numerischen Flussterm (siehe(2.16)) ersetzt wird. Zum einen muss damit der Austausch von Impuls undMasse (beziehungsweise Wasserhohe) zwischen den Elementen approximativbeschrieben werden. Denn durch das Fehlen von globalen Forderungen wieder Stetigkeit der Losung uber das gesamte Berechnungsgebiet ist das Ober-flachenintegral das einzige mit den Nachbardreiecken interagierende Elementin (2.15). Außerdem entgeht man dadurch der Entscheidung, welcher Wertder Zustandsvariablen, die ja beim diskreten Verfahren auf dem Rand liegen(siehe Abbildung (2.5)), zu verwenden ist.Zur Analyse dieser Unstetigkeit zwischen zwei Zellen T+ und T− betrach-tet man dabei ein sogenanntes Riemann-Problem. Das Riemann-Problem istdefiniert als ein Anfangswertproblem, in welchem zum Zeitpunkt t = 0 zweiStartzustande durch eine Unstetigkeit voneinander getrennt sind [7] (sie-he Abbildung (4.1)). Man kann dies als Dammbruchproblem zum Zeitpunkt

38

KAPITEL 4. NUMERISCHE FLUSSE 39

x

q

q−

q+

Abbildung 4.1: Das Riemann-Problem mit den zwei Initialzustanden q− undq+

t = 0 interpretieren (der Damm ist unendlich dunn). Fur Riemann-Problemewurden auch zu nichtlinearen Gleichungen, wie den Euler-Gleichungen, zubestimmten Anfangswerten exakte Losungen konstruiert (siehe Leveque [7]).

Zur Analyse der Unstetigkeit bei numerischen Verfahren wurden diese zuerstvon Godunov eingesetzt. Da allerdings die exakte Losung nicht fur alle An-fangswerte moglich oder zumindest bekannt ist, wurden bald Riemann-Losereingesetzt, die dies approximativ losen. Denn auch bei Startzustanden, beidenen die exakte Losung des Riemann-Problems bekannt ist, ist die Berech-nung dieser mit sehr großem Aufwand verbunden.

Seien nun q+ und q− die beiden Zustande links und rechts der Trennwand(siehe Abbildung (4.1)). Der numerische Fluss F ∗

e uber die Unstetigkeit e

(spater ist dies die Kante zwischen zwei Elementen) ist im Eindimensionalendefiniert als Funktion F ∗

e (q+, q−). Er soll den Ubergang einer Große q (beider Transportgleichung die Masse, bei den Flachwassergleichungen Masseoder Impuls) uber die Unstetigkeit angemessen beschreiben. Dazu werdennach Cockburn [4] folgende Eigenschaften gefordert:

1. Der numerische Fluss F ∗e (q+, q−) soll fur q+ −→ q− gegen den exakten

Fluss F (q) konvegieren (Lipschitz-Bedingung).

2. F ∗e (q+, q−) soll nicht wachsend sein in q+. Wenn q+ also großer wird

(und damit der Unterschied zwischen q+ und q− kleiner), dann sollF ∗

e (q+, q−) nicht anwachsen.

3. F ∗e (q+, q−) soll entprechend nicht sinkend sein in q−.


q−

q+

q

x

q

x

q−

q+

q*

Abbildung 4.2: Ausbreitung der Schockfronten fur die eindimensionalenFlachwasser-Gleichungen nach Leveque [7]

Die letzen beiden Bedingungen werden Monotonie-Bedingungen genannt.

Fur die in dieser Arbeit zu losenden Flachwasser-Gleichungen bedeutet dies,dass der numerische Fluss uber eine Kante e die ”Fortpflanzung“ der jewei-ligen Variable ξ, ξvx, oder ξvy an das nachste Element beschreibt. An dieserStelle ist zu erwahnen, dass bei Berechnungen in mehreren Raumdimensio-nen als numerischer Fluss nur der Fluss in Normalenrichtung zur Kante indie Rechnung einfließt. Es ist also immer die Anderung senkrecht zur Kantee gemeint, wenn im folgenden der numerische Fluss erwahnt wird.

4.2 Schockwellen

Zum besseren Verstandnis der in den nachsten beiden Unterkapiteln folgen-den Herleitung zweier spezieller numerischer Flusse werden hier kurz einigeInterpretationsmoglichkeiten aus der Theorie der hyperbolischen Differenti-algleichungen zusammengefasst:

Wird ein Riemann-Problem betrachtet, so sind zum Beginn der Berechnungdie zwei Zustande q+ und q− unstetig voneinander getrennt. Zum Zeitpunktt = 0 wird diese Trennwand beseitigt, und die beiden Levels beginnen sicheinander anzugleichen (siehe Bild (4.2)). In der Theorie der linearen hyper-bolischen Erhaltungsgleichungen spricht man von Schockwellen. Nach Leve-que [7] existieren fur die eindimensionalen Flachwassergleichungen bei linea-risierter Betrachtungsweise zwei verschiedene Arten von Schockwellen, dieals Wellenfront durch das Wasser laufen:

• Die echten Schockwellen

Diese breiten sich in Bild (4.2) von links nach rechts aus. Anschaulich


x

t

λ1λ2

q− q

+

q*

Abbildung 4.3: Die Zustande q−, q∗ und q+ getrennt durch Unstetigkeitenuber die Charakteristiken λ1 und λ2.

schiebt sich eine Wasserwand, die Schockfront, von links nach rechts.

• Die Verdunnungswelle

Im Bild wandert die Verdunnungswelle von der q-Achse nach links undverringert die Wassermenge auf der linken Seite.

Die Eigenwerte λ1 und λ2 der Jacobi-Matrix des Flussterms F (q) bestim-men die Geschwindigkeiten, mit der sich die beiden Wellen durchs Wasserfortpflanzen (sehe Leveque [7]).Bei den eindimensionalen Flachwasser-Gleichungen sind dies:

λ1 = u +√

gh, λ2 = u−√

gh.

Die beiden Schockwellen konnen im Phasenraum als charakteristische Funk-tionen dargestellt werden (siehe Abbildung (4.3)). Charakteristische Funk-tionen oder kurz Charakteristiken einer hyperbolischen Differentialgleichungsind Funktionen, entlang welchen die Losung der Differentialgleichung kon-stant ist. Die Zustande links beziehungsweise rechts der beiden Charak-teristiken λ1/2 sind die Ausgangszustande, die mit Fortschreiten der Zeitimmer weiter vom Ort der ursprunglichen Unstetigkeit verdrangt werden.Der Zustand q∗ verbindet die beiden anderen Zustande, ist aber uber dieSchockfronten unstetig.

4.3 Der Lax-Friedrichs-Fluss

In Gleichung (2.16) wurde das Oberflachenintegral aus (2.15) durch einennumerischen Fluss ersetzt. Fur die nachsten beiden Unterkapitel wird nun


Tk kT’

e

q− q+

Abbildung 4.4: Die zur Berechnung des Lax-Friedrichs-Flusses notwendigenKnoten q− und q+ auf der Kante e.

der numerische Fluss uber eine Kante e eines Dreiecks Tk betrachtet. Dabeiseien q− Zustandsvariablen des Knotens auf der Kante, welcher zum aktu-ellen Dreieck Tk gehort, und q+ die Zustandsvariablen des NachbardreiecksTk′ auf eben dieser Kante.

Der bekannteste und wegen seiner Einfachheit wohl meistverwendete nu-merische Fluss, der die vorangegangenen Bedingungen auch fur nichtlineareErhaltungsgleichungen erfullt, ist der Lax-Friedrichs-Fluss. Seine Gleichunglautet nach Cockburn [4]:

(F · n)∗ =12(F(q−) + F(q+)) +

12λ(q− − q+). (4.1)

Dabei werden der lokale und der globale Lax-Friedrichs-Fluss unterschieden:Wahrend beim lokalen Lax-Friedrichs-Fluss

λlLF = maxi|λi| |λi ist Eigenwert von ∂F (q−) oder ∂F (q+) (4.2)

der betragsmaßig großte Eigenwert entlang der Kante e ist, so ist λgLF beimglobalen Lax-Friedrich-Fluss der betragsmaßig großte Eigenwert entlang ei-ner Kante im ganzen Berechnungsgebiet:

λgLF = maxi|λi| |λi ist Eigenwert von ∂F (q)entlang aller Kanten e .

(4.3)Anschaulich wird also die Ausbreitung der durch die Unstetigkeit entste-henden Schockwelle (echter Schock oder Verdunnungswelle) mit der großt-


Abbildung 4.5: Verschmieren der Lax-Friedrichs-Losung (blaue, gestrichelteLinie) im Vergleich zur exakten Losung (rote, durchgezogene Linie).

moglichen Ausbreitungsgeschwindigkeit entlang dieser Kante oder mit dergroßten Ausbreitungsgeschwindigkeit auf dem gesamten Berechnungsgebietabgeschatzt. Dadurch entsteht beim Lax-Friedrichs-Fluss eine kunstliche,zu hohe Dispersion. Es werden folglich durch die zu schnelle Ausbreitungscharfe Schockkanten ”verschmiert“ (skizziert in Abbildung (4.5)), was zueiner ”Verwasserung“ der Losung fuhrt. Dies kann durch eine Verringerungdes Faktors λ gemindert werden, allerdings geht dies zulasten der Stabilitatdes Losers, die bei Anwendung von Gleichung (4.2) oder (4.3) gegeben ware(siehe Schwanenberg [14]).

Ausgeschrieben lautet der Lax-Friedrichs-Fluss:

(F · n)∗e =12

(n ·

(ξv)−(ξvxv + 1

2gξ2ex)−(ξvyv + 1

2gξ2ey)−

+ n ·

(ξv)+(ξvxv + 1

2gξ2ex)+(ξvyv + 1

2gξ2ey)+

)

+12λ( ξ−

(ξvx)−(ξvy)−

− ξ+

(ξvx)+(ξvy)+

) (4.4)

Dabei sollen wieder die Indizes (·)− die Werte im Dreieck Tk und (·)+ dieWerte im Nachbardreieck T

′k bezeichnen. λ wird dabei nach Hesthaven et al.

[6] folgendermaßen berechnet:

λlLF = max(x,y)∈e

∥∥∥∥∥(

vx(x, y)vy(x, y)

)∥∥∥∥∥+√

gξ(x, y), (4.5)


x

t

λ1

λ2

Abbildung 4.6: Nichtlineare Charakteristiken λi fur nichtlineare hyperboli-sche Erhaltungsgleichungen

λgLF = maxe∈S

Tk

max(x,y)∈e

∥∥∥∥∥(

vx(x, y)vy(x, y)

)∥∥∥∥∥+√

gξ(x, y). (4.6)

Fur die Transportgleichung (2.1) sind die Eigenwerte und damit der Faktorλ in Gleichung (4.1) konstant, es gilt:

λlLF = λgLF =

∥∥∥∥∥(

vx

vy

)∥∥∥∥∥ , (4.7)

und damit fur den Lax-Friedrichs-Fluss:

(F·n)∗e =12(n·

((ξvx)−(ξvy)−

)+n·

((ξvx)+(ξvy)+

))+

12

∥∥∥∥∥(

u

v

)∥∥∥∥∥ (ξ−−ξ+). (4.8)

4.4 Der Fluss nach Roe

Fur nichtlineare hyperbolische Erhaltungsgleichungen wie die Flachwasser-gleichungen sind die Charakteristiken nicht mehr wie in Abbildung (4.2)lineare, sondern eben nichtlineare Funktionen. Im Eindimensionalen werdensie durch die Eigenwerte λi = v(x, t) ±

√gξ(x, t) der Jacobi-Matrix des

Flussterms F (q) beschrieben, hangen also von der Geschwindigkeit und derWellenhohe am jeweiligen Ort und zum jeweiligen Zeitpunkt ab. In der (x, t)-Ebene sind sie daher keine Geraden mehr, sondern nichtlineare Funktionen(siehe Abbildung (4.6)).


4.4.1 Herleitung

Analytische Losungen fur solche nichtlinearen Probleme sind nur fur be-stimmte Problemklassen (wie bereits erwahnt die Euler-Gleichungen) be-kannt und außerdem kompliziert zu berechnen, was sie ungeeignet fur einnumerisches Verfahren macht. Deshalb wird das Problem bei approximativenRiemann-Losern durch ein vereinfachtes Problem ersetzt, das man mathe-matisch gut im Griff hat.

Als solche Verfeinfachung fuhrte Roe in [10] die nach ihm benannte Roe-Linearisierung ein. Dabei wird zur Bestimmung des numerischen Flusses(des Kernstucks eines Riemann-Losers) angenommen, dass der linearisiertenumerische Fluss F∗ durch eine lineare Abbildung der Ausgangsgroßen q−

und q+ dargestellt werden kann.

F∗ = A(q−, q+)q∗ (4.9)

Dabei ist q∗ ein Zwischenzustand im Bereich zwischen q− und q+. Es giltnun eine solche Matrix A zu konstruieren, die nach Roe [10] außerdem nochverschiedenen Anforderungen genugen muß:

1. A(q−, q+) −→ ∂F∂q fur q− −→ q+,

2. A(q−, q+)(q− − q+) = F(q−)− F(q+) fur alle q−, q+,

3. A soll eine reell diagonalisierbare Matrix mit paarweise verschiedenenEigenwerten sein (und damit eine Basis aus Eigenvektoren besitzen).

Bedingung (1) garantiert dabei, dass Riemann-Probleme ohne Unstetigkeitexakt gelost werden, wahrend Bedingung (2) die Erhaltungseigenschaft desVerfahrens garantiert. Es werden also die Großen ξ, ξvx und ξvy, die alleder Massenerhaltung beziehungsweise der Impulserhaltung genugen, kon-serviert. Um das linearisierte Ersatzproblem losen zu konnen, ist Bedingung(3) notig. Des weiteren geht mit Bedingung (2) folgende Eigenschaft einher:Es wird angenommen, dass fur die Schocklosung mit Schockgeschwindigkeits, also einem Eigenwert von A, der exakten Gleichung zu den Anfangs-zustanden q− und q+ gilt:

F(q−)− F(q+) = s(q− − q+), s ∈ R (4.10)


Dies ist fur lineare Gleichungen erfullt, fur nichtlineare Gleichungen bedeu-tet diese Annahme, dass der Schock von einer einzigen Welle dominiert wird,und nicht durch mehrere interagierende Wellen erzeugt wird. Bedingung (2)garantiert, dass (q− − q+) ein Eigenvektor von A zum Eigenwert s von A

ist. Es fuhren also das exakte Modell zur selben Schocklosung wie das linea-risierte Modell.

Da A nun nach Bedingung (3) eine Basis aus den Eigenvektoren ri besitzt,kann der Startzustand als Linearkombination dieser Eigenvektoren geschrie-ben werden:

q− − q+ =∑

i

αiri.

Nach Hesthaven te al. [6] lost

q∗ = q− +∑λi≤0

αiri = q+ −∑λi≥0

αiri (4.11)

das linearisierte Riemann-Problem. Damit gewinnt man nun aus den Glei-chungen (4.9) und (4.11) die vollstandige Formel fur den Roe-Fluss:

Aq∗ =12ne ·

(F(q−) + F(q+)

)+

12|A| (q− − q+). (4.12)

Dabei ist|A| = S |Λ|S−1, wo A = SΛS−1.

4.4.2 Berechnung der Matrix A

Es bleibt jetzt noch zu klaren, wie genau die Matrix A bestimmt wird. LautBedingung (2) muss A den Sprung in den Zustandsvariablen q− − q+ aufden Sprung in den exakten Flusstermen F(q−) + F(q+) abbilden.

Als erste Herangehenweise zur Gewinnung von A setzt man den Sprung imFlussterm als Integral an. Man nimmt dabei an, dass q(β) linear von denbeiden Initialzustanden q− und q+ abhangt:

F(q−)− F(q+) =∫ 1

0

dF(q(β))dβ

dβ =∫ 1

0

dF(q(β))dq

dq

dβdβ (4.13)


Mit der linearen Abhangigkeit

q(β) = q+ + β(q− − q+), β ∈ [0; 1]

folgt daraus fur Gleichung (4.13):

F(q−)− F(q+) = A(q− − q+), mit A =∫ 1

0

dF(q(β))dq

dβ. (4.14)

Das Integral zur Bestimmung von A analytisch zu berechnen stellt sich alsschwierige Aufgabe heraus, da F bei den Flachwasser-Gleichungen nichtline-ar von q abhangt. Als Alternative konnte man die Integrale durch Quadraturberechnen. Allerdings fuhrt dies zu einem hohen Rechenaufwand. Zum an-deren konnte dies negative Auswirkungen auf die Stabilitat des Verfahrenshaben.

Aus diesem Grund wird die Matrix A im Folgenden nach einer Idee kon-struiert, die P. L. Roe in [10] exemplarisch fur die Euler-Gas-Gleichungeneingefuhrt hat:Nach Bedingung (2) sucht man sich eine Matrix, die nur abhangig von denInitialzustanden q− und q+ ist und q− − q+ auf F(q−) − F(q+) abbildet.Vereinfachend sei nun ∆(·) = (·)− − (·)+ definiert. Nun zerlegt man dieFlussterm in den Flachwassergleichungen

divF(q) = div

ξv

ξvxv + 12gξ2ex

ξvyv + 12gξ2ey

=∂

∂xf(q) +

∂

∂yg(q) :=

:=∂

∂x

ξvx

ξv2x + 1

2gξ2

ξvxvy

+∂

∂y

ξvy

ξvxvy

ξv2y + 1

2gξ2

(4.15)

Gesucht ist nun eine Matrix A1, die ∆q auf ∆f(q) beziehungsweise A2, die∆q auf ∆g(q) abbildet. Um diese Abbildungen zu beschreiben wird nun einParameter-Vektor

w =

w1

w2

w3

=

ξ12

ξ12 vx

ξ12 vy

(4.16)


eingefuhrt, mit dessen Hilfe sich die Matrizen Af und Ag leicht konstruierenlassen. Dazu werden noch Rechenregeln fur die ∆-Schreibweise benotigt, dieleicht von Hand verifiziert werden konnen:

• ∆(p + q) = ∆p + ∆q,

• ∆(pq) = p∆q + q∆p,

• ∆(1q ) = ∆q

(q∗)2.

Dabei ist q = 12(q− + q+) das arithmetische und q∗ = √q−q+ das geometri-

sche Mittel.

Mit diesen Rechenregeln lasst sich leicht verifizieren, dass

∆q = B∆w, ∆f = C∆w, und ∆g = D∆w (4.17)

ist. Hieraus ergeben sich folgende Matrizen:

B =

2w1 0 0w2 w1 0w3 0 w1

, C =

w2 w1 02w2

1w1 2w2 00 w3 w2

, (4.18)

D =

w3 0 w1

0 w3 w2

2gw21w1 0 2w3

(4.19)

Mit ihnen lasst sich nun eine Abbildung Af von ∆q nach ∆f beziehungsweiseAg von ∆q nach ∆g konstruieren:

∆f = CB−1∆q, ∆g = DB−1∆q. (4.20)

Dabei ist

CB−1 =

w22w1

12 0

− w22

2w12 + 2gw2

13w22w1

0

−w2w3

2w12

w32w1

w2w1

, (4.21)

DB−1 =

w32w1

0 12

w2w3

2w12

w3w1− w2

2w1

−w32

2w12 + 2gw2

1 0 3w32w1

(4.22)


Diese Matrizen mussen nun laut (4.12) in ihre Eigenwerte und Eigenvektorenzerlegt werden, um |A| und damit den numerischen Fluß (F ·n)∗ berechnenzu konnen. Diese Eigenwerte λi und die zugehorigen Eigenvektoren ri derMatrix Af = CB−1 lauten:

λ1 =w2

w1, λ2 =

w2

w1+√

gw21, λ3 =

w2

w1−√

gw21, (4.23)

r1 =

001

, r2 =

1

w2w1

+ 2√

gw21

w3w1

, r3 =

1

w2w1− 2√

gw21

w3w1

. (4.24)

Fur die Matrix Ag = DB−1 lauten die Eigenwerte µi und die zugehorigenEigenvektoren si:

µ1 =w3

w1, µ2 =

w3

w1+√

gw21, µ3 =

w3

w1−√

gw21, (4.25)

s1 =

010

, s2 =

1−w2

w1

w3w1

+ 2√

gw21

, s3 =

1−w3

w1

w3w1− 2√

gw21

. (4.26)

Aus diesen Komponenten kann nun der vollstandige Fluss nach der Vor-schrift (4.12) berechnet werden.

Um den letzten Summanden aus Formel (4.12) leichter berechnen zu konnenwird nun darauf zuruckgegriffen, dass die Matrizen Af = CB−1 beziehungs-weise Ag = DB−1 nach Forderung 3 linear unabhangige Eigenvektoren ri be-ziehungsweise si und folglich eine Basis aus Eigenvektoren besitzen. Deshalbkann der Sprung q− − q+ in den Zustandsvariablen als Linearkombinationder Eigenvektoren dargestellt werden. Es gilt

q− − q+ =3∑

i=1

αiri, beziehungsweise q− − q+ =3∑

i=1

βisi, (4.27)

fur Sprunge entlang der x- beziehungsweise y-Achse.


Seien nun u := w2w1

, v := w3w1

, a :=√

gw21 und ∆q =

∆q1

∆q2

∆q3

der Sprung in

den Zustandsvariablen q.Dann lassen sich die Koeffizienten αi bestimmen als

α1 = ∆q3 −∆q1v,

α2 =14a

∆q2 +(12− u

4a

)∆q1,

α3 = ∆q1 − α2. (4.28)

Die Koeffizienten βi ergeben sich analog zu

β1 = ∆q2 + ∆q1u,

β2 = ∆14a

∆q3 +(12− v

4a

)∆q1,

β3 = ∆q1 − β2. (4.29)

Mit diesen Koeffizienten kann nun Formel (4.12) fur einen Fluss entlangder x-Achse geschrieben werden als

Aq∗ =12ne ·

(F(q−) + F(q+)

)+

12|A| (q− − q+) =

=12ne ·

(F(q−) + F(q+)

)+

12

3∑i=1

|λi|αiri. (4.30)

Dabei sind mit λi die Eigenwerte und mit ri die Eigenvektoren der MatrixAf bezeichnet.

Kapitel 5

Details zur Implementierung

In diesem Kapitel soll nun geklart werden, wie die oben aus mathematischerSicht eingefuhrten neuen Flussterme implementiert wurden. Zuerst wird be-schrieben, welche Anderungen sich fur das ganze Programmgerust ergaben.Anschließend wird die gewahlte Implementierung der Flussfunktionen unddabei insbesondere des Faktors λ aus Gleichung (4.5) beziehungsweise (4.6)beim Lax-Friedrichs-Fluss, sowie der Matrix A (siehe Gleichung (4.12)) beimRoe-Fluss, und dabei zu beachtende Feinheiten geklart.

5.1 Grundlegende Annahmen

Tk

e1

n−

e2

n−

e3

n−

Abbildung 5.1: Flusse in Richtung der außeren Normalen nei .

Der numerische Fluss (F · n)∗ beschreibt den Austausch der Zustands-großen q uber die Grenzen der Dreieckselemente hinweg. Ein Element Tk

interagiert dabei mit samtlichen Elementen Tk′ , die eine gemeinsame Kante

51

KAPITEL 5. DETAILS ZUR IMPLEMENTIERUNG 52

mit Tk besitzen.

Wie in Gleichung (4.15) beschrieben, kann der numerische Fluss in zweiFlusse in die jeweilige Koordinatenrichtung zerlegt werden. In unserem Ver-fahren wird als Fluss uber die Elementgrenzen immer der Fluss in die Nor-malenrichtung der jeweiligen Kante betrachtet (siehe Abbildung (5.1)). Umdies zu gewahrleisten, soll im Folgenden daher bei der Berechnung des nu-merischen Flusses (F · n)∗ der Flussterm F(q) immer mit der Normalen deraktuell verwendeten Kante multipliziert werden:

(F · n)∗ =12(n · F(q−) + n · F(q+)) +

12α(q− − q+) =

=12

(nx

(f(q−) + g(q−)

)+ ny

(f(q+) + g(q+)

))+

12α(q− − q+).

Dabei ist n =

(nx

ny

)der außere Normalenvektor der aktuell verwendeten

Kante.

5.2 Anpassung des Algorithmus und der Daten-

struktur

5.2.1 Der modifizierte Lax-Friedrichs-Fluss

Ein Riemann-Loser wird hauptsachlich durch die Wahl des numerischenFlusses (F · n)∗ bestimmt. So wurde im bisherigen Verfahren der soge-nannte modifizierte Lax-Friedrichs-Fluss verwendet. Dieser modifizierte Lax-Friedrichs-Fluss berechnet sich nach der Formel

(F · n)∗ =12n · (F(q−) + F(q+)) +

12α(q− − q+). (5.1)

Der Faktor α wurde dabei als konstant angenommen und je nach Notwendig-keit experimentell so bestimmt, dass die Stabilitat des Losers gewahrleistetwar. So wurde zum Beispiel zur Losung der Flachwasser-Gleichungen beikonstanten Ansatzfunktionen α = 2.0 gesetzt.

Dieser Fluss erfullt, wie leicht zu uberprufen ist, die Bedingungen, die nachKapitel (4.1) an einen numerischen Fluss gestellt werden. So ist dieser Fluss


nicht wachsend im ersten Argument (Bedingung (2)), nicht fallend im zwei-ten Argument (Bedingung(3)) und erfullt die Lipschitz-Bedingung (1). Fallsnun q− gegen q+ strebt, so konvergiert dieser numerische Fluss F∗(q−, q+)auch gegen den exakten Fluss F(q). Angewendet auf das hier vorgestellteVerfahren hat dieser Fluss den Vorteil, dass die vollstandige Berechnung desneuen Flusses in zwei Schritten vollzogen werden kann. So wird bei innerenKanten beim ersten Besuch der Kante (im Dreieck Tk) der Anteil des Flussesberechnet, der durch dieses Dreieck Tk in den Flussterm eingeht:

flownew←− 12(n−e · F(q−) + αq−).

Dabei ist n−e die außere Normale des Dreiecks Tk an der Kante e. Wird dieKante nun im Dreieck Tk

′ zum zweiten Mal besucht, so wird der erste Teildes Flusses, der im Dreieck Tk berechnet wurde, ubergeben und der Anteilaus dem Element Tk′ addiert:

flownew←− flownew− 12(n+

e · F(q+) + αq+).

Auch hier ist n+e die außere Normale des Dreiecks Tk′ an der Kante e und

es gilt n+e = −n−e . Damit kann der vollstandige numerische Fluss uber eine

Kante gemaß (5.1) also nach einer einfachen Regel berechnet werden.

5.2.2 Der lokale Lax-Friedrichs-Fluss

Im Gegensatz zur bisherigen Implementierung wird nun der lokale Lax-Fried-richs-Fluss als numerischer Fluss verwendet. Er berechnet sich nach derFormel (4.1):

(F · n)∗ =12n ·(F(q−) + F(q+)) +

12λ(q− − q+).

Auch hier ist n wieder der außere Normalenvektor der verwendeten Kante.Der Faktor λ ist nun allerdings keine Konstante mehr, sondern berechnetsich wie in Formel (4.5) gegeben als großter Eigenwert der Jacobi-Matrix∂F(q)

∂q .

Die eben vorgestellte Berechnungsvorschrift fur den modifizierten Lax-Fried-richs-Fluss kann nun nicht mehr verwendet werden, da λ erst nach demzweiten Besuch der Kante berechnet werden kann. Denn λ ist das Maximum


der großten Eigenwerte der Jacobi-Matrix beider Seiten der Kante e:

λ = maxλ1, λ2,

wobei λ1 der großte Eigenwert der Kante e im Element Tk und λ2 der großteEigenwert der Kante e im Element Tk′ ist.

Aus diesem Grund musste die Datenstruktur geandert werden:Der Wert von λ1, also der Eigenwert, der beim ersten Besuch der Kante be-rechnet wurde, muss gespeichert werden und beim zweiten Besuch ubergebenwerden, um den Faktor λ = maxλ1, λ2 vollstandig zu berechnen. Aller-dings muss nur die Datenstruktur der Farbstapel geandert werden. Fur dieDaten auf dem Ein- oder Ausgabestapel ist es nicht notwendig, die Eigenwer-te zu speichern, da diese nur zur Flussberechnung benotigt werden, welchenach dem zweiten Besuch abgeschlossen ist. Auch werden die Eigenwerte injedem Zeitschritt neu berechnet, so dass es nicht notig ist, diese langer alsbis zur vollstandigen Flussberechnung zu speichern.

Da der maximale Eigenwert λ beim ersten Besuch einer Kante noch nichtbekannt ist, kann vom zu berechnenden neuen Fluss aber noch nicht derkomplette erste Anteil berechnet werden. Deshalb wird nun beim erstenBesuch einer Kante nur der Flusswert

flownew←− 12n−e · F(q−)

berechnet. Beim zweiten Besuch muss dann zuerst λ vollstandig berechnetwerden und anschließend die vollstandige Flussberechnung erfolgen:

flownew←− flownew +12n+

e · F(q+) +12λ(q− − q+).

Auch hier sind n−e und n+e wieder die außeren Normalenvektoren auf der

Kante e der Dreiecke ”-“ und ”+“. Die Zustandsvariablen q− wird dabeierst beim zweiten Teil der Berechnung benotigt, obwohl sie die Werte derZustandvariablen des ersten Elements Tk beschreiben. Also mussen die Wer-te q− ebenso wie λ1 mit auf den Farbstapel gelegt werden, damit sie beimzweiten Besuch der Kante zur vollstandigen Berechnung zur Verfugung ste-hen. Ein schematischer Ablauf der Flussberechnung fur den lokalen Lax-


Friedrichs-Fluss ist in Algorithmus 5.1 abgebildet.

Algorithmus 5.1: Berechnung des lokalen Lax-Friedrichs-Flusses

% Erster Besuch der Kante e im Dreieck Tk

Berechne: λ1

Berechne: flownew = 12n

−e · F(q−)

Lege q−,flownew, λ1 auf Farbstapel...% Zweiter Besuch der Kante e im Dreieck Tk′

Hole q−,flownew, λ1 vom FarbstapelBerechne: λ2

Berechne: λ = max(λ1, λ2)Berechne: flownew = −flownew + 1

2n+e · F(q+) + 1

2λ(q+ − q−)

5.2.3 Der globale Lax-Friedrichs-Fluss

Wahrend beim lokalen Lax-Friedrichs-Fluss als Parameter λ der großte Ei-genwert auf der Kante des zu berechnenden Flusses verwendet wird, wirdbeim globalen Lax-Friedrichs-Fluss der großte Eigenwert entlang einer Kanteim ganzen Berechnungsgebiet wahrend eines Zeitschrittes als λ benutzt. Esist also fur jede Kante einer Traversierung das zugehorige λ zu berechnen.Am Ende des Durchlaufs muss dann das Maximum uber diese λ gebildetwerden. Dabei tritt das Problem auf, dass der Parameter erst am Ende ei-ner Traversierung feststeht und folglich erst jetzt der Fluss uber die Kantenberechnet werden konnte. Man konnte dies durch eine zusatzliche Traversie-rung beheben, in der dann, da nun der Lax-Friedrichs-Parameter feststeht,der Fluss berechnet wird.

Um den hoheren Aufwand einer zusatzlichen Traversierung zu vermeiden,wurde jedoch folgene Vorschrift gewahlt:Bei jedem Kurvendurchlauf wird der globale Lax-Friedrichs-Parameter λ be-rechnet. Dieser wird anschließend in der nachsten Traversierung verwendet,um die Flussterme zu berechnen. Es existieren im Algorithmus also zweiLax-Friedrichs-Parameter:

• Ein aktueller Parameter λactual, mit dem in der aktuellen Traversie-


rung die numerischen Flusse uber die Kanten berechnet werden, und

• ein Parameter λnext, in dem das Maximum uber alle neu berechne-ten λ gebildet wird. Dieser wird im nachsten Zeitschritt als λactual

verwendet.

Der Vorteil dieser Strategie ist, dass hierbei in jedem Zeitschritt der Para-meter λactual a priori bekannt ist. Die Informationen einer bereits besuchtenKante mussen also, wie beim modifizierten Lax-Friedrichs-Fluss, nicht wei-ter mitgeschleppt werden, um den vollstandigen Flussterm zu berechnen.Algorithmus 5.2 beschreibt das Vorgehen bei der Berechnung des globalenLax-Friedrichs-Flusses.

Algorithmus 5.2: Berechnung des globalen Lax-Friedrichs-Flusses


Berechne: λ

Berechne: λnext = max(λ, λnext)Berechne: flownew = 1

2n−e · F(q−) + 1

2λactualq−

Lege flownew auf den Farbstapel...% Zweiter Besuch der Kante e im Dreieck Tk′

Hole flownew vom FarbstapelBerechne: λ

Berechne: λnext = max(λ, λnext)Berechne: flownew = −flownew + 1

2n+e · F(q+) + 1

2λactualq+

...% Am Ende einer Traversierungλactual = λnext

λnext = 0.0


5.3 Details zur Berechnung des Lax-Friedrichs -

Parameters λ

Wie wird aber nun der Lax-Friedrichs-Parameter λ genau berechnet?Zur Klarung dieser Frage muss unterschieden werden nach den verschiedenenGleichungen und dem Raum der Ansatzfunktionen.

5.3.1 Die Transportgleichung

Fur die Transportgleichung (2.1) ist der Parameter λ sehr leicht zu be-rechnen. Wie bereits erwahnt soll nach Cockburn [4] fur λ der betragsmaßiggroßte Eigenwert der Jacobi-Matrix ∂F

∂q verwendet werden. Im Fall der Trans-portgleichung bei konstanter Wassertiefe ergibt sich als Eigenwert:

λ =

(vx

vy

)= const. (5.2)

Da bei der Transportgleichung diese Geschwindigkeiten konstant sind, bleibthier nichts mehr zu tun. Der Lax-Friedrichs-Parameter λ ist, wie in Glei-chung (5.2) zu sehen, folglich sowohl fur den lokalen, als auch fur den glo-balen Lax-Friedrichs-Fluss, unabhangig vom Raum der Ansatzfunktionen,konstant.

5.3.2 Die Flachwasser-Gleichungen mit konstanten Ansatz-

funktionen

Dies ist bei den Flachwasser-Gleichungen nicht mehr der Fall. Nach Hestha-ven et al. [6] ist der Lax-Friedrichs-Parameter als

λ = max∥∥∥∥∥(

vx

vy

)∥∥∥∥∥+√

gξ

(5.3)

zu wahlen. Je nach dem, ob man den globalen oder den lokalen Lax-Friedrichs-Fluss verwendet, muss dabei das Maximum uber alle Kanten oder nur uberdie aktuelle Kante gebildet werden.

Wird das Discontinuous-Galerkin-Verfahren mit konstanten Ansatzfunktio-nen verwendet, so sind auf jedem Element die Werte fur die Zustandsvaria-


blen konstant.

q =

ξ

ξvx

ξvy

= const (5.4)

Fur die Berechnung des Parameters λ bedeutet dies, dass λeTk

= const. DieIndizes Tk und e sollen dabei verdeutlichen, dass es sich um den λ-Wertdes Elements Tk zur Kante e handelt. Ebenso ist der Parameterwert λe

Tk′

des Nachbarelements Tk′ zur selben Kante e konstant. Zur vollstandigenBerechnung des λ-Werts der Kante e muss nun das Maximum aus den λ derbeiden Seiten der Kante e gebildet werden:

λe = maxλeTk

, λeT

k′ . (5.5)

Um den vollstandigen lokalen oder globalen Lax-Friedrichs-Fluss zu erhaltenmuss nun nur noch λ in den Algorithmus 5.1 beziehungsweise 5.2 eingesetztwerden.

5.3.3 Die Flachwasser-Gleichungen mit linearen Ansatzfunk-

tionen

Deutlich schwieriger gestaltet sich die Berechung von λ fur die Flachwasser-Gleichung mit linearen Ansatzfunktionen. Grundsatzlich mussen dabei eben-so wie bei den konstanten Ansatzfunktionen die Maxima auf beiden Seiteneiner Kante berechnet werden (siehe Gleichung (4.4)). Das Maximum ausden erhaltenen beiden Werten wird dann als λ fur die Flussberechnung ver-wendet.

Allerdings sind die Parameter λTkund λT

k′ sehr viel schwerer zu berechnen.

Denn auf jedem Element Tk sind die Zustandsvariablen (5.4) nun lineareFunktionen:

q =

ξ(x, y)ξvx(x, y)ξvy(x, y)

∈ P31. (5.6)


Damit ergibt sich als Problemstellung folgende Aufgabe:

λeTkbzw.T

k′ = max

(x,y)∈e

∥∥∥∥∥(

vx(x, y)vy(x, y)

)∥∥∥∥∥+√

gξ(x, y), e ∈ Tkbzw.Tk′ .

(5.7)Im Folgenden werden die Argumente (x, y) sowie die Indizes e und Tk/Tk′

der Ubersicht halber weggelassen. Es handelt sich jedoch, solange nicht an-ders erwahnt, um die Berechnung eines maximalen λ auf der Kante e ineinem Dreieck Tk.

Der Versuch, analytisch ein Maximum von (5.7) entlang einer Kante e zuberechnen, ist nicht sinnvoll, da sich einerseits die Losung als kompliziert er-weist, und andererseits damit wertvolle Rechenzeit vergeudet wurde. Dennlaut Cockburn [4] ist die Qualitat der Losung mit steigendem Polynom-grad der Ansatzfunktionen immer weniger verwendeten numerischen Flussabhangig. Aus diesen Grunden wurden drei verschiedene Losungsansatzeimplementiert, um das Maximum aus (5.7) abzuschatzen.

Abschatzung mit der Maximumsnorm

Im dieser ersten Abschatzung wurde fur Gleichung (5.7) als Norm die Ma-ximumsnorm gewahlt. Damit ergibt sich fur (5.7):

λ = max(x,y)∈e

∥∥∥∥∥(

vx(x, y)vy(x, y)

)∥∥∥∥∥max

+√

gξ(x, y)

=

= max(x,y)∈e

max|vx(x, y)|, |vy(x, y)|+√

gξ(x, y). (5.8)

Im implementierten Verfahren werden bei der Abschatzung mit der Maxi-mumsnorm dabei einfach die λ-Werte an den beiden Endpunkten der Kantee berechnet, und λ mit deren Maximum als Wert belegt. Sei nun vxl

derGeschwindigkeits-Wert in x-Richtung am linken Ende und vxr der Geschwin-digkeits-Wert in x-Richtung am rechten Ende der Kante e, ebenso vyl

, vyr , ξl

und ξr. Weiter seienvxm , vym , ξm die entsprechenden Werte in der Mitte derKante. Dann wird im Algorithmus λ vereinfachend folgendermaßen berech-net:

λ = maxmaxvxi , vyi+√

gξi, i ∈ l, m, r (5.9)


Abschatzung mit der Dreiecksungleichung

Als zweite Moglichkeit der Abschatzung von (5.7) wird die Dreiecksunglei-chung verwendet:

λ = max(x,y)∈e

∥∥∥∥∥(

vx(x, y)vy(x, y)

)∥∥∥∥∥2

+√

gξ(x, y) ≤

≤ max(x,y)∈e

∥∥∥∥∥(

vx(x, y)vy(x, y)

)∥∥∥∥∥2

+ max(x,y)∈e

√gξ(x, y) (5.10)

Zur Berechnung des Maximums des ersten Summanden von Formel (5.10)ist nun lediglich die Auswertung von

√v2xl

+ v2yl

und√

v2xr

+ v2yr

notwendigund das Maximum aus diesen beiden Werten zu bilden.

Dies kann wie folgt begrundet werden:Da sowohl vx(x) als auch vy(x) lineare Funktionen in der Variablen x sind,werden die Maxima von v2

x(x) und v2y(x) auf dem Rand des betrachteten

Gebiets, also fur x = l oder x = r angenommen. Dies ist der Fall, weilv2x(x) und v2

y(x) Teilstucke einer nach oben geoffneten Parabel sind. Alsoliegt auch das Maximum von v2

x(x) + v2y(x) auf dem Rand. Dies lasst sich

schnell durch eine Fallunterscheidung zeigen:

• 1. Fall: v2xl

> v2xr∨ v2

yl> v2

yr

Da die beiden Werte an x = l großer sind als die beiden Werte anx = r, muss auch dass Maximum der Summenfunktion v2

x > v2y an der

Stelle x = l liegen.

• 2. Fall: v2xl

< v2xr∨ v2

yl< v2

yr

Analog zu Fall 1 muss hier das Maximum an der Stelle x = r liegen.

• 3. Fall: sonstIn diesem Fall liegt je ein Maximum der beiden Funktionen v2

x(x)beziehungsweise v2

y(x) an x = l und eines an x = r. Man betrachtetnun die auf Bild (5.1) eingezeichneten Sekanten v∗x und v∗y . Da v2

x undv2y konvex sind, gilt v∗x(x) ≥ v2

x(x) und v∗y(x) ≥ v2y(x) fur x ∈ [l, r].

Falls nun das Maximum von v2x(x) + v2

y(x) im Inneren des Intervalles[l, r] lage, musste das Maximum von v∗x(x)+v∗y(x) ebenfalls im innerenvon [l, r] liegen. Da v∗x(x) + v∗y(x) eine lineare Funktion ist, fuhrt dieszum Widerspruch. Also liegt das Maximum von v2

x(x)+v2y(x) auf auch


xv*

yv*

yv2

xv2

rl

x

y

Abbildung 5.2: Abschatzung von v2x und v2

y durch die Sekanten v∗x und v∗y .

in diesem Fall dem Rand.

Da die Quadratwurzel eine streng monotone Funktion ist, liegt nun auchdas Maximum von

√v2x(x) + v2

y(x) an x = l oder x = r.

Im Programmcode wird λ also folgendermaßen berechnet:

λ = max√

v2xl

(x) + v2yl

(x),√

v2xr

(x) + v2yr

(x)+√

g maxξl, ξr. (5.11)

Abschatzung durch quadratische Interpolation

Als letzte Moglichkeit zur naherungsweisen Berechnung des maximalen λ

wurde die quadratische Interpolation der Funktion λ(ξ(x), vx(x), vy(x)) ge-testet. Dabei wird λ(ξ(x), vx(x), vy(x)) an den Stellen x = l, x = m undx = r ausgewertet. Die sich ergebenden Werte werden mit λl, λm und λr

bezeichnet. Mit diesen Werten wird dann ein quadratisches Interpolations-polynom mit den Stutzstellen (1, λl), (3

2 , λm) und (2, λr) bestimmt.

Berechnet man nun die Extremstelle dieses Polynoms so erhalt man als x-Koordinate

xextr =3l − 4m + r

4(l − 2m + r)

Um den richtigen Extremwert im Intervall [l, r] zu erhalten, muss nun noch


eine Fallunterscheidung eingebaut werden:

λ =

maxλl, λextr, λr, 1 < xextr < 2,

maxλl, λr, sonst.(5.12)

5.4 Behandlung der Randkanten

Noch nicht geklart wurde bisher, wie die Randkanten im Rahmen des Lax-Friedrichs-Flusses behandelt werden. Im Gegensatz zu inneren Kanten wirdeine Randkante pro Gittertraversierung nur ein einziges mal besucht. Des-halb ist hier ein gesondertes Vorgehen notwendig, damit nicht nur der exakteFluss F(q−) berechnet wird.

Zwei verschiedene Randbedingungen, namlich die Auslaufrandbedingungund die Wandrandbedingung wurden in der Arbeit von Schwaiger [13] ein-gefuhrt. Dabei werden fur die nicht vorhandenen Flusswerte und Zustands-variablen Werte an imaginaren Punkten außerhalb des Berechnungsgebietsverwendet, die nach bestimmten Kriterien die entsprechende Randbedin-gung simulieren. So ist die Berechnung des kompletten Flusses an der Rand-kante wahrend eines einzigen Besuches moglich (und auch notig).

So werden zum Beispiel die Koordinaten des imaginaren Punkts q+, derdem regularen Punkt q− auf der Randkante gegenuber liegt, bei konstantenAnsatzfunktionen fur die Auslaufrandbedingung folgendermaßen bestimmt:

q− =

ξ−

(ξvx)−(ξvy)−

=

ξ+

(ξvx)+(ξvy)+

= q+.

Fur die Auslaufrandbedingung bei linearen Ansatzfunktionen, die Wand-randbedingung bei konstanten und linearen Ansatzfunktionen sowie einegenauere Herleitung siehe Schwaiger [13] und Aizinger [5].

Fur die Berechnung des lokalen und des globalen Lax-Friedrichs-Fluss wer-den diese imaginaren Werte, die sich aus den Randbedingungen ergeben,ebenfalls benotigt. Zur Berechnung des Parameters λ einer Randkante wirdje nach gewahltem Fluss ein λTk

nach der entsprechenden Vorschrift be-rechnet. Dann wird mit den imaginaren Werten wie bei inneren Kanten ein


Parameter λTk′ berechnet und das Maximum λ = maxλTk

, λTk′ gebildet.

Mit diesem Parameter λ wird nun der gesamte Fluss berechnet.

5.5 Details zur Berechnung des Roe-Flusses

Ahnlich wie bei der Bestimmung des Lax-Friedrichs-Flusses taucht bei derBerechnung des Roe-Flusses das Problem auf, dass der vollstandige nume-rische Fluss uber eine Kante e erst beim zweiten Besuch dieser Kante be-rechnet werden kann. Laut Formel (4.30) mussen zur Berechnung des Roe-Flusses die beiden Zustande q− und q+ bekannt sein. Aus diesen werdendann die Eigenwerte und Eigenvektoren der Roe-Matrix A berechnet, unddamit dann der zweite Summand aus (4.30) bestimmt.

Im implementierten Verfahren wird dabei folgendermaßen vorgegangen:Da wie beim Lax-Friedrichs-Fluss der Vektor der Zustandsvariablen q+ desElements auf der gegenuberliegenden Seite der Kante e benotigt wird, kannder vollstandige Roe-Fluss erst beim zweiten Besuch einer Kante berech-net werden. Deshalb wird beim ersten Besuch einer Kante lediglich der Zu-standsvariablenvektor q+ gespeichert und ubergeben. Wird die Kante e imVerlauf der Traversierung das zweite Mal besucht, so wird der Fluss gemaßden Formeln (4.23), (4.24), (4.28) und (4.30) berechnet.

Wie in Kapitel 5.1 bereits beschrieben, wird bei der Auswertung des nu-merischen Flusses nur der Fluss in Richtung der außeren Normale n derbetreffenden Kante betrachtet. Da in dem hier implementieren Verfahrenalle Elemente aus rechtwinklig-gleichschenkligen Dreiecken bestehen, sindbis auf das Vorzeichen nur 4 verschiedene Normalenrichtungen moglich. Istdie Kante e parallel zur y-Achse, so ist der Normalenvektor n = (1, 0)T . DerRoe-Fluss kann dann einfach durch Formel (4.30) berechnet werden. Ist dieKante e parallel zur y-Achse, so ist in Formel (4.30) lediglich αi durch βi,λi durch µi und ri durch si zu ersetzen.

Um Fallunterscheidungen im Programmcode zu vermeiden, wird bei der Im-plementierung ein eleganter Trick, verwendet, den Hesthaven et al. in [6]vorgeschlagen haben:Man nutzt die Rotationsinvarianz der Impulsgleichungen, also der zwei-


ten und dritten Zeile der Flachwasser-Gleichungen. Die Geschwindigkeits-impulse ξvx und ξvy werden dabei in Normalen- und Tangentenrichtungtransformiert. Dies geschieht durch eine Multiplikation des Geschwindikeit-simpulsvektors mit der Rotations-Matrix R(

ξvrotx

ξvroty

)=

(nx ny

−ny nx

)(ξvx

ξvy

). (5.13)

Dabei bezeichnet ξvrotx den Geschwindigkeitsimpuls in Normalenrichtung,

ξvroty den Geschwindigkeitsimpuls in Tangentialrichtung und n = (nx, ny)T

den außeren Normalenvektor der Kante e. Da bei der numerischen Flussbe-rechnung lediglich der Geschwindikeitsimpuls in Normalenrichtung beruck-sichtigt wird, fuhrt dies zur einer deutlichen Vereinfachung der Berechnungs-vorschrift. Außerdem kann diese Vorschrift fur jeden der 4 verschiedenenKanten-Normalenrichtungen verwendet werden. Fur die in Normalen- undTangentialrichtung gedrehten Vektoren kann der Roe-Fluss dann gemaß For-mel (4.30) berechnet werden.

Algorithmus 5.3: Berechnung des Roe-Flusses


Lege q− auf Farbstapel...% Zweiter Besuch der Kante e im Dreieck Tk′

Hole q− vom FarbstapelTransformiere die Vektoren des Geschwindigkeitsimpulses mittels derDreh-

matrix R von q− und q+ in Normal- und Tangential-KoordinatenBerechne: die Koeffizienten w− und w+ aus Gleichung (4.16)Berechne: die Eigenvektoren ri und die Eigenwerte λi der Roe-Matrix A

gemaß den Vorschriften (4.23) und (4.24)Berechne: die Koeffizienten αi der Linearkombination der Eigenvektoren

ri wie in Gleichung (4.28)Berechne: mit diesen Werten den Roe-Fluss nach Formel (4.30)Drehe die erhaltenen Fluss-Vektoren mittels R−1 zuruck in kartesische

Koordinaten

Kapitel 6

Numerische Ergebnisse

In diesem Kapitel soll nun das implementierte Verfahren anhand einigerModellrechnungen getestet werden. Dazu wird zuerst eine kurze Rechenzeit-und Speicheranalyse durchgefuhrt. Anschließend wird ein neu implementier-tes Testszenario, das radiale Dammbruch-Problem, vorgestellt. Anhand die-ses Problems werden die oben hergeleiteten verschiedenen Flussfunktionenbezuglich Stabilitat und numerischer Dispersion verglichen. Hauptaugen-merk liegt dabei auf dem Einfluss des Lax-Friedrichs-Parameters λ bezuglichdieser Kriterien.

6.1 Rechenzeit- und Speicheranalyse

Wie eingangs erwahnt, sind die Rechenzeit und der Speicherbedarf nebender Stabilitat die Hauptkriterien eines konvergenten numerischen Verfah-rens. Deshalb ist es von großer Wichtigkeit, einordnen zu konnen, inwiefernsich die Modifizierung eines Algorithmus auf diese beiden Großen auswirkt.

Fur die in dieser Arbeit am implementierten Programm vorgenommenenModifikationen wird daher eine kurze Zusammenfassung der Anderung be-zuglich dieser Kriterien gegeben. Als Bezugsgroße dient dabei der Algorith-mus mit dem bisher angewandten modifizierten Lax-Friedrichs-Fluss. Eineprazise Analyse von dessen Laufzeit findet sich in der Arbeit von Schwaiger[13].

65

KAPITEL 6. NUMERISCHE ERGEBNISSE 66

6.1.1 Speicheranalyse

Grundsatzlich hangt der Speicherbedarf des Verfahrens von der Feinheit desGitters und damit von der Berechnungsgenauigkeit ab. Fur die durch dasDiscontinuous-Galerkin-Verfahren mit konstanten Ansatzfunktionen gelostenFlachwassergleichungen mussen fur jedes Dreieckselement 3 Werte fur dieZustandvariablen q und fur die Flusse uber die 3 Kanten eines Elementsje 6 Werte (flow new und flow old) fur den numerischen Fluss gespeichertwerden, also insgesamt 21 real - oder double-Werte. Fur den linearen Ansatzsind ebensoviele Flusswerte, aber nun 3 × 3 = 9 Werte fur q zu speichern,insgesamt also 27 Werte pro Element.

Nun mussen aufgrund der veranderten Berechnungsvorschriften bei der Fluss-berechnung fur jede Kante eines Dreiecks zusatzlich noch die unvollstandigberechneten Lax-Friedrichs-Parameter λ und die Zustandsvariablen q desElements gegenuber der jeweiligen Kante gespeichert werden. Allerdings hatdies keinen signifikanten Einfluss auf den Speicherbedarf, da diese Werte nurzwischen dem ersten und dem zweiten Besuch einer Kante, also nur auf denFlusstapeln, gespeichert werden muss. Weil die Flussstapel im Vergleich zuden Ein- und Ausgabestapeln aber sehr klein sind, wird sich der Speicher-bedarf fur das Verfahren kaum andern.

6.1.2 Rechenzeitanalyse

Veranderungen in der Rechenzeit ergeben sich lediglich durch den Mehrauf-wand bei der Berechnung des Lax-Friedrichs-Parameters λ. Denn alle ande-ren Operationen, wie die Berechnung des neuen ξ-Wertes oder der Flussbe-rechnung laufen im Wesentlichen genauso ab wie im bisherigen Verfahren.

Je nach gewahltem Ansatz fur die Flussberechnung sind die Kosten fur dieBerechnung des Parameters λ unterschiedlich. Fur jedes Dreieckselement desGitters mussen sowohl beim lokalen als auch beim globalen Lax-Friedrichs-Verfahren drei Parameter λ berechnet werden. In Tabelle (5.1) sind die Kos-ten fur die verschiedenen Verfahren aufgelistet. Die Anzahl der Operationenfur die Berechnung der λ bei Gittertiefe 15 ergeben sich abhangig von derGittertiefe k nach folgender Formel:

#Operationen(k) = #Kanten(k)(2 ∗Op(λ) + 1). (6.1)


Verfahren OperationenOp(λ) zur

λ-Berechnung

Operationen beiGittertiefe 13

konst. LLF 9 1872640konst. GLF 9 1872640

lin. LLF-MaxNorm 32 6406400lin. LLF-Dreieck 34 6800640

Tabelle 6.1: Zusatzlicher Rechenaufwand bei Verwendung des Lax-Friedrichs-Flusses.

Dabei ergeben sich die Anzahl der Kanten abhangig von k nach der Glei-chung

#Kanten(k) =12(3 ∗#Dreiecke(k) + #Randkanten(k)

)+ #Randkanten(k),

wobei die Anzahl der Dreiecke und die Anzahl der Randkanten in Abhangig-keit von k folgendermaßen berechnet werden:

#Dreiecke(k) = 2k+1, #Randkanten(k) = 22+bn2 c.

6.2 Testrechnungen anhand des radialsymmetri-

schen Dammbruch-Problems

Fur die Testrechnungen wurde ein neues Startszenario implementiert, das ra-dialsymmetrische Dammbruch-Problem. Dazu wird die Wassertiefe auf demBerechnungsgebiet Ω = [0, 1] × [0, 1] auf einen konstanten Wert gesetzt.Um den Mittelpunkt des Einheitsquadrats wird ein Kreis Radius r = 1

8

mit erhohtem Wasserspiegel festgelegt, Zu sehen beispielsweise in Abbil-dung (6.1). Die Komponenten ξvx und ξvy sind zum Startzeitpunkt auf nullgesetzt.

Gewahlt wurde dieses Szenario hauptsachlich wegen der ihm eigenen Radi-alsymmetrie. So pflanzt sich die Welle idealerweise in beide Raumrichtungenuber die Zeit in konzentrischen Kreisen fort. Die Restriktion der Losung aufeine Ebene, die senkrecht auf der xy-Ebene steht und durch den Mittel-punkt des Berechnungsgebiets verlauft, kann als Losung eines eindimensio-


nalen Problems betrachtet werden. Diese ist sehr anschaulich und daher gutgeeignet, verschiedene Verfahren zu vergleichen.

6.2.1 Konstante Ansatzfunktionen

In Abbildung (6.1) ist der Losungsverlauf der Flachwasser-Gleichungen zusehen. Als Flussterm wurde hier der lokale Lax-Friedrichs-Fluss verwen-det. Auf der z-Achse ist im ersten Bild die Wasserhohe ξ zum Zeitpunktt = 0.001s abgebildet (ξ = 5 auf dem Turm, ξ = 3 sonst). Die weiteren Bilderfolgen im Abstand von ∆t = 0.015s. Gerechnet wurde dabei mit konstantenAnsatzfunktionen auf Gittertiefe 15 (65536 Elemente) ohne Adaption. AlsZeitschrittweite wurde τ = 0.00001 verwendet. Fur die Randwerte wurde dieWandbedingung benutzt, so dass die sich ausbreitende Welle ab Bild 6 vomRand des Berechnungsgebiets reflektiert wird. Zur Losung der gewohnlichenDifferentialgleichung in der Zeit wurde dabei das Euler-Verfahren verwendet(siehe Anhang).

Zum Vergleich der verschiedenen Flussfunktionen fur den konstanten Ansatzsind in Abbildung (6.2) eindimensionale Schnitte durch die Ebene mit y =0, 5 durch den Losungsverlauf zu den Zeitpunkten t = 0.005s, 0.01s0.015s

und 0.02s abgebildet. Hier wurde die Auslaufrandbedingung verwendet, an-sonsten sind die Parameter des Verfahrens die gleichen wie in Abbildung(6.1). Dabei ist auffallig, dass der globale Lax-Friedrichs-Fluss (blaue Li-nie) fur ein sehr viel starkeres ”zerfließen“ der Losung sorgt als der lokaleLax-Friedrichs-Fluss. Dagegen wird von der schwarzen Kurve, berechnet mitdem modifizierten Lax-Friedrichs-Fluss mit α = 1, sehr viel starker die sichausbreitende Schockwelle mit steileren Flanken der Wellenfront erhalten.

Der lokale Lax-Friedrichs-Fluss (grun) und der modifizierte Lax-Friedrichs-Fluss mit α = 5 (rot) nehmen eine Mittelstellung ein. Sie reproduzierendie Schockwelle nicht so genau wie die schwarze Linie, zerfließen aber auchnicht so schnell wie die blaue Linie. Allerdings ist im dritten Bild in Ab-bildung (6.2) zu erkennen, dass die grune Linie die Absenkung in der Bild-mitte darstellen kann, die von der roten Linie ignoriert wird. Der lokaleLax-Friedrichs-Fluss reagiert also sehr viel genauer auf lokale Schwankun-gen, was durch die Lokalitat des berechneten λ auch zu erwarten ist.


Abbildung 6.1: Lokaler Lax-Friedrichs-Fluss (Wassertiefe 3, Turmhohe 5)mit Wandrandbedingung. Gerechnet auf Tiefe 15 (65536 Elemente) mit kon-stanten Ansatzfunktionen.


Abbildung 6.2: Schnitt durch den Losungsverlauf bei konstanten Ansatz-funktionen. Der Parameter λ wurde nach folgenden Vorschriften berech-net: Modifizierter Lax-Friedrichs-Fluss mit α = 1.0 (schwarz), modifizierterLax-Friedrichs-Fluss mit α = 5.0 (rot), lokaler (grun) und globaler Lax-Friedrichs-Fluss (blau).


Abbildung 6.3: Nicht konvergentes Verfahren bei lokalem Lax-Friedrich-Fluss mit linearen Ansatzfunktionen. λ wird durch quadratische Interpolati-on abgeschatzt. Im unteren Bild ist eine mogliche schlechte Approximationder Funktion λ(x) durch eine Parabel dargestellt.

6.2.2 Lineare Ansatzfunktionen

Fur die Simulation, die Abbildung (6.5) zugrunde liegt, wurde das Dis-continuous-Galerkin-Verfahren mit linearen Ansatzfunktionen und Auslau-frandbedingung verwendet. Als Flussterm wurde der lokale Lax-Friedrichs-Fluss verwendet, der Lax-Friedrichs-Parameter λ wurde mit der Maximums-norm abgeschatzt. Die Wasserhohen betragen am Startzustand ξ = 51 furden Turm in der Mitte und ξ = 50 fur die restliche Wasserflache. Die Ge-schwindigkeiten sind zu Beginn der Berechnungen wieder auf null gesetzt.Als Zeitdiskretisierung wurde auch hier wieder τ = 0.00001s verwendet unddie Zeitintegration mit dem Euler-Verfahren durchgefuhrt. Auch hier wur-de mit Gittertiefe 15 (65536 Elemente) ohne Adaption gerechnet. Das ersteBild wurde zum Zeitpunkt t = 0.0005s gemacht, die weiteren Bilder folgenim Abstand von 0.004s.

Auch das Verfahren mit dem lokalen Lax-Friedrichs-Fluss und λ-Abschatzung


nach der Dreiecksungleichung zeigt einen sehr ahnlichen Verlauf. Wird λ

dagegen mit Hilfe der quadratischen Interpolation approximiert, so konver-giert das Verfahren nicht. Dies lasst sich nur dadurch erklaren, dass dieAbschatzung von

λ = max√

v2x + v2

y +√

gξ,

wie im unteren Bil von Abbildung (6.3) angedeutet, durch ein quadratischesPolynom qualitativ zu schlecht ist und deshalb zur Divergenz des Verfahrensfuhrt. Ein ahnliches, nicht konvergentes Verhalten ist festzustellen, falls derglobale Lax-Friedrichs-Fluss verwendet wird.

In Abbildung (6.5) sind zum Vergleich der verschiedenen konvergenten Ver-fahren wiederum eindimensionale Schnitte entlang der Ebene y = 0.5 ab-gebildet. Zu Beobachten ist, dass sich die Verlaufe fur den lokalen Lax-Friedrichs-Fluss mit λ-Abschatzung nach der Maximumsnorm (grun) undnach der Dreiecksungleichung (blau) auf der einen Seite sowie nach demmodifizierte Lax-Friedrichs-Fluss mit α = 50 (rot) auf der anderen Seitesehr ahnlich verhalten, und zwar in sehr viel starkerem Maße als noch beiden konstanten Ansatzfunktionen in Abbildung (6.2). Diese Beobachtung,dass fur steigenden Polynomgrad bei den Ansatzfunktionen die Wahl desnumerischen Flusses an Einfluß auf die Qualitat der Losung verliert, decktsich mit der Aussage, die Cockburn in [4] bereits veroffentlicht hat.

6.2.3 Der Parameter λ

Was nun noch bezuglich des Lax-Friedrichs-Flusses zu klaren bleibt: Wiehangt die Stabilitat der Losung von dem Parameter λ ab?Um den Verlauf der λ-Werte zu beobachten, wurde fur Abbildung (6.6)der im aktuellen Zeitschritt maximale λ-Wert des jeweiligen Verfahrens furkonstante ind lineare Ansatzfunktionen im Verlauf der Berechnungszeit dar-gestellt. Es ist zu beobachten, das der globale Lax-Friedrichs-Fluss mit sehrviel großeren Werten arbeitet als der lokale Lax-Friedrichs-Fluss. Da die Ei-genwerte der Jacobi-Matrix des Flussterms F, deren Maximum λ ja ist, dieAusbreitungsgeschwindigkeit von Schockwellen beschreiben, ist damit auchdas starkere Zerfließen der durch den globalen Lax-Friedrichs-Fluss berech-neten Losung in Abbildung (6.2) erklart.


Abbildung 6.4: Lokaler Lax-Friedrichs-Fluss (Wassertiefe 50, Turmhohe 51)mit Auslaufrandbedingung. λ wird nach der Maximumsnorm abgeschatzt.Gerechnet auf Tiefe 15 (65536 Elemente) mit linearen Ansatzfunktionen.


Abbildung 6.5: Schnitt durch den Losungsverlauf bei linearen Ansatzfunk-tionen. Der Parameter λ wurde nach folgenden Vorschriften berechnet:Modifizierter Lax-Friedrichs-Fluss mit α = 10.0 (schwarz), modifizierterLax-Friedrichs-Fluss mit α = 50.0 (rot), lokaler Lax-Friedrichs-Fluss λabgeschatzt mit der Maximumsnorm (grun) und der Dreiecksungleichung(blau).


Abbildung 6.6: Verlauf der maximalen λ-Werte pro Zeitschritt. Im linkenBild die Verlaufe der Eigenwerte bei konstanten Ansatzfunktionen, im rech-ten Bild bei linearen Ansatzfunktionen. Es wurden die selben Farben ver-wendet wie in Bild (6.2) beziehungsweise (6.5).

Des weiteren ist zu erkennen, dass die rote Linie in den beiden Bildern,die den modifizierten Lax-Friedrichs-Fluss mit α = 5 bei konstanten, be-ziehungsweise α = 50 bei linearen Ansatzfunktionen, darstellt, etwa einenmittleren Wert der Verlaufe bei den lokalen Lax-Friedrichs-Flussen darstellt.Dies erklart auch das ahnliche Verhalten der roten und grunen beziehungs-weise der roten, grunen und blauen Linie in Bild (6.2) bei konstanten undin Bild (6.5) bei linearen Ansatzfunktionen. Allerdings kann, wie bereitsgesagt, der lokale Lax-Friedrichs-Fluss lokale Phanomene besser darstellen,wie im dritten Bild von Abbildung (6.2) zu sehen ist. Die schwarzen Linien,welche den modifizierten Lax-Friedrichs-Fluss mit α = 1 beziehungsweiseα = 10 darstellen, erzeugen dagegen qualitativ schlechte Losungen.

Wie sich die Wahl des Parameters λ im modifizierten Lax-Friedrichs-Flussauswirkt, ist in Tabelle (6.2) und Tabelle (6.3) sehen. Speziell bei linearenAnsatzfunktionen ist es dabei interessant, wie tief beziehungsweise seicht dasWasserniveau sein darf, um noch ein konvergents Verfahren zu erhalten. Beilinearen Ansatzfunktionen besteht namlich das Problem, dass bei starkenAnderungen der Zustandsvariablen q Anomalien zu beobachten sind. Dieseentstehen dadurch, dass die ”linearen Platten“, aus denen ein Dreiecksele-ment besteht, bis an den Wassergrund oder sogar noch darunter reichen, wasdie Konvergenz des Verfahrens beeintrachtigt. In Tabelle (6.2) und Tabelle(6.3) sind einige dieser Falle getestet.


λ = 1 λ = 5 λ = 20 λ = 100H = 1,5/3,5 konvergent konvergent konvergent konvergent

H = 3/5 konvergent konvergent konvergent konvergentH = 23/25 nicht

konvergentkonvergent konvergent konvergent

H = 53/55 nichtkonvergent

nichtkonvergent

konvergent konvergent

Tabelle 6.2: Einfluss des Parameters λ und der Wassertiefe auf die Konver-genz bei konstanten Ansatzfunktionen.

λ = 5 λ = 50 λ = 100 λ = 200H = 10/11 nicht

konvergentnicht

konvergentnicht

konvergentnicht

konvergentH = 20/21 nicht

konvergentkonvergent nicht

konvergentnicht

konvergentH = 50/51 nicht

konvergentkonvergent konvergent nicht

konvergent

Tabelle 6.3: Einfluss des Parameters λ und der Wassertiefe auf die Konver-genz bei linearen Ansatzfunktionen.

Um das Verhalten des Lax-Friedrichs-Parameters λ zu studieren ist in Ab-bildung (6.7) auf den Bildern in der linken Spalte der Verlauf der Was-seroberflache der Flachwassergleichung mit Auslaufrandbedingung zu denZeitpunkten t = 0, 0005s und weiterhin in 0, 004s Schritten abgebildet. ZumVergleich ist der Parameter-Wert λ zu den gleichen Zeitpunkten in der rech-ten Spalte visualisiert. Es ist zu erkennen, dass der λ-Wert genau an den-jenigen Orten groß ist, an denen zum gleichen Zeitpunkt eine Wellenfront,also ein steiler Anstieg in der Wasserhohe erkennbar ist. So ist im drittenλ-Bild sehr gut zu erkennen, dass sich hier der Wasserstand wieder starkverandern wird, wahrend im dritten Bild der Wasseroberflache dies nochnicht zu erkennen ist. Die Front der λ-Werte scheint also ein kleines Stuckvor der Wellenfront zu verlaufen.

Dieser Eindruck wird in Abbildung (6.8) bestatigt. Die Bilder sind eben-so wie in Abbildung (6.2) und (6.5) eindimensionale Schnitte durch das


Abbildung 6.7: Lokaler Lax-Friedrichs-Fluss (links) bei linearen Ansatzfunk-tionen, λ abgeschatzt mit der Dreiecksungleichung. Auf der rechten Seite istder Lax-Friedrichs-Parameter λ zum selben Zeitpunkt dargestellt.


Abbildung 6.8: Verlauf der Losung (blau) und des Lax-Friedrichs-Parameters(grun) λ. Berechnet mit linearen Ansatzfunktionen, λ wurde nach der Drei-ecksungleichung abgeschatzt.


Berechnungsgebiet. Hier ist ganz deutlich zu erkennen, dass die Front derLax-Friedrichs-Parameter λ etwas vor der Wellenfront verlauft, sozusagenein ”Fruhanzeiger“ fur starkere Veranderungen des Losungsverlaufs ist. Dieslegt nahe, diesen Parameter λ oder die Anderung des Parameters ∂λ

∂t als Ad-aptionskriterium zu wahlen. Als bisheriges Adaptionskriterium wurde dieAnderung der Wasserhohe verwendet. Die Adaption hinkte damit der Wel-lenfront immer etwas hinterher (siehe Kapitel (3.4)). Durch die Verwendungvon λ konnte damit durchaus eine Verbesserung des adaptiven Verfahrenserzielt werden.

Kapitel 7

Abschließende Bemerkung

Wie zu Beginn der Arbeit bereits angesprochen, ist es von großter Wichtig-keit, Verfahren mit moglichst vor allem geringem Speicherbedarf aber auchgeringer Rechenzeit zu konstruieren. Werden solche Verfahren dann quali-tativ weiterentwickelt, so ist besonders darauf zu achten, dass diese beidenGroßen durch die Verbesserungen nicht gefahrdet werden. Ich hoffe, durchmeine Arbeit einen Teil dazu beigetragen zu haben, das Verfahren robusterund genauer und damit bei kaum erhohter Rechenzeit und praktisch glei-chem Speicherbedarf verbessert zu haben.

Als ersten Schritt zur weiteren Verbesserung des Verfahrens hat sich, wiein Kapitel 6.2.3 angedeutet, durch die Beobachtung des Lax-Friedrichs-Parameters λ ergeben, diesen als Adaptionskriterium zu wahlen. Dadurchsollte eine genauere Losung zu realisieren sein, weil das Gitter schon verfei-nert sein kann, wenn die Wellenfront diesen Ort passiert.

Als weitere Verbesserungsmoglichkeit ware eine Zeitschrittweitensteuerungzu nennen, die durch eine Schatzung des maximalen lokalen Fehlers denZeitschritt vergroßert oder verringert. So konnte auch entlang der ZeitachseAdaptivitat erzielt werden, und damit die Berechnungszeit verringert wer-den.

Abschließend ware es wunschenswert, eine vollstandigere Form der Flach-wassergleichungen (2.1) - (2.2) losen zu konnen. Dazu mussen in die Berech-nungsvoschriften weitere Parameter, wie zum Beispiel veranderlicher Mee-resgrund inklusive starker Steigungen am Meeresgrund (die fur die Entste-

80

KAPITEL 7. ABSCHLIESSENDE BEMERKUNG 81

hung großer Brandungswellen verantwortlich sind), eingearbeitet werden.

Anhang A

Verwendete mathematische

Verfahren

A.1 Explizites Euler-Verfahren

Das explizite Eulerverfahren ist ein Standardverfahren zur Integration ge-wohnlicher Differentialgleichungen

dy(t)dt

= f(t, y(t)), y(t0) = y0. (A.1)

Fur diskrete Zeitpunkte tk berechnet es eine Punktfolge (yk)k=1,...,n, welchedie Funktionswerte y(tk) approximiert. Es ist gegeben durch die Vorschrift:

yk+1 = yk + (tk−1 − tk)f(tk, yk), y0 = y(0). (A.2)

Das Euler-Verfahren hat Konvergenzorndung 1, bei Halbierung der Schritt-weite verringert sich der Fehler yk− y(tk) folglich auch um den Faktor 1

2 . InAbbildung (A.1) ist das Euler-Verfahren grafisch dargestellt.

A.2 Runge-Kutta-Verfahren zweiter Ordnung

Das Runge-Kutta-Verfahren zweiter Ordnung ist ein verbessertes Verfah-ren zur Losung von Gleichung (A.1). Im Gegensatz zum expliziten Euler-Verfahren hat es Konvergenzordnung zwei. Wird die Schrittweite um einenFaktor α verringert, so verringert sich der Fehler yk − y(tk) um den Faktorα2.

82

ANHANG A. VERWENDETE MATHEMATISCHE VERFAHREN 83

y

t

yk+1

yk

tk tk+1

Abbildung A.1: Berechnung eines Zeitschrittes beim expliziten Euler-Verfahren.

yk+0,5

yk

yk+1

tk tk+0,5 tk+1

t

y

Abbildung A.2: Berechnung eines Zeitschrittes beim Runge-Kutta-Verfahrenzweiter Ordnung.Der Fehler ist aufgrund des Zwischenschrittes deutlich klei-ner als beim expliziten Euler-Verfahren.

ANHANG A. VERWENDETE MATHEMATISCHE VERFAHREN 84

Es gehorcht der Vorschrift

yk+1 = yk + (τ)f(tk +

τ

2, yk +

τ

2f(tk, yk)

)(A.3)

fur k = 1, . . . , n. Dabei ist τ = tk+1 − tk. Eine Veranschaulichung der Vor-gehensweise ist in Abbildung (A.2) zu sehen.

Anhang B

matlab-Funktionen

B.1 matlab-Funktion fur Daten in Matrixform

function visualMatrix(A)

number = A(1,2);low = 2; high = number + 1;sizeA = size (A);

aviobj = avifile(’filename.avi’);

while(high ¡ sizeA)x = A(low:high,1); y = A(low:high,2); z = A(low:high,3);

tri = delaunay(x,y,’QJ’,’Pp’);%fur Gittervisualisierungtri = reshape([1:number],[3:number/3]); tri = tri’;

trisurf(tri,x,y,z);%fur Gittervisualisierungtrimesh(tri,x,y,z);

view(-37,60);colormap(winter)lighting phong

85

ANHANG B. MATLAB-FUNKTIONEN 86

shading interp

%konstante Ansatzfunktionen, Transportgleichung, Gittervisualisierungaxis([0,1,0,1,1,5]); caxis([2.5,3.5]);% lineare Ansatzfunktionen (nicht Gittervisualisierung)axis([0,1,0,1,49,51]); caxis([49.5,50.5]);axis offframe = getframe;aviobj = addframe(aviobj,frame);

number = A(high + 1,2);low = high + 2;high = low + number - 1;

end

aviobj = close(aviobj);

end


B.2 matlab-Funktion fur Daten in Vektorform

function visualVektor(A)

number = A(2);low = 4;high = 3*number + 3;sizeA = size(A);B = A(low:high);x = B(1:number);y = B(number+1:2*number);z = B(2*number+1:3*number);tri = delaunay(x,y,’QJ’,’Pp’);

aviobj = avifile(’filename.avi’);

while (high ¡ sizeA(1))trisurf(tri,x,y,z);

view(-37,60);colormap(winter)lighting phongshading interp

%konstante Ansatzfunktionen, Transportgleichung, Gittervisualisierungaxis([0,1,0,1,1,5]); caxis([2.5,3.5]);% lineare Ansatzfunktionen (nicht Gittervisualisierung)axis([0,1,0,1,49,51]); caxis([49.5,50.5]);axis off

frame = A(high + 2);aviobj = addframe(aviobj,frame);

number = A(high + 2);low = high + 4;high = low + number - 1;


z = A(low:high);endaviobj = close(aviobj);end


B.3 matlab-Funktion fur Bild-Plots

function plotSWE(A,frame)

number = A(2);low = 4;high = 3*number+3;

B=A(low:high);x=B(1:number);y=B(number+1:2*number);z=B(2*number+1:3*number);

%fur Gittervisualisierungtri = delaunay(x,y);

% Hole das angegebene framelow = (3*number+4) + (frame-1)*(number+3) + 3;high = low + number - 1;z=A(low:high);

%Baue die Oberflachetrisurf(tri,x,y,z);

%fur konstante Ansatzfunktionenaxis([0,1,0,1,0,5]); caxis([1.5,4.5]);%fur lineare Ansatzfunktionenaxis([0,1,0,1,49,51]); caxis([49.5,50.5]);axis off;

view(-37,60);lighting phongshading interpcolormap(winter)

%Ausgabe von Bildern


print -djpeg100 name.jpg%print -deps name.eps

end


B.4 matlab-Funktion fur zweidimensionalen Schnitt

function schnitt2d(A,frame)

number = A(2);low = 4;high = 3*number+3;

B=A(low:high);x=B(1:number);y=B(number+1:2*number);

% Hole das angegebene framelow = (3*number+4) + (frame-1)*(number+3) + 3;high = low + number - 1;z=A(low:high);

xd = 0:0.01:1;[XI,YI] = meshgrid(xd,0.5);

ZI = griddata(x,y,z,XI,YI);

%verschiedene Farben fur die Verfahrenplot(XI,ZI,’LineWidth’,2,’Color’,[0,0,0]); %Modlf 1 / 10%plot(XI,ZI,’LineWidth’,2,’Color’,[1,0,0]); %Modlf 5 / 50%plot(XI,ZI,’LineWidth’,2,’Color’,[0,1,0]); %Llf / Maxnorm%plot(XI,ZI,’LineWidth’,2,’Color’,[0,0,1]); % Glf /Dreieck

axis([0, 1, 49.0, 51.5]); % fur lineare Ansatzfunktionen%axis([0, 1, 2.0, 6.0]); % fur konstante Ansatzfunktionen

hold on;

%Ausgabe von Bildern:print -djpeg100 test.png


end

Literaturverzeichnis

[1] BADER, M.: Raumfullende Kurven. Vorlesungsskript, TUMunchen, 2005.

[2] BADER, M., SCHRAUFSTETTER, S., VIGH, C. und BEHRENS,J.: Efficient Storage and Processing of Adaptive Triangular Gridsusing Sierpinski Curves. 2005.

[3] BOCK, C.: Discontinuous-Galerkin-Verfahren zur Losung derFlachwasser-Gleichungen auf adaptiven Dreiecksgittern. Diplomar-beit, TU Munchen, April 2008.

[4] COCKBURN, B.: Discontinuous-Galerkin Methods for convection-dominated Problems. Journal of Scientific Computing, 16(3): S. 173-261, 2002.

[5] DAWSON, C. und AIZINGER, V.: A Discontinuous Galerkin Me-thod for Three-Dimensional Shallow Water Equations. Journal ofScientific Computing, 22(1):S. 245-267, 2005.

[6] HESTHAVEN, J., WARBURTON, T.: Nodal Discontinuous Galer-kin Methods. ISBN: 978-0-387-72065-4, Springer Verlag, 2008.

[7] LEVEQUE, R.: Finite Volume Methods for Hyperbolic Problems.Cambridge University Press, Cambridge, 2002.

[8] RADZIEOWSKI, C.: Numerische Simulation zeitabhangiger Pro-bleme auf dynamisch-adaptiven Dreiecksgittern. Diplomarbeit, TUMunchen, November 2007.

[9] REMACLE, J.-F., FRANZAO, S., LI, X. und SHEPARD, M.: Ad-aptive Discontinuous Galerkin Method for the Shallow Water Equa-tions. International Journal for Numerical Methods in Fluids, 2003.

93

LITERATURVERZEICHNIS 94

[10] ROE, P. L.: Approximate Riemann Solvers, Parameter Vectors, andDifference Schemes. Journal of Computational Physics, 43, S.357-372, 1981.

[11] ROE, P. L.: in Proceedings, Seventh International Conferenceof Numerical Methods in Fluid Dynamics. Springer Verlag, NewYork/berlin, 1981.

[12] SCHRAUFSTETTER, S.: Speichereffiziente Algorithmen zumLosen partieller Differentialgleichungen auf adaptiven Dreiecksgit-tern. Diplomarbeit, TU Munchen, Juli 2006.

[13] SCHWAIGER, J.: Adaptive Discontinuous-Galerkin-Verfahren zumLosen der Flachwassergleichungen mit verschiedenen Randbedin-gungen. Diplomarbeit, TU Munchen, September 2008.

[14] SCHWANENBERG, D.: Die Runge-Kutta-Discontinuous-Galerkin-Methode zur Losung konvektionsdominierter tiefen-gemittelter Flachwasserprobleme. Dissertation, RWTH Aachen,Marz 2003.

Adaptive Discontinuous-Galerkin-Verfahren zur L¨osung der ... · 1 Hiermit erkl¨are ich, dass ich...

Documents

Transcript of Adaptive Discontinuous-Galerkin-Verfahren zur L¨osung der ... · 1 Hiermit erkl¨are ich, dass ich...