Home

On this site we want to give a short introduction on how to implement a new set of equations with the ExaHyPE Teaching code.

The Shallow Water Equations

The system of hyperbolic PDEs we regard are the Shallow Water Equations (wiki) The Shallow Water Equations (SWE) model the flow of flat water. The equations neglect any vertical velocities, which is as sufficient approximation as long as the the height of the water is significantly smaller than the regarded domain. We thus only need to consider the water height h and velocities in x any y dimension, u and v.

In conservative form the SWE are defined as:

$$\frac{\delta}{\delta t} \vec{Q} + \frac{\delta}{\delta x} \vec{F_1}\left(\vec{Q}\right) + \frac{\delta}{\delta y} \vec{F_2}\left(\vec{Q}\right) = 0,$$

where

$$\vec{Q}=\begin{pmatrix} h\\\ hu\\\ hv \end{pmatrix} \textrm{, } \left(\vec{F_1}\left(\vec{Q}\right)\textrm{ }\vec{F_2}\left(\vec{Q}\right)\right) =\begin{pmatrix} hu & hv \\\ \frac{1}{2}gh^2 + hu^2 & huv \\\ huv & \frac{1}{2}gh^2 + hv^2 \end{pmatrix}$$

The flux models the conservation of mass and momentum and considers hydrostatic pressure. The gravitational constant g will be set to 1.

The eigenvalues of the PDE system are needed later on for the computation of an admissible time step size. For the SWE in two dimension they are:

$$u\textrm{, }u-\sqrt{g\cdot h}\textrm{ and } u+\sqrt{g\cdot h}$$

Setting up the specification file

After the equation is defined we need to set up the specification file which for all applications carries the suffix *.exahype. In the folder Demonstrators/TrivialProject you can find a template TrivialProject.exahype for an application that is solved by the ADER-DG method.

Set up a new directory SWE and copy the template there (We name it SWE_ADERDG.exahype). In this file all place holders are marked by ## and need to be filled to get a working application. (All other options are only used for optimizations and don't need to be touched or understood):

`exahype-project`

The name of the project. We use SWE.

exahype-project SWE

`output-directory const`

Defines where the the code body will be generated relative to the root directory
In our case we use ./Demonstrators/SWE

output-directory const  = ./Demonstrators/SWE

`computational-domain`

Here we define the number of dimensions, the simulated domain and the physical end time of the application
The SWE are defined in two dimensions we want to do our simulation on the cube [-0.5,0.5]^2 and for one second:

  computational-domain
    dimension const          = 2
    width                    = 1.0, 1.0
    offset                   = -0.5, -0.5
    end-time                 = 1.0
  end computational-domain

`solver ADER-DG` ...

This region defines the characteristics of the simulated PDE and the used solver.
variables const By looking at the formula you see that the SWE consist of 3 quantities 'h,hu,hv' which we will denote with h and a vector p of length 2.
order const The polynomial degreee of the high order representation. We'll use order 3.
maximum-mesh-size The maximum size of a single cell. ExaHyPE always has to hold a power of 3 cells in one dimension¹. By setting this value you define an lower limit for the number of cells which is the rounded up to the next power of 3. For our example setting maximum-mesh-size to 0.05 gives us a theoretical number of 20 cells in one dimension. The next power of three is 27, so our mesh will end up with 27 cells in each dimension and a mesh size of 1/27
type const Here you can choose between a linear or nonlinear PDE. SWE are nonlinear (as the flux is).
terms const The type of terms that occur in the PDE system. SWE only hold a flux.

  solver ADER-DG SWE_ADERDG
    variables const    = h:1,p:2
    order const        = 3
    /* 27 points: 0.05, 9 points: 0.15 */
    maximum-mesh-size  = 0.05
    time-stepping      = globalfixed
    type const         = nonlinear
    terms const        = flux
    optimisation const = generic
    language const     = C
  end solver

¹ The explanation for this takes too long at this point, but believe me there is a reason.

Adding a plotter

Plotting is part of the solver, so we place the plotter into the solver block.
In this example we add a simple vtk output plotter which is sufficient for our needs.
Place these lines after language const and before end solver:

    plot vtk::Cartesian::vertices::ascii VtkWriter
      variables const = 3
      time            = 0.0
      repeat          = 1e-2
      output          = ./swe
    end plot

This results in a writer plotting all 3 variables every 0.02 seconds starting at time 0. The files will have a prefix swe and be placed in the project directory.

Run the Toolkit

To generate the code body apply the toolkit to the specification file we just created. Go back to the ExaHyPE root directory and hit :

./Toolkit/toolkit.sh Demonstrators/SWE/SWE_ADERDG.exahype

The Toolkit will end without output if it terminated successful, otherwise you'll receive a detailed error message.

Go back to Demonstrators/SWE/. This folder should contain several automatically generated files:

ls SWE_ADERDG
> AbstractSWE_ADERDG.cpp	AbstractSWE_ADERDG.h  KernelCalls.cpp  Makefile  README_generated.md  SWE_ADERDG.cpp  SWE_ADERDG.h  SWE_ADERDG_Variables.h  VtkWriter.cpp  VtkWriter.h

Implement the SWE

The PDE system is implemented in SWE_ADERDG.cpp. If you take look at the corresponding header (SWE_ADERDG.h) you see the functions you'll need to specify the PDE.

   void adjustPointSolution(const double* const x,const double t,const double dt,double* Q) final override;
   void eigenvalues(const double* const Q,const int d,double* lambda) final override;
       void boundaryValues(const double* const x,const double t,const double dt,const int faceIndex,const int direction,const double * const fluxIn,const double* const stateIn, const double* const gradStateIn,double *fluxOut,double* stateOut) final override;
   void flux(const double* const Q,double** F) final override;

Flux implementation and Variables class

We'll start with the implementation of the flux. To make things easier ExaHyPE comes with a nice way to directly work with symbolic variables. You can access all arrays that hold state or flux variables with the mathematical notation you defined in the specification.

We evolve the last two terms of the flux to matrix vector products as it simplifies the implementation.

$$\begin{pmatrix} \frac{1}{2}gh^2 + \frac{(hu)^2}{h} & \frac{hu \cdot hv}{h} \\\ \frac{hu \cdot hv}{h}& \frac{1}{2}gh^2 + \frac{(hv)^2}{h} \end{pmatrix} = \frac{1}{2}gh^2 \cdot \mathbf{I}_2 + \begin{pmatrix} hu \\hv \end{pmatrix} \cdot \begin{pmatrix} hu & hv \end{pmatrix} \cdot \frac{1}{h}$$

By using symbolic variables and upper definition of the flux we can implement the method:

void SWE::SWE_ADERDG::flux(const double* const Q,double** F) {                                                                                                                                                                  
  ReadOnlyVariables vars(Q);
  Fluxes f(F);

  tarch::la::Matrix<2,2,double> I;
  I = 1, 0,
      0, 1;

  double h_inv = 1.0 / vars.h();
  f.h(vars.p());
  f.p(outerDot(vars.p(),vars.p()) * h_inv + 0.5 * vars.h() * vars.h() * I);
}

Eigenvalues

The eigenvalues can directly be taken from upper mathematical formulation. They need to passed with respect to the dimension defined by the argument const int d. For d=0 in x and d=1 in y (in 3D we additionally get d=2 in z).

In our case d defines the considered velocity field.

void SWE::SWE_ADERDG::eigenvalues(const double* const Q,const int direction,double* lambda) {                                                                                                                                                       
  ReadOnlyVariables vars(Q);
  Variables eigs(lambda);

  double u = vars.p(direction)/vars.h();
  double c = std::sqrt(vars.h());

  eigs.h() = u - c;
  eigs.p(0)= u;
  eigs.p(1)= u + c;
}

Initial Values

void SWE::SWE_ADERDG::adjustPointSolution(const double* const x,const double t,const double dt,double* Q) {                                                                                                                                                          
  Variables vars(Q) ;
  if (tarch::la::equals(t,0.0)) {
    vars.h() = std::exp(-0.5*(x[0]*x[0]+x[1]*x[1])/0.01) * 0.01 + 4.0;
    vars.p(0) = 0.0;
    vars.p(1) = 0.0;
  }
}

Boundary Conditions

This method implements conditions that are applied along the boundary of the domain. The engine passes the fields stateIn and fluxIn which hold the respective state and flux of the simulated domain. In this short example we implement reflecting boundary conditions. For SWE this means we have to invert the velocity vector along the normal of the interface. As our mesh is Cartesian, the normal can be defined by the single index where it is one const int normalNonZero.

The state can then be inverted by:

  // Dimensions                        = 2                                                                                                                                                                  
  // Number of variables + parameters  = 3 + 0                                                                                                                                                              

  int normalZero = (direction+1) % 2;

  ReadOnlyVariables in(stateIn) ;
  Variables out(stateOut);

  out.h() = in.h();
  out.p(normalZero)   =  in.p(normalZero);
  out.p(direction)= -in.p(direction);
  //...

To get the new flux we apply the inverted velocity vector to the flux term (in this example for normalNonZero=0):

$$\vec{F_{out}}\cdot \vec{n}=F(stateOut)\cdot \vec{n}= \begin{pmatrix} -hu \\\ \frac{1}{2}gh^2 + (-hu)^2 \\\ (-hu)v \end{pmatrix}$$

which results in the implementation

  //...
  ReadOnlyVariables f_in (fluxIn) ;
  Variables f_out(fluxOut);

  f_out.h() = -f_in.h();
  f_out.p(normalNonZero) =  f_in.p(normalNonZero);
  f_out.p(normalZero)    = -f_in.p(normalZero);
}

Shared and Distributed Memory

Using shared and distributed memory to run is optional.
To use multiple cores add the shared or distributed memory simply add the following to the .exahype file:

  shared-memory
    identifier               = dummy
    configure                = {}
    cores                    = 1
    properties-file          = sharedmemory.properties
  end shared-memory  
  
  distributed-memory
    identifier       = static_load_balancing
    configure        = {hotspot,fair,ranks-per-node:1}
    buffer-size      = 64
    timeout          = 60
  end distributed-memory

Then in the terminal export the correct settings:

export SHAREDMEM=TBB
export DISTRIBUTEDMEM=MPI

Navigation Menu

Repository navigation

Home

The Shallow Water Equations

Setting up the specification file

`exahype-project`

`output-directory const`

`computational-domain`

`solver ADER-DG` ...

Adding a plotter

Run the Toolkit

Implement the SWE

Flux implementation and Variables class

Eigenvalues

Initial Values

Boundary Conditions

Shared and Distributed Memory

Clone this wiki locally

Home

The Shallow Water Equations

Setting up the specification file

exahype-project

output-directory const

computational-domain

solver ADER-DG ...

Adding a plotter

Run the Toolkit

Implement the SWE

Flux implementation and Variables class

Eigenvalues

Initial Values

Boundary Conditions

Shared and Distributed Memory

Clone this wiki locally

`exahype-project`

`output-directory const`

`computational-domain`

`solver ADER-DG` ...