path: root/src/SetupPGH.c

 /*@@
   @file      SetupPGH.c
   @date      Fri Feb 21 10:18:17 1997
   @author    Paul Walker
   @desc 
   Initializes the pgh. You should really consult the
   pGH description (@seeheader pGH.h) in order to clearly
   understand what is going on.
   @enddesc 
   @version $Header$
 @@*/


#include <stdio.h>
#include <stdlib.h>
#include <math.h>

#include "cctk.h"
#include "pugh.h"
#include "cctk_Parameters.h"

static char *rcsid = "$Header$";

void BoundingBox(pGH *GH, int *nproc, int staggertype);
void Neighbours1D(pGH *GH, int *nproc, int *maxpoints, int *minpoints, double *avgpoints);
void Neighbours3D(pGH *GH, int *nproc, int *maxpoints, int *minpoints, double *avgpoints);
void DestroyPGF(pGH *GH, pGF **GFin);
void pGH_DumpInfo(pGH *GH);
int SetupDefaultTopology(int dim, int *nprocs);
int pGH_SetupnProcs(pGH *GH,int dim);

 /*@@
   @routine    PUGH_SetupPGH
   @date       Fri Feb 21 10:21:36 1997
   @author     Gabrielle Allen, Tom Goodale, Paul Walker
   @desc 
   Sets up the PGH distribution based on processor number and
   problem size. This includes setting an optimal domain
   decomposition for each dimension, and setting ghosts and overlaps.
   @enddesc 
   @hdate 4 June 1999
   @hauthor Thomas Radke
   @hdesc Check for unique precision of CCTK_INT and CCTK_REAL types
          in metacomputing environment
   @hdate 30 June 1999
@@*/

pGH *PUGH_SetupPGH(void *callerid, 
                   int dim, 
                   int *nsize, 
                   int *nghostzones,
                   int staggertype,
		   int periodic) 
{
  pGH *GH;
  int idim;

  void pGH_SetupBasics(pGH *, int, int *, int *, int);
  void pGH_SetupBounds(pGH *GH);
  void pGH_SetupOwnership(pGH *,int dim);
  void pGH_SetupStaggering(pGH *GH, int dim);

  pConnectivity *pugh_SetupConnectivity
    (
     int dim, 
     int total_procs, 
     int *nprocs,
     int periodic
    );

  pGExtras *pugh_SetupPGExtras
    (
     int dim, 
     int periodic,
     int stagger,
     int *sh,
     int *nghosts,
     int total_procs,
     int *nprocs,
     int this_proc
    );

  /* Allocate for myself */
  GH = (pGH*) malloc(sizeof(pGH));

  /* Set an identifier for my parent */
  GH->callerid = callerid ;

  GH->GFExtras     = (pGExtras **)malloc(dim*sizeof(pGExtras*));
  GH->Connectivity = (pConnectivity **)malloc(dim*sizeof(pConnectivity*));

  /* Set the total number of processors */
  pGH_SetupnProcs(GH,dim);

  /* Set up connectivity and extras for each dimension */
  for (idim=1;idim<=dim;idim++)
  {
    int *nprocs;

    nprocs = (int *) malloc(idim*sizeof(int));

    SetupDefaultTopology(idim,nprocs);

    GH->Connectivity[idim-1]=pugh_SetupConnectivity(idim, 
                                                    GH->nprocs, 
                                                    nprocs,
                                                    periodic);
      
    /* FIXME : Should not be on GH */
    GH->nproc = ((pConnectivity **)GH->Connectivity)[idim-1]->nprocs;

    GH->GFExtras[idim-1]=pugh_SetupPGExtras(idim, 
                                            periodic,
                                            staggertype,
                                            nsize,
                                            nghostzones,
                                            GH->nprocs,
                                            ((pConnectivity **)GH->Connectivity)[idim-1]->nprocs,
                                            GH->myproc);

      free(nprocs);

  }

  pGH_SetupBasics(GH,dim,nsize,nghostzones,staggertype);
  pGH_SetupBounds(GH);
  pGH_SetupOwnership(GH,dim);
  pGH_SetupStaggering(GH,dim);

  return GH;

}


int pGH_SetupnProcs(pGH *GH, int dim)
{

  int i,j,k,idir;
  int *nproc;

#ifdef MPI
  CCTK_REAL4 sum_sizes [2], compiled_sizes [2];
#endif 

  DECLARE_CCTK_PARAMETERS;

  GH->dim = dim;

#ifdef MPI
  /* Set up my communicator. This would allow us to run pugh on 
     a subset of processors at a later date if, for instance, we
     were using panda or what not.
   */
  CACTUS_MPI_ERROR(MPI_Comm_dup(MPI_COMM_WORLD, &(GH->PUGH_COMM_WORLD)));

  CACTUS_MPI_ERROR(MPI_Comm_size(GH->PUGH_COMM_WORLD, &(GH->nprocs)));
  CACTUS_MPI_ERROR(MPI_Comm_rank(GH->PUGH_COMM_WORLD, &(GH->myproc)));

  /* check that all executables uses the same integer and fp precision
     within a metacomputing environment
     NOTE: We cannot use CCTK_INT4s in MPI_Allreduce() since
           (although they have the same size) they might be correspond to
           different MPI datatypes.
           CCTK_REAL4 should always refer to MPI_FLOAT. */
  compiled_sizes [0] = sizeof (CCTK_INT);
  compiled_sizes [1] = sizeof (CCTK_REAL);
  CACTUS_MPI_ERROR (MPI_Allreduce (compiled_sizes, sum_sizes, 2, PUGH_MPI_REAL4,
                    MPI_SUM, GH->PUGH_COMM_WORLD));
  if (compiled_sizes [0] * GH->nprocs != sum_sizes [0])
  {
    CCTK_WARN (0, "Cannot run executables with different precision for "
                  "CCTK_INTs within a metacomputing environment !\n"
                  "Please configure with unique CCTK_INTEGER_PRECISION !");
  }

  if (compiled_sizes [1] * GH->nprocs != sum_sizes [1])
  {
    CCTK_WARN (0, "Cannot run executables with different precision for "
                  "CCTK_REALs within a metacomputing environment !\n"
                  "Please configure with unique CCTK_REAL_PRECISION !");
  }

  /* define the complex datatype as a concatanation of 2 PUGH_MPI_REALs */
  CACTUS_MPI_ERROR(MPI_Type_contiguous (2, PUGH_MPI_REAL, &GH->pugh_mpi_complex));
#else
  GH->nprocs = 1;
  GH->myproc = 0;
#endif

}

void pGH_SetupBasics(pGH *GH, 
                     int dim,
                     int *nsize,        
                     int *nghostzones,
                     int staggertype) 
{

  int i,j,k,idir;

  DECLARE_CCTK_PARAMETERS;

  /* Activate GH (conv.grids can be inactivated by nancheck 
     so the evolution will carry on) */
  GH->active=1;

  /* Set up arrays for ownership,overlap and ghosts */
  for (i=0;i<PUGH_NSTAGGER;i++)
    for (j=0;j<2;j++)
    {
      GH->ownership[i][j] = (int *)malloc(dim*sizeof(int));
      GH->ghosts[i][j] = (int **)malloc(2*dim*sizeof(int *));
      for (k=0;k<2*dim;k++)
      {
        GH->ghosts[i][j][k] = (int *)malloc(dim*sizeof(int));
      }
      GH->overlap[i][j] = (int **)malloc(2*dim*sizeof(int *));
      for (k=0;k<2*dim;k++)
      {
        GH->overlap[i][j][k] = (int *)malloc(dim*sizeof(int));
      }
    }

  GH->nsize = dim > 0 ? GH->GFExtras[dim-1]->nsize : NULL;
  GH->nghostzones = dim > 0 ? GH->GFExtras[dim-1]->nghostzones : NULL;

  /*
  if (Contains("commmodel","loc") != 0) {
    GH->commmodel = PUGH_ALLOCATEDBUFFERS;
  } else {
    GH->commmodel = PUGH_DERIVEDTYPES;
  }
  */
  GH->commmodel = PUGH_ALLOCATEDBUFFERS;

  /* get the 'verbose output during synchronization' flag */
  GH->sync_verbose = sync_verbose;

  /* initialize the timer variable */
#if 0
  CactusResetTimer (&GH->comm_time);
#endif
  /* set staggering flag */
  GH->stagger = staggertype;

  if (periodic) 
  {
    GH->periodic = 1;
  } 
  else 
  {
    GH->periodic = 0;
  }

  /* Variable list */
  GH->nvariables = 0;
  GH->variables = NULL;

  /* Allocate the bounds memory */
  GH->lb        = (int **)malloc(GH->nprocs*sizeof(int *));
  GH->ub        = (int **)malloc(GH->nprocs*sizeof(int *));
  GH->neighbors = (int **)malloc(GH->nprocs*sizeof(int *));

  GH->lb[0] = (int *)malloc(GH->nprocs*dim*sizeof(int));
  for (i=1;i<GH->nprocs;i++)
  {
    GH->lb[i] = GH->lb[0] + i*dim;
  }

  GH->ub[0] = (int *)malloc(GH->nprocs*dim*sizeof(int));
  for (i=1;i<GH->nprocs;i++)
  {
    GH->ub[i] = GH->ub[0] + i*dim;
  }

  GH->neighbors[0] = (int *)malloc(GH->nprocs*2*dim*sizeof(int));
  for (i=1;i<GH->nprocs;i++)
  {
    GH->neighbors[i] = GH->neighbors[0] + i*2*dim;
  }

  /* Remote neighbor points */
  GH->rnpoints = (int *)malloc(GH->nprocs*sizeof(int));

  /* rnsize[nproc][dir] */
  GH->rnsize = (int **)malloc(GH->nprocs*sizeof(int *));
  for (i=0;i<GH->nprocs;i++)
  {
    GH->rnsize[i] = (int *)malloc(dim*sizeof(int));
  }

  /* In general - we will noodle this later.*/
  GH->level = 0;
  GH->mglevel = 0;
  GH->convlevel = 0;
  GH->forceSync = 0;

  

}


void pGH_SetupBounds(pGH *GH) 
{
  
  /* Set up bounding box information */
    BoundingBox(GH,GH->nproc,GH->stagger);

}



void pGH_SetupOwnership(pGH *GH,int dim) 
{
  int idir,i,j,k,dir,staglcv;
  int *one;
#ifdef MPI
  int dirs[3];
#endif
  int maxpoints, minpoints;
  double avgpoints;
  int *istart;
  int *iend;
  int *nproc;

  istart = (int *)malloc(dim*sizeof(int));
  iend = (int *)malloc(dim*sizeof(int));

  maxpoints = 0; 
  minpoints = 1;
  for (i=0;i<dim;i++)
    minpoints = minpoints*GH->nsize[i]; 
  avgpoints = 0;

  nproc = (int *)malloc(dim*sizeof(int));
  for (i=0;i<dim;i++)
  {
    nproc[i] = GH->nproc[i];
  }

  /* Set up neighbour information */
  switch (GH->dim) 
  {
    case 1:
      Neighbours1D(GH, nproc, &maxpoints, &minpoints, &avgpoints);
      break;
    case 3:
      Neighbours3D(GH, nproc, &maxpoints, &minpoints, &avgpoints);
      break;
case 0:
break;
    default:
      CCTK_WARN(0,"PUGH can't cope with this dimension grid yet");
      return;
  }

  avgpoints = (1.0 * avgpoints) / (GH->nprocs);
  if (maxpoints > minpoints)
  {
    GH->maxskew = (maxpoints-minpoints)/avgpoints * 100.0;
  }
  else
  {
    GH->maxskew = 0;
  }

  /* Local information */
  GH->npoints = dim > 0 ?GH->GFExtras[dim-1]->npoints : 0;
  GH->lnsize  = dim > 0 ?GH->GFExtras[dim-1]->lnsize : 0;



  /* So now set up the ownership, ghosts, and overlap */
  
  /* Ownership is pretty easy. Remember ownership is indexed as
     [stagger][ijk][min/max]. See pGF_Reduction for a use of this,
     among others.
     Note: Ownership is same for staggered and non-staggered grids
   */
  for (idir=0;idir<dim;idir++)
  {
    GH->ownership[PUGH_VERTEXCTR][0][idir] = 
      (GH->lb[GH->myproc][idir] == 0 ? 0 : GH->nghostzones[idir]);
    GH->ownership[PUGH_VERTEXCTR][1][idir] = 
      (GH->ub[GH->myproc][idir] == GH->nsize[idir]-1 ? 
       GH->lnsize[idir] : GH->lnsize[idir] - GH->nghostzones[idir]);
  }

  /* correct for periodic identification : Tue Jun 17 08:40:15 CDT 1997 */
  if(GH->periodic)
  {
    for (idir=0;idir<dim;idir++)
    {
      GH->ownership[PUGH_VERTEXCTR][0][idir] 
	= GH->nghostzones[idir];
      GH->ownership[PUGH_VERTEXCTR][1][idir] 
	= GH->lnsize[idir] - GH->nghostzones[idir];
    }
  } 

  for (idir=0;idir<dim;idir++)
  {
    if (GH->nsize[idir] == 1) 
    {
      GH->ownership[PUGH_VERTEXCTR][0][idir] = 0;
      GH->ownership[PUGH_VERTEXCTR][1][idir] = 1;
    }
  }

  /* Great now set up the overlap region. This is the region we own,
     but which is covered by our neighbors ghost. So this is what we
     want to SEND
   */
  for (dir=0;dir<2*dim;dir++) 
  {
    int tmp;

    for (idir=0;idir<dim;idir++)
    {

      istart[idir] = 0; 
      iend[idir] = GH->lnsize[idir];

      /* We want to send at the dir - 1 */

      /* Correct to direction speicif sweeps ... */
      if (dir == 2*idir) 
      {
        /* Want to go from sw to sw + sw (eg, 1 to < 2, 2 to < 4) */
        istart[idir] = GH->nghostzones[idir];
        iend[idir]   = istart[idir] + GH->nghostzones[idir];
      }
      if (dir == 2*idir+1) 
      {
        /* Want to go from nx - 2 sw to nx - sw
           (eg nx-2 to < nx-1, nx-4 to < nx-2)
           Luckily iend is aready set to nx ... */
        tmp = iend[idir];
        istart[idir] = tmp - 2 * GH->nghostzones[idir];
        iend[idir]   = tmp - GH->nghostzones[idir];
      }

      GH->overlap[PUGH_VERTEXCTR][0][dir][idir] = istart[idir];
      GH->overlap[PUGH_VERTEXCTR][1][dir][idir] = iend[idir];
    }
  }

  /* And finally the ghost indices. This used to be in
     pGF_FinishRecv.c, but now it is here.
   */
  for (dir=0;dir<2*dim;dir++) 
  {

    int tmp;

    /* We want to send at the dir - 1 */
    for (idir=0;idir<dim;idir++)
    {
      istart[idir] = 0; 
      iend[idir] = GH->lnsize[idir];

      /* Correct to direction speicif sweeps ... */
      /* Remember this is SW less than the other way... */
      if (dir == 2*idir) 
      {
      /* Want to go from sw to sw + sw (eg, 1 to < 2, 2 to < 4) */
        istart[idir] = 0;
        iend[idir]   = GH->nghostzones[idir];
      }
      if (dir == 2*idir+1) 
      {
      /* Want to go from nx - 2 sw to nx - sw
         (eg nx-2 to < nx-1, nx-4 to < nx-2)
         Luckily iend is aready set to nx ... */
        tmp = iend[idir];
        istart[idir] = tmp - GH->nghostzones[idir];
        iend[idir]   = tmp;
      }

      GH->ghosts[PUGH_VERTEXCTR][0][dir][idir] = istart[idir];
      GH->ghosts[PUGH_VERTEXCTR][1][dir][idir] = iend[idir];

    }

  }

  /* if staggering, set up ownership, overlaps, ghosts */
  if(GH->stagger == PUGH_STAGGER) 
  {
    /* copy ownership from PUGH_VERTEXCTR */
    for(i=1;i<PUGH_NSTAGGER;i++)
    {
      for(j=0;j<dim;j++)
      {
        for(k=0;k<2;k++)
        {
          GH->ownership[i][k][j] = GH->ownership[PUGH_VERTEXCTR][k][j];
        }
      }
    }

    /* correct ownership for staggered grids */
    for(i=0;i<PUGH_NSTAGGER;i++)
    {
      for (idir=0;idir<dim;idir++)
      { 
        if(GH->ub[GH->myproc][idir] != GH->nsize[idir]-1)
        {
          GH->ownership[i][1][idir] --;
        }
        else
        {
          if(i == PUGH_FACECTR(idir))
          {
            GH->ownership[i][1][idir] --;
          }
        }
      }
    }

    /* set overlaps for staggered grids */
    for(staglcv=0;staglcv<PUGH_NSTAGGER;staglcv++)
    for(dir=0;dir<2*dim;dir++)
    {
      for (idir=0;idir<dim;idir++)
      {
        istart[idir] = 0;  iend[idir] = GH->lnsize[idir];
        if(staglcv == PUGH_FACECTR(idir)) iend[idir] --;
        if(dir == 2*idir)
        {
          istart[idir] = GH->nghostzones[idir];
          iend[idir] = 2 * GH->nghostzones[idir];
          if(staglcv != PUGH_FACECTR(idir))
          {
            istart[idir] ++;
            iend[idir] ++;
          }
        }
        if(dir == 2*idir+1)
        {
          istart[idir] = GH->lnsize[idir] - 2 * GH->nghostzones[idir] -1;
          iend[idir] = istart[idir] + GH->nghostzones[idir];
        }

        GH->overlap[staglcv][0][dir][idir] = istart[idir];
        GH->overlap[staglcv][1][dir][idir] = iend[idir];
      }
    }

       
    /* set ghosts for staggered grids */
    for(staglcv=0;staglcv<PUGH_NSTAGGER;staglcv++)
    {
      for(dir=0;dir<2*dim;dir++)
      {
        for (idir=0;idir<dim;idir++)
        {
          istart[idir] = 0;  
          iend[idir] = GH->lnsize[idir];
          if(staglcv == PUGH_FACECTR(idir)) iend[idir] --;
          if(dir == 2*idir)
          {
            istart[idir] = 0;
            iend[idir] = GH->nghostzones[idir];
          }
          if(dir == 2*idir+1)
          {
            istart[idir] = GH->lnsize[idir] - GH->nghostzones[idir];
            iend[idir] = GH->lnsize[idir];
            if(staglcv == PUGH_FACECTR(idir))
            {
              istart[idir] --;
              iend[idir] --;
            }
          }
          GH->ghosts[staglcv][0][dir][idir] = istart[idir];
          GH->ghosts[staglcv][1][dir][idir] = iend[idir];
        }
      }
    }
  }


#ifdef MPI
  /* Great, now we have made these, lets set up the derived data types */
  for (i=0;GH->commmodel==PUGH_DERIVEDTYPES && i<PUGH_NSTAGGER;i++) 
  {
    int *dtindices, *one, dir,maxs;
    int ii,jj,kk,ll;
    int maxghosts;

    /* OK so now lets set up the communications datatypes */
    dirs[0] =  GH->lnsize[1] * GH->lnsize[2];
    dirs[1] =  GH->lnsize[0] * GH->lnsize[2];
    dirs[2] =  GH->lnsize[0] * GH->lnsize[1];
    
    maxs = 0;
    maxghosts = 0;
    for (idir=0;idir<dim;idir++)
    {
      maxs = (dirs[idir] > maxs ? dirs[idir] : maxs);
      maxghosts = (GH->nghostzones[idir] > maxghosts ? GH->nghostzones[idir] : maxghosts);
    }
    
    maxs *= maxghosts;

    dtindices = (int *)malloc(maxs*sizeof(int));
    one       = (int *)malloc(maxs*sizeof(int));

    for (dir=0;dir<6;dir++) {
      /* Great so now calculate. Sends first */
      ll = 0;
      for (kk = GH->overlap[i][0][dir][2];
           kk < GH->overlap[i][1][dir][2];
           kk ++)
      {
        for (jj = GH->overlap[i][0][dir][1];
             jj < GH->overlap[i][1][dir][1];
             jj ++)
        {
          for (ii = GH->overlap[i][0][dir][0];
               ii < GH->overlap[i][1][dir][0];
               ii ++) 
          {
            dtindices[ll] = DATINDEX(GH,ii,jj,kk);
            one[ll] = 1;
            ll++;
          }
        }
      }

      CACTUS_MPI_ERROR (MPI_Type_indexed (ll, one, dtindices, PUGH_MPI_CHAR,
                        &(GH->send_char_dt [i][dir])));
      CACTUS_MPI_ERROR (MPI_Type_commit (&(GH->send_char_dt [i][dir]))); 
      CACTUS_MPI_ERROR (MPI_Type_indexed (ll, one, dtindices, PUGH_MPI_INT,
                        &(GH->send_int_dt [i][dir])));
      CACTUS_MPI_ERROR (MPI_Type_commit (&(GH->send_int_dt [i][dir]))); 
      CACTUS_MPI_ERROR (MPI_Type_indexed (ll, one, dtindices, PUGH_MPI_REAL,
                        &(GH->send_real_dt [i][dir])));
      CACTUS_MPI_ERROR (MPI_Type_commit (&(GH->send_real_dt [i][dir]))); 
      CACTUS_MPI_ERROR (MPI_Type_indexed (ll, one, dtindices, GH->pugh_mpi_complex,
                        &(GH->send_complex_dt [i][dir])));
      CACTUS_MPI_ERROR (MPI_Type_commit (&(GH->send_complex_dt [i][dir]))); 

      /* And recieves */
      ll = 0;
      for (kk = GH->ghosts[i][0][dir][2];
           kk < GH->ghosts[i][1][dir][2];
           kk ++)
      {
        for (jj = GH->ghosts[i][0][dir][1];
             jj < GH->ghosts[i][1][dir][1];
             jj ++)
        {
          for (ii = GH->ghosts[i][0][dir][0];
               ii < GH->ghosts[i][1][dir][0];
               ii ++) 
          {
            dtindices[ll] = DATINDEX(GH,ii,jj,kk);
            one[ll] = 1;
            ll++;
          }
        }
      }

      CACTUS_MPI_ERROR (MPI_Type_indexed (ll, one, dtindices, PUGH_MPI_CHAR,
                        &(GH->recv_char_dt [i][dir])));
      CACTUS_MPI_ERROR (MPI_Type_commit (&(GH->recv_char_dt [i][dir]))); 
      CACTUS_MPI_ERROR (MPI_Type_indexed (ll, one, dtindices, PUGH_MPI_INT,
                        &(GH->recv_int_dt [i][dir])));
      CACTUS_MPI_ERROR (MPI_Type_commit (&(GH->recv_int_dt [i][dir]))); 
      CACTUS_MPI_ERROR (MPI_Type_indexed (ll, one, dtindices, PUGH_MPI_REAL,
                        &(GH->recv_real_dt [i][dir])));
      CACTUS_MPI_ERROR (MPI_Type_commit (&(GH->recv_real_dt [i][dir]))); 
      CACTUS_MPI_ERROR (MPI_Type_indexed (ll, one, dtindices, GH->pugh_mpi_complex,
                        &(GH->recv_complex_dt [i][dir])));
      CACTUS_MPI_ERROR (MPI_Type_commit (&(GH->recv_complex_dt [i][dir]))); 
    }

    free(dtindices);
    free(one);
  }
#endif

  free(nproc);
  free(istart);
  free(iend);
}

void pGH_SetupStaggering(pGH *GH, int dim) 
{
  int s,d,k,i,dir;
  int *upperbnd;

  upperbnd = (int*) malloc(dim*sizeof(int));
  
  for (d=0;d<dim;d++)
  {
    upperbnd[d] = (GH->ub[GH->myproc][d] == GH->nsize[d]-1) ? 1 : 0;
  }
  /* copy stuff from PUGH_M */
  for(s=1;s<PUGH_NSTAGGER;s++)
  {
    for(d=0;d<dim;d++)
    {
      for(k=0;k<2;k++)
      {
        GH->ownership[s][k][d]    = GH->ownership[0][k][d];
        for (dir=0;dir<2*dim;dir++)
        {
          GH->overlap[s][k][dir][d] = GH->overlap[0][k][dir][d];
          GH->ghosts [s][k][dir][d] = GH->ghosts [0][k][dir][d]; 
        }
      }
    }
  }

  /* correct ownership for staggered grids for the C,P cases */
  for(s=1;s<PUGH_NSTAGGER;s++)
  {
    for (d=0;d<dim;d++)
    {   
      /* if we are at a upper physical bnd, decrease by one */
      if(upperbnd[d]) 
      {
        GH->ownership[s][1][d] --;
      }
    }
  }

  /* correct overlaps for staggered grids */
  for(s=1; s<PUGH_NSTAGGER; s++) 
  {
    for(dir=0;dir<2*dim;dir++)
    {
      for (d=0;d<dim;d++)
      {
        /* if we are at a upper physical bnd && at a upperface (1,3,5) 
           && we are in the dir direction decrease by one */
        if((upperbnd[d])&&((dir+1)%2==0)&&(dir==2*d+1))
        {
          GH->overlap[s][0][dir][d]--;
          GH->overlap[s][1][dir][d]--;
        }
      }
    }
  }
       
  /* correct ghosts for staggered grids */
  for(s=1; s<PUGH_NSTAGGER; s++)
  { 
    for(dir=0;dir<2*dim;dir++)
    {
      for (d=0;d<dim;d++)
      { 
        /* if we are at a upper physical bnd && at a upperface (1,3,5) 
           && we are in the dir direction decrease by one */
        if(upperbnd[d]&&((dir+1)%2==0)&&(dir==2*d+1))
        {
          GH->ghosts[s][0][dir][d]--;
          GH->ghosts[s][1][dir][d]--;
        }
      }
    }
  }
  free(upperbnd);
}

 /*@@
   @routine    DestroyPGH
   @date       Thu Aug 21 11:38:10 1997
   @author     Paul Walker
   @desc 
   Destroys a GH object.
   @enddesc 
@@*/
void DestroyPGH(pGH **GHin) 
{
  /* First remove me from the list. */
  pGH    *GH;
  cGroup pgroup;
  int i;
  int didit;
  int variable;
  int group;
  int this_var;

  GH = *GHin;

#ifdef MPI
  CACTUS_MPI_ERROR(MPI_Comm_free(&(GH->PUGH_COMM_WORLD)));
#endif

  didit = 0;

  /* Great. Now go about the work of destroying me. */

  variable = 0;

  for(group = 0; group < CCTK_NumGroups(); group++)
  {

#ifdef DEBUG_PUGH
    printf("Calling Destroying Group %s \n",CCTK_GroupName(group));
#endif 

    CCTK_GroupData(group,&pgroup);

    for (this_var = 0; this_var < pgroup.numvars; this_var++, variable++) 
    {
      switch(pgroup.grouptype)
      {
        case GROUP_GF:
          for(i = 0; i < pgroup.numtimelevels; i++)
          {
            DestroyPGF(GH, &(((pGF ***)GH->variables)[variable][i]));
          }
          break;
        case CCTK_ARRAY:
          CCTK_WARN(1,"Need to add code for Destroying an Pugh array");
          break;
        case GROUP_SCALAR:
          for(i = 0; i < pgroup.numtimelevels; i++)
          {
            free(GH->variables[variable][i]);
          }
          break;
      }
      free(GH->variables[variable]);
    }
  }

  free(GH->variables);

  free(GH->lb[0]);
  free(GH->lb);

  free(GH->ub[0]);
  free(GH->ub);

  free(GH->rnpoints);
  free(GH->rnsize);

  free(GH->neighbors[0]);
  free(GH->neighbors);

  free(GH);
  *GHin = NULL;

  fflush(stdout);
}

void pGH_DumpInfo(pGH *GH) 
{
  int i,j,k;
  printf("INFO: \n");
  printf("   myproc: %d/%d \n",GH->myproc,GH->nprocs);
  printf("   ownership: \n");
  for (i=0;i<PUGH_NSTAGGER;i++)
  {
    for (j=0;j<GH->dim;j++)
    {
      printf("     GH->ownership[%d][0/1][%d]: %d %d \n",
             i,j,GH->ownership[i][0][j],GH->ownership[i][1][j]);
    }
  }

  for (i=0;i<PUGH_NSTAGGER;i++)
  {
    for (j=0;j<GH->dim;j++)
    {
      for (k=0;k<2*GH->dim;k++)
      {
        printf("     GH->ghosts[%d][0/1][%d][%d]: %d %d \n",
               i,k,j,GH->ghosts[i][0][k][j],GH->ghosts[i][1][k][j]);
      }
    }
  }

  for (i=0;i<PUGH_NSTAGGER;i++)
  {
    for (j=0;j<GH->dim;j++)
    {
      for (k=0;k<2*GH->dim;k++)
      {
        printf("     GH->overlap[%d][0/1][%d][%d]: %d %d \n",
               i,k,j,GH->overlap[i][0][j][k],GH->overlap[i][1][k][j]);
      }
    }
  }
}


 /*@@
   @routine    SetupDefaultTopology
   @date       Sun 23 Jan
   @author     Gabrielle Allen
   @desc 
   Use PUGH parameters to set the processor decomposition 
   is required
   @enddesc 
   @hdate 
   @hauthor 
   @hdesc 
   @hdate 
@@*/

int SetupDefaultTopology(int dim, int *nprocs)
{

  DECLARE_CCTK_PARAMETERS

  int retval=PUGH_SUCCESS;

  switch(dim)
  {
    case 1:
      nprocs[0] = processor_topology_1d_x;
      break;
    case 2:
      nprocs[0] = processor_topology_2d_x;
      nprocs[1] = processor_topology_2d_y;
      break;
    case 3:
      nprocs[0] = processor_topology_3d_x;
      nprocs[1] = processor_topology_3d_y;
      nprocs[2] = processor_topology_3d_z;
      break;
    default:
      nprocs=NULL;
      retval=PUGH_ERROR;
  }

  return retval;
}