/** # 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. ~~~gnuplot Topography 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 '' ~~~ 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. ~~~gnuplot Evolution of the free surface 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' ~~~ The flow is supercritical at the entrance to the lake. ~~~gnuplot Evolution of the Froude number 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)' ~~~ We can check the performance for all schemes. ~~~bash cat lake-tr/out cat lake-tr-explicit/out cat lake-tr-ml/out ~~~ On my computer, this gives for the implicit scheme: ~~~bash # Multigrid, 4460 steps, 1.42008 CPU, 1.42 real, 3.22e+06 points.step/s, 14 var ~~~ for the multilayer implicit scheme: ~~~bash # Multigrid, 4524 steps, 1.46437 CPU, 1.464 real, 3.16e+06 points.step/s, 14 var ~~~ and for the explicit scheme : ~~~bash # 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. */