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)

    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)

    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)

    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.