# sandbox/wind_pool.c

## Constant wind over constant depth pool

Example of a constant wind over a square constant depth pool

A constant wind of 1 m/s blows in the x direction over a square pool of 1000 m side length with constant depth 2 m.

We assume quadratic bottom friction with coefficient ${C}_{f}=2.5×{10}^{-3}$. Wind stress $\tau$ is defined by the formula $\frac{\tau }{\rho }={C}_{10}{U}_{10}^{2},$ where ${C}_{10}$ is taken from Wu (1982) JGRC ${C}_{10}=\left(0.8+0.065{U}_{10}\right)×{10}^{-3}.$ So for a ${U}_{10}=10m/s$ we have a wind stress of $\frac{\tau }{{\rho }_{air}}=0.145{m}^{2}/{s}^{2}$ and thus $\frac{\tau }{{\rho }_{water}}=1.7×{10}^{-4}{m}^{2}/{s}^{2}$.

``````#include "grid/cartesian1D.h"
#include "saint-venant_ss.h"

#define MAXLEVEL 8
#define MINLEVEL 4
#define ETAE 1e-8``````

Here we can set standard parameters in Basilisk

``````int main()
{
init_grid(1 << MINLEVEL);
#else // Cartesian
// 1024^2 grid points
init_grid(1 << MAXLEVEL);
#endif``````

Here we setup the domain geometry. For the moment Basilisk only supports square domains. This case uses metres east and north. We set the size of the box `L0` is 1000 m and the coordinates of the lower-left corner `(X0,Y0)` are (-500,-500). In this case we are assuming a square ‘pool’ of length 1000 m

``````  // the domain is
size (1000.);
origin(-L0/2.,-L0/2.);``````
`````` `G` is the acceleration of gravity required by the Saint-Venant
solver. This is the only dimensional parameter..``````
``````  // acceleration of gravity in m/s^2
G = 9.81;
run();
}``````

Here we define an auxilliary function which we will use several times in what follows. Again we have two `#if...#else` branches selecting whether the simulation is being run on an (adaptive) quadtree or a (static) Cartesian grid.

We want to adapt according to two criteria: an estimate of the error on the free surface position – to track the wave in time – and an estimate of the error on the maximum wave height `hmax` – to make sure that the final maximum wave height field is properly resolved.

We first define a temporary field (in the automatic variable `η`) which we set to $h+{z}_{b}$ but only for “wet” cells. If we used $h+{z}_{b}$ everywhere (i.e. the default $\eta$ provided by the Saint-Venant solver) we would also refine the dry topography, which is not useful.

``````int adapt() {
scalar η[];
foreach()
η[] = h[] > dry ? h[] + zb[] : 0;
boundary ({η});``````

We can now use wavelet adaptation on the list of scalars `{η,hmax}` with thresholds `{ETAE,HMAXE}`. The compiler is not clever enough yet and needs to be told explicitly that this is a list of `double`s, hence the `(double[])` (type casting).

The function then returns the number of cells refined.

``````  astats s = adapt_wavelet ({eta}, (double[]){ETAE},
MAXLEVEL, MINLEVEL);
fprintf (stderr, "# refined %d cells, coarsened %d cells\n", s.nf, s.nc);
return s.nf;
#else // Cartesian
return 0;
#endif
}``````

## Initiation

``````event initiate(i=0)
{``````

The initial still water surface is at $z=0$ so that the water depth $h$ is…

``````  foreach() {
zb[] = -2;
h[] = max(0., - zb[]);
ts.x[] = 0;
ts.y[] = 0;
}
boundary ({h,zb});
}``````

## Stopping condition

We want the simulation to stop when we are close to steady state. To do this we store the `h` field of the previous timestep in an auxilliary variable `hn`.

``````scalar hn[];

event init_hn (i = 0) {
foreach() {
hn[] = h[];
}
}``````

Every 10 timesteps we check whether $u$ has changed by more than 10-8. If it has not, the event returns 1 which stops the simulation. We also output running statistics to the standard error.

``````event logfile (i+= 10; i <= 10000000) {
double dh = change (h, hn);
if ( (i > 100 && dh < 1e-8) || i==10000000) {
foreach()
fprintf (stderr, "%g %g\n", x, h[]);
return 1; /* stop */
}
stats s = statsf (h);
norm n = normf (u.x);
if (i == 0)
fprintf (stderr, "t i h.min h.max h.sum u.x.rms u.x.max dh  dt\n");
fprintf (stderr, "%12.8g %d %12.8g %12.8g %12.8g %12.8g %12.8g %12.8g %12.8g\n", t, i,
s.min, s.max, s.sum, n.rms, n.max, dh, dt);
}``````

## Bottom Friction and Wind Forcing

We also use a simple implicit scheme to implement quadratic bottom friction i.e. $\frac{d\mathbf{\text{u}}}{dt}=-{C}_{f}\mid \mathbf{\text{u}}\mid \frac{\mathbf{\text{u}}}{h}$ with ${C}_{f}=2.5×{10}^{-3}$.

Also assume that we have a constant wind blowing in the x direction of

``````event source (i++) {
double ramp = t < 12000. ? t/12000. : 1.;
struct { double x, y; } ts = {1.7e-4,0};

foreach() {
ts.x[] = ramp*1.7e-4;
ts.y[] = 0;
double a_inv = h[] < dry ? 0. : h[]/(h[] + 2.5e-3*dt*norm(u));
//double a = h[] < dry ? HUGE : 1. + 2.5e-3*dt*norm(u)/h[];
foreach_dimension()
u.x[] = u.x[] * a_inv;
}
boundary ({h,u});
}``````

## Output

Every 5 minutes, the $h$, ${z}_{b}$ and `hmax` fields are interpolated bilinearly onto a `n x n` regular grid and written on standard output.

``````event snapshots (t += 300) {
printf ("# file: t-%g\n", t);
//  output_field ({h, zb}, stdout, n = 1 << MAXLEVEL, linear = true);
foreach()
printf ("%10.4g %10.4g %10.8g %10.8g\n", x, y, h[], zb[]);
printf ("\n");
}

event movies (t+=60) {
static FILE * fp = NULL;
if (!fp) fp = popen ("ppm2mpeg > eta.mpg", "w");
scalar m[], etam[];
foreach() {
etam[] = η[]*(h[] > dry);
m[] = etam[] - zb[];
}
boundary ({m, etam});
output_ppm (etam, fp, min = -0.005, max = 0.005 , n = 512, linear = true);
//  output_ppm (etam, fp, n = 512, linear = true);
#if 0
static FILE * fp1 = NULL;
if (!fp1) fp1 = popen ("ppm2mpeg > level.mpg", "w");
scalar l = etam;
foreach()
l[] = level;
output_ppm (l, fp1, min = MINLEVEL, max = MAXLEVEL, n = 512);
#endif
}``````

We apply our `adapt()` function at every timestep.
``event do_adapt (i++) adapt();``