# Lake flowing into itself

We consider a periodic domain, 10 km long, with a 50-metres-deep central lake. The domain is tilted with a slope of 1/1000 and a Darcy-Weissbach friction is imposed on the bottom.

set xlabel 'x (m)'
set ylabel 'z (m)'
set size ratio -1
zb(x) = - 50.*exp(-(x/1000.)**2) - 0.001*x
plot [-5000:5000] zb(x) w l t ''

We use either the explicit or the implicit Saint-Venant solver.

#include "grid/multigrid1D.h"
#if EXPLICIT
#include "saint-venant.h"
#else
#include "saint-venant-implicit.h"
#endif

## Domain setup and initial conditions

#define LEVEL 10

double slope = 0.001, eta0 = 0.5, f = 0.05;

int main() {
L0 = 10000.;
X0 = -L0/2.;
G = 9.81;
N = 1 << LEVEL;
periodic (right);
DT = 10;
run();
}

event init (i = 0) {
foreach() {
zb[] = - 50.*exp(-sq(x/1000.));
h[] = eta0 - zb[];
}
}

## Source terms

We add the slope explicitly and the Darcyâ€“Weissbach friction implicitly.

event source (i++) {
#if EXPLICIT
foreach()
u.x[] = (u.x[] + G*slope*dt)/(1. + dt*f/8.*u.x[]/sq(h[]));
boundary ((scalar *){u});
#else
foreach()
q.x[] = (q.x[] + h[]*G*slope*dt)/(1. + dt*f/8.*q.x[]/cube(h[]));
boundary ((scalar *){q});
#endif
}

## Outputs

scalar herr[], uerr[];
event init_herr(i = 0) {
foreach()
herr[] = uerr[] = 0;
}

event logfile (i++) {
double dh = change (h, herr),
#if EXPLICIT
du = change (u.x, uerr);
#else
du = change (q.x, uerr);
#endif
fprintf (stderr, "%g %g %g\n", t, dh, du);
#if 0
if (i > 100 && dh < 1e-7 && du < 1e-7) {
return 1;
}
#endif
}

We output the free-surface, Froude number etc..

event printprofile (t = {600, 3600, 7200}) {
char name[80];
sprintf (name, "prof-%g", t);
FILE * fp = fopen (name, "w");
foreach() {
#if EXPLICIT
fprintf (fp, "%g %g %g %g %g %g\n", x, h[], u.x[], zb[],
u.x[]/sqrt(G*h[]), u.x[]*h[]);
#else
fprintf (fp, "%g %g %g %g %g %g\n", x, h[], q.x[]/h[], zb[],
q.x[]/(h[]*sqrt(G*h[])), q.x[]);
#endif
}
fclose (fp);
}

## Results

After two hours a steady-state is reached. Both schemes give comparable results, even for the transient solution at 10 minutes.

set xlabel 'x (m)'
set ylabel 'z (m)'
set key top right
plot [][-6:]								\
'prof-600' u 1:(zb($1)+$2) w l t 't = 10 min (implicit)',		\
'../lake-tr-explicit/prof-600' u 1:(zb($1)+$2) w l t		\
't = 10 min (explicit)',						\
'prof-7200' u 1:(zb($1)+$2) w l t 't = 2 hours (implicit)',	\
'../lake-tr-explicit/prof-7200' u 1:(zb($1)+$2) w l t		\
't = 2 hours (explicit)',						\
zb(x) t 'topography'

The flow is supercritical at the entrance to the lake.

set ylabel 'Froude'
set key top right
plot 'prof-600' u 1:5 w l t 't = 10 min (implicit)', \
'../lake-tr-explicit/prof-600' u 1:5 w l t 't = 10 min (explicit)', \
'prof-7200' u 1:5 w l t 't = 2 hours (implicit)',			\
'../lake-tr-explicit/prof-7200' u 1:5 w l t 't = 2 hours (explicit)'

We can check the performance for both schemes.

cat lake-tr/out
cat lake-tr-explicit/out

On my computer, this gives for the implicit scheme:

# Multigrid, 4421 steps, 2.63557 CPU, 2.638 real, 1.72e+06 points.step/s, 14 var

and for the explicit scheme :

# Multigrid, 32519 steps, 14.4444 CPU, 14.73 real, 2.26e+06 points.step/s, 16 var

The gain in number of timesteps is a factor of ~7.3 and the gain in runtime a factor of ~5.6, which reflects the fact that the implicit scheme is slightly slower (per grid point) than the explicit scheme.