sandbox/fcor_test.c

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
/*
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](http://en.wikipedia.org/wiki/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`. */

scalar hn[];

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

/**
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.
$$
\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();