1. Riemann Solver

This chapter implements and tests the most basic part of our project: The f-wave solver for the one-dimensional shallow water equations. The shallow water equations are a system of nonlinear hyperbolic conservations laws with an optional source term:

(1.1)\[\begin{split}\begin{bmatrix} h \\ hu \end{bmatrix}_t + \begin{bmatrix} hu \\ hu^2 + \frac{1}{2}gh^2 \end{bmatrix}_x = S(x,t).\end{split}\]

\(x \in \mathbb{R}\) is the location in space and \(t \in \mathbb{R}^+\) time. The quantities \(q=[h, hu]^T\) are given by the height of the water column \(h(x,t)\) and the momentum \(hu(x,t)\), \(u\) is the particle velocity. \(g\) is the used gravity, typically we use \(g:=9.80665\text{m}/\text{s}^2\). \(f := [hu, hu^2 + \frac{1}{2}gh^2]^T\) is the flux function.

../_images/tsunami_quantities.svg

Fig. 1.1 Sketch of the quantities appearing in the one-dimensional shallow water equations.

Our source term \(S(x,t)\) will include the effect of space-dependent bathymetry (topography of the ocean) in future parts of the project. The introduction of additional forces, e.g., friction or the coriolis effect, as part of the term is possible. Fig. 1.1 illustrates all involved variables.

Note: As units we use meters (m) and seconds (s) for all computations.

1.1. Literature

We discuss the basics with respect to the numerical formulation, software, and development strategies in our meetings. Nevertheless, not all details can be covered in such a short time. For now, the following list of books, papers and guides are recommended for your personal studies:

1.2. Getting Started

Our project starts with a given code to jump-start your developments. First, we make sure that you can use the code and work on it collaboratively.

Tasks

  1. The project’s initial code is available from Github. Fork the git-repository for your work. You may do this in your preferred way, e.g., by using your own resources or by using a provider of your choice. Make sure that all team members have access to your fork.

  2. The initial code ships with an implementation of the Roe solver. The Roe solver is similar to the targeted f-wave solver which you will develop as part of Ch. 1.3. Make yourself familiar with the initial code. The software uses SCons as a build tool. Build the code, run the unit tests, run the solver, and visualize the results.

  3. Generate a Doxygen documentation from the code. In future, comment your code using Doxygen syntax, especially functions and function parameters.

  4. Use git’s version control features for all changes and write meaningful commit messages. Consider going through pull requests and code reviews for changes. Make yourself familiar with git.

  5. It is paramount to test new features of our software before going to production. For this, we’ll introduce unit tests through Catch2 whenever possible. Make yourself familiar with Catch2.

Warning

The last three tasks are ongoing and have to be fulfilled throughout the entire project. This means that you have to write meaningful documentation, meaningful commit messages and unit tests for every single part of the project!

1.3. F-wave Solver

The f-wave solver approximately solves the following Initial Value Problem (IVP) for the shallow water equations Eq. 1.1 over time:

(1.3.1)\[\begin{split}q(x,0) = \begin{cases} q_{l} \quad &\text{if } x < 0 \\ q_{r} \quad &\text{if } x > 0 \end{cases} \qquad q_l, q_r \in \mathbb{R}^+ \times \mathbb{R}.\end{split}\]

Theory shows that the solution, arising from the discontinuity at \(x=0\), consist of two waves. Each wave is either a shock or a rarefaction wave. The f-wave solver uses two shock waves to approximate the true solution.

First, we use the Roe eigenvalues \(\lambda^{\text{Roe}}_{1/2}\) in terms of the left and right quantities \(q_l\) and \(q_r\) with respect to position \(x=0\) to approximate the true wave speeds:

(1.3.2)\[\begin{split}\begin{aligned} \lambda^{\text{Roe}}_{1}(q_l, q_r) &= u^{\text{Roe}}(q_l, q_r) - \sqrt{gh^{\text{Roe}}(q_l, q_r)}, \\ \lambda^{\text{Roe}}_{2}(q_l, q_r) &= u^{\text{Roe}}(q_l, q_r) + \sqrt{gh^{\text{Roe}}(q_l, q_r)}, \end{aligned}\end{split}\]

where the height \(h^{\text{Roe}}\) and particle velocity \(u^{\text{Roe}}\) are given as:

\[\begin{split}\begin{aligned} h^{\text{Roe}}(q_l, q_r) &= \frac{1}{2} (h_l + h_r), \\ u^{\text{Roe}}(q_l, q_r) &= \frac{u_l \sqrt{h_l} + u_r \sqrt{h_r}}{\sqrt{h_l}+\sqrt{h_r}}. \end{aligned}\end{split}\]

Using the Roe eigenvalues we can define corresponding eigenvectors \(r_{1/2}^{\text{Roe}}\):

\[\begin{split}\begin{aligned} r_1^{\text{Roe}} &= \begin{bmatrix} 1 \\ \lambda^{\text{Roe}}_1 \end{bmatrix}, \\ r_2^{\text{Roe}} &= \begin{bmatrix} 1 \\ \lambda^{\text{Roe}}_2 \end{bmatrix}. \end{aligned}\end{split}\]

The decomposition of the jump in the flux function \(f\), \(\Delta f := f(q_r) - f(q_l)\), into the eigenvectors gives the waves \(Z_{1/2}\):

(1.3.3)\[\Delta f = \sum_{p=1}^2 \alpha_p r_p \equiv \sum_{p=1}^2 Z_p, \qquad \alpha_p \in \mathbb{R}.\]

The left “cell” \(\mathcal{C}_{l}\) is influenced by the left-going waves (\(\lambda_p < 0\)) and the right “cell” \(\mathcal{C}_r\) by the right-going waves (\(\lambda_p > 0\)). This leads to the definition of net updates which summarize the net effect of the waves to the left and right “cell”:

(1.3.4)\[\begin{split}\begin{split} A^- \Delta Q := \sum_{p:\{ \lambda_p^\text{Roe} < 0 \}} Z_p \\ A^+ \Delta Q := \sum_{p:\{ \lambda_p^\text{Roe} > 0 \}} Z_p \end{split}\end{split}\]

Note: The eigencoefficients \(\alpha_p\) in Equation (Eq. 1.3.3) are obtained by multiplying the inverse of the matrix of right eigenvectors \(R=[r_1^\text{Roe}, r_2^\text{Roe}]\) with the jump in fluxes:

\[\begin{split}\begin{bmatrix} \alpha_1 \\ \alpha_2 \end{bmatrix} = \begin{bmatrix} 1 & 1 \\ \lambda^{\text{Roe}}_1 & \lambda^{\text{Roe}}_2 \end{bmatrix}^{-1} \Delta f.\end{split}\]

Tasks

  1. Implement the f-wave solver for the homogeneous, i.e., \(S(x,t)=0\), shallow water equations. This implementation should be independent of the main code. The already implemented Roe solver in the files Roe.h, Roe.cpp and Roe.test.cpp might serve as a guideline.

    • Input values are the left state \(q_l = [h_l, (hu)_l]^T\) and the right state \(q_r = [h_r, (hu)_r]^T\).

    • Output values are the left- and right-going net updates: \(A^- \Delta Q\) and \(A^+ \Delta Q\).

  2. As discussed above, write meaningful unit-tests using Catch2 for all implemented features. Examples for the basic f-wave solver are:

    • Verification of the eigenvalue computation: You may derive a basic set of eigenvalues for given input values \(q_l\) and \(q_r\) manually.

    • Zero – with respect to machine precision – net updates in the case of steady states, e.g., \(q_l = q_r\).

    • Correctness tests for supersonic problems \(\lambda_{1/2}^\text{Roe} < 0 \vee \lambda_{1/2}^\text{Roe} > 0\): You may derive requirements using Eq. 1.3.2. Remark: This implies that one of the net updates is zero as stated in Eq. 1.3.4.