.. _ch01: ############### Riemann Solver ############### .. contents:: Inhaltsverzeichnis :local: :depth: 2 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 Mathematische Grundlagen ======================== Shallow Water Equations ----------------------- Unser Ausgangspunkt sind die eindimensionalen Shallow Water Equations — ein System aus zwei Erhaltungsgleichungen, das die Dynamik von Flachwasser beschreibt: .. math:: \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) Die erste Zeile beschreibt die **Massenerhaltung**: die Wasserhöhe :math:`h` ändert sich, wenn Wasser mit dem Impuls :math:`hu` zu- oder abfliesst. Die zweite Zeile beschreibt die **Impulserhaltung**: der Impuls :math:`hu` ändert sich durch den konvektiven Transport :math:`hu^2/h` und den hydrostatischen Druck :math:`\frac{1}{2}g h^2`. Wir definieren dazu zwei zentrale Grössen: - **Erhaltungsgrössen:** :math:`q = [h,\; hu]^T` — Wasserhöhe und Impuls - **Flussfunktion:** :math:`f(q) = [hu,\; hu^2/h + \frac{1}{2}g h^2]^T` Für Assignment 1 betrachten wir den homogenen Fall :math:`S(x,t) = 0` (keine Bathymetrie, keine Reibung). Das Riemann-Problem ------------------- An jeder Zellkante treffen zwei Zustände aufeinander — links :math:`q_l` und rechts :math:`q_r`. Diese Unstetigkeit nennt man ein **Riemann-Problem**: .. math:: q(x, 0) = \begin{cases} q_l & \text{wenn } x < 0 \\ q_r & \text{wenn } x > 0 \end{cases} 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. F-wave Solver — Theorie ======================== Der f-wave Solver löst das Riemann-Problem in drei Schritten: 1. Wellengeschwindigkeiten bestimmen (Eigenwerte) 2. Flusssprung in Eigenvektoren zerlegen (Eigenkoeffizienten) 3. Wellen den Zellen zuordnen (Net-Updates) 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: .. math:: 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}} Daraus ergeben sich die **Eigenwerte** (Wellengeschwindigkeiten): .. math:: \lambda_1 &= u^\text{Roe} - \sqrt{g \cdot h^\text{Roe}} \lambda_2 &= u^\text{Roe} + \sqrt{g \cdot h^\text{Roe}} Die physikalische Bedeutung: :math:`\lambda_1` ist die Geschwindigkeit der nach links laufenden Welle, :math:`\lambda_2` die der nach rechts laufenden. Der Term :math:`\sqrt{g \cdot h^\text{Roe}}` ist die lokale Wellenausbreitungsgeschwindigkeit. .. admonition:: Im Code :class: tip Die Eigenwerte werden in der Funktion ``eigenvalues()`` berechnet. Die Wurzel :math:`\sqrt{g}` ist als Konstante ``gSqrt`` in ``constants.h`` vordefiniert, sodass nur noch :math:`\sqrt{h^\text{Roe}}` berechnet werden muss. Schritt 2: Eigenkoeffizienten (Flusssprung-Zerlegung) ------------------------------------------------------ Dies ist der **Kernschritt** und gleichzeitig der wesentliche Unterschied zum Roe-Solver. Wir berechnen zunächst die Flussfunktion :math:`f(q)` für beide Seiten: .. math:: 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} und daraus den **Flusssprung**: .. math:: \Delta f = f(q_r) - f(q_l) .. important:: **Unterschied zum Roe-Solver:** Der Roe-Solver zerlegt den Sprung in den *Zustandsgrössen* :math:`\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* :math:`\Delta f`. Diesen Flusssprung zerlegen wir nun in die Eigenvektoren :math:`r_1, r_2` der Roe-Matrix: .. math:: \Delta f = \alpha_1 \cdot r_1 + \alpha_2 \cdot r_2 mit den Eigenvektoren: .. math:: r_1 = \begin{bmatrix} 1 \\ \lambda_1 \end{bmatrix}, \qquad r_2 = \begin{bmatrix} 1 \\ \lambda_2 \end{bmatrix} In Matrixform geschrieben ist das ein lineares Gleichungssystem: .. math:: \underbrace{ \begin{bmatrix} 1 & 1 \\ \lambda_1 & \lambda_2 \end{bmatrix} }_{R} \begin{bmatrix} \alpha_1 \\ \alpha_2 \end{bmatrix} = \Delta f Die Lösung erhalten wir durch Invertieren der 2×2-Matrix :math:`R`: .. math:: R^{-1} = \frac{1}{\lambda_2 - \lambda_1} \begin{bmatrix} \lambda_2 & -1 \\ -\lambda_1 & 1 \end{bmatrix} Damit ergibt sich: .. math:: \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) .. admonition:: Im Code :class: tip In der Funktion ``eigencoefficients()`` werden die drei Teilschritte ausgeführt: 1. **Flussfunktion auswerten** 2. **Flusssprung bilden** 3. **R⁻¹ · Δf berechnen** Schritt 3: Net-Updates ----------------------- Aus den Eigenkoeffizienten :math:`\alpha_p` und Eigenvektoren :math:`r_p` bilden wir die **f-Wellen**: .. math:: 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} Jede f-Welle :math:`Z_p` wird der linken oder rechten Zelle zugeordnet — abhängig vom Vorzeichen des zugehörigen Eigenwertes :math:`\lambda_p`: .. math:: 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)} Anschaulich: Wellen, die sich nach links ausbreiten (:math:`\lambda_p < 0`), verändern die linke Zelle. Wellen, die sich nach rechts ausbreiten (:math:`\lambda_p > 0`), verändern die rechte Zelle. .. admonition:: Im Code :class: tip ``netUpdates()`` berechnet mit den anderen Funktionen die Updates: 1. ``eigenvalues()`` berechnen 2. ``eigencoefficients()`` berechnen 3. f-Wellen berechnen 4. basierend auf Vorzeichen auf Update für rechte o. linke Zelle addieren Sonderfälle ----------- **Steady State** (:math:`q_l = q_r`): Wenn beide Seiten identisch sind, ist :math:`\Delta f = 0`, also :math:`\alpha_1 = \alpha_2 = 0` und alle Net-Updates verschwinden. Die Simulation bleibt stabil — es passiert nichts. **Supersonische Strömung** (:math:`\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 :math:`\lambda_1 > 0` und :math:`\lambda_2 > 0` laufen beide Wellen nach rechts, d.h. :math:`A^- \Delta Q = 0` (kein Update links). Unterschied zum Roe-Solver =========================== Beide Solver verwenden dieselben Eigenwerte und Eigenvektoren. Der entscheidende Unterschied liegt in der Zerlegung: .. list-table:: :header-rows: 1 :widths: 20 40 40 * - - **Roe-Solver** - **F-wave Solver** * - Zerlegt wird - :math:`\Delta q = q_r - q_l` - :math:`\Delta f = f(q_r) - f(q_l)` * - Gleichungssystem - :math:`R \cdot \alpha = \Delta q` - :math:`R \cdot \alpha = \Delta f` * - Wellen - :math:`\text{wave}_p = \lambda_p \cdot \alpha_p \cdot r_p` - :math:`Z_p = \alpha_p \cdot r_p` * - Skalierung - Zusätzlich mit :math:`\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. 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 :math:`\lambda_{1,2}`) - ``eigencoefficients()`` → Schritt 2 (Flusssprung-Zerlegung :math:`\alpha_{1,2}`) - ``netUpdates()`` → Schritt 3 (f-Wellen bilden und zuordnen) ``FWave.test.cpp`` Unit-Tests mit Catch2 (siehe nächster Abschnitt). Unit-Tests ========== Unsere Tests decken die wesentlichen Szenarien ab: Eigenwert-Verifikation ----------------------- Test mit :math:`h_l = 10,\; h_r = 9,\; u_l = -3,\; u_r = 3`: Die Roe-Mittelwerte ergeben sich zu :math:`h^\text{Roe} = 9.5` und :math:`u^\text{Roe} \approx -0.079`. Daraus folgen :math:`\lambda_1 \approx -9.731` und :math:`\lambda_2 \approx 9.573`. Diese Werte werden gegen die analytische Lösung geprüft. Eigenkoeffizienten-Verifikation -------------------------------- Für denselben Testfall verifizieren wir, dass die Zerlegung konsistent ist: :math:`\alpha_1 + \alpha_2 = \Delta f_0` und :math:`\alpha_1 \lambda_1 + \alpha_2 \lambda_2 = \Delta f_1`. Das stellt sicher, dass :math:`R \cdot \alpha = \Delta f` tatsächlich erfüllt ist. Steady State ------------ Bei :math:`q_l = q_r` (z.B. :math:`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 (:math:`hu = 0`) als auch für einen nicht-trivialen Gleichgewichtszustand (:math:`h = 5,\; hu = 3` auf beiden Seiten). Dambreak -------- Symmetrischer Dammbruch (:math:`h_l = 10,\; h_r = 8,\; hu = 0`): Da keine Strömung vorliegt, sind die Eigenwerte betragsmässig gleich (:math:`\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. Supersonisch ------------ Wir konstruieren Fälle, in denen beide Eigenwerte dasselbe Vorzeichen haben: - **Alle Wellen nach rechts** (:math:`h = 0.1,\; u = 10` → :math:`\lambda_{1,2} > 0`): Das linke Net-Update muss exakt null sein, das rechte muss dem gesamten Flusssprung entsprechen. - **Alle Wellen nach links** (:math:`h = 0.1,\; u = -10` → :math:`\lambda_{1,2} < 0`): Spiegelbildlicher Fall — das rechte Net-Update ist null. Tests Ausführen --------------- .. code-block:: bash # Alle FWave Tests ./build/tests "[FWave]" # einzelnen Testcases ./build/tests "[FWaveEigenvalues]" ./build/tests "[FWaveEigencoefficients]" ./build/tests "[FWaveSteadyState]" ./build/tests "[FWaveDamBreak]" ./build/tests "[FWaveSupersonic]" Beiträge ======== .. list-table:: :header-rows: 1 :widths: 30 70 * - Person - Beitrag * - Moritz Arnhold - Build-Setup, Dokumentation, Code Review * - Moritz Martin - FWave-Implementierung, FWave-Tests, Dokumentation