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
| /* tangential interpolation on face vector fields (in parallel)
see also: interpu.c */
#include "refine_unbalanced.h"
scalar h[];
face vector u[];
double R0 = 0.1;
u.n[right] = dirichlet(exp(-(x*x + y*y)/(R0*R0)));
u.t[top] = dirichlet(exp(-(x*x + y*y)/(R0*R0)));
int main (int argc, char ** argv)
{
int maxlevel = argc > 1 ? atoi(argv[1]) : 6;
int minlevel = argc > 2 ? atoi(argv[2]) : 5;
origin (-1, -1);
init_grid (1);
foreach_cell() {
cell.pid = pid();
cell.flags |= active;
}
tree->dirty = true;
refine_unbalanced (level <= minlevel*(1. - sqrt(sq((x + 0.5) - 0.1) +
sq((y + 0.5)- 0.1))), NULL);
mpi_partitioning();
refine_unbalanced (level <= maxlevel*(1. - sqrt(sq((x + 0.5) - 0.1) +
sq((y + 0.5) - 0.1))), NULL);
foreach()
h[] = exp(-(x*x + y*y)/(R0*R0));
foreach_face(x) u.x[] = exp(-(x*x + y*y)/(R0*R0));
foreach_face(y) u.y[] = exp(-(x*x + y*y)/(R0*R0));
// output_cells (stdout);
double max = 0;
foreach_face(x, reduction(max:max))
for (int i = -1; i <= 1; i++) {
double xu = x, yu = y + i*Delta;
double e = exp(-(xu*xu+yu*yu)/(R0*R0)) - u.x[0,i];
if (fabs(e) > max)
max = fabs(e);
if (fabs(e) > 1e-6)
fprintf (stdout, "u %g %g %d %d %g %g\n",
xu, yu, level, cell.neighbors, u.x[0,i], e);
}
double maxv = 0;
foreach_face (y, reduction(max:maxv))
for (int i = -1; i <= 1; i++) {
double xv = x + i*Delta, yv = y;
double e = exp(-(xv*xv+yv*yv)/(R0*R0)) - u.y[i,0];
if (fabs(e) > maxv)
maxv = fabs(e);
if (fabs(e) > 1e-6)
fprintf (stdout, "v %g %g %d %d %g %g\n",
xv, yv, level, cell.neighbors, u.y[i,0], e);
}
mpi_all_reduce (max, MPI_DOUBLE, MPI_MAX);
mpi_all_reduce (maxv, MPI_DOUBLE, MPI_MAX);
fprintf (stderr, "maximum error on halos: %g %g\n", max, maxv);
}
|
