src/test/lake-tr.c
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)'
zb(x) = - 50.*exp(-(x/1000.)**2) - 0.001*x
plot [-5000:5000] zb(x) w l t ''
Topography (script)
We use either the explicit or the implicit Saint-Venant solver.
#include "grid/multigrid1D.h"
#if EXPLICIT
# include "saint-venant.h"
#elif ML
# include "layered/hydro.h"
# include "layered/nh.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.025;
int main() {
L0 = 10000.;
X0 = -L0/2.;
G = 9.81;
N = 1 << LEVEL;
periodic (right);
#if ML
We add some damping by off-centering the implicit scheme.
theta_H = 0.6;
#endif
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 || ML
foreach()
u.x[] = (u.x[] + G*slope*dt)/(1. + dt*f/8.*u.x[]/h[]);
#else
foreach()
q.x[] = (q.x[] + h[]*G*slope*dt)/(1. + dt*f/8.*q.x[]/sq(h[]));
#endif
}
Outputs
We check for steady-state.
scalar herr[], uerr[];
event init_herr(i = 0) {
foreach()
herr[] = uerr[] = 0;
}
event logfile (i++) {
double dh = change (h, herr),
#if EXPLICIT || ML
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 || ML
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);
}
#if 0
event gnuplot (i += 10)
{
static FILE * fp = popen ("gnuplot", "w");
if (i == 0)
fputs ("set term x11\nset yrange [-0.2:2]\n", fp);
fprintf (fp, "plot '-' u 1:2 w l\n");
foreach()
fprintf (fp, "%g %g\n", x, u.x[]/sqrt(G*h[]));
fprintf (fp, "e\n");
fflush (fp);
}
#endif
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-ml/prof-600' u 1:(zb($1)+$2) w l t \
't = 10 min (implicit ML)', \
'../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-ml/prof-7200' u 1:(zb($1)+$2) w l t \
't = 2 hours (implicit ML)', \
'../lake-tr-explicit/prof-7200' u 1:(zb($1)+$2) w l t \
't = 2 hours (explicit)', \
zb(x) t 'topography'
Evolution of the free surface (script)
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-ml/prof-600' u 1:5 w l t 't = 10 min (implicit ML)', \
'../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-ml/prof-7200' u 1:5 w l t 't = 2 hours (implicit ML)', \
'../lake-tr-explicit/prof-7200' u 1:5 w l t 't = 2 hours (explicit)'
Evolution of the Froude number (script)
We can check the performance for all schemes.
cat lake-tr/out
cat lake-tr-explicit/out
cat lake-tr-ml/out
On my computer, this gives for the implicit scheme:
# Multigrid, 4460 steps, 1.42008 CPU, 1.42 real, 3.22e+06 points.step/s, 14 var
for the multilayer implicit scheme:
# Multigrid, 4524 steps, 1.46437 CPU, 1.464 real, 3.16e+06 points.step/s, 14 var
and for the explicit scheme :
# Multigrid, 32519 steps, 8.59387 CPU, 8.685 real, 3.83e+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.