# Transcritical flow over a bump with multiple layers

We want to reproduce the transcritical test case of Audusse et al, 2011, section 5.6.2.

``````#include "grid/cartesian1D.h"
#include "saint-venant.h"``````

We need a field to store the variable bottom friction.

``````scalar λ[];

int main() {
X0 = 0.;
L0 = 21.;
G = 9.81;
N = 256;``````

The viscosity is set to $\nu =0.01{m}^{2}/s$ and the bottom friction is variable.

``````  ν = 0.01;
lambda_b = λ;``````

We vary the number of layers.

``````  nl = 2;  run();
nl = 5;  run();
nl = 15; run();
}``````

We impose the outlet water level.

``````h[right]   = dirichlet(0.6);
η[right] = dirichlet(0.6);``````

## Initialisation

We initialise the topography, the initial water depth h and we create a field hc to check convergence on h.

``````scalar hc[];

event init (i = 0) {
foreach() {
zb[] = max(0., 0.2*(1. - 1./sq(5.75/2.)*sq(x - 10.)));
hc[] = h[]  = 0.6 - zb[];
}``````

We call the friction event (below) to initialize the bottom friction.

``  event ("friction");``

## Boundary conditions on velocity

We impose a constant inflow of 1 m^2/s at the inlet and a Neumann condition at the outlet.

``````  for (vector u in ul) {
u.n[left] = dirichlet(h[left] ? 1./h[left] : 0.);
u.n[right] = neumann(0.);
}
}``````

## Bottom friction

We use the Strickler relation: $k\left(h,\mathbf{\text{U}}\right)=\frac{g}{{S}^{2}{h}^{1/3}}\mid \mathbf{\text{U}}\mid$ with $S=25{m}^{1/3}/s$ the Strickler coefficient, $h$ the water depth and $\mathbf{\text{U}}$ the depth-averaged velocity. Note that we have to use a lower Strickler coefficient (i.e. larger friction) to get results comparable to those of Audusse et al, 2011.

``````event friction (i++) {
foreach() {
double U = 0.;
int l = 0;
for (vector u in ul)
U += u.x[]*layer[l++];
double S = 25., k = G/(sq(S)*pow(h[],1./3.))*fabs(U);
λ[] = k > 0. ? ν/k : 0.;
}
boundary ({λ});
}``````

We check for convergence.

``````event logfile (t += 0.1; i <= 100000) {
double dh = change (h, hc);
if (i > 0 && dh < 1e-4)
return 1;
}``````

Uncomment this part if you want on-the-fly animation.

``````#if 0
event output (i++) {
static FILE * fp = popen ("gnuplot", "w");
fprintf (fp,
"set title 'nl=%d, t=%f'\n"
"set xl 'x'\nset yl 'h'\n"
"plot [0:21][] '-' u 1:2 w l t 'eta', '-' u 1:3 w l t 'zb'\n",
nl, t);
foreach()
fprintf (fp, "%g %g %g\n", x, η[], zb[]);
fprintf (fp, "e\n");
fflush (fp);
}
#endif``````

## Outputs

At the end of the simulation we save the profiles.

``````event output (t = end) {
char name[80];
sprintf (name, "end-%d", nl);
FILE * fp = nl == 15 ? stderr : fopen (name, "w");
foreach() {
fprintf (fp, "%g %g %g\n", x, η[], zb[]);
if (nl == 15) {
double z = zb[];
int l = 0;
printf ("%g %g %g\n", x, z, u.x[]);
for (vector u in ul) {
z += layer[l++]*h[];
printf ("%g %g %g\n", x, z, u.x[]);
}
printf ("\n");
}
}
}``````