.. _ch03: ###################################### Bathymetry & Boundary Conditions ###################################### .. contents:: Inhaltsverzeichnis :local: :depth: 2 Einleitung ========== Erweiterung des 1D-Lösers aus :ref:`ch02` um **Bathymetrie**, **Wet-Dry-Randbedingungen** und einen ersten **Tsunami-Lauf** entlang der japanischen Pazifikküste. Bisher haben wir nur Riemann-Probleme auf flachem Boden gelöst — jetzt kommen Topographie, Küste und ein realer Erdbeben-Auslöser dazu. Aufgabenstellung: https://scalable.uni-jena.de/opt/tsunami/chapters/assignment_3.html Bathymetrie & Quellterm ======================= Die Bathymetrie :math:`b(x)` beschreibt den Meeresboden auf Meereshöhe als Referenz: :math:`b \geq 0` ist Land, :math:`b < 0` ist Ozean. Bei variabler Bathymetrie kommt zur SWE ein **Quellterm** hinzu: .. math:: \partial_t q + \partial_x f(q) = \Psi(q, b), \qquad \Psi = \begin{bmatrix} 0 \\ -g\,h\,\partial_x b \end{bmatrix} Diskret an einer Zellkante :math:`(i-1/2)`: .. math:: \Delta x\,\Psi_{i-1/2} = \begin{bmatrix} 0 \\ -\tfrac{1}{2}\,g\,(b_r - b_l)\,(h_l + h_r) \end{bmatrix} F-Wave mit Quellterm -------------------- Der f-wave-Solver zerlegt **Flusssprung minus Quellterm** in die Eigenrichtungen: .. math:: \Delta f - \Delta x\,\Psi = \alpha_1 \, r_1 + \alpha_2 \, r_2 In ``src/solvers/FWave.cpp::eigencoefficients``: .. code-block:: cpp t_real l_deltaXPsi[2] = {0}; l_deltaXPsi[1] = -t_real(0.5) * g * (i_bmR - i_bmL) * (i_hL + i_hR); t_real l_diff[2] = {l_fJump[0] - l_deltaXPsi[0], l_fJump[1] - l_deltaXPsi[1]}; o_eigencoefficient1 = l_detInv * ( i_eigenvalue2 * l_diff[0] - l_diff[1] ); o_eigencoefficient2 = l_detInv * ( -i_eigenvalue1 * l_diff[0] + l_diff[1] ); Bathymetrie im Patch -------------------- ``WavePropagation1d`` hält ein zusätzliches Array ``m_bm`` mit :math:`N+2` Einträgen (zwei Ghost-Cells), value-initialisiert auf 0. Setter und Getter folgen derselben :math:`+1`-Offset-Konvention wie ``m_h``/``m_hu``. Die Ghost-Methoden verlängern die Bathymetrie an beiden Rändern (:math:`b_0 := b_1`, :math:`b_{n+1} := b_n`). CSV-Output ---------- Sobald Bathymetrie im Spiel ist, ist die **freie Oberfläche** :math:`\eta = h + b` interessanter als ``h`` allein. Der Writer gibt deshalb in der ``height``-Spalte :math:`h+b` aus und schreibt zusätzlich die Spalte ``bathymetry``: .. code-block:: text x,y,height,momentum_x,bathymetry 0.05,0.05,-2.22921e-05,4.41986,-2 Plot-Skripte rekonstruieren :math:`h = \eta - b` bei Bedarf. Lake-at-rest (well-balanced) ---------------------------- Eine Schlüssel-Eigenschaft: Bei :math:`h+b = \mathrm{const}` und :math:`hu = 0` gilt :math:`\Delta f = \Delta x\,\Psi` exakt; die :math:`\alpha_p` verschwinden, das Wasser bleibt bit-stationär. Zwei Tests decken das ab: ``[FWaveBathy]`` direkt am Solver, ``[WaveProp1dLakeAtRest]`` als Integrationstest über ein welliges Bett unter ruhendem Wasser. Reflektierende Randbedingungen ============================== Wet-Dry-Übergänge werden per Spec als reflektierende Wand modelliert (Spiegel-Regel an Kante :math:`(i-1/2)`): .. math:: h_i := h_{i-1}, \qquad (hu)_i := -\,(hu)_{i-1}, \qquad b_i := b_{i-1} Per-Seiten-Auswahl (Enum-API) ----------------------------- Domain-Ränder werden über eine kombinierbare API gesteuert: .. code-block:: cpp enum class BoundaryCondition { OUTFLOW, //!< zero-gradient: Ghost = innere Zelle REFLECTING //!< Wand: h und b kopieren, hu negieren }; wp.setGhostOutflow(); // beide Seiten outflow wp.setGhostReflection(); // beide Seiten reflektierend wp.setGhostBoundary(BoundaryCondition::OUTFLOW, BoundaryCondition::REFLECTING); // gemischt Wet-Dry im Inneren ------------------ Vor dem Riemann-Aufruf wird im Zeitschritt-Loop geprüft, ob beide benachbarten Zellen nass sind. Andernfalls werden die Zustände gespiegelt oder die Edge übersprungen: .. code-block:: cpp const bool l_dryL = l_hL_ref <= 0; const bool l_dryR = l_hR_ref <= 0; if( l_dryL && l_dryR ) continue; // dry-dry: Riemann undefiniert if( l_dryL ) { // wet-dry: dry-Seite spiegeln l_hL_ref = l_hR_ref; l_huL_ref = -l_huR_ref; l_bmL_ref = l_bmR_ref; } else if( l_dryR ) { l_hR_ref = l_hL_ref; l_huR_ref = -l_huL_ref; l_bmR_ref = l_bmL_ref; } /* solver.netUpdates(...) */ if( !l_dryL ) { /* update L */ } // dry stays dry if( !l_dryR ) { /* update R */ } Zwei Subtilitäten: 1. **Wet/Dry per Wassersäule**, nicht per Bathy-Vorzeichen. Klassische Riemann-Setups (DamBreak/ShockShock) haben :math:`b=0`, wären sonst "onshore". Der Test ``h > 0`` löst das. 2. **Updates nur auf nasse Zellen.** Sonst würde die Wet-Dry-Spiegelung netto Masse zur dry-Seite transportieren — die Forderung "dry stays dry" wäre verletzt. Demo: Shock-shock aus q_l + Wand -------------------------------- Setzt man :math:`q_l` überall, links Outflow und rechts Reflexion, so wirkt die Wand wie ein gespiegelter Riemann-Partner :math:`(h_l, -hu_l)` — das Gebiet "sieht" dasselbe Shock-Shock-Problem wie aus :ref:`ch02`. In der Mitte entsteht ein Plateau mit :math:`h^* > h_l` und :math:`u \approx 0`. ``[WaveProp1dSetGhostBoundary]`` führt den Demo-Lauf durch und verifiziert die erwarteten Eigenschaften. Hydraulische Sprünge ==================== Übergänge zwischen subkritischer (:math:`F < 1`) und superkritischer (:math:`F > 1`) Strömung über einem Bump testen die Quellterm-Implementierung. Die **Froude-Zahl** .. math:: F = \frac{u}{\sqrt{g\,h}} vergleicht Strömungs- mit Wellengeschwindigkeit. Bei :math:`F<1` können Druckwellen stromauf laufen; bei :math:`F>1` ist der Strom schneller als die Wellen — ein hydraulischer Sprung wird möglich. Setups ------ ``Subcritical`` und ``Supercritical`` initialisieren je eine Lake-at-rest- Konfiguration über einem identisch geformten parabolischen Bump in :math:`x \in (8, 12)`, mit unterschiedlichen Tiefen und Strömungen: .. list-table:: :header-rows: 1 :widths: 25 35 35 * - - **Subcritical** - **Supercritical** * - Bathy am Bump - :math:`-1{,}8 - 0{,}05\,(x-10)^2` - :math:`-0{,}13 - 0{,}05\,(x-10)^2` * - Bathy Flat-Shelf - :math:`-2` - :math:`-0{,}33` * - :math:`h(x,0)` - :math:`-b(x)` - :math:`-b(x)` * - :math:`hu(x,0)` - :math:`4{,}42` - :math:`0{,}18` * - Domain / Endzeit - :math:`(0, 25)` m / 200 s - :math:`(0, 25)` m / 200 s Maximale Froude-Zahl analytisch (am Bump-Peak, :math:`t=0`): - **Subcritical** (:math:`h_{\text{Peak}}=1{,}8`): :math:`F \approx 0{,}585` - **Supercritical** (:math:`h_{\text{Peak}}=0{,}13`): :math:`F \approx 1{,}227` Subkritischer Lauf ------------------ .. figure:: ../_static/ch03/subcritical_states.png :alt: Subkritisch: Initial vs. Final η/hu/b :align: center :width: 90% Bei :math:`t=200` s ist :math:`\eta` über dem Bump leicht abgesenkt (~10 cm), :math:`hu` nahezu konstant, das Bett unverändert. .. figure:: ../_static/ch03/subcritical_froude.png :alt: Subkritisches Froude-Profil :align: center :width: 85% :math:`F` bleibt überall unter 1 — wie analytisch vorhergesagt. .. figure:: ../_static/ch03/subcritical_evolution.gif :alt: Animation Subcritical :align: center :width: 85% Aus dem Initial-Zustand bildet sich rasch ein stationäres Profil mit kleiner :math:`\eta`-Senke über dem Bump. Keine Wellen, keine Sprünge. Superkritischer Lauf -------------------- .. figure:: ../_static/ch03/supercritical_states.png :alt: Superkritisch: Initial vs. Final mit hydraulischem Sprung :align: center :width: 90% :math:`\eta` fällt am Bumpaufstieg deutlich (Strömung beschleunigt), stromabwärts der Bump-Spitze (:math:`x \approx 11{,}7` m) springt :math:`\eta` abrupt nach oben — der erwartete **hydraulische Sprung**. .. figure:: ../_static/ch03/supercritical_froude.png :alt: Superkritisches Froude-Profil :align: center :width: 85% :math:`F` startet bei :math:`\approx 0{,}3` (subkritisch im Flat-Shelf), schiesst über dem Bump auf :math:`\approx 2{,}6` (stark superkritisch) und kollabiert beim hydraulischen Sprung wieder unter 1. .. figure:: ../_static/ch03/supercritical_evolution.gif :alt: Animation Supercritical :align: center :width: 85% Einschwingvorgang: aus dem Initial-Zustand entsteht über die ersten ~10 s das stationäre Profil mit hydraulischem Sprung. Momentum-Drift -------------- Bei einer stationären Strömung sollte :math:`hu(x)` analytisch konstant sein. Der f-wave-Solver hält das **nicht exakt** ein: .. figure:: ../_static/ch03/supercritical_momentum_drift.png :alt: Momentum-Drift bei x = 20 m über die Zeit :align: center :width: 85% Momentum :math:`hu` an einer Zelle stromabwärts des Bumps. Initial 0,18; nach Einschwingen pendelt es sich bei :math:`\approx 0{,}125` ein (~30 % Drift). .. admonition:: Warum? :class: tip Der Drift stammt aus der grundlegenden Schwäche des f-wave-Solvers für nicht-Lake-at-rest-stationäre Lösungen: er ist well-balanced **nur** für :math:`hu \equiv 0`. Für allgemeinere stationäre Profile bräuchte es Augmented-Roe oder HLLE-Plus. Tsunami-Simulation ================== Erster echter Tsunami-Lauf entlang einer realen Bathymetrie aus dem **GEBCO_2025-Grid**. Datenextraktion (GMT ``grdtrack``) ---------------------------------- Per Aufgabenstellung verwenden wir ein 1D-Profil zwischen :math:`p_1 = (141{,}024949,\;37{,}316569)` (Ostküste Fukushima) und :math:`p_2 = (146{,}0,\;37{,}316569)` (offener Pazifik), mit 250 m Sampling. Extraktion mit GMT, abgelegt unter ``resources/GEBCO_2026_sub_ice_bathy.csv`` (1761 Datenzeilen, ``distance ∈ [0, 439{,}9] km``). ``setups/TsunamiEvent1d.cpp`` erkennt automatisch, ob die ``distance``-Spalte Meter oder Kilometer enthält, und konvertiert nach Meter. .. figure:: ../_static/ch03/tsunami1d_bathymetry.png :alt: Bathymetrie p1 → p2 :align: center :width: 95% Bathymetrie entlang der 440 km. Die Küste startet bei +14 m, der Schelf fällt sanft auf :math:`-1\,000` m, dann bricht das Profil in den **Japan-Trench** ab (~:math:`-7\,800` m bei :math:`x \approx 240` km). Rote Zone: Bereich des Erdbeben-Displacements. Setup-Gleichungen ----------------- Per Spec wird die Roh-Bathymetrie :math:`b_{\text{in}}(x)` mit einer Regularisierung :math:`\delta = 20` m und einem künstlichen Erdbeben- Displacement :math:`d(x)` kombiniert: .. math:: h(x, 0) &= \begin{cases} \max(-b_{\text{in}},\,\delta) & b_{\text{in}} < 0 \\ 0 & \text{sonst} \end{cases} b(x) &= \begin{cases} \min(b_{\text{in}},\,-\delta) + d(x) & b_{\text{in}} < 0 \\ \max(b_{\text{in}},\,+\delta) + d(x) & \text{sonst} \end{cases} d(x) &= \begin{cases} A\,\sin\!\Bigl(\tfrac{x - x_0}{(x_1-x_0)/2}\,\pi + \pi\Bigr) & x_0 < x < x_1 \\ 0 & \text{sonst} \end{cases} Standardwerte: :math:`A = 10` m, :math:`x_0 = 175\,000` m, :math:`x_1 = 250\,000` m. Das :math:`\delta`-Clamp verhindert, dass der Solver an der Küste durch :math:`u = hu/h`-Divisionen instabil wird. CLI: parametrisiertes Displacement ---------------------------------- ``TsunamiEvent1d`` nimmt :math:`A`, :math:`x_0`, :math:`x_1` optional als Konstruktor-Argumente; ``main.cpp`` reicht sie über die ``ARG``-Form durch: .. code-block:: bash ./build/tsunami_lab 1000 tsunami1d # Defaults ./build/tsunami_lab 1000 tsunami1d 20 175000 250000 # doppelte A ./build/tsunami_lab 1000 tsunami1d 10 250000 325000 # weiter offshore Anfangszustand -------------- .. figure:: ../_static/ch03/tsunami1d_initial.png :alt: Initial state :align: center :width: 95% Anfangszustand bei :math:`t=0`. Außerhalb der roten Zone ist :math:`\eta = 0` (Lake-at-rest auf Meeresspiegel); innerhalb hebt/senkt der Sinus-Displacement das Wasser um :math:`\pm 10` m. :math:`hu \equiv 0`. Wellenausbreitung & Runup ------------------------- .. figure:: ../_static/ch03/tsunami1d_evolution.gif :alt: Tsunami-Animation :align: center :width: 95% 1 h Simulation. Der Erdbeben-Bump zerfällt in zwei Wellen: eine **west-gehende** Richtung Küste, eine **ost-gehende** Richtung Pazifik. .. figure:: ../_static/ch03/tsunami1d_evolution.png :alt: Tsunami: Mehrere Zeitscheiben :align: center :width: 95% Sechs Zeitscheiben über 60 min. Bei :math:`t \approx 12` min sind die Wellen separiert; bei :math:`t \approx 36` min wird die west-gehende Welle durch **Shoaling** stark amplifiziert; bei :math:`t \approx 48` min trifft sie die Küste; danach Reflexion im Küstennahbereich. .. figure:: ../_static/ch03/tsunami1d_runup.png :alt: Runup-Zoom :align: center :width: 90% Zoom auf die ersten 30 km vor der Küste. Maximale Amplitude :math:`\eta \approx +3` m kurz vor der Küste; nach Reflexion schwingt das Wasser mit :math:`\approx 1{,}5` m im Küstenstreifen weiter. .. admonition:: Was wir beobachten — und was nicht :class: important - **Sehen wir.** Aufspaltung in zwei Wellen → Shoaling → Amplifikation auf :math:`\approx 3` m → Reflexion an der Küstenwand. - **Sehen wir nicht.** Echtes Auflaufen über Land. Per Spec ist die Küste eine perfekt reflektierende Wand-Cell. Auch das :math:`\delta = 20` m-Clamp glättet die Wassertiefe nahe Null mathematisch und dämpft die Amplitude unmittelbar an der Küste. Variation des Displacements --------------------------- .. figure:: ../_static/ch03/tsunami1d_variants_initial.png :alt: Vier Displacement-Varianten — Initial :align: center :width: 95% Vier Varianten: Default :math:`(A, x_0, x_1)=(10, 175, 250)` km, halbe und doppelte Amplitude, sowie eine weiter östlich verschobene Quellzone. .. figure:: ../_static/ch03/tsunami1d_variants_t30min.png :alt: Vergleich nach 30 min :align: center :width: 95% Wellenfeld nach :math:`\approx 30` min Simulation. **Erkenntnisse.** - **Lineare Amplituden-Skalierung.** Halbe Amplitude → ~halbe Welle an der Küste; doppelte → ~doppelte. Im flachen Bereich ist die SWE näherungsweise linear. - **Verspätete Ankunft bei verschobener Quelle.** Die "weiter offshore"- Variante liegt bei :math:`t = 30` min noch deutlich östlicher — 75 km zusätzliche Laufstrecke. - **Nichtlineares Shoaling.** Die :math:`A=20`-Variante hat an der Küste überproportional starke Spitzen — Shoaling geht über abnehmende Wassertiefe quadratisch in die Amplitude ein. Reproduzierbarkeit ================== Alle Plots werden von ``scripts/plot_ch03.py`` erzeugt. Voraussetzungen: - vorgebautes ``build/tsunami_lab`` (``scons`` einmal ausführen) - Python mit ``matplotlib`` und ``Pillow`` - ``resources/GEBCO_2026_sub_ice_bathy.csv`` muss vorhanden sein .. code-block:: bash python3 scripts/plot_ch03.py # alles python3 scripts/plot_ch03.py --only=subcritical python3 scripts/plot_ch03.py --only=supercritical python3 scripts/plot_ch03.py --only=tsunami1d python3 scripts/plot_ch03.py --only=tsunami_variants Unit-Tests ========== .. list-table:: :header-rows: 1 :widths: 30 70 * - Tag - Inhalt * - ``[FWaveBathy]`` - Solver-Ebene: Δb≠0 mit identischem q liefert Source-Term-Update; Lake-at-rest mit nicht-trivialem b liefert exakt null. * - ``[WaveProp1dLakeAtRest]`` - Welliges Bett unter ruhendem Wasser, 50 Zeitschritte. Tolerance :math:`10^{-3}` für :math:`\eta` und :math:`hu`. * - ``[WaveProp1dSetGhostBoundary]`` - Demo aus dem Boundary-Conditions-Abschnitt: :math:`q_l` überall, outflow links + reflektierend rechts → Shock-Shock-Plateau. * - ``[WaveProp1dWetDryInterface]`` - Insel in einer Strömung. Dry stays dry, Wet-Side bremst, Bathy unverändert. * - ``[Subcritical]`` / ``[Supercritical]`` - Setup-Werte gegen Spec-Formeln + analytische Froude-Zahl am Bump. * - ``[TsunamiEvent1d]`` - Spec-Eigenschaften: Momentum=0, :math:`\delta`-Clamps, freie Oberfläche bei z=0 außerhalb der Displacement-Zone, Küstenzelle dry. Insgesamt **25 Test-Cases / 779 Assertions**, alle grün im ``mode=debug+san``-Build. Beiträge ======== .. list-table:: :header-rows: 1 :widths: 30 70 * - Person - Beitrag * - Moritz Arnhold - F-Wave-Quellterm, Bathymetrie-Storage in ``WavePropagation1d``, Wet-Dry-Reflexion im timeStep-Loop, Subcritical / Supercritical / TsunamiEvent1d-Setups, GEBCO-CSV mit ``grdtrack``, SConscript * - Moritz Martin - ``BoundaryCondition``, Bug-Fixes in ``WavePropagation1d``, TsunamiEvent1d-Parametrisierung, Plot-Skript ``plot_ch03.py``, Test-Suite-Erweiterung, Kapitel-Dokumentation