2. Finite Volume Discretization

2.1. Einleitung

Dieses Kapitel dokumentiert die Implementierung der Finite-Volumen-Diskretisierung für die eindimensionalen Shallow Water Equations sowie die Einbettung unseres f-wave-Solvers in den Zeit-Integrator. Nachdem Kapitel 01 einen einzelnen Riemann-Zustandssprung aufgelöst hat, bauen wir hier den kompletten 1D-Löser: eine Zellstruktur mit Ghost-Cells, eine Zeitschleife mit CFL-stabilem Zeitschritt und eine Validierung gegen ~1 000 000 vorgerechneter Riemann-Lösungen.

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

2.2. Finite-Volumen-Diskretisierung

2.2.1. Zellstruktur und Erhaltungsform

Das Rechengebiet wird in \(N\) Zellen gleicher Breite \(\Delta x\) zerlegt. Jede Zelle \(i\) trägt einen Zellmittelwert

\[Q_i^{n} \approx \frac{1}{\Delta x} \int_{x_{i-1/2}}^{x_{i+1/2}} q(x, t_n)\,dx\]

der Erhaltungsgrössen \(q = [h,\; hu]^T\). Integriert man die SWE über eine Zelle und ein Zeitintervall \([t_n, t_{n+1}]\), ergibt sich die konservative Update-Formel:

\[Q_i^{n+1} = Q_i^{n} - \frac{\Delta t}{\Delta x}\,\bigl( F_{i+1/2} - F_{i-1/2} \bigr)\]

mit den Flüssen \(F_{i\pm 1/2}\) an den Zellkanten. An jeder Kante steht ein Riemann-Problem zwischen den Zellmittelwerten links und rechts der Kante — genau die Situation, die unser f-wave-Solver löst.

2.2.2. Net-Updates statt Fluss-Differenzen

Der f-wave-Solver liefert direkt die Net-Updates \(A^\pm\Delta Q\) (siehe Riemann Solver), also den Anteil des Flusssprungs, der auf die linke bzw. rechte Zelle wirkt. Damit wird aus der obigen Update-Formel:

\[Q_i^{n+1} = Q_i^{n} - \frac{\Delta t}{\Delta x}\,\bigl( \underbrace{A^{+}\Delta Q_{i-1/2}}_{\text{von Kante links}} + \underbrace{A^{-}\Delta Q_{i+1/2}}_{\text{von Kante rechts}} \bigr)\]

Für jede Kante werden \(A^-\Delta Q\) und \(A^+\Delta Q\) einmal berechnet und beiden Nachbarzellen zugeordnet. Die innere Schleife in src/patches/WavePropagation1d.h iteriert daher über Kanten, nicht über Zellen.

2.2.3. CFL-Bedingung

Damit Wellen in einem Zeitschritt höchstens eine halbe Zelle weit wandern (Courant-Friedrichs-Lewy-Bedingung), wählen wir

\[\Delta t = \tfrac{1}{2}\, \frac{\Delta x}{\max_i\bigl( |u_i| + \sqrt{g\,h_i}\bigr)}\]

In main.cpp wird \(\Delta t\) einmalig aus der Anfangsbedingung bestimmt und für die gesamte Simulation festgehalten. Die Skalierung \(\lambda = \Delta t / \Delta x\) geht direkt in WavePropagation1d::timeStep(lambda).

2.3. Implementierung

2.3.1. WavePropagation1d

Die Klasse hält zwei Arrays pro Erhaltungsgrösse (m_h[2], m_hu[2]): eines für den aktuellen und eines für den neuen Zeitschritt. Statt jeden Schritt zu kopieren, wechseln wir den Lese-/Schreibpuffer über m_step ^= 1 (Double-Buffering). Jeweils am Rand des Gebiets liegt eine Ghost-Zelle, insgesamt also \(N+2\) Einträge pro Array.

2.3.2. Solver-Umschaltung per Template

Der konkrete Riemann-Solver wird als Template-Parameter übergeben:

template<class T>
class WavePropagation1d : public WavePropagation<T> { ... };

// in main.cpp
WavePropagation1d<FWave>  l_waveProp(l_nx);
// alternativ:
// WavePropagation1d<Roe> l_waveProp(l_nx);

Damit teilen sich Roe- und f-wave-Solver exakt denselben Patch-Code, ohne virtuelle Funktionsaufrufe in der inneren Schleife. Ein Wechsel des Solvers ist eine einzige Zeile plus Rebuild.

2.3.3. Outflow-Randbedingung

setGhostOutflow() kopiert den Zustand der jeweiligen Randzelle in die Ghost-Zelle:

l_h[0]         = l_h[1];              // linke Ghost-Zelle
l_h[m_nCells+1] = l_h[m_nCells];      // rechte Ghost-Zelle

An der Ghost-Kante ist damit \(q_l = q_r\), der Flusssprung \(\Delta f = 0\) — die Welle verlässt das Gebiet spurlos. Das ist die offene Randbedingung und vermeidet künstliche Reflexionen.

2.4. Middle-States-Validierung

Für einen Sanity-Check liegt uns resources/middle_states.csv mit ≈ 1 Mio. exakt gelösten Riemann-Problemen \((h_l, h_r, hu_l, hu_r, h^*)\) vor.

2.4.1. Zweistufige Prüfung

Die Validierung läuft auf zwei Ebenen:

Unit-Test (schnell, CI-tauglich). In WavePropagation1d.test.cpp prüft der Catch2-Testcase [WaveProp1dMiddleStates] die ersten 50 Zeilen der CSV. Für jede Zeile wird das Gebiet in linke/rechte Hälfte aufgeteilt, bis zu einer CFL-sicheren Plateauzeit simuliert und die mittlere Zelle mit \(h^*\) verglichen. Es dürfen höchstens 5 % der Zeilen einen relativen Fehler > 1 % haben.

Full-Sweep-Binary (komplett). build/validate_middle_states [nCells] [maxRows] [relTol] geht über beliebig viele Zeilen der CSV und schreibt eine Detail-CSV mit abs_err, rel_err und pass-Flag pro Problem:

./build/validate_middle_states 100 2000 1e-2
# → validation_middle_states.csv, Summary auf stdout

2.4.2. Ergebnisse

Relative Fehler des f-wave-Solvers gegen die exakten Middle-States. Bei 2000 Zeilen (100 Zellen, 0,4 · Halbdomänzeit) ist die Durchfallquote 0 %, der Fehler liegt im Bereich \(10^{-7}\), also tief unter unserer Toleranz von 1 %.

Erkenntnis

Der Fehler liegt nicht bei 1e-7 wegen der Finite-Volumen-Diskretisierung, sondern weil unser t_real float ist — das ist die Rundungsgrenze einfacher Gleitkommazahlen. Mit double würde der Fehler um Grössenordnungen fallen, bei 100 Zellen aber auch der Gitterfehler sichtbar werden. Für die Sanity-Prüfung ist diese Auflösung völlig ausreichend.

2.5. 2.1 Shock-Shock und Rare-Rare

Die shock-shock- und rare-rare-Setups sind zwei symmetrische Riemann-Probleme, in denen zwei Ströme gleicher Höhe mit entgegengesetztem Impuls aufeinandertreffen bzw. auseinanderdriften. Beide sind unter src/setups/ implementiert und werden aus main.cpp über die Szenarien shockshock bzw. rarerare ausgewählt.

2.5.1. Setup

Szenario	shock-shock	rare-rare
\(h_l,\,h_r\)	\(h\) (beidseitig gleich)	\(h\) (beidseitig gleich)
\(hu_l,\,hu_r\)	\(+hu,\;-hu\) (kollidierend)	\(-hu,\;+hu\) (divergierend)
Erwartete Lösung	zwei Schocks nach aussen, hohes Plateau	zwei Verdünnungswellen nach aussen, tiefes Plateau
CLI-Aufruf	./build/tsunami_lab N shockshock h hu loc	./build/tsunami_lab N rarerare h hu loc

2.5.2. Shock-shock: kollidierende Ströme

Zeitliche Entwicklung. Zwei Schocks wandern mit \(\pm\sqrt{g\,h^*}\) nach aussen; das Mittelplateau mit \(u=0\) wächst.

Fünf Zeitscheiben desselben Laufs: oben Wasserhöhe, unten Teilchengeschwindigkeit \(u=hu/h\). Zum Endzeitpunkt sind beide Schocks bereits durch den offenen Rand gewandert; im gesamten Gebiet steht das Plateau mit \(u=0\).

\(|hu|\)-Sweep bei festem \(h=10\). Je grösser der aufeinander zulaufende Impuls, desto höher das Plateau \(h^*\).

\(h\)-Sweep bei festem \(|hu|=5\). Aufgetragen ist \(h - h_0\) — also der Plateauaufbau. Bei \(h=20\) bleibt der Aufbau klein (~0,36 m), bei \(h=2\) wächst er auf über 1 m an. Dünneres Wasser reagiert also deutlich empfindlicher auf denselben Impuls.

Erkenntnisse.

• Ausbreitungsrichtung. Die Eigenwerte \(\lambda_{1,2} = u \mp \sqrt{g\,h}\) (Kapitel 1.3) geben die Schockgeschwindigkeiten direkt vor. Im Plateau gilt \(u=0\), also \(\lambda_{1,2} = \mp\sqrt{g\,h^*}\) — genau die beiden Schock-Geschwindigkeiten in den Plots.

• Plateauhöhe aus Massenerhaltung. Die Rankine-Hugoniot-Bedingung über die zwei symmetrischen Schocks liefert:

  \[hu = (h^* - h)\,\sqrt{\tfrac{1}{2} g\,h^*\,(h^* + h)/h}\]

  Für \(|hu|=5,\,h=10\) folgt \(h^* \approx 10{,}51\) — exakt das, was die Simulation zeigt (\(\Delta h \approx 0{,}51\), grüne Kurve im h-Sweep).

• Nichtlinearität. Eine Verdopplung des Impulses führt zu deutlich mehr als doppeltem Plateauanstieg (\(|hu|=5 \to \Delta h^* \approx 0{,}51\); \(|hu|=10 \to \Delta h^* \approx 1{,}03\)). Genauso hat der h-Sweep den umgekehrten Trend: dünneres Wasser verstärkt den Aufbau nichtlinear. Das ist die Energieaufstauung durch den Schock — kinetische Energie wird in potentielle Energie umgewandelt, und bei kleinem \(h\) ist mehr Headroom pro Einheit Masse verfügbar.

2.5.3. Rare-rare: divergierende Ströme

Zwei Verdünnungswellen wandern nach aussen. In der Mitte entsteht eine Mulde mit reduziertem Wasserstand.

Fünf Zeitscheiben. Die Geschwindigkeit \(u\) zeigt den divergierenden Strom am Anfang (Sprung bei \(x=5\)) und wandert allmählich gegen den Plateauwert \(u \approx 0\) an beiden Rändern.

\(|hu|\)-Sweep bei festem \(h=10\). Je grösser der divergierende Impuls, desto tiefer die Mulde. Für \(|hu| \to 99\) entsteht ein Dry-Bed.

\(h\)-Sweep bei festem \(|hu|=5\) (\(h - h_0\), negativ heisst Mulde). Bei \(h=2\) ist der Trog mit knapp 1 m am tiefsten — der Solver läuft hier bereits sehr nahe an die Dry-Bed-Schranke \(|u| < \sqrt{g h} \approx 4{,}4\) m/s.

Erkenntnisse.

• Dry-Bed-Schranke. Aus den Riemann-Invarianten der Verdünnungswelle folgt

  \[|u_r - u_l| < 2\sqrt{g\,h}\]

  Für \(h=10\) ergibt das \(|u| < \sqrt{g\,h} \approx 9{,}9\) m/s, also \(|hu| < 99\). In unserem \(|hu|\)-Sweep ist \(|hu|=80\) schon nah an der Grenze — das Plateau fällt auf \(h^* \approx 3{,}5\) m. Bei \(|hu|=100\) wird die Simulation instabil (negative Höhe in der Rarefaktion).

• h-Empfindlichkeit. Im h-Sweep zeigt sich, dass bei gleichem \(|hu|\) ein kleineres \(h\) eine tiefere Mulde erzeugt. Das ist das symmetrische Gegenstück zum shock-shock-Befund: dünneres Wasser reagiert stärker auf denselben Impuls, weil \(u = hu/h\) wächst und damit näher an die Dry-Bed-Schranke rückt.

• Symmetrie. Shock-shock und rare-rare sind Spiegelprobleme: das Vorzeichen des Impulses entscheidet, ob ein Schock (Druckaufbau, \(h^* > h\)) oder eine Rarefaktion (Druckabbau, \(h^* < h\)) entsteht.

2.6. 2.2 Dam-Break

Das klassische Dam-Break-Problem ist unser erstes asymmetrisches Setup. Bei \(t=0\) fällt der Damm zwischen einem hohen Reservoir (\(h_l\)) und einem niedrigen Flussbett (\(h_r\)), beide zunächst in Ruhe.

2.6.1. Klassischer Dam-Break

Nach links wandert eine Rarefaktionswelle (glatter Abfall von 10 auf \(h^*\)), nach rechts eine Schockwelle (steile Front von \(h^*\) auf 5). Die Geschwindigkeit im Mittelplateau ist positiv und konstant.

Zeitscheiben. Das charakteristische Profil aus Rarefaktion + Plateau + Schock ist bei \(t=0{,}32\) s gut erkennbar; bei \(t \geq 0{,}63\) s haben die Wellen beide Ränder bereits erreicht.

Endzustand für \(h_l = 10\) bei unterschiedlichen \(h_r\). Je grösser das Höhenverhältnis \(h_l/h_r\), desto stärker die Wellen — und desto schneller die Schockfront.

2.6.2. Einfluss der Flussströmung \(u_r\)

Die Aufgabenstellung fragt explizit nach dem Einfluss der Flussgeschwindigkeit \(u_r\). Da unser DamBreak1d-Setup mit \(hu=0\) arbeitet (der Sonderfall für die Höhenverhältnis-Analyse), beantworten wir die Frage mit dem analytischen wet-wet-Riemann-Solver aus scripts/plot_ch02.py — derselbe, der auch die Evakuierungszeit bestimmt. Als Basisfall dient die Reservoir-Konfiguration \(h_l=14,\,h_r=3{,}5,\,u_l=0\).

Links: Schockgeschwindigkeit \(s(u_r)\). Rechts: Mittelhöhe \(h^*(u_r)\). Rote Marker: unser Reservoir-Fall \(u_r = 0{,}2\) m/s.

Erkenntnisse.

• Struktur (Dry-Bed ausgeschlossen): Für \(h_l > h_r\) bildet sich immer links eine Rarefaktion und rechts ein Schock. Das folgt aus \(\lambda_1 < 0 < \lambda_2\) zusammen mit der Entropiebedingung.

• Einfluss von \(h_l/h_r\): Je grösser das Höhenverhältnis, desto grösser die Sprunghöhe \(h^* - h_r\) und damit die Schockgeschwindigkeit. Für \(h_r=1,\,2\) reichen die Wellen bis \(t=1{,}25\) s noch nicht ganz durch — der Rarefaktions-Fächer ist sichtbar.

• Einfluss von \(u_r\): Der Einfluss ist nahezu linear und vergleichsweise klein: über das Intervall \(u_r \in [-2, 3]\) m/s variiert die Schockgeschwindigkeit nur zwischen etwa 10 und 12,6 m/s. Analytisch:

  \[s \approx u_r + \sqrt{g\,h^*(h^*+h_r)/(2\,h_r)}\]

  Die Dominanz der Wurzel (bei \(h_l=14,\,h_r=3{,}5\) ist sie \(\approx 10{,}9\) m/s) erklärt, warum die geringe Reservoir-Strömung \(u_r = 0{,}2\) m/s die Evakuierungszeit nur um Sekunden verschiebt (rechnerisch: bei \(u_r = 0\) wäre \(t \approx 2293\) s, bei \(u_r = 0{,}2\) sind es \(2248\) s).

• h* sinkt mit \(u_r\): eine Strömung mit dem Schock (positives \(u_r\)) zieht Wasser nach rechts, das Plateau wird niedriger; eine Strömung gegen den Schock (negatives \(u_r\)) staut Wasser auf. Das ist physikalisch genau das, was man erwartet.

2.6.3. Reservoir → Dorf

Aufgabe. Ein unendliches Reservoir bricht in einen Fluss mit Dorf 25 km stromabwärts. Anfangsbedingungen aus der Aufgabenstellung:

\[\begin{split}q_l = \begin{bmatrix} 14 \\ 0 \end{bmatrix},\qquad q_r = \begin{bmatrix} 3{,}5 \\ 0{,}7 \end{bmatrix}\end{split}\]

2.6.3.1. Umsetzung

• Setup src/setups/Reservoir.cpp: hardcoded Werte, Dammposition \(x = 25\,000\) m.

• Main-Loop src/main.cpp: Domain \(50\,000\) m (Damm in der Mitte, Dorf am rechten Rand). Die Zeitschleife bricht ab, sobald sich die Höhe in der letzten Zelle (= Dorfposition) ändert — dann gibt main die abgelaufene Simulationszeit als Evakuierungszeit aus.

• Aufruf:

  ./build/tsunami_lab 2000 reservoir
  

2.6.3.2. Analytische Lösung (Referenzwert)

Das ist ein wet-wet Dam-Break mit Rarefaktion links + Schock rechts. Der Mittelstatus \((h^*, u^*)\) ergibt sich aus der Kopplung

\[u_l - 2\bigl(\sqrt{g\,h^*} - \sqrt{g\,h_l}\bigr) = u_r + (h^* - h_r)\sqrt{\tfrac{g(h^*+h_r)}{2\,h^*\,h_r}}\]

Numerisches Lösen (Bisektion in scripts/plot_ch02.py) liefert:

Grösse	Wert
Mittelhöhe \(h^*\)	7,64 m
Mittelgeschwindigkeit \(u^*\)	6,12 m/s
Schockgeschwindigkeit \(s\)	11,12 m/s
Evakuierungszeit (25 km / s)	≈ 2248 s ≈ 37 min 28 s

2.6.3.3. Visualisierung der Wellenfront

Analytische Dam-Break-Lösung auf dem 50-km-Gebiet, Damm bei 25 km (rot), Dorf bei 50 km (grün). Oben \(h(x,t)\) mit dem Rarefaktionsfächer nach links und der Schockfront nach rechts; unten die dazugehörige Teilchengeschwindigkeit. Die Schockfront erreicht das Dorf nach ≈ 37,5 min — das gleiche Ergebnis liefert die Simulation (Abschnitt Simulation unten).

Sechs Zeitscheiben derselben Lösung. Gut zu sehen: das Plateau \(h^* \approx 7{,}64\) m bewegt sich konstant mit \(s \approx 11{,}12\) m/s nach rechts, der Rarefaktionsfächer verbreitert sich linear in der Zeit.

2.6.3.4. Simulation (nach Rebuild)

Der Simulations-Ouput wird über main.cpp erzeugt, nachdem das Domain auf 50 km gesetzt und der Evakuierungszeit-Print eingebaut wurde:

scons
./build/tsunami_lab 2000 reservoir
# → Ausgabe enthält:
# wave arrived at village (rightmost cell) after simulation time 22XX s (~37.X min)

Das Plot-Skript kann den gleichen Lauf fahren und eine Simulations-GIF erzeugen; die analytische Lösung dient dabei als Referenz.

Erkenntnis — Evakuierungszeit

Nach dem Dammbruch hat das Dorf ca. 37 Minuten Vorwarnzeit, bevor der Schock eintrifft. Die Strömung \(u_r=0{,}2\) m/s (entspricht \(hu_r=0{,}7\)) beschleunigt die Welle nur geringfügig; dominant ist \(\sqrt{g\,h^*(h^*+h_r)/(2\,h_r)} \approx 10{,}9\) m/s.

Note

Nach Änderung von main.cpp (Domain \(50\,000\) m, Print der Evakuierungszeit) muss einmal gebaut werden:

scons
python3 scripts/plot_ch02.py --only=reservoir

Danach ist das reservoir-GIF aktualisiert und der von main.cpp gemessene Wert (Simulation) stimmt mit der analytischen Lösung überein.

2.7. Reproduzierbarkeit

Alle Plots und GIFs in diesem Kapitel werden von scripts/plot_ch02.py erzeugt. Voraussetzungen:

• vorgebautes build/tsunami_lab und build/validate_middle_states (scons einmal ausführen)

• Python mit matplotlib und Pillow (zum Beispiel in einem venv: python3 -m venv .venv; .venv/bin/pip install matplotlib)

Regenerieren:

python3 scripts/plot_ch02.py                    # alles
python3 scripts/plot_ch02.py --only=riemann     # nur shock-shock + rare-rare
python3 scripts/plot_ch02.py --only=dambreak    # nur dam-break
python3 scripts/plot_ch02.py --only=reservoir   # nur Reservoir + Evak-Zeit
python3 scripts/plot_ch02.py --only=middle_states \
        --validation-rows 2000 --validation-cells 100

2.8. Unit-Tests

Das WavePropagation-Patch hat zwei Tests in src/patches/WavePropagation1d.test.cpp:

[WaveProp1d]

Ein Zeitschritt auf einem 100-Zellen-Dam-Break (\(h_l=10,\,h_r=8\)). Die Net-Updates an der Dammkante sind aus dem Roe-Solver analytisch bekannt, alle anderen Zellen müssen im Steady-State bleiben.

[WaveProp1dMiddleStates]

Sanity-Sweep über die ersten 50 Zeilen von resources/middle_states.csv mit dem f-wave-Solver. Akzeptiert bis zu 5 % Ausreisser über 1 % relativem Fehler — in der Praxis sind es 0 %.

Ausführen:

./build/tests "[WaveProp1d]"
./build/tests "[WaveProp1dMiddleStates]"

2.9. Continuous Integration

Zwei Pipelines laufen auf jeden Push bzw. Merge-Request:

GitHub Actions (.github/workflows/main.yml) — pro Push auf main und als tägliche Cron-Job. Jobs:

• cppcheck Statische Analyse

• Release-Build + ./build/tests

• mode=debug+san + Tests (AddressSanitizer + UBSan)

• valgrind ./build/tests (Leak-Check im Debug-Build)

GitLab CI (.gitlab-ci.yml) — Build-Matrix, plus Doku-Pipeline: build:tests (scons + Catch2), build:sphinx, build:doxygen, pages (kombiniert Sphinx + Doxygen in ein Pages-Deployment).

Damit ist jedes Commit in main automatisch grün geprüft und die Online-Doku ist aktuell.

2.10. Beiträge

Person	Beitrag
Moritz Arnhold	csv parser / Solver-agnostisches WavePropagation1d (Template), ShockShock1d / RareRare1d / Reservoir-Setups, erste Middle-States-Prüfung
Moritz Martin	validate_middle_states-Binary, Middle-States-Sanity-Unit-Test, Fix der Ghost-Zellen-Initialisierung, Plot-Skript & Evakuierungszeit-Analytik, Kapitel-Dokumentation