src/grid/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;

foreach_dimension()
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 (scalar s 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 (scalar s 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;
}

foreach_dimension()
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 (scalar s 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 (scalar s 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 (scalar s 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};
  foreach_dimension() {
    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;
  foreach_dimension()
    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++;
  }
  foreach_dimension()
    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};
  foreach_dimension() {
    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()*(1 << dimension*depth());
  else
    i = pid()*((1 << dimension*(depth() + 1)) - 1)/((1 << dimension) - 1);
  foreach_cell() {
    if (!leaves || is_leaf(cell))
      index[] = i++;
    if (is_leaf(cell))
      continue;
  }
  boundary ({index});
  return pid() == 0 ? i*npe() - 1 : -1;
}