Finite Volume Discretization
============================

This chapter assumes a finite one-dimensional domain :math:`\Omega := [a,b]; \; a,b \in \mathbb{R}; \; a < b` which is discretized by :math:`n` non-overlapping cells :math:`\mathcal{C}_i`: :math:`\Omega = \cup_{i=1}^n \mathcal{C}_i`.
In each cell :math:`i` we specify the quantities :math:`Q_i = (h, hu)_i` and the bathymetry :math:`b_i`: Again, :math:`h_i \in \mathbb{R}^+` is the water height, :math:`(hu)_i \in \mathbb{R}` the momentum, and :math:`b_i \in \mathbb{R}` the bathymetry relative to sea level.
Additionally, we have to specify proper boundary conditions which is done by the two ghost cells :math:`\mathcal{C}_0` and :math:`\mathcal{C}_{n+1}` neighboring the left and right boundary.
The obtained spatial discretization is called finite volume discretization.

With a given finite volume discretization we can define a set of :math:`n+1` edge-local Riemann problems:

.. math::

   q(x,t^n) =
         \begin{cases}
           Q_{i-1} \quad   &\text{if } x < x_{i-1/2} \\
           Q_{i} \quad &\text{if }     x > x_{i-1/2}
         \end{cases} \quad \forall i \in \{ 1, \ldots, n+1 \}.

The solution of these Riemann problems is obtained analogue to the IVP in :eq:`eq:rp`:
The f-wave solver determines two net-updates :math:`A^\mp \Delta Q_{i-1/2}` per edge :math:`x_{i-{1/2}}`.
The net-update :math:`A^- \Delta Q_{i-1/2}` is the net-effect of left-going waves and influences the quantities of :math:`\mathrm{C}_{i-1}`.
:math:`A^+ \Delta Q_{i-1/2}` summarizes the net-effect of the right-going waves and influences the quantities in :math:`\mathrm{C}_i`.

The final update formula from time :math:`t^n` to time :math:`t^{n+1} = t^n + \Delta t` using the time step :math:`\Delta t \in \mathbb{R}^+` is given as:

.. math:: Q_i^{n+1} = Q_i^n - \frac{\Delta t}{\Delta x} \left( A^+ \Delta Q_{i-1/2} + A^- \Delta Q_{i+1/2} \right)  \quad \forall i \in \{ 1, \ldots, n \}.

Thus a cell :math:`\mathrm{C}_i` is influenced by the right going waves -- summarized by the net-update :math:`A^+ \Delta Q_{i-1/2}` – of its left edge at position :math:`x_{i-1/2}` and the left going waves -- summarized by the net-update :math:`A^- \Delta Q_{i+1/2}` -- of its right edge at position :math:`x_{i+1/2}`. :math:`\Delta x \in \mathbb{R}^+` is the cell size which is constant in our implementations.

:math:`\Delta t` is chosen in a way such that the waves do not interact with each other.
This is implied by waves which do not cross the cell centers :math:`x_i`.
In that case stability of the method is guaranteed if the time step :math:`\Delta t` is chosen to satisfy the CFL-criterion:

.. math:: \Delta t < \frac{1}{2} \cdot \frac{\Delta x}{\lambda^\text{max}},

where :math:`\lambda^\text{max}` is maximum over all absolute values of the eigenvalues computed for the :math:`n+1` edges:

.. math:: \lambda^\text{max} = \max_{i=1}^{n+1} \left( \left|\lambda_{1/2}^{\text{Roe}}\right| \right)_{i-1/2}

**Note:** We do not use bathymetry in this part of the project, you can simply set a dummy
value in all cells: :math:`b_i = 0 \quad \forall i \in 0, \ldots, n+1`.
The initial code comes with outflow boundary conditions, which set
appropriate values in the ghost cells :math:`\mathcal{C}_0` and
:math:`\mathcal{C}_{n+1}` before every time step (see method
``setGhostOutflow()``). We will add bathymetry to our f-wave solver and
implement new boundary conditions soon.

.. admonition:: Tasks

    #. Integrate your f-wave solver into the one-dimensional implementation of a wave propagation patch given in files ``WavePropagation1d.h``, ``WavePropagation1d.cpp``, ``WavePropagation1d.test.cpp``.
       Find a good solution to switch between the provided Roe solver and your f-wave solver.

    #. The CSV-file `middle_states\.csv <https://scalable.uni-jena.de/assets/tsunami_lab/middle_states.csv>`_ contains a collection of constant middle states which arise immediately in the Riemann solution at the initial discontinuity.
       Use these middle states as a sanity check.

    #. Embed your solver into a continuous integration tool.
       Ensure to run at least the solver's unit tests after every commit to your git repository.
       You may use any suitable tool for this task including `GitHub Actions <https://github.com/features/actions>`__, `GitLab Runner <https://docs.gitlab.com/runner/>`__, `Travis CI <https://www.travis-ci.com/>`__, `Buildkite <https://buildkite.com/>`__, `GoCD <https://www.gocd.org/>`__ or `Jenkins <https://www.jenkins.io/>`__.

.. _ch:shock_rare:

Shock and Rarefaction Waves
---------------------------

**Shock-Shock Problem**:
Let's use our solver to solve shock-shock Riemann problems.
Imagine two streams of water which move in opposite directions and smash into each other at some position :math:`x_\text{dis}`.
The scenario is given by the following setup:

.. math::

   \begin{cases}
           Q_i = q_{l} \quad &\text{if } x_i \le x_\text{dis} \\
           Q_i = q_{r} \quad &\text{if }   x_i > x_\text{dis}
         \end{cases} \qquad q_l \in \mathbb{R}^+ \times \mathbb{R}^+, \; q_r \in \mathbb{R}^+ \times \mathbb{R}^-,

with initial conditions:

.. math::
   :label: eq:shockcond

   q_l=
           \begin{bmatrix}
             h_l \\ (hu)_l
           \end{bmatrix}, \quad
         q_r =
           \begin{bmatrix}
             h_r \\ (hu)_r
           \end{bmatrix} =
           \begin{bmatrix}
             h_l \\ -(hu)_l
           \end{bmatrix}.



.. _ch:rarerare:

**Rare-Rare Problem:**
We can setup rare-rare Riemann problems by two streams of water, which move away from each other at some position :math:`x_\text{dis}`.
The scenario is defined as:

.. math::
   \begin{cases}
           Q_i = q_{r} \quad &\text{if } x_i \le x_\text{dis} \\
           Q_i = q_{l} \quad &\text{if }   x_i > x_\text{dis}
         \end{cases} \qquad q_l \in \mathbb{R}^+ \times \mathbb{R}^+, \; q_r \in \mathbb{R}^+ \times \mathbb{R}^-,

with initial conditions identical to :eq:`eq:shockcond`.

.. admonition:: Tasks

    #. Implement the shock-shock and rare-rare problems as setups.
    #. Play
       around with different sets of initial water heights :math:`h_l` and
       particles velocities :math:`u_l`. What do you observe? Is there a
       connection to the wave speeds :math:`\lambda_{1/2} = u \mp \sqrt{gh}`
       in :numref:`ch:fwave`?


.. _ch:dambreak:


Dam-Break
---------

.. figure:: data_1/dam_break_initial_eng.svg
  :name: fig:dambreak_initial

  Illustration of the initial setting for the dam break problem.

Now, let's solve the dam-break problem.
You can imagine a water reservoir which is separated, as illustrated in :numref:`fig:dambreak_initial`, from a river by a dam initially.
We assume a total failure of the dam, thus nothing keeps the water from moving downstream:

.. math::

   \begin{cases}
           Q_i = q_{l} \quad &\text{if } x_i \le x_\text{dis} \\
           Q_i = q_{r} \quad &\text{if }   x_i > x_\text{dis}
         \end{cases} \qquad q_l \in \mathbb{R}^+ \times 0, \; q_r \in \mathbb{R}^+ \times \mathbb{R}^+, \; h_l > h_r

:numref:`fig:dambreak_solution` shows the analytical solution of
the dam-break problem, which consists of a rarefaction and a shock wave.

.. figure:: data_1/dam_break_solution_eng.svg
  :name: fig:dambreak_solution


  Illustration of the analytical solution of the dam break problem.


.. admonition:: Tasks

  #. Apply your solver to the dam-break setup and play around with different sets of initial water heights :math:`h_l` and :math:`h_r`.
     What do you observe? How large is the impact of the particle velocity :math:`u_r` in the river?
  #. Assume a water reservoir of unlimited size and a village 25 km downstream with initial values :math:`q_l=[14, 0]^T` and :math:`q_r=[3.5, 0.7]^T`.
     How much time do you have to evacuate the village in our model before the shock wave arrives?