src/test/mpi-interpu.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
  
/* 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));
    boundary ({h});

    foreach_face(x) u.x[] = exp(-(x*x + y*y)/(R0*R0));
    foreach_face(y) u.y[] = exp(-(x*x + y*y)/(R0*R0));
    boundary ((scalar *){u});
    //  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 (ferr, "maximum error on halos: %g %g\n", max, maxv);  
  }