... | ... | @@ -39,323 +39,3 @@ export SHAREDMEM=None |
|
|
and in the file ```./ExaHyPE/Makefile``` comment out all lines containing a ```-lrt```
|
|
|
|
|
|
There currently is no support for Windows, please use a virtual machine like [VirtualBox](https://www.virtualbox.org/) |
|
|
|
|
|
|
|
|
## SWE Tutorial ##
|
|
|
|
|
|
# The Shallow Water Equations
|
|
|
The system of hyperbolic PDEs we regard are the Shallow Water Equations ([wiki](https://en.wikipedia.org/wiki/Shallow_water_equations))
|
|
|
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:
|
|
|
```math
|
|
|
\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
|
|
|
```math
|
|
|
\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:
|
|
|
```math
|
|
|
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. Start from the following template for an application that is solved by the ADER-DG method.
|
|
|
|
|
|
```
|
|
|
exahype-project SWE
|
|
|
|
|
|
peano-kernel-path const = ./Peano
|
|
|
exahype-path const = ./ExaHyPE
|
|
|
output-directory const = ./ApplicationExamples/SWE_Test
|
|
|
architecture const = noarch
|
|
|
log-file = mylogfile.log
|
|
|
plotter-subdirectory const = Writers
|
|
|
|
|
|
|
|
|
/*Configure the domain*/
|
|
|
computational-domain
|
|
|
dimension const = 2
|
|
|
width = ##width
|
|
|
offset = ##offset
|
|
|
end-time = ##final time
|
|
|
end computational-domain
|
|
|
|
|
|
/*Configure shared memory settings*/
|
|
|
/*shared-memory
|
|
|
identifier = dummy
|
|
|
configure = {background-tasks:8}
|
|
|
cores = 8
|
|
|
properties-file = sharedmemory.properties
|
|
|
end shared-memory*/
|
|
|
|
|
|
/*Optimiziation: DO N0T TOUCH*/
|
|
|
global-optimisation
|
|
|
fuse-algorithmic-steps = all
|
|
|
fuse-algorithmic-steps-rerun-factor = 0.99
|
|
|
fuse-algorithmic-steps-diffusion-factor = 0.99
|
|
|
spawn-predictor-as-background-thread = off
|
|
|
spawn-amr-background-threads = off
|
|
|
/* 0.0 und 0.8 sind schon mal zwei Faktoren */
|
|
|
disable-vertex-exchange-in-time-steps = on
|
|
|
time-step-batch-factor = 0.0
|
|
|
disable-metadata-exchange-in-batched-time-steps = off
|
|
|
double-compression = 0.0
|
|
|
spawn-double-compression-as-background-thread = off
|
|
|
end global-optimisation
|
|
|
|
|
|
/*Defines the hyperbolic PDE*/
|
|
|
solver ADER-DG ADERSolver
|
|
|
variables const = h:1,p:2
|
|
|
order const = ##polynomial order of the simulation
|
|
|
/* 27 points: 0.05, 9 points: 0.15 */
|
|
|
maximum-mesh-size = ##mesh size
|
|
|
time-stepping = globalfixed
|
|
|
type const = nonlinear
|
|
|
terms const = ##PDE terms, ie fluxes, sources etc
|
|
|
optimisation const = generic
|
|
|
language const = C
|
|
|
|
|
|
/*##Add plotters here */
|
|
|
/*plot vtk::Cartesian::vertices::ascii VtkWriter
|
|
|
variables const = ##Number of variables to be plotted
|
|
|
time = ##Start time of the plotter
|
|
|
repeat = ##Plotting interval
|
|
|
output = ##File suffix
|
|
|
end plot */
|
|
|
end solver
|
|
|
|
|
|
end exahype-project
|
|
|
```
|
|
|
|
|
|
Set up a new directory *SWE_Test* in ApplicationExamples and copy the template there and 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.
|
|
|
|
|
|
```javascript
|
|
|
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 *./ApplicationExamples/SWE_Test*
|
|
|
|
|
|
```javascript
|
|
|
output-directory const = ./ApplicationExamples/SWE_Test
|
|
|
```
|
|
|
|
|
|
## `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:
|
|
|
|
|
|
```javascript
|
|
|
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 lokking 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<sup>1</sup>. 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 $`\frac{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.
|
|
|
|
|
|
```javascript
|
|
|
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
|
|
|
```
|
|
|
<sup>1</sup> The explanation for this takes too long at this point, but believe me there is a reason.
|
|
|
|
|
|
## Adding a plotter
|
|
|
* In this example we add a simple vtk output plotter which is sufficient for our needs.
|
|
|
* Uncomment the plotter lines after `language const` and before `end solver` and fill in the missing information:
|
|
|
|
|
|
```javascript
|
|
|
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 :
|
|
|
|
|
|
```sh
|
|
|
./Toolkit/toolkit.sh ApplicationExamples/SWE_Test/SWE_ADERDG.exahype
|
|
|
```
|
|
|
The Toolit will end with `setup build environment ... ok` if it terminated successful, else you'll receive a detailed error message.
|
|
|
|
|
|
Go back to ApplicationExamples/SWE_Test/ you should see several automatically generated files:
|
|
|
|
|
|
```sh
|
|
|
ls SWE_Test
|
|
|
> 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.
|
|
|
|
|
|
```c++
|
|
|
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 normalNonZero,const double * const fluxIn,const double* const stateIn,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.
|
|
|
```math
|
|
|
\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:
|
|
|
|
|
|
```c++
|
|
|
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.
|
|
|
|
|
|
```c++
|
|
|
void SWE::SWE_ADERDG::eigenvalues(const double* const Q,const int d,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
|
|
|
|
|
|
```++c
|
|
|
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:
|
|
|
|
|
|
```c++
|
|
|
void SWE::SWE_ADERDG::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, double *fluxOut,double* stateOut) {
|
|
|
// Dimensions = 2
|
|
|
// Number of variables + parameters = 3 + 0
|
|
|
ReadOnlyVariables in(stateIn) ;
|
|
|
Variables out(stateOut);
|
|
|
|
|
|
out.h() = in.h();
|
|
|
out.p(0) = in.p(0);
|
|
|
out.p(1) = in.p(1);
|
|
|
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):
|
|
|
|
|
|
```math
|
|
|
\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
|
|
|
|
|
|
```c++
|
|
|
//...
|
|
|
ReadOnlyVariables f_in (fluxIn) ;
|
|
|
Variables f_out(fluxOut);
|
|
|
|
|
|
f_out.h() = f_in.h();
|
|
|
f_out.p(0) = f_in.p(0);
|
|
|
f_out.p(1) = f_in.p(1);
|
|
|
}
|
|
|
```
|
|
|
|
|
|
## Compile & Run the program
|
|
|
After all methods are filled you can compile the program by running
|
|
|
|
|
|
```bash
|
|
|
make -j8
|
|
|
```
|
|
|
|
|
|
After compilation is finished run
|
|
|
|
|
|
```bash
|
|
|
./ExaHyPE_SWE SWE_ADER.exahype
|
|
|
```
|
|
|
|
|
|
You can look at the result files with Paraview. |
|
|
\ No newline at end of file |