Wind-driven lake

This is a simple test case of a wind-driven lake where we can compare results with an analytical solution. For the bottom of the domain we impose no-slip condition (that is the default), for the top we impose a Neumann condition (see viscous friction between layers for details).

#include "grid/cartesian1D.h"
#include "saint-venant.h"

int main() {
L0 = 1.;
G = 100.;
N = 64;
nu = 1.;
dut = unity;

We vary the number of layers.

  for (nl = 4; nl <= 32; nl *= 2)
run();
}

We set the initial water level to 1 and we allocate a scalar field uc to check the convergence of the velocity in the first layer.

scalar uc[];

event init (i = 0) {
vector u0 = ul[0];
foreach() {
h[] = 1.;
uc[] = u0.x[];
}
}

We check for convergence.

event logfile (t += 0.1; i <= 100000) {
vector u0 = ul[0];
double du = change (u0.x, uc);
if (i > 0 && du < 1e-5)
return 1;
}

We compute the error between the numerical solution and the analytical solution.

#define uan(z)  ((z)/4.*(3.*(z) - 2.))

event error (t = end) {
int i = 0;
foreach() {
if (i++ == N/2) {
double sumh = zb[], emax = 0.;
int l = 0;
for (vector u in ul) {
double z = sumh + h[]*layer[l]/2.;
double e = fabs(u.x[] - uan (z));
if (e > emax)
emax = e;
sumh += h[]*layer[l++];
}
fprintf (stderr, "%d %g\n", nl, emax);
}
}
}

We save the horizontal velocity profile at the center of the domain and the two components of the velocity field for the case with 32 layers.

event output (t = end) {
char name[80];
sprintf (name, "uprof-%d", nl);
FILE * fp = fopen (name, "w");
int i = 0;
foreach() {
if (i++ == N/2) {
int l = 0;
double z = zb[] + h[]*layer[l]/2.;
for (vector u in ul)
fprintf (fp, "%g %g\n", z, u.x[]), z += h[]*layer[l];
}
if (nl == 32) {
double sumh = zb[];
int l = 0;
scalar w;
vector u;
for (w,u in wl,ul) {
double z = sumh + h[]*layer[l]/2.;
printf ("%g %g %g %g\n", x, z, u.x[], w[]);
sumh += layer[l++]*h[];
}
printf ("\n");
}
}
fclose (fp);
}

Results

set xr [0:1]
set xl 'z'
set yl 'u'
set key left top
plot [0:1]x/4.*(3.*x-2.) t 'analytical', \
'uprof-4' t '4 layers', \
'uprof-8' t '8 layers', \
'uprof-16' t '16 layers', \
'uprof-32' t '32 layers'
reset
set cbrange [1:2]
set logscale
set xlabel 'Number of layers'
set ylabel 'max|e|'
set xtics 4,2,32
set grid
fit a*x+b 'log' u (log($1)):(log($2)) via a,b
plot [3:36]'log' u 1:2 pt 7 t '', \
exp(b)*x**a t sprintf("%.2f/N^{%4.2f}", exp(b), -a)
reset
unset key
set xlabel 'x'
set ylabel 'z'
plot [0:1][0:1]'out' u 1:2:($3/5.):($4/5.) w vectors