sandbox/fcor_test.c
/*
Western Boundary Current intensification
This is using Stommel's western boundary current as a test case
Stommel, Henry (April 1948). "The Westward Intensification of Wind-Driven Ocean Currents". Transactions, American Geophysical Union 29 (2): 202–206. doi:10.1029/tr029i002p00202
See http://ido.at.fcen.uba.ar/index_archivos/Stommel_1948.pdf
A sinusoidal wind with wind stress in the x direction of $1.7 \times 10^{-4} \cos(\pi y / L0)$ m/s blows in the x direction over a square pool of L0 = 1000 km side length with constant depth 2 km. There is linear friction with coefficient 2.5e-3.
A beta-plane approximation is made with
f=-9.5e-5+1.7e-11*y
Which corresponds to the origin being at 41 degrees South.
$$ */
#include "grid/cartesian.h"
#include "storm_surge.h"
#define MAXLEVEL 8
#define MINLEVEL 4
#define ETAE 1e-8
double t1=8*604800.;
Here we can set standard parameters in Basilisk
int main()
{
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
and the origin to be the lower-left corner (X0,Y0)
. In this case we are assuming a square ‘pool’ of length 1000 km.
// the domain is
L0 = 1000000.;
G
is the acceleration of gravity required by the Saint-Venant solver. This is the only dimensional parameter..
G = 9.81;
#if QUADTREE
// 32^2 grid points to start with
init_grid( 1 << MINLEVEL );
#else // Cartesian
// 1024^2 grid points
init_grid( 1 << MAXLEVEL );
#endif
run();
}
Adaptation
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() {
#if QUADTREE
scalar eta[];
foreach()
eta[] = h[] > dry ? h[] + zb[] : 0;
boundary ({eta});
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
}
event initiate(i=0)
{
foreach() {
fcor[]=-9.5e-5+1.7e-11*y;
zb[] = -2000.;
h[] = max(0., - zb[]);
ts.x[] = -1.7e-4*cos(pi*y/L0);
ts.y[] = 0.;
}
boundary ({h,zb});
}
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
.
We output running statistics to the standard error.
event logfile (i+=10; t<t1) {
double dh = change (h, hn);
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);
}
We also use a simple implicit scheme to implement linear bottom friction i.e. \displaystyle \frac{d\mathbf{u}}{dt} = - C_f\frac{\mathbf{u}}{h} with C_f=2.5 \times 10^{-3}.
event source (i++) {
foreach() {
fcor[]=-9.5e-5+1.7e-11*y;
ts.x[] = - 1.7e-4*cos(pi*y/L0);
ts.y[] = 0.;
double a_inv = (h[] < dry ? 0. : h[]/(h[] + 2.5e-3*dt));
foreach_dimension()
u.x[] = u.x[] * a_inv;
}
boundary ({h,u,fcor,ts});
}
Every 12 hours, various fields are interpolated bilinearly onto a n x n
regular grid and written on standard output.
event snapshots (t += 86400.) {
scalar ux[], uy[], tsx[], tsy[];
int day = (int) t/86400.;
foreach(){
ux[]=u.x[];
uy[]=u.y[];
tsx[]=ts.x[];
tsy[]=ts.y[];
}
printf ("%% file: t-%d\n", day);
output_field ({h, eta, zb, ux, uy, tsx, tsy, fcor}, stdout, n = 1 << MAXLEVEL, linear = true);
}
event movies (t+=60) {
static FILE * fp = NULL;
if (!fp) fp = popen ("ppm2mpeg > eta.mpg", "w");
scalar m[], etam[];
foreach() {
etam[] = eta[]*(h[] > dry);
m[] = etam[] - zb[];
}
boundary ({m, etam});
output_ppm (etam, fp, min = -0.01, max = 0.01 , n = 512, linear = true);
#if QUADTREE
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
}
Adaptivity
We apply our adapt()
function at every timestep.
event do_adapt (i++) adapt();