# Bouncing Saint-Venant bump

This test case is identical to bump2D.c but using the generic solver for systems of conservation laws rather than the Saint-Venant solver.

``#include "conservation.h"``

The only conserved scalar is the water depth `h` and the only conserved vector is the flow rate `q`.

``````scalar h[];
vector q[];
scalar * scalars = {h};
vector * vectors = {q};``````

Using these notations, the Saint-Venant system of conservation laws (assuming a flat topography) can be written ${\partial }_{t}\left(\begin{array}{c}\hfill h\hfill \\ \hfill {q}_{x}\hfill \\ \hfill {q}_{y}\hfill \end{array}\right)+\nabla \cdot \left(\begin{array}{cc}\hfill {q}_{x}\hfill & {q}_{y}\\ \hfill \frac{{q}_{x}^{2}}{h}+\frac{G{h}^{2}}{2}\hfill & \frac{{q}_{x}{q}_{y}}{h}\\ \hfill \frac{{q}_{y}{q}_{x}}{h}\hfill & \frac{{q}_{y}^{2}}{h}+\frac{G{h}^{2}}{2}\end{array}\right)=0$ with $G$ the acceleration of gravity.

This system is entirely defined by the `flux()` function called by the generic solver for conservation laws. The parameter passed to the function is the array `s` which contains the state variables for each conserved field, in the order of their definition above (i.e. scalars then vectors). In the function below, we first recover each value (`h`, `qx` and `qy`) and then compute the corresponding fluxes (`f[0]`, `f[1]` and `f[2]`). The minimum and maximum eigenvalues for the Saint-Venant system are the characteristic speeds $u±\sqrt{\left(}Gh\right)$.

``````double G = 1.;

void flux (const double * s, double * f, double e[2])
{
double h = s[0], qx = s[1], u = qx/h, qy = s[2];
f[0] = qx;
f[1] = qx*u + G*h*h/2.;
f[2] = qy*u;
// min/max eigenvalues
double c = sqrt(G*h);
e[0] = u - c; // min
e[1] = u + c; // max
}``````

The solver is now fully defined and we proceed with initial conditions etc… as when using the standard Saint-Venant solver (see bump2D.c for details).

``````#define LEVEL 7

int main()
{
origin (-0.5, -0.5);
init_grid (1 << LEVEL);
run();
}

event init (i = 0)
{
foreach()
h[] = 0.1 + 1.*exp(-200.*(x*x + y*y));
}

event logfile (i++) {
stats s = statsf (h);
fprintf (stderr, "%g %d %g %g %.8f\n", t, i, s.min, s.max, s.sum);
}

event outputfile (t <= 2.5; t += 2.5/8) {
static int nf = 0;
printf ("file: eta-%d\n", nf);
output_field ({h}, stdout, N, linear = true);

scalar l[];
foreach()
l[] = level;
printf ("file: level-%d\n", nf++);
output_field ({l}, stdout, N);

/* check symmetry */
foreach() {
double h0 = h[];
point = locate (-x, -y);
//    printf ("%g %g %g %g %g\n", x, y, h0, h[], h0 - h[]);
assert (fabs(h0 - h[]) < 1e-12);
point = locate (-x, y);
assert (fabs(h0 - h[]) < 1e-12);
point = locate (x, -y);
assert (fabs(h0 - h[]) < 1e-12);
}
}

#if TREE