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
| #include "saint-venant.h"
#include "terrain.h"
#define MAXLEVEL 9
#define MINLEVEL 5
// metres to degrees
const double mtd = 360./40075e3;
int main()
{
// 512^2 grid points
init_grid (1 << MAXLEVEL);
// the domain is 5 degrees squared
size (5.);
// centered on 174,-40.8 longitude,latitude
origin (174 - L0/2., -40.8 - L0/2.);
/* rescale G so that time is in minutes, horizontal length scales in
degrees and vertical length scales in metres */
G = 9.81*sq(mtd)*sq(60.);
run();
}
// M2 tidal frequency. The period is 12h25 (745 minutes).
const double u0 = 1., M2F = 2.*pi/745.;
// "radiation" boundary conditions on left,right,top,bottom
u.n[left] = - radiation ( u0*sin(M2F*t));
u.n[top] = radiation ( u0*sin(M2F*t));
u.n[right] = + radiation (- u0*sin(M2F*t));
u.n[bottom] = - radiation (- u0*sin(M2F*t));
event init (i = 0)
{
// use several databases for topography zb
terrain (zb,
"/home/popinet/terrain/niwa",
"/home/popinet/terrain/geo",
"/home/popinet/terrain/etopo2",
NULL);
// ensure that water level is conserved during refinement/coarsening
// the default is to conserve volume
conserve_elevation();
// sealevel at z = 0
foreach()
h[] = max(0., - zb[]);
}
// every timestep
event logfile (i++) {
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 dt\n");
fprintf (stderr, "%g %d %g %g %g %g %g %g\n", t, i, s.min, s.max, s.sum,
n.rms, n.max, dt);
double C_f = 1e-4;
foreach() {
// quadratic bottom friction, coefficient 1e-4 (dimensionless)
double a = h[] < dry ? HUGE : 1. + C_f*dt*norm(u)/(h[]*mtd);
foreach_dimension()
u.x[] /= a;
}
}
// snapshots every hour
event snapshots (t += 60; t <= 6000) {
printf ("file: t-%g\n", t);
output_field ({h, zb, eta, u}, stdout, n = N, linear = true);
}
// movies every 10 minutes
event movies (t += 10) {
scalar m[], etam[];
foreach() {
etam[] = eta[]*(h[] > dry);
m[] = etam[] - zb[];
}
output_ppm (etam, mask = m, min = -1, max = 1, n = 512, linear = true,
file = "eta.mp4");
scalar vort = etam;
vorticity (u, vort);
output_ppm (vort, mask = m, // min = -1e-2, max = 1e-2,
n = 512, linear = true, file = "vort.mp4");
scalar l = etam;
foreach()
l[] = level;
output_ppm (l, min = MINLEVEL, max = MAXLEVEL, n = 512, file = "level.mp4");
}
// tide gauges
Gauge gauges[] = {
{"well", 96.88, -12.13, "Wellington, New Zealand"},
{NULL}
};
event gauges1 (i++) output_gauges (gauges, {eta});
event adapt (i++) {
scalar eta[];
foreach()
eta[] = h[] > dry ? h[] + zb[] : 0;
double cmax = 1e-2;
astats s = adapt_wavelet ({eta}, (double[]){cmax}, MAXLEVEL, MINLEVEL);
fprintf (stderr, "# refined %d cells, coarsened %d cells\n", s.nf, s.nc);
}
|
