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
| #include "advection.h"
#include "vof.h"
scalar c[];
scalar * interfaces = {c}, * tracers = NULL;
int MAXLEVEL;
int main()
{
// coordinates of lower-left corner
origin (-0.5, -0.5);
for (MAXLEVEL = 5; MAXLEVEL <= 7; MAXLEVEL++) {
init_grid (1 << MAXLEVEL);
run ();
}
}
#define circle(x,y) (sq(0.1) - (sq(x-0.25) + sq(y)))
event init (i = 0)
{
fraction (c, circle(x,y));
}
#define end 0.785398
event velocity (i++) {
#if TREE
double cmax = 1e-2;
adapt_wavelet ({c}, &cmax, MAXLEVEL, list = {c});
#endif
double a = -8.;
trash ({u});
foreach_face(x) u.x[] = - a*y;
foreach_face(y) u.y[] = a*x;
}
event logfile (t = {0,end}) {
stats s = statsf (c);
fprintf (stderr, "# %f %.12f %f %g\n", t, s.sum, s.min, s.max);
}
event interface (t += end/10.) {
static FILE * fp = fopen ("interface", "w");
if (N == 64) {
output_facets (c, fp);
if (t == end)
output_cells (fp);
}
}
event field (t = end) {
scalar e[];
fraction (e, circle(x,y));
foreach()
e[] -= c[];
norm n = normf (e);
fprintf (stderr, "%d %g %g %g\n", N, n.avg, n.rms, n.max);
if (N == 64)
output_field ({e}, stdout, N, linear = false);
}
|
