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 (stderr, "maximum error on halos: %g %g\n", max, maxv);  
    }