src/test/galilean_invariance.c
Advection of a rippled interface
We test the ability of the multilayer solver to transport a corrugated interface with a constant velocity profile, i.e. we check that the solver preserves the property of Galilean invariance.
In the absence of body forces, the interface is initialised with a sinusoidal shape which is then swept away with a uniform velocity profile U (this moving interface is therefore not to be confused with a inertial-gravity propagating wave).
The resulting solid-body like motion is an exact solution of the multilayer set of equations with the non-hydrostatic corrections, with \phi = 0.
#include "grid/multigrid1D.h"
#include "layered/hydro.h"
#include "layered/nh.h"
#include "layered/remap.h"We non-dimensionalise the problem with the advection velocity U and the wavelength \lambda, and we set the amplitude of the wave to one tenth of the channel depth.
#define wavelength 1.
#define U 1.
#define wavenumber (2.*pi/wavelength)
#define depth 0.5
#define amplitude (depth/10.)
#define endoftime (1.*wavelength/U)
double maxwaveerror;Main function
The boundary conditions are set to periodic. The test is done in a zero-gravity environment. The default minmod slope limiter is turned off to avoid blunting off the wave crests.
int main()
{
size (1 [0]); // space is dimensionless
periodic (right);
G = 0.;
gradient = NULL;
breaking = 0.; // fixme: without this the computation diverges for N = 256
for (N = 64; N <= 256; N *= 2)
for (nl = 1; nl <= 16; nl *= 2)
run();
}Initialisation
The initial shape of the interface is a simple sine, and the velocity field is set to a constant velocity U = 1 and w = 0 everywhere.
event init (i = 0)
{
foreach() {
double H = depth + amplitude*cos(wavenumber*x);
foreach_layer() {
h[] = H/nl;
u.x[] = U;
}
}
maxwaveerror = 0.;
}Outputs
We keep track of the amplitude of the wave vs time, and for a typical case (N = 128 and \text{nl} = 4) we export deeper checks on the non-hydrostatic pressure \phi and velocity u.
event monitoring (t += endoftime/100; t = 0.; t <= endoftime)
{
stats s = statsf(eta);
double etaamp = fabs((s.max - s.min)/(2*amplitude) - 1);
if (N == 128 && nl == 4) {
char name[80];
sprintf (name, "output-N-%d-nl-%d", N, nl);
static FILE * fpout = fopen (name, "w");
double absdevU = 0., absPhi = 0.;
foreach()
foreach_layer() {
if (fabs(u.x[] - U) > absdevU)
absdevU = fabs(u.x[] - U);
if (fabs(phi[]) > absPhi)
absPhi = fabs(phi[]);
}
fprintf (fpout, "%g %g %g %g\n", t, absPhi, absdevU/U, etaamp);
assert (absPhi < 2e-13 && absdevU/U < 1e-14);
}
maxwaveerror = max(maxwaveerror, etaamp);
}The reference file of this test is just based on the overall maximum wave amplitude error.
Movie
This is how we export the headline animated gif.
event movie (i += 5; t <= (wavelength/U))
{
if (N == 128 && nl == 4) {
static FILE * fp = popen ("gnuplot 2> /dev/null", "w");
if (i == 0)
fprintf (fp,
"set term pngcairo enhanced font ',10' size 800,250\n"
"set size ratio -1\n"
"unset key\n");
fprintf (fp,
"set output 'plot%04d.png'\n"
"set title 't = %.2f'\n"
"p [%g:%g][0:]'-' u 1:3:2 w filledcu lt rgb \"#993498DB\", "
" '' u 1:2 w l lw 4 lt rgb \"#3498DB\" ",
i/5, t, X0, X0 + L0);
int i = 4;
foreach_layer()
fprintf (fp, ", '' u 1:%d w l dt 2 lt rgb \"gray50\"", i++);
fprintf (fp, "\n");
foreach (serial) {
double H = 0.;
foreach_layer()
H += h[];
fprintf (fp, "%g %g %g", x, zb[] + H, zb[]);
double z = zb[];
foreach_layer() {
fprintf (fp, " %g", z);
z += h[];
}
fprintf (fp, "\n");
}
fprintf (fp, "e\n\n");
fflush (fp);
}
}
event moviemaker (t = end)
{
if (N == 128 && nl == 4)
system ("mogrify -format gif plot*.png && "
"gifsicle --delay 10 --loop plot*.gif > wave.gif && "
"rm -f plot*.*");
}Results
For the typical case N = 128 and \text{nl} = 4 the time history of the Poisson error \|\phi\|_\infty and of any deviation of the velocity to the setpoint \|u-U\|_\infty – which would be responsible for phase errors – are represented and are observed to be of the order of the machine precision.
set term @SVG size 960,240 font ',10'
set size ratio 2./(1. + sqrt(5.))
set logscale y
set grid
set style line 1 \
linecolor rgb '#774F38' \
linetype 1 linewidth 2 \
pointtype 7 pointsize 1.
set style line 2 \
linecolor rgb '#E08E79' \
linetype 1 linewidth 2 \
pointtype 7 pointsize 1.
set style line 3 \
linecolor rgb '#F1D4AF' \
linetype 1 linewidth 2 \
pointtype 7 pointsize 1.
set style line 4 \
linecolor rgb '#ECE5CE' \
linetype 1 linewidth 2 \
pointtype 7 pointsize 1.
set style line 5 \
linecolor rgb '#C5E0DC' \
linetype 1 linewidth 2 \
pointtype 7 pointsize 1.
set ylabel 'Error (phi)'
set xlabel "t"
set key left
set multiplot layout 1, 2
plot [0:1][1e-16:1e-10]'./output-N-128-nl-4' u 1:2 t 'Error on phi' with lp ls 1
set ylabel 'Error (u)'
plot [0:1][1e-16:1e-10]'./output-N-128-nl-4' u 1:3 t 'Error on u' with lp ls 2
unset multiplot
For the same case (N = 128 and \text{nl} = 4) the relative error in wave amplitude (dissipation error) is monitored through time, and is seen to be of the order of 10^{-4}.
set term @SVG size 640,320 font ',10'
set ylabel 'Error (amplitude)'
plot './output-N-128-nl-4' u 1:4 t 'Relative error on amplitude' with lp ls 3
Convergence
We track the value of the relative error on wave amplitude for various number of gridpoints N and number of layers \text{nl}. For up to 16 stacked layers the results are undistinguishable, and the error is observed to decrease as N^{-2}, indicating a second-order precision.
set xlabel 'Number of grid points'
set ylabel 'Error'
set key bottom left
set logscale x 2
set label 1 "N^{-2}" at 180,5e1*180**-2 font ',12' textcolor rgb '#E08E79' offset 0,1
filter1 = '< awk ''{ if ($2 == 1) print $1, $2, $3}'' '
filter2 = '< awk ''{ if ($2 == 2) print $1, $2, $3}'' '
filter4 = '< awk ''{ if ($2 == 4) print $1, $2, $3}'' '
filter8 = '< awk ''{ if ($2 == 8) print $1, $2, $3}'' '
filter16 = '< awk ''{ if ($2 == 16) print $1, $2, $3}'' '
plot filter1.'log' u 1:3 t "nl = 1" with linespoints linestyle 1, \
filter2.'log' u 1:3 t "nl = 2" with linespoints linestyle 2, \
filter4.'log' u 1:3 t "nl = 4" with linespoints linestyle 3, \
filter8.'log' u 1:3 t "nl = 8" with linespoints linestyle 4, \
filter16.'log' u 1:3 t "nl = 16" with linespoints linestyle 5, \
[96:96<<2]5e1*x**-2 t '' with lines linestyle 2 lw 3
