src/test/parabola.c

    Oscillations in a parabolic container

    We solve the Saint-Venant equations (explicitly or semi-implicitly) for the damped oscillations of a free-surface in a parabolic container. An analytic solution can be obtained due to the fact that the velocity is independent of space (for a linear interface intersecting a parabola). This was first observed by Thacker and later generalised to the (linearly) damped case by Sampson, Easton, Singh, 2006.

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

double e1 = 0., e2 = 0., emax = 0.;
int ne = 0;

double A = 3000.;
double h0 = 10.;
double tau = 1e-3;
double Bf = 5.;

int main()
{
  origin (-5000);
  size (10000);
  G = 9.81;
  DT = 10.;
#if ML
  dry = 1e-6;
#endif
#if !EXPLICIT && !ML
  dry = 1e-6;
  CFLa = 0.25;
#endif
  for (N = 32; N <= 512; N *= 2) {
    e1 = e2 = emax = 0.;
    ne = 0;
    run();
    fprintf (stderr, "%d %g %g %g\n", N, e1/ne/h0, sqrt(e2/ne)/h0, emax/h0);
  }
}

double Psi (double x, double t)
{
  // Analytical solution, see Sampson, Easton, Singh, 2006
  double p = sqrt (8.*G*h0)/A;
  double s = sqrt (p*p - tau*tau)/2.;
  return A*A*Bf*Bf*exp (-tau*t)/(8.*G*G*h0)*(- s*tau*sin (2.*s*t) + 
	    (tau*tau/4. - s*s)*cos (2.*s*t)) - Bf*Bf*exp(-tau*t)/(4.*G) -
    exp (-tau*t/2.)/G*(Bf*s*cos (s*t) + tau*Bf/2.*sin (s*t))*x;
}

event init (i = 0)
{
  foreach() {
    zb[] = h0*sq(x/A);
    h[] = max(h0 + Psi(x,0) - zb[], 0.);
  }
}

event friction (i++) {
#if EXPLICIT || ML
  // linear friction (implicit scheme)
  foreach()
    u.x[] /= 1. + tau*dt;
#else
  // linear friction (implicit scheme)
  foreach()
    q.x[] /= 1. + tau*dt;
#endif
}

scalar e[];

event error (i++) {
  foreach()
    e[] = h[] - max(h0 + Psi(x,t) - zb[], 0.);
  norm n = normf (e);
  e1 += n.avg;
  e2 += n.rms*n.rms;
  ne++;
  if (n.max > emax)
    emax = n.max;
  printf ("e %g %g %g %g %g\n", t, n.avg, n.rms, n.max, dt);
}

event field (t = 1500) {
  if (N == 64) {
#if EXPLICIT || ML
    foreach()
      printf ("p %g %g %g %g %g\n", x, h[], u.x[], zb[], e[]);
#else
    foreach()
      printf ("p %g %g %g %g %g\n", x, h[], q.x[], zb[], e[]);
#endif
    printf ("p\n");
  }
}

event umean (t += 50; t <= 6000) {
  if (N == 128) {
    double sq = 0., sh = 0.;
    foreach() {
#if EXPLICIT || ML
      sq += Delta*h[]*u.x[];
#else
      sq += Delta*q.x[];
#endif
      sh += Delta*h[];
    }
    printf ("s %g %g %f\n", t, sq/sh, sh);
  }
}

#if 0
event gnuplot (i += 10) {
  static FILE * fp = popen ("gnuplot 2> /dev/null", "w");
  fprintf (fp,
	   "set title 't = %.2f'\n"
	   "p [-5000:5000]'-' u 1:3:2 w filledcu lc 3 t '',"
	   " '' u 1:(-1):3 t '' w filledcu lc -1\n", t);
  foreach()
    fprintf (fp, "%g %g %g\n", x, zb[] + h[], zb[]);
  fprintf (fp, "e\n\n");
  //  fprintf (fp, "pause 0.02\n");
}
#endif
    h0 = 10.
a = 3000.
tau = 1e-3
B = 5.
G = 9.81
p = sqrt (8.*G*h0)/a
s = sqrt (p*p - tau*tau)/2.
u0(t) = B*exp (-tau*t/2.)*sin (s*t)

set xlabel 'x (m)'
set ylabel 'z (m)'
t = 1500
psi(x) = a*a*B*B*exp (-tau*t)/(8.*G*G*h0)*(- s*tau*sin (2.*s*t) + \
      (tau*tau/4. - s*s)*cos (2.*s*t)) - B*B*exp(-tau*t)/(4.*G) - \
      exp (-tau*t/2.)/G*(B*s*cos (s*t) + tau*B/2.*sin (s*t))*x + h0
bed(x) = h0*(x/a)**2
set key top center
plot [-5000:5000] \
      '< grep ^p out' u 2:5:($5+$3) w filledcu lc 3 t 'Implicit', \
      psi(x) > bed(x) ? psi(x) : bed(x) lc 2 t 'Analytical', \
      bed(x) lw 3 lc 1 lt 1 t 'Bed profile'
    Free surface and topography at t=1500 for N=64 grid points. (script)

    Free surface and topography at t=1500 for N=64 grid points. (script)

    reset
set key top right
set ylabel 'u0'
set xlabel 'Time'
plot u0(x) t 'analytical', \
     '< grep ^s out' u 2:3 every 2 w p t 'Implicit', \
     '< grep ^s ../parabola-explicit/out' u 2:3 every 2 w p t 'Explicit', \
     '< grep ^s ../parabola-ml/out' u 2:3 every 2 w p t 'Multilayer'
    Evolution of the axial velocity with time (script)

    Evolution of the axial velocity with time (script)

    reset
set xlabel 'Resolution'
set ylabel 'Relative error norms'
set key bottom left
set logscale
set cbrange [1:2]
set xtics 32,2,512
set grid
ftitle(a,b) = sprintf("order %4.2f", -b)
f1(x)=a1+b1*x
fit f1(x) 'log' u (log($1)):(log($2)) via a1,b1
f2(x)=a2+b2*x
fit f2(x) 'log' u (log($1)):(log($3)) via a2,b2
fm(x)=am+bm*x
fit fm(x) 'log' u (log($1)):(log($4)) via am,bm
plot exp (f1(log(x))) t ftitle(a1,b1), \
     exp (f2(log(x))) t ftitle(a2,b2), \
     exp (fm(log(x))) t ftitle(am,bm),  \
     'log' u 1:2 t '|h|_1', \
     'log' u 1:3 t '|h|_2', \
     'log' u 1:4 t '|h|_{max}' lc 0, \
     '../parabola-explicit/log' u 1:2 t '|h|_1 (explicit)', \
     '../parabola-explicit/log' u 1:3 t '|h|_2 (explicit)', \
     '../parabola-explicit/log' u 1:4 t '|h|_{max} (explicit)', \
     '../parabola-ml/log' u 1:2 t '|h|_1 (multilayer)', \
     '../parabola-ml/log' u 1:3 t '|h|_2 (multilayer)', \
     '../parabola-ml/log' u 1:4 t '|h|_{max} (multilayer)'
    Convergence of the error on the free surface position (script)

    Convergence of the error on the free surface position (script)

    See also