Peano
The Shallow Water Equations

Table of Contents

First Exercise

Now that we have an example of how to implement an equation, we can try implementing our own equation, let's try out the shallow-water equations, or SWE for short:

\( \begin{equation*} \label{SWE_wo} \frac{\partial}{\partial t}\left( \begin{array}{lr} h \\ h u_1 \\ h u_2 \\ b \end{array} \right) + \nabla \underbrace{\begin{pmatrix} h u_1 & h u_2 \\ h u_1^2 + \frac{g*h^2}{2} & h u_1 u_2 \\ h u_2 u_1 & h u_2^2 + \frac{g*h^2}{2} \\ 0 & 0 \end{pmatrix}} = \vec{0} \end{equation*} \)

With eigenvalues:

\( \begin{equation*} \left( \begin{array}{lr} \lambda_1 \\ \lambda_2 \\ \lambda_3 \end{array} \right) = \left( \begin{array}{lr} u_n + \sqrt{g (h+b)} \\ u_n \\ u_n - \sqrt{g (h+b)} \\ \end{array} \right) \end{equation*} \)

These are an approximation for the movement of shallow non-viscous fluids, i.e., they assume that there is no pressure differential within the fluid, and no loss of velocity due to friction. Their components are the water height h, the momentum in x- and y- direction hu_1 and hu_2 respectively, and the bathymetry b, which is the height of the ground below the liquid. They also use the gravity g=9.81.

We want to simulate a so-called dam-break scenario, that is a scenario where a large mass of fluid spreads rapidly from an initial resting position, such as might happen when you jump into a pool, or when you jump into a bathtub, or when you jump into a puddle. A configuration file has already been created for you at tutorials/exahype2/swe/SWE.ipynb. Try running it, and then look into the generated implementation file SWESolver.cpp.

  • Exercises
    • First, let's implement some initial conditions. Let's start simple and set all values to 0, except for the initial water height, which should be 1.0 on the left half of the domain, and \bb 2.0 on the right half of the domain. You can access the position of the point you are at via the parameter x. (Hint: The size and position of your domain can be found in the configuration file.)
    • Second, let's do the same thing for the boundary conditions. Here you will need to set the values outside of the domain through the parameter $Qoutside$, set the velocities and the bathymetry to be equal to $0$, and the water height to $1.0$ everywhere.
    • Thirdly, fill out the \texttt{maxEigenvalue} function. Take care to use the proper velocity for whatever direction is being computed.\ (Hint: remember that the base unknowns are not $u$, but $hu$)
    • Next, let's implement the flux according to SWE_wo}, in addition to this, zero out the NCP, that is to say set \texttt{BTimesDeltaQ} to $0$ for all its components, after all we do not need an NCP here, do we?
    • Once all of this is done, re-compile the project using by running the Jupyterlab cell again, then run the executable \texttt{DamBreak} and go do some stretches while it runs. Maybe get a glass of water as well. Once you are finished stretching and the code is finished running, convert the output and visualize it, does it look like what you'd expect?

The initial conditions should look like this:

Solution to the First Exercise

If you get stuck or have a typo you cannot quite find, you can find an example solution at:

Extending SWE

So, now that you hopefully have a working simulation we have to admit something. When we described the shallow-water equations to you, we sort of lied. Well not really, it's just that the equation above is the form of the SWE without bathymetry. With bathymetry, it looks as follows (note that the gravity term was moved from the flux to the NCP):

\( \begin{equation*} \label{SWE_w} \frac{\partial}{\partial t}\left( \begin{array}{lr} h \\ h u_1 \\ h u_2 \\ b \end{array} \right) + \nabla \underbrace{\begin{pmatrix} h u_1 & h u_2 \\ h u_1^2 & h u_1 u_2 \\ h u_2 u_1 & h u_2^2 \\ 0 & 0 \end{pmatrix}}_{ flux } + \underbrace{\begin{pmatrix} 0 & 0\\ h g*\partial_x(b+h) & 0\\ 0 & h g*\partial_y(b+h) \\ 0 & 0 \end{pmatrix}}_{ ncp } = \vec{0} \end{equation*} \)

Let's modify our implementation to take this into account, and while we're at it change our scenario a little, and heck why not learn about mesh refinement too? Mesh refinement in ExaHyPE 2 is controlled via a so-called refinement criterion. Cells can be sub-divided into 3^d cells so long as this would not cause the sub-cells to be smaller than the lowest cell size as defined in the configuration.

  • Exercises
    • First, let's modify our flux to remove the term that shouldn't be in it, and implement our NCP. The NCP is implemented into the parameter BTimesDeltaQ, and depends on the gradients of the values, which are stored in the parameter deltaQ.
    • Then, let's modify our starting conditions. Let's try running a radial dam-break, that is a circular one. And let's have it be dependent on the bathymetry. For this, set the water height to be constant at 1.0, and the bathymetry to be 0.2 within a distance of 0.5 of the center of the domain, and 0.0 otherwise. The velocities can remain 0. (Hint: You can use the function tarch::la::norm2() to compute the norm of a vector)
    • With these parameters, re-compile and re-run the simulation.
    • Finally, let's try some adaptive mesh refinement. There is a function refinementCriterion within your implementation file which currently always returns ::exahype2::RefinementCommand::Keep. Conditionally set it to return ::exahype2::RefinementCommand::Refine if within a distance of 0.4 to 0.6 of the center of the domain. With this, re-compile the project, and run it one final time.

The initial conditions should look like this:

And once you are finished with these tasks, you can move on to The Euler Equations and Air Flow Over a Symmetric Airfoil, where you will build your own application from scratch.

Solution to the Second Exercise

If you get stuck or have a typo you cannot quite find, you can find an example solution at: