multigrid-mpi.h

Go to the documentation of this file.

/** @file multigrid-mpi.h
 */
#define MULTIGRID_MPI 1
 
#if dimension == 1
 
macro2 foreach_slice_x (int start, int end, int l) {
  {
    int ig = 0; NOT_UNUSED(ig);
    Point point = {0};
    point.level = l; SET_DIMENSIONS();
    for (point.i = start; point.i < end; point.i++)
      {...}
  }
}
 
#elif dimension == 2
 
macro2 foreach_slice_x (int start, int end, int l) {
  {
    int ig = 0, jg = 0; NOT_UNUSED(ig); NOT_UNUSED(jg);
    Point point = {0};
    point.level = l; SET_DIMENSIONS();
    for (point.i = start; point.i < end; point.i++)
      for (point.j = 0; point.j < point.n.y + 2*GHOSTS; point.j++)
  {...}
  }
}
 
macro2 foreach_slice_y (int start, int end, int l) {
  {
    int ig = 0, jg = 0; NOT_UNUSED(ig); NOT_UNUSED(jg);
    Point point = {0};
    point.level = l; SET_DIMENSIONS();
    for (point.i = 0; point.i < point.n.x + 2*GHOSTS; point.i++)
      for (point.j = start; point.j < end; point.j++)
  {...}
  }
}
 
#elif dimension == 3
 
macro2 foreach_slice_x (int start, int end, int l) {
  {
    int ig = 0, jg = 0, kg = 0; NOT_UNUSED(ig); NOT_UNUSED(jg); NOT_UNUSED(kg);
    Point point = {0};
    point.level = l; SET_DIMENSIONS();
    for (point.i = start; point.i < end; point.i++)
      for (point.j = 0; point.j < point.n.y + 2*GHOSTS; point.j++)
  for (point.k = 0; point.k < point.n.z + 2*GHOSTS; point.k++)
    {...}
  }
}
 
macro2 foreach_slice_y (int start, int end, int l) {
  {
    int ig = 0, jg = 0, kg = 0; NOT_UNUSED(ig); NOT_UNUSED(jg); NOT_UNUSED(kg);
    Point point = {0};
    point.level = l; SET_DIMENSIONS();
    for (point.i = 0; point.i < point.n.x + 2*GHOSTS; point.i++)
      for (point.j = start; point.j < end; point.j++)
  for (point.k = 0; point.k < point.n.z + 2*GHOSTS; point.k++)
    {...}
  }
}
 
macro2 foreach_slice_z (int start, int end, int l) {
  {
    int ig = 0, jg = 0, kg = 0; NOT_UNUSED(ig); NOT_UNUSED(jg); NOT_UNUSED(kg);
    Point point = {0};
    point.level = l; SET_DIMENSIONS();
    for (point.i = 0; point.i < point.n.x + 2*GHOSTS; point.i++)
      for (point.j = 0; point.j < point.n.y + 2*GHOSTS; point.j++)
  for (point.k = start; point.k < end; point.k++)
    {...}
  }
}
 
#endif // dimension == 3
 
typedef struct {
  Boundary b;
  MPI_Comm cartcomm;
} MpiBoundary;
 
for (int _d = 0; _d < dimension; _d++)
static void * snd_x (int i, int dst, int tag, int level, scalar * list,
         MPI_Request * req)
{
  if (dst == MPI_PROC_NULL)
    return NULL;
  size_t size = 0;
  for (int _s = 0; _s < 1; _s++) /* scalar in list */
    size += s.block;
  size *= pow((1 << level) + 2*GHOSTS, dimension - 1)*GHOSTS*sizeof(real);
  double * buf = (double *) malloc (size), * b = buf;
  foreach_slice_x (i, i + GHOSTS, level)
    for (int _s = 0; _s < 1; _s++) /* scalar in list */
      for (scalar sb = s; sb.i < s.i + s.block; sb.i++, b++)
  memcpy (b, &sb[], sizeof(real));
  MPI_Isend (buf, size, MPI_BYTE, dst, tag, MPI_COMM_WORLD, req);
  return buf;
}
 
for (int _d = 0; _d < dimension; _d++)
static void rcv_x (int i, int src, int tag, int level, scalar * list)
{
  if (src == MPI_PROC_NULL)
    return;
  size_t size = 0;
  for (int _s = 0; _s < 1; _s++) /* scalar in list */
    size += s.block;
  size *= pow((1 << level) + 2*GHOSTS, dimension - 1)*GHOSTS*sizeof(real);
  double * buf = (double *) malloc (size), * b = buf;
  MPI_Status s;
  MPI_Recv (buf, size, MPI_BYTE, src, tag, MPI_COMM_WORLD, &s);
  foreach_slice_x (i, i + GHOSTS, level)
    for (int _s = 0; _s < 1; _s++) /* scalar in list */
      for (scalar sb = s; sb.i < s.i + s.block; sb.i++, b++)
  memcpy (&sb[], b, sizeof(real));
  free (buf);
}
 
trace
static void mpi_boundary_level (const Boundary * b, scalar * list, int level)
{
  scalar * list1 = NULL;
  for (int _s = 0; _s < 1; _s++) /* scalar in list */
    if (!is_constant(s) && s.block > 0)
      list1 = list_add (list1, s);
  if (!list1)
    return;
 
  prof_start ("mpi_boundary_level");
  
  if (level < 0) level = depth();  
  MpiBoundary * mpi = (MpiBoundary *) b;
  struct { int x, y, z; } dir = {0,1,2};
  for (int _d = 0; _d < dimension; _d++) {
    int left, right;
    MPI_Cart_shift (mpi->cartcomm, dir.x, 1, &left, &right);  
    MPI_Request reqs[2];
    void * buf[2];
    int npl = (1 << level) + 2*GHOSTS, nr = 0;
    if ((buf[0] = snd_x (npl - 2*GHOSTS, right, 0, level, list1, &reqs[nr])))
      nr++;
    if ((buf[1] = snd_x (2, left,  1, level, list1, &reqs[nr])))
      nr++;
    rcv_x (0, left,  0, level, list1);
    rcv_x (npl - GHOSTS,   right, 1, level, list1);
    MPI_Status stats[nr];
    MPI_Waitall (nr, reqs, stats);
    free (buf[0]); free (buf[1]);
  }
 
  free (list1);
 
  prof_stop();
}
 
static void mpi_boundary_destroy (Boundary * b)
{
  MpiBoundary * m = (MpiBoundary *) b;
  MPI_Comm_free (&m->cartcomm);
  free (m);
}
 
static void mpi_dimensions_error (int n)
{
  fprintf (stderr,
     "%s:%d: error: the number of MPI processes must be equal to ",
     __FILE__, LINENO);
  if (n > 1)
    fprintf (stderr, "%dx", n);
  fprintf (stderr, "%d^i\n", 1 << dimension);
  exit (1);  
}
 
Boundary * mpi_boundary_new()
{
  MpiBoundary * m = qcalloc (1, MpiBoundary);
  int n = 1;
  for (int _d = 0; _d < dimension; _d++)
    n *= Dimensions.x;
  if (npe() % n)
    mpi_dimensions_error (n);
  int j = npe()/n, i = 0;
  while (j > 1) {
    if (j % (1 << dimension))
      mpi_dimensions_error (n);
    j /= 1 << dimension;
    i++;
  }
  for (int _d = 0; _d < dimension; _d++)
    Dimensions.x *= 1 << i;
  MPI_Dims_create (npe(), dimension, &Dimensions.x);
  MPI_Cart_create (MPI_COMM_WORLD, dimension,
       &Dimensions.x, &Period.x, 0, &m->cartcomm);
  MPI_Cart_coords (m->cartcomm, pid(), dimension, mpi_coords);
 
  // make sure other boundary conditions are not applied
  struct { int x, y, z; } dir = {0,1,2};
  for (int _d = 0; _d < dimension; _d++) {
    int l, r;
    MPI_Cart_shift (m->cartcomm, dir.x, 1, &l, &r);
    if (l != MPI_PROC_NULL)
      periodic_boundary (left);
    if (r != MPI_PROC_NULL)
      periodic_boundary (right);
  }
 
  // rescale the resolution
  Dimensions_scale = Dimensions.x;
  N /= Dimensions.x;
  int r = 0;
  while (N > 1)
    N /= 2, r++;
  grid->depth = grid->maxdepth = r;
  N = Dimensions.x*(1 << r);
  grid->n = 1 << dimension*depth();
  grid->tn = npe()*grid->n;
  
  // setup boundary methods and add to list of boundary conditions
  Boundary * b = (Boundary *) m;
  b->level = mpi_boundary_level;
  b->destroy = mpi_boundary_destroy;
  add_boundary (b);
 
  return b;
}
 
trace
double z_indexing (scalar index, bool leaves)
{
  long i;
  if (leaves)
    i = pid()*(1L << (dimension*depth()));
  else
    i = pid()*((1L << (dimension*(depth() + 1))) - 1)/((1L << dimension) - 1);
  for (int _i = 0; _i < _N; _i++) /* foreach_cell */ {
    if (!leaves || is_leaf(cell))
      index[] = i++;
    if (is_leaf(cell))
      continue;
  }
  boundary ({index});
  return pid() == 0 ? i*npe() - 1 : -1;
}