.. _ch02:

#############################
Finite Volume Discretization
#############################

.. contents:: Inhaltsverzeichnis
   :local:
   :depth: 2

Einleitung
==========

Dieses Kapitel dokumentiert die Implementierung der **Finite-Volumen-Diskretisierung**
für die eindimensionalen Shallow Water Equations sowie die Einbettung unseres
:ref:`f-wave-Solvers <ch01>` 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


Finite-Volumen-Diskretisierung
==============================

Zellstruktur und Erhaltungsform
-------------------------------

Das Rechengebiet wird in :math:`N` Zellen gleicher Breite :math:`\Delta x` zerlegt.
Jede Zelle :math:`i` trägt einen **Zellmittelwert**

.. math::

   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 :math:`q = [h,\; hu]^T`.
Integriert man die SWE über eine Zelle und ein Zeitintervall :math:`[t_n, t_{n+1}]`,
ergibt sich die **konservative Update-Formel**:

.. math::

   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 :math:`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.

Net-Updates statt Fluss-Differenzen
-----------------------------------

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

.. math::

   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 :math:`A^-\Delta Q` und :math:`A^+\Delta Q` **einmal**
berechnet und beiden Nachbarzellen zugeordnet. Die innere Schleife in
``src/patches/WavePropagation1d.h`` iteriert daher über **Kanten**,
nicht über Zellen.

CFL-Bedingung
-------------

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

.. math::

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

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


Implementierung
===============

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
:math:`N+2` Einträge pro Array.

Solver-Umschaltung per Template
-------------------------------

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

.. code-block:: cpp

   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.

Outflow-Randbedingung
---------------------

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

.. code-block:: cpp

   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 :math:`q_l = q_r`, der Flusssprung
:math:`\Delta f = 0` — die Welle verlässt das Gebiet spurlos.
Das ist die **offene Randbedingung** und vermeidet künstliche Reflexionen.


Middle-States-Validierung
=========================

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

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 :math:`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:

.. code-block:: bash

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

Ergebnisse
----------

.. figure:: ../_static/ch02/middle_states_errors.png
   :alt: Histogramm der relativen Fehler über 2000 Riemann-Probleme
   :align: center
   :width: 90%

   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 :math:`10^{-7}`, also tief unter unserer
   Toleranz von 1 %.

.. admonition:: Erkenntnis
   :class: tip

   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.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.

Setup
-----

.. list-table::
   :header-rows: 1
   :widths: 20 40 40

   * - Szenario
     - **shock-shock**
     - **rare-rare**
   * - :math:`h_l,\,h_r`
     - :math:`h` (beidseitig gleich)
     - :math:`h` (beidseitig gleich)
   * - :math:`hu_l,\,hu_r`
     - :math:`+hu,\;-hu`  (kollidierend)
     - :math:`-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``

Shock-shock: kollidierende Ströme
---------------------------------

.. figure:: ../_static/ch02/shock_shock_evolution.gif
   :alt: Animation: shock-shock, h=10, |hu|=5
   :align: center
   :width: 85%

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

.. figure:: ../_static/ch02/shock_shock_evolution.png
   :alt: Shock-shock: fünf Zeitscheiben
   :align: center
   :width: 85%

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

.. figure:: ../_static/ch02/shock_shock_momentum_sweep.png
   :alt: Shock-shock: Plateauhöhe für unterschiedliche |hu|
   :align: center
   :width: 85%

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

.. figure:: ../_static/ch02/shock_shock_height_sweep.png
   :alt: Shock-shock: h-Sweep
   :align: center
   :width: 85%

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

**Erkenntnisse.**

- **Ausbreitungsrichtung.** Die Eigenwerte :math:`\lambda_{1,2} = u \mp \sqrt{g\,h}`
  (Kapitel 1.3) geben die Schockgeschwindigkeiten direkt vor.
  Im Plateau gilt :math:`u=0`, also :math:`\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:

  .. math::

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

  Für :math:`|hu|=5,\,h=10` folgt :math:`h^* \approx 10{,}51` — exakt das,
  was die Simulation zeigt (:math:`\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 (:math:`|hu|=5 \to \Delta h^* \approx 0{,}51`;
  :math:`|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
  :math:`h` ist mehr Headroom pro Einheit Masse verfügbar.

Rare-rare: divergierende Ströme
-------------------------------

.. figure:: ../_static/ch02/rare_rare_evolution.gif
   :alt: Animation: rare-rare, h=10, |hu|=5
   :align: center
   :width: 85%

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

.. figure:: ../_static/ch02/rare_rare_evolution.png
   :alt: Rare-rare: fünf Zeitscheiben
   :align: center
   :width: 85%

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

.. figure:: ../_static/ch02/rare_rare_momentum_sweep.png
   :alt: Rare-rare: Plateauhöhe für unterschiedliche |hu|
   :align: center
   :width: 85%

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

.. figure:: ../_static/ch02/rare_rare_height_sweep.png
   :alt: Rare-rare: h-Sweep
   :align: center
   :width: 85%

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

**Erkenntnisse.**

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

  .. math::

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

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

- **h-Empfindlichkeit.** Im h-Sweep zeigt sich, dass bei gleichem :math:`|hu|`
  ein kleineres :math:`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 :math:`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,
  :math:`h^* > h`) oder eine Rarefaktion (Druckabbau, :math:`h^* < h`)
  entsteht.


2.2 Dam-Break
=============

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

Klassischer Dam-Break
---------------------

.. figure:: ../_static/ch02/dam_break_evolution.gif
   :alt: Animation: dam-break, h_l=10, h_r=5
   :align: center
   :width: 85%

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

.. figure:: ../_static/ch02/dam_break_evolution.png
   :alt: Dam-break: fünf Zeitscheiben
   :align: center
   :width: 85%

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

.. figure:: ../_static/ch02/dam_break_hr_sweep.png
   :alt: Dam-break: h_r-Sweep
   :align: center
   :width: 85%

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

Einfluss der Flussströmung :math:`u_r`
--------------------------------------

Die Aufgabenstellung fragt explizit nach dem Einfluss der
Flussgeschwindigkeit :math:`u_r`. Da unser ``DamBreak1d``-Setup mit
:math:`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 :math:`h_l=14,\,h_r=3{,}5,\,u_l=0`.

.. figure:: ../_static/ch02/dam_break_ur_impact.png
   :alt: Dam-break: analytischer Einfluss der Flussgeschwindigkeit u_r
   :align: center
   :width: 95%

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

**Erkenntnisse.**

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

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

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

  .. math::

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

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

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

Reservoir → Dorf
----------------

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

.. math::

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

Umsetzung
~~~~~~~~~

- **Setup** ``src/setups/Reservoir.cpp``: hardcoded Werte, Dammposition
  :math:`x = 25\,000` m.
- **Main-Loop** ``src/main.cpp``: Domain :math:`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:

  .. code-block:: bash

     ./build/tsunami_lab 2000 reservoir

Analytische Lösung (Referenzwert)
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

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

.. math::

   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:

.. list-table::
   :header-rows: 1
   :widths: 30 70

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

Visualisierung der Wellenfront
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

.. figure:: ../_static/ch02/reservoir_evolution.gif
   :alt: Animation: Reservoir bricht, Welle läuft 25 km zum Dorf
   :align: center
   :width: 92%

   **Analytische Dam-Break-Lösung** auf dem 50-km-Gebiet, Damm bei 25 km (rot),
   Dorf bei 50 km (grün). Oben :math:`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).

.. figure:: ../_static/ch02/reservoir_evolution.png
   :alt: Analytische Lösung, mehrere Zeitpunkte
   :align: center
   :width: 92%

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

Simulation (nach Rebuild)
~~~~~~~~~~~~~~~~~~~~~~~~~

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

.. code-block:: bash

   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.

.. admonition:: Erkenntnis — Evakuierungszeit
   :class: important

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

.. note::

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

   .. code-block:: bash

      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.


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:

.. code-block:: bash

   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


Unit-Tests
==========

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

``[WaveProp1d]``
   Ein Zeitschritt auf einem 100-Zellen-Dam-Break (:math:`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:

.. code-block:: bash

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


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.


Beiträge
========

.. list-table::
   :header-rows: 1
   :widths: 30 70

   * - 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