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");
}
#endifh0 = 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'
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'
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)'
