.. _ch:dimensional_splitting:
Two-Dimensional Solver
======================
In this part of out project, we'll cover the extension of our one-dimensional discretization to two spatial dimensions.
Then, we'll compare our two-dimensional solver to the one-dimensional one of the previous parts.
After this chapter, we have all numerics at hand and will be able to conduct tsunami simulations in the following parts.
Dimensional Splitting
---------------------
In their differential form, the two-dimensional shallow water equations are given as:
.. math::
:label: eq:swe2d
\begin{bmatrix}
h \\ hu \\ hv
\end{bmatrix}_t
+
\begin{bmatrix}
hu \\ hu^2 + \frac{1}{2}gh^2 \\ huv
\end{bmatrix}_x
+
\begin{bmatrix}
hv \\ huv \\ hv^2 + \frac{1}{2} g h^2
\end{bmatrix}_y
=
S(x,y,t).
:math:`t \in \mathbb{R}^+` is time. :math:`x \in \mathbb{R}` and :math:`y \in \mathbb{R}` are the Cartesian coordinates.
The quantities :math:`q = [h, hu, hv]^T` are given through the space-time dependent water height :math:`h(x,y,t)`, the momentum in :math:`x`-direction :math:`hu(x,y,t)` and the momentum in :math:`y`-direction :math:`hv(x,y,t)`. :math:`F := [ hu, hu^2 + \frac{1}{2}gh^2, huv ]^T` is the flux function in :math:`x`-direction, and :math:`G := [ hv, huv, hv^2 + \frac{1}{2} g h^2]^T` the one in :math:`y`-direction.
:math:`g` is the used gravity, usually we use :math:`g:=9.80665 \text{m}/\text{s}^2`.
As for the one-dimensional shallow water equations, we only consider the space-dependent effect of bathymetry :math:`B(x,y)` in the source term :math:`S(x,y) = [0, -ghB_x, -ghB_y]^T`.
The method of dimensional splitting allows us to apply our one-dimensional f-wave solver to the shallow water equations in two spatial dimensions.
For this we split the two-dimensional problem into two one-dimensional ones which we solve after each other in so-called sweeps.
The :math:`x`-sweep is given as:
.. math:: Q_{i,j}^* = Q_{i,j}^n - \frac{\Delta t}{\Delta x} \left( A^+ \Delta Q_{i-1/2,j} + A^- \Delta Q_{i+1/2,j} \right) \quad \forall i \in \{ 1, .., n \}, \; j \in \{ 0, .., n+1 \}.
The :math:`y`-sweep is applied to the intermediate values :math:`Q_{i,j}^*` and defined as:
.. math:: Q_{i,j}^{n+1} = Q_{i,j}^* - \frac{\Delta t}{\Delta y} \left( B^+ \Delta Q^*_{i,j-1/2} + B^- \Delta Q^*_{i,j+1/2} \right) \quad \forall i,j \in \{ 1, .., n \}.
As in the one-dimensional case, a cell :math:`\mathcal{C}_i` is influenced by all neighboring cells which are connected through the net-updates :math:`A^\pm \Delta Q_{i\mp1/2,j}` of the vertical edges :math:`x_{i\mp1/2,j}` and :math:`B^\pm \Delta Q^*_{i,j\mp1/2}` of the horizontal edges :math:`y_{i, j\mp1/2}`.
The time step restriction due to the CFL criterion for :math:`\Delta t` is similar to the one-dimensional case but has to be true in the :math:`x`- and :math:`y`-sweep:
.. math::
:label: eq:cfl2d
\Delta t < \frac{1}{2} \cdot \frac{\Delta x}{\lambda_x^\text{max}}\quad \wedge \quad \Delta t < \frac{1}{2} \cdot \frac{\Delta y}{\lambda_y^\text{max}},
where :math:`\lambda_x^\text{max}` is the maximum eigenvalue appearing in the :math:`x`-sweep and :math:`\lambda_y^\text{max}` the maximum eigenvalue in the :math:`y`-sweep.
**Note:** We continue using our constant time step derived in the setup phase.
This is valid as long as the CFL-criterion in :eq:`eq:cfl2d` is satisfied for all time steps.
.. _tasks-2-2:
.. admonition:: Tasks
#. Add support for two-dimensional problems to the solver. Implement dimensional splitting through a new class ``patches::WavePropagation2d``. Change other parts of the software as required but keep supporting one-dimensional settings.
#. Implement a circular dam break setup in the computational domain :math:`[-50, 50]^2` by using the following initial values:
.. math::
\begin{cases}
[h, hu, hv]^T = [10, 0, 0]^T &\text{if } \sqrt{x^2+y^2} < 10 \\
[h, hu, hv]^T = [5, 0, 0]^T \quad &\text{else}
\end{cases}
#. Illustrate your support for bathymetry in two dimensions by adding an obstacle to the computational domain of the circular dam break setup.
Stations
--------
When solving wave propagation problems, we are often times interested in output at specific points (or stations) of the computational domain.
For a given station :math:`s=(x,y)`, a quantity of interest, e.g., the deviation of the water colum from sea level :math:`\eta_s(t) = h_s(t) + b_s(t)` is then typically visualized, compared or analyzed.
An example workload, often requiring per-station output, is given through verification and validation.
In verification, we might run two solvers with the same inputs and compare the solution at a set of important stations.
In validation, we would run our solver, and compare the resulting synthetic data to observed data.
.. figure:: data_4/tohoku_dart_stations.png
:name: fig:tohoku_dart_stations
:align: center
:width: 100%
Map of some DART stations in the Pacific. The stations are shown as yellow diamonds with their respective ids.
:numref:`fig:tohoku_dart_stations` shows a set of stations which are of interest for the 2011 M 9.1 Tohoku tsunami event.
Specifically, the locations of some `Deep-ocean Assessment and Reporting of Tsunamis (DART) <https://nctr.pmel.noaa.gov/Dart/>`__ stations are given.
For these stations, we can `access <https://www.ncei.noaa.gov/maps/hazards/>`__ data of past tsunami events and thus do validation by comparing our numerical solutions to the measured data.
An example of such a comparison is given in :numref:`fig:tohoku_21419`.
.. figure:: data_4/tohoku_21419.svg
:name: fig:tohoku_21419
:align: center
:width: 100%
Comparison of synthetic data ("Computed Solution") to observed data ("Measured Data") at DART station 21419 for the 2011 M 9.1 Tohoku tsunami event.
.. admonition:: Tasks
#. Add a new class ``tsunami_lab::io::Stations`` which summarizes a collection of user-defined stations.
Each station has a name and a location in space.
All stations share the output frequency in seconds.
The output for each station should be written in a separate ASCII-CSV file.
#. Find a suitable way to provide the names and locations of each stations to your solver.
Further, implement a time step-independent output frequency for the stations.
Keep in mind, that this runtime-configuration should be extensible and usable in later parts of the project.
One possible solution are XML-based runtime configurations through the library `pugixml <https://pugixml.org/>`__.
#. Use a symmetric problem setup, e.g., a circular dam break problem, to compare your two-dimensional solver to your one-dimensional one at a set of stations.
#. (optional) Run a "convergence study", i.e., use a synthetic setup and decrease the size of your grid cells.
For example, for the computational domain :math:`\Omega = [0, 100\text{m}]^2` you could test mesh-spacings :math:`h \in \{ 50, 25, 10, 7.5, 5, 4, 3, 2, 1 \}`.
Compare your solution at a set of points, what do you observe?