1. Riemann Solver

1.1. Einleitung 

Dieses Kapitel dokumentiert die Implementierung und das Testen des f-wave Solvers für die eindimensionalen Shallow Water Equations (SWE). Der f-wave Solver ist der zentrale Baustein unserer Tsunami-Simulation: Er berechnet, wie sich eine Unstetigkeit zwischen zwei benachbarten Zellen in Wellen auflöst, die sich nach links und rechts ausbreiten.

Die zugehörige Aufgabenstellung findet sich unter: https://scalable.uni-jena.de/opt/tsunami/chapters/assignment_1.html

1.2. Mathematische Grundlagen 

1.2.1. Shallow Water Equations 

Unser Ausgangspunkt sind die eindimensionalen Shallow Water Equations — ein System aus zwei Erhaltungsgleichungen, das die Dynamik von Flachwasser beschreibt:

\[\begin{split}\underbrace{ \begin{bmatrix} h \\ hu \end{bmatrix}_{t} }_{\text{zeitliche Änderung}} + \underbrace{ \begin{bmatrix} hu \\ hu^2 + \frac{1}{2}g h^2 \end{bmatrix}_{x} }_{\text{räumlicher Fluss}} = S(x,t)\end{split}\]

Die erste Zeile beschreibt die Massenerhaltung: die Wasserhöhe \(h\) ändert sich, wenn Wasser mit dem Impuls \(hu\) zu- oder abfliesst. Die zweite Zeile beschreibt die Impulserhaltung: der Impuls \(hu\) ändert sich durch den konvektiven Transport \(hu^2/h\) und den hydrostatischen Druck \(\frac{1}{2}g h^2\).

Wir definieren dazu zwei zentrale Grössen:

Erhaltungsgrössen: \(q = [h,\; hu]^T\) — Wasserhöhe und Impuls
Flussfunktion: \(f(q) = [hu,\; hu^2/h + \frac{1}{2}g h^2]^T\)

Für Assignment 1 betrachten wir den homogenen Fall \(S(x,t) = 0\) (keine Bathymetrie, keine Reibung).

1.2.2. Das Riemann-Problem 

An jeder Zellkante treffen zwei Zustände aufeinander — links \(q_l\) und rechts \(q_r\). Diese Unstetigkeit nennt man ein Riemann-Problem:

\[\begin{split}q(x, 0) = \begin{cases} q_l & \text{wenn } x < 0 \\ q_r & \text{wenn } x > 0 \end{cases}\end{split}\]

Die Lösung besteht aus zwei Wellen, die sich vom Punkt der Unstetigkeit nach links und rechts ausbreiten. Der f-wave Solver approximiert beide Wellen als Schockwellen.

1.3. F-wave Solver — Theorie 

Der f-wave Solver löst das Riemann-Problem in drei Schritten:

Wellengeschwindigkeiten bestimmen (Eigenwerte)
Flusssprung in Eigenvektoren zerlegen (Eigenkoeffizienten)
Wellen den Zellen zuordnen (Net-Updates)

1.3.1. Schritt 1: Eigenwerte (Wellengeschwindigkeiten)

Die Geschwindigkeiten der beiden Wellen werden über Roe-gemittelte Grössen approximiert. Aus den linken und rechten Zuständen berechnen wir zunächst:

\[ \begin{align}\begin{aligned}h^\text{Roe} &= \frac{1}{2}(h_l + h_r)\\u^\text{Roe} &= \frac{u_l \sqrt{h_l} + u_r \sqrt{h_r}}{\sqrt{h_l} + \sqrt{h_r}}\end{aligned}\end{align} \]

Daraus ergeben sich die Eigenwerte (Wellengeschwindigkeiten):

\[ \begin{align}\begin{aligned}\lambda_1 &= u^\text{Roe} - \sqrt{g \cdot h^\text{Roe}}\\\lambda_2 &= u^\text{Roe} + \sqrt{g \cdot h^\text{Roe}}\end{aligned}\end{align} \]

Die physikalische Bedeutung: \(\lambda_1\) ist die Geschwindigkeit der nach links laufenden Welle, \(\lambda_2\) die der nach rechts laufenden. Der Term \(\sqrt{g \cdot h^\text{Roe}}\) ist die lokale Wellenausbreitungsgeschwindigkeit.

Im Code

Die Eigenwerte werden in der Funktion eigenvalues() berechnet. Die Wurzel \(\sqrt{g}\) ist als Konstante gSqrt in constants.h vordefiniert, sodass nur noch \(\sqrt{h^\text{Roe}}\) berechnet werden muss.

1.3.2. Schritt 2: Eigenkoeffizienten (Flusssprung-Zerlegung)

Dies ist der Kernschritt und gleichzeitig der wesentliche Unterschied zum Roe-Solver.

Wir berechnen zunächst die Flussfunktion \(f(q)\) für beide Seiten:

\[\begin{split}f(q_l) = \begin{bmatrix} (hu)_l \\ \frac{(hu)_l^2}{h_l} + \frac{1}{2} g\, h_l^2 \end{bmatrix}, \qquad f(q_r) = \begin{bmatrix} (hu)_r \\ \frac{(hu)_r^2}{h_r} + \frac{1}{2} g\, h_r^2 \end{bmatrix}\end{split}\]

und daraus den Flusssprung:

\[\Delta f = f(q_r) - f(q_l)\]

Important

Unterschied zum Roe-Solver: Der Roe-Solver zerlegt den Sprung in den Zustandsgrössen \(\Delta q = q_r - q_l = [h_r - h_l,\; (hu)_r - (hu)_l]^T\). Der f-wave Solver zerlegt stattdessen den Sprung in den Flüssen \(\Delta f\).

Diesen Flusssprung zerlegen wir nun in die Eigenvektoren \(r_1, r_2\) der Roe-Matrix:

\[\Delta f = \alpha_1 \cdot r_1 + \alpha_2 \cdot r_2\]

mit den Eigenvektoren:

\[\begin{split}r_1 = \begin{bmatrix} 1 \\ \lambda_1 \end{bmatrix}, \qquad r_2 = \begin{bmatrix} 1 \\ \lambda_2 \end{bmatrix}\end{split}\]

In Matrixform geschrieben ist das ein lineares Gleichungssystem:

\[\begin{split}\underbrace{ \begin{bmatrix} 1 & 1 \\ \lambda_1 & \lambda_2 \end{bmatrix} }_{R} \begin{bmatrix} \alpha_1 \\ \alpha_2 \end{bmatrix} = \Delta f\end{split}\]

Die Lösung erhalten wir durch Invertieren der 2×2-Matrix \(R\):

\[\begin{split}R^{-1} = \frac{1}{\lambda_2 - \lambda_1} \begin{bmatrix} \lambda_2 & -1 \\ -\lambda_1 & 1 \end{bmatrix}\end{split}\]

Damit ergibt sich:

\[ \begin{align}\begin{aligned}\alpha_1 &= \frac{1}{\lambda_2 - \lambda_1} \left( \lambda_2 \cdot \Delta f_0 - \Delta f_1 \right)\\\alpha_2 &= \frac{1}{\lambda_2 - \lambda_1} \left( -\lambda_1 \cdot \Delta f_0 + \Delta f_1 \right)\end{aligned}\end{align} \]

Im Code

In der Funktion eigencoefficients() werden die drei Teilschritte ausgeführt:

Flussfunktion auswerten
Flusssprung bilden
R⁻¹ · Δf berechnen

1.3.3. Schritt 3: Net-Updates 

Aus den Eigenkoeffizienten \(\alpha_p\) und Eigenvektoren \(r_p\) bilden wir die f-Wellen:

\[\begin{split}Z_p = \alpha_p \cdot r_p = \alpha_p \cdot \begin{bmatrix} 1 \\ \lambda_p \end{bmatrix} = \begin{bmatrix} \alpha_p \\ \alpha_p \cdot \lambda_p \end{bmatrix}\end{split}\]

Jede f-Welle \(Z_p\) wird der linken oder rechten Zelle zugeordnet — abhängig vom Vorzeichen des zugehörigen Eigenwertes \(\lambda_p\):

\[ \begin{align}\begin{aligned}A^- \Delta Q &= \sum_{p:\; \lambda_p < 0} Z_p \qquad \text{(Net-Update für die linke Zelle)}\\A^+ \Delta Q &= \sum_{p:\; \lambda_p > 0} Z_p \qquad \text{(Net-Update für die rechte Zelle)}\end{aligned}\end{align} \]

Anschaulich: Wellen, die sich nach links ausbreiten (\(\lambda_p < 0\)), verändern die linke Zelle. Wellen, die sich nach rechts ausbreiten (\(\lambda_p > 0\)), verändern die rechte Zelle.

Im Code

netUpdates() berechnet mit den anderen Funktionen die Updates:

eigenvalues() berechnen
eigencoefficients() berechnen
f-Wellen berechnen
basierend auf Vorzeichen auf Update für rechte o. linke Zelle addieren

1.3.4. Sonderfälle 

Steady State (\(q_l = q_r\)): Wenn beide Seiten identisch sind, ist \(\Delta f = 0\), also \(\alpha_1 = \alpha_2 = 0\) und alle Net-Updates verschwinden. Die Simulation bleibt stabil — es passiert nichts.

Supersonische Strömung (\(\lambda_1, \lambda_2\) gleiches Vorzeichen): Wenn die Strömungsgeschwindigkeit so gross ist, dass beide Wellen in dieselbe Richtung laufen, erhält eine Seite das gesamte Update, die andere bleibt unverändert. Beispiel: Bei \(\lambda_1 > 0\) und \(\lambda_2 > 0\) laufen beide Wellen nach rechts, d.h. \(A^- \Delta Q = 0\) (kein Update links).

1.4. Unterschied zum Roe-Solver 

Beide Solver verwenden dieselben Eigenwerte und Eigenvektoren. Der entscheidende Unterschied liegt in der Zerlegung:

	Roe-Solver	F-wave Solver
Zerlegt wird	\(\Delta q = q_r - q_l\)	\(\Delta f = f(q_r) - f(q_l)\)
Gleichungssystem	\(R \cdot \alpha = \Delta q\)	\(R \cdot \alpha = \Delta f\)
Wellen	\(\text{wave}_p = \lambda_p \cdot \alpha_p \cdot r_p\)	\(Z_p = \alpha_p \cdot r_p\)
Skalierung	Zusätzlich mit \(\lambda_p\) skaliert	Keine weitere Skalierung nötig

Beim Roe-Solver wird der Zustandssprung zerlegt und die resultierenden Wellen müssen anschliessend noch mit der Wellengeschwindigkeit multipliziert werden, um die korrekte Flussdifferenz zu erhalten. Beim f-wave Solver ist diese Skalierung bereits implizit enthalten, weil von Anfang an der Flusssprung zerlegt wird.

1.5. Codestruktur 

Die Implementierung besteht aus drei Dateien in src/solvers/:

FWave.h

Header mit der Klasse tsunami_lab::solvers::FWave. Die öffentliche Schnittstelle ist identisch zum Roe-Solver: eine statische Methode netUpdates() mit den Höhen und Impulsen beider Seiten als Eingabe und den Net-Updates als Ausgabe.

FWave.cpp

Implementierung mit drei Funktionen, die den drei mathematischen Schritten entsprechen:

eigenvalues() → Schritt 1 (Roe-Eigenwerte \(\lambda_{1,2}\))
eigencoefficients() → Schritt 2 (Flusssprung-Zerlegung \(\alpha_{1,2}\))
netUpdates() → Schritt 3 (f-Wellen bilden und zuordnen)

FWave.test.cpp

Unit-Tests mit Catch2 (siehe nächster Abschnitt).

1.6. Unit-Tests 

Unsere Tests decken die wesentlichen Szenarien ab:

1.6.1. Eigenwert-Verifikation 

Test mit \(h_l = 10,\; h_r = 9,\; u_l = -3,\; u_r = 3\):

Die Roe-Mittelwerte ergeben sich zu \(h^\text{Roe} = 9.5\) und \(u^\text{Roe} \approx -0.079\). Daraus folgen \(\lambda_1 \approx -9.731\) und \(\lambda_2 \approx 9.573\). Diese Werte werden gegen die analytische Lösung geprüft.

1.6.2. Eigenkoeffizienten-Verifikation 

Für denselben Testfall verifizieren wir, dass die Zerlegung konsistent ist: \(\alpha_1 + \alpha_2 = \Delta f_0\) und \(\alpha_1 \lambda_1 + \alpha_2 \lambda_2 = \Delta f_1\). Das stellt sicher, dass \(R \cdot \alpha = \Delta f\) tatsächlich erfüllt ist.

1.6.3. Steady State 

Bei \(q_l = q_r\) (z.B. \(h = 10,\; hu = 0\) auf beiden Seiten) müssen alle Net-Updates im Rahmen der Maschinengenauigkeit null sein. Wir testen dies sowohl für den trivialen Fall (\(hu = 0\)) als auch für einen nicht-trivialen Gleichgewichtszustand (\(h = 5,\; hu = 3\) auf beiden Seiten).

1.6.4. Dambreak 

Symmetrischer Dammbruch (\(h_l = 10,\; h_r = 8,\; hu = 0\)): Da keine Strömung vorliegt, sind die Eigenwerte betragsmässig gleich (\(\lambda_1 = -\lambda_2\)). Daraus folgt, dass die Höhen-Updates entgegengesetztes Vorzeichen haben und die Impuls-Updates gleich sind. Zusätzlich prüfen wir, dass die Summe der Net-Updates den gesamten Flusssprung ergibt.

1.6.5. Supersonisch 

Wir konstruieren Fälle, in denen beide Eigenwerte dasselbe Vorzeichen haben:

Alle Wellen nach rechts (\(h = 0.1,\; u = 10\) → \(\lambda_{1,2} > 0\)): Das linke Net-Update muss exakt null sein, das rechte muss dem gesamten Flusssprung entsprechen.
Alle Wellen nach links (\(h = 0.1,\; u = -10\) → \(\lambda_{1,2} < 0\)): Spiegelbildlicher Fall — das rechte Net-Update ist null.

1.6.6. Tests Ausführen 

# Alle FWave Tests
./build/tests "[FWave]"

# einzelnen Testcases
./build/tests "[FWaveEigenvalues]"
./build/tests "[FWaveEigencoefficients]"
./build/tests "[FWaveSteadyState]"
./build/tests "[FWaveDamBreak]"
./build/tests "[FWaveSupersonic]"

1.7. Beiträge 

Person	Beitrag
Moritz Arnhold	Build-Setup, Dokumentation, Code Review
Moritz Martin	FWave-Implementierung, FWave-Tests, Dokumentation