path: root/src/SetupPGF.c

 /*@@
   @file      SetupPGF.c
   @date      Fri Feb 21 11:04:20 1997
   @author    
   @desc 
   Sets up (eg allocates) the pugh Grid Function. This is
   crucial to understaiding the grid function struct.
   @enddesc 
   @version $Id$
 @@*/

#include <stdio.h>
#include <stdlib.h>
#include <assert.h>
#include <string.h>

#include "cctk.h"
#include "declare_parameters.h"
#include "cGH.h"
#include "Groups.h"

#include "pugh.h"

static char *rcsid = "$Id$";

void SetGFComm(pGH *GH, pGF *res, int docomm);

 /*@@
   @routine    SetupPGF
   @date       Fri Feb 21 11:05:04 1997
   @author     
   @desc 
   Allocates the Grid Func structure (I wanted to say object ...)
   There are several steps to this
   <ul>
     <li> Set up the name
     <li> Set up the docomm flag and IO flags
     <li> Allocate output file pointer arrays
     <li> Allocate the data buffer
     <li> Allocate the send buffer
     <li> Add myself to the Grid Hierarchy
   </ul>
   But this is pretty straightforward code nonetheless.
   <p>
   Note that the input variable "storage" determines whether to
   allocate the 3-D data or not. You can toggle this later
   with @seeroutine DisableGFDataStorage and @seeroutine
   EnableGFDataStorage
   @enddesc 
   @calls DisableGFDataStorage, EnableGFDataStorage
@@*/


pGF *SetupPGF(pGH *GH, const char *name, int dim, int varsize, int vtype) 
{
  DECLARE_PARAMETERS

  pGF *res;			/* The result */
  pGF **rlist;			/* A temporary to add to the GH */
  int i, dir, dirp1, dirp2, sz;	/* Counter thingies */

  /* Fudge for CCTK. */
  int stagger = PUGH_VERTEXCTR;
  int storage = PUGH_NOSTORAGE;
  int docomm = PUGH_NOCOMM;

  /* Space for the struct */
  res = (pGF *)malloc(sizeof(pGF));

  /* Set my parent GH */
  res->parentGH = GH;

  /* Set up the name */
  res->name = (char *)malloc((1+strlen(name))*sizeof(char));
  strcpy(res->name,name);

  res->stagger = stagger;
  if( (GH->stagger == PUGH_NO_STAGGER) &&
      (res->stagger != PUGH_VERTEXCTR) )
  {  printf ("FATAL ERROR! Cannot have staggered grids inside a GH with \n");
     printf ("             designation PUGH_NO_STAGGER. \n");
     STOP;
  }

  res->varsize = varsize;
  res->vtype = vtype;

#if 0
  
   /* Finally add this to the list of grid funcs in the gh */
  GH->ngridFuncs ++;
  rlist = (pGF **)malloc(GH->ngridFuncs * sizeof(pGF *));
  if (GH->ngridFuncs > 1) {
    for (i=0;i<GH->ngridFuncs-1;i++)
      rlist[i] = GH->gridFuncs[i];
    free(GH->gridFuncs);
  }
  res->gfno = GH->ngridFuncs - 1;
  rlist[GH->ngridFuncs-1] = res;
  GH->gridFuncs = rlist;

#endif

  /* FIXME: special pad in enable GF storage. */
  res->gfno = GH->nvariables-1;


 SetGFComm(GH, res, docomm);

  /* Set up the comm layer */
  res->padddata = NULL;
  res->data = NULL;
  res->send_buffer = NULL;
  res->recv_buffer = NULL;

  if (storage == PUGH_NOSTORAGE) {
    res->storage = PUGH_STORAGE; /* Fake it out... */
    DisableGFDataStorage(GH, res);
  } else if (storage == PUGH_STORAGE) {
    res->storage = PUGH_NOSTORAGE; /* Fake it out... */
    EnableGFDataStorage(GH, res);
    printf("Enabling storage in pugh\n");
  } else {
    printf ("WARNING: Storage keyword set incorrectly for %s [%d]\n",res->name,storage);
    printf ("         Perhaps misspelled in GridFuncList? Defaulting to storage active.\n");
    EnableGFDataStorage(GH, res);
  }

  /* FIXME I don't understand if the above needs to be there */
  /* It can be useful for debugging to have memory on for all GFs */
  if (enable_all_storage) EnableGFDataStorage(GH,res);

#ifdef MPI
  /* Null my send and recieve requests */
  for (i=0;i<6;i++) {
    res->sreq[i] = MPI_REQUEST_NULL;
    res->rreq[i] = MPI_REQUEST_NULL;
  }
#endif

  /* And return myself... */
  return res;
}

 /*@@
   @routine    DestroyPGF
   @date       Thu Aug 21 11:44:22 1997
   @author     Paul Walker
   @desc 
   Destroys a GF object.
   @enddesc 
@@*/

void DestroyPGF(pGH *GH, pGF **GFin) {
  pGF *GF;
  GF = *GFin;

  if (GF->storage)
    DisableGFDataStorage(GH, GF);
  
  free(GF->name);
  free(GF->padddata);

  GF->data = NULL;

  free(GF);
  *GFin = NULL;

}


 /*@@
   @routine    EnableGFDataStorage
   @date       Thu Apr  3 12:44:38 1997
   @author     Paul Walker
   @desc 
   This routine toggles the data storage to the "on"
   position. That means that data points to a 3D
   array, rather than a single Double, and the
   storage flag is set to one. This is used quite
   a lot by thorns which need to toggle memory. for
   instance, see the trK Driver in the util thron.
   @enddesc 
@@*/

int EnableGFDataStorage(pGH *GH, pGF *GF) {
  int i, dir, dirp1, dirp2, sz;
  int special_pad, cache_size, start;

  char *char_temp;
  Double *Double_temp;
  int *int_temp;
  Complex *Complex_temp;
  Complex Czero = {0.0,0.0};

  static int ive_blathered = 0;

  DECLARE_PARAMETERS;

  if (zero_memory && !ive_blathered) {
    printf ("Zeroing memory for allocated GFs at alloc time\n");
  }
  
  if (padding_active && !ive_blathered) {
    printf ("PUGH Memory padding active. Stats:\n");
    printf ("  cacheline_bits   : %d\n",
	    padding_cacheline_bits);
    printf ("  size             : %d\n",
	    padding_size);
    printf ("  spacing          : %d\n",
	    padding_address_spacing);
  }
  
  ive_blathered = 1;

  if (GF->storage == PUGH_STORAGE) {
    printf ("WARNING: Tried to enable %s when already enabled\n",
	    GF->name);
    return;
  }

  /* Set up the send buffers. Note we only do this if
     we are not using the derived types comm style.
     */
  assert(!GF->send_buffer);
  GF->send_buffer = (Double **)malloc(6*sizeof(Double *));

  assert(!GF->recv_buffer);
  GF->recv_buffer = (Double **)malloc(6*sizeof(Double *));
  
  for (i=0;i<6;i++) {
    dir = i/2;
    dirp1 = (i/2+1)%3;
    dirp2 = (i/2+2)%3;
    
    if (dir == 0) {
      sz = GH->stencil_width * GH->lny * GH->lnz;
    } else if (dir == 1) {
      sz = GH->stencil_width * GH->lnx * GH->lnz;
    } else {
      sz = GH->stencil_width * GH->lnx * GH->lny;
    }
    
    if (GH->neighbors[GH->myproc][i] >= 0) {
      GF->buffer_sz[i] = sz;
      GF->send_buffer[i] = (Double *)malloc(sz*sizeof(Double));
      GF->recv_buffer[i] = (Double *)malloc(sz*sizeof(Double));
      assert(GF->recv_buffer[i]);
      assert(GF->send_buffer[i]);
    } else {
      GF->buffer_sz[i] = 0;
      GF->send_buffer[i] = NULL;
      GF->recv_buffer[i] = NULL;
    }
  }

  /* Set up the storage */
  if (GF->padddata) free (GF->padddata);
  if (!padding_active) {
    GF->padddata = (void *)malloc((GH->npoints)*GF->varsize);
    GF->data     = GF->padddata;	/* No padding case */

    if (zero_memory) {
      switch(GF->vtype)
      {
	case VARIABLE_CHAR: char_temp = GF->padddata;
	  for (i=0;i<GH->npoints;i++) char_temp[i] = 0;
	  break;
	case VARIABLE_INTEGER: int_temp = GF->padddata;
	  for (i=0;i<GH->npoints;i++) int_temp[i] = 0;
	  break;
	case VARIABLE_REAL: Double_temp = GF->padddata;
	  for (i=0;i<GH->npoints;i++) Double_temp[i] = 0.0;
	  break;
	case VARIABLE_COMPLEX: Complex_temp = GF->padddata;
	  for (i=0;i<GH->npoints;i++) Complex_temp[i] = Czero;
	  break;
	default:
	  fprintf(stderr, "Unknown variable size in PUGH enable GF storage\n");
      }

    }

  } else {
    GF->padddata  = (void *)malloc((GH->npoints + padding_size)
				    *GF->varsize);
   
    /* Align the actual starting address on a cache line.
       More specifically align on a unique cache line for
       each field. 
       */
    cache_size = 2 << (padding_cacheline_bits - 1);
    start = ((long) (GF->padddata) / GF->varsize)%cache_size;

    special_pad = ((GF->gfno) * padding_address_spacing + cache_size - start)%cache_size;

    char_temp = (char *)(GF->padddata) + special_pad*GF->varsize;
    GF->data = char_temp;

    if (special_pad > padding_size){
      printf("FATAL ERROR: padding size not large enough for cache type specified %d %d %d \n",special_pad,padding_size,cache_size);
      STOP;
    }

    if (zero_memory) {
      switch(GF->vtype)
      {
	case VARIABLE_CHAR: char_temp = GF->padddata;
	  for (i=0;i<GH->npoints+padding_size;i++) char_temp[i] = 0;
	  break;
	case VARIABLE_INTEGER: int_temp = GF->padddata;
	  for (i=0;i<GH->npoints+padding_size;i++) int_temp[i] = 0;
	  break;
	case VARIABLE_REAL: Double_temp = GF->padddata;
	  for (i=0;i<GH->npoints+padding_size;i++) Double_temp[i] = 0.0;
	  break;
	case VARIABLE_COMPLEX: Complex_temp = GF->padddata;
	  for (i=0;i<GH->npoints+padding_size;i++) Complex_temp[i] = Czero;
	  break;
	default:
	  fprintf(stderr, "Unknown variable size in PUGH enable GF storage\n");
      }
    }
  }

  if (!GF->padddata) {
    printf ("FATAL ERROR: Cannot allocate data for %s [%d]\n",
	    GF->name, GF->gfno);
    STOP;
  }
  GF->storage = PUGH_STORAGE;

  return 1;

}

 /*@@
   @routine    DisableGFDataStorage
   @date       Thu Apr  3 12:45:34 1997
   @author     Paul Walker
   @desc 
   This routine disables the grid function storage.
   That is, it un-allocates the 3D array and in its
   place allocates a single Double with the value
   0.0. This allows us to still have something to
   pass around (an array of size (1,1,1) in fortran
   speak) but also to not need all our 3D arrays 
   "on" all the time.
   @enddesc 
@@*/


int DisableGFDataStorage(pGH *GH, pGF *GF) {
  int i;

  if (GF->storage == PUGH_NOSTORAGE) {
    printf ("Warning: Tried to disable %s when already disabled\n",
	    GF->name);
    return;
  }

  if (GF->padddata) {
    free(GF->padddata);
    GF->padddata = NULL;
    GF->data = NULL;
  }

  if (GF->send_buffer) {
    for (i=0;i<6;i++)
      if (GF->send_buffer[i]) {
#ifdef MPI
	if(GF->sreq[i] != MPI_REQUEST_NULL){
	  CACTUS_MPI_ERROR(MPI_Request_free(&(GF->sreq[i])));
	}
#endif
    	free(GF->send_buffer[i]);
	GF->send_buffer[i] = NULL;
      }
    free(GF->send_buffer);
    GF->send_buffer = NULL;
  }


  if (GF->recv_buffer) {
    for (i=0;i<6;i++)
      if (GF->recv_buffer[i]) {
	free(GF->recv_buffer[i]);
	GF->recv_buffer[i] = NULL;
      }
    free(GF->recv_buffer);
    GF->recv_buffer = NULL;
  }

  GF->padddata = (Double *)malloc(sizeof(Double));
  GF->data = GF->padddata;
#if 0
  GF->data[0] = 0.0;		/* Very important! */
#endif
  GF->storage = PUGH_NOSTORAGE;
}


 /*@@
   @routine    SetGFComm
   @date       Thu Apr  3 12:46:23 1997
   @author     Paul Walker
   @desc 
   This sets the docomm[3] array of the GF based
   on the setting of the comm flag. Note the
   comm flag constants are defined in
   @seeheader pughDriver.h. As well as SetupPGF
   this is called by thorns (eg, elliptic) which
   need to toggle this state.
   @enddesc 
   @calledby   SetupPGF
@@*/


void SetGFComm(pGH *GH, pGF *res, int docomm) {
  /* Copy docomm */
  int i;
  res->commflag = docomm;

  for (i=0;i<6;i++)
    res->docomm[i] = 0;

  if (docomm == PUGH_NOCOMM) {
    for (i=0;i<6;i++)
      res->docomm[i] = 0;
  } else if (docomm == PUGH_ALLCOMM) {
    for (i=0;i<6;i++)
      res->docomm[i] = 1;
  } else if (docomm == PUGH_XCOMM) {
    res->docomm[0] = 1;
    res->docomm[1] = 1;
  } else if (docomm == PUGH_YCOMM) {
    res->docomm[2] = 1;
    res->docomm[3] = 1;
  } else if (docomm == PUGH_ZCOMM) {
    res->docomm[4] = 1;
    res->docomm[5] = 1;
  } else if (docomm == PUGH_PLUSFACESCOMM) {
    res->docomm[1] = 1;
    res->docomm[3] = 1;
    res->docomm[5] = 1;
  } else if (docomm == PUGH_MINUSFACESCOMM) {
    res->docomm[0] = 1;
    res->docomm[2] = 1;
    res->docomm[4] = 1;
  } else {
    printf("WARNING: Unsupported comm. in grid function %s .\n",res->name);
    printf("         Perhaps misspelled in GridFuncList? Defaulting to allcomm.\n");
    for (i=0;i<6;i++)
      res->docomm[i] = 1;
  }

  /* Handle nx = 1 type cases. This is only important for one
     processor MPI periodic boundaries.
   */
  /* This is a hack to get a GH. Unfortunately I decided not
     to pass a GH to this routine, and I should go back and
     change it, but that is a lot of changes so for now I'll
     do the following.
    */
  if (GH->nx == 1) {
    res->docomm[0] = 0;
    res->docomm[1] = 0;
  }

  if (GH->ny == 1) {
    res->docomm[2] = 0;
    res->docomm[3] = 0;
  }

  if (GH->nz == 1) {
    res->docomm[4] = 0;
    res->docomm[5] = 0;
  }
}