# 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).

We run the test case for three different solvers: multilayer Saint-Venant, layered hydrostatic, layered non-hydrostatic.

#include "grid/multigrid1D.h"
#if STVT
#include "saint-venant.h"
#else
#include "layered/hydro.h"
#if NH
#include "layered/nh.h"
#endif
#include "layered/remap.h"
#endif

There are five parameters L, h, g, \dot{u}, \nu. Only three are independent e.g. \displaystyle a = \frac{L}{h}, \displaystyle \text{Re} = \frac{\dot{u} h^2}{\nu} = s \frac{gh^3}{\nu^2}, \displaystyle s = \frac{\dot{u} \nu}{gh} = \frac{\dot{u}^2 h}{\text{Re} g} where a is the aspect ratio, Re the Reynolds number and s the slope of the free-surface. A characteristic time scale is \displaystyle t_{\nu} = \frac{h^2}{\nu} We choose a small Reynolds number and a small slope.

double Re = 10.;
double s = 1./1000.;

double du0;
scalar duv[];

int main()
{
L0 = 10.;
X0 = -L0/2.;
G = 9.81;
N = 64;
nu = sqrt(s*G/Re);
du0 = sqrt(s*Re*G);
dut = duv;

We vary the number of layers.

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

We set the initial water level to 1 and set the surface stress.

event init (i = 0) {
foreach() {
#if STVT
h[] = 1.;
#else
for (scalar h in hl)
h[] = 1./nl;
#endif
duv[] = du0*(1. - pow(2.*x/L0,10));
}
}

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

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

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

Uncomment this part if you want on-the-fly animation.

#if 0
event gnuplot (i += 20) {
static FILE * fp = popen ("gnuplot 2> /dev/null", "w");
if (i == 0)
fprintf (fp, "set term x11\n");
fprintf (fp,
"set title 'nl = %d, t = %.2f'\n"
"p [%g:%g][0:]'-' u 1:3:2 w filledcu lc 3 t '',"
" '' u 1:(-1):3 t '' w filledcu lc -1", nl, t,
X0, X0 + L0);
#if STVT
fprintf (fp, "\n");
foreach()
fprintf (fp, "%g %g %g\n", x, zb[] + h[], zb[]);
#else
int i = 4;
for (scalar h in hl)
fprintf (fp, ", '' u 1:%d w l lw 2 t ''", i++);
fprintf (fp, "\n");
foreach() {
double H = 0.;
for (scalar h in hl)
H += h[];
fprintf (fp, "%g %g %g", x, zb[] + H, zb[]);
double z = zb[];
for (scalar h in hl)
fprintf (fp, " %g", z += h[]);
fprintf (fp, "\n");
}
#endif
fprintf (fp, "e\n\n");
//  fprintf (fp, "pause 0.05\n");
fflush (fp);
}
#endif

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;
#if !NH && !STVT
scalar * wl = NULL;
for (scalar h in hl) {
scalar w = new scalar;
wl = list_append (wl, w);
}
vertical_velocity (wl);
foreach() {
double wm = 0.;
scalar h, w;
for (h,w in hl,wl) {
double w1 = w[];
w[] = (w1 + wm)/2.;
wm = w1;
}
}
boundary (wl);
#endif // !NH && !STVT
foreach() {
if (i++ == N/2) {
#if STVT
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++];
#else
double z = zb[];
scalar h;
vector u;
for (u,h in ul,hl)
fprintf (fp, "%g %g\n", z + h[]/2., u.x[]), z += h[];
#endif
}
if (nl == 32) {
double z = zb[];
#if STVT
int l = 0;
scalar w;
vector u;
for (w,u in wl,ul) {
printf ("%g %g %g %g\n", x, z + h[]*layer[l]/2., u.x[], w[]);
z += layer[l++]*h[];
}
#else
scalar w, h;
vector u;
for (h,w,u in hl,wl,ul)
printf ("%g %g %g %g\n", x, z + h[]/2., u.x[], w[]), z += h[];
#endif
printf ("\n");
}
}
fclose (fp);
#if !NH && !STVT
delete (wl), free (wl);
#endif
}

## Results

set xr [0:1]
set xl 'z'
set yl 'u'
set key left top
G = 9.81
s = 1./1000.
Re = 10.
du0 = sqrt(s*Re*G)
plot [0:1]du0*x/4.*(3.*x-2.) t 'analytical', \
'uprof-4' pt 5 t '4 layers', \
'uprof-8' pt 6 t '8 layers', \
'uprof-16' pt 9 t '16 layers', \
'uprof-32' pt 10 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) lt 1
reset
unset key
set xlabel 'x'
set ylabel 'z'
scale = 10.
plot [-5:5][0:1]'out' u 1:2:($3*scale):($4*scale) w vectors