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 EXPLICIT
# include "saint-venant.h"
#else
# include "saint-venant-implicit.h"
#endif

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

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

int main()
{
  origin (-5000);
  size (10000);
  G = 9.81;
  DT = 10.;
#if !EXPLICIT
  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 - τ*τ)/2.;
  return A*A*Bf*Bf*exp (-τ*t)/(8.*G*G*h0)*(- s*τ*sin (2.*s*t) + 
	    (τ*τ/4. - s*s)*cos (2.*s*t)) - Bf*Bf*exp(-τ*t)/(4.*G) -
    exp (-τ*t/2.)/G*(Bf*s*cos (s*t) + τ*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
  // linear friction (implicit scheme)
  foreach()
    u.x[] /= 1. + τ*dt;
  boundary ({u.x});
#else
  // linear friction (implicit scheme)
  foreach()
    q.x[] /= 1. + τ*dt;
  boundary ({q.x});
#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
    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
      sq += Δ*h[]*u.x[];
#else
      sq += Δ*q.x[];
#endif
      sh += Δ*h[];
    }
    printf ("s %g %g %f\n", t, sq/sh, sh);
  }
}
Free surface and topography at t=1500 for N=64 grid points.

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

Evolution of the axial velocity with time

Evolution of the axial velocity with time

Convergence of the error on the free surface position

Convergence of the error on the free surface position

See also