path: root/src/DumpVar.c

/*@@
   @file      DumpVar.c
   @date      Fri May 21 1999
   @author    Thomas Radke
   @desc 
   Do the actual writing of a variable, for output or for checkpointing.
   @enddesc 
   @history
   @hauthor Thomas Radke @hdate May 21 1999
   @hdesc Just copied from thorn FlexIO.
   @hendhistory
 @@*/


#include "cctk.h"
#include "cctk_Parameters.h"
#include "CactusPUGH/PUGH/src/include/pugh.h"
#include "CactusPUGH/PUGHSlab/src/PUGHSlab.h" 
#include "CactusBase/IOUtil/src/ioGH.h"
#include "ioHDF5GH.h"

#include <stdio.h>
#include <stdlib.h>
#ifdef SGI
#include <time.h>
#endif

#define IOTAGBASE 20000 /* This may break on more than 2000 processors */

#define IOHDF5_DEBUG

/* info structure describing how to dump a given CCTK variable type */
typedef struct {
  int iohdf5_type;              /* the HDF5 type to use */
  int element_size;             /* size of a single data element */
#ifdef CCTK_MPI
  MPI_Datatype mpi_type;        /* the data type for MPI communication */
#endif
} dumpInfo_t;

static char *char_time_date = NULL;

/* local function prototypes */
static void IOHDF5_DumpGS (cGH *GH, int index, int timelevel, hid_t iof,
                           int iohdf5_type);

static int  IOHDF5_DumpGA (cGH *GH, int index, int timelevel, ioHDF5Geo_t *slab_geo, 
			   hid_t iof,  dumpInfo_t *info);

static int  IOHDF5_getDumpData_ND (cGH *GH, int index, int timelevel, ioHDF5Geo_t *slab_geo, 
				   void **outme, int *free_outme, 
				   CCTK_INT *geom, int element_size);

static void IOHDF5_procDump (cGH *GH, int index, int timelevel, ioHDF5Geo_t *slab_geo, 
			     void *outme, CCTK_INT *geom, int proc, 
			     int iohdf5_type, hid_t iof);

static void IOHDF5_AddCommonAttributes (cGH *GH, int index, int timelevel, 
					CCTK_INT *global_shape, 
					hid_t dataset);

static void IOHDF5_AddCommonAttributes_ND (cGH *GH, int index, ioHDF5Geo_t *slab_geo, 
					   CCTK_INT *global_shape, hid_t dataset);

#ifdef CCTK_MPI
#ifdef HAVE_PARALLEL
static void IOHDF5_collectiveDump (cGH *GH, int index, int timelevel, ioHDF5Geo_t *slab_geo,
				   void *outme, CCTK_INT *geom, 
				   int hdf5io_type, hid_t iof);
#endif /* HAVE_PARALLEL */
#endif /* MPI */



/*@@
  @routine    IOHDF5_DumpVar
  @date       16 Apr 1999
  @author     Thomas Radke
  @desc
              Generic dump routine, just calls the appropriate dump routine
  @enddesc
  @calledby   IOHDF5_DumpGH IOHDF5_Write3D
  @history
  @endhistory
  @var        GH
  @vdesc      Pointer to CCTK grid hierarchy
  @vtype      cGH
  @vio        in
  @endvar
  @var        index
  @vdesc      global index of the variable to be dumped
  @vtype      int
  @vio        in
  @endvar
  @var        timelevel
  @vdesc      the timelevel to store
  @vtype      int
  @vio        in
  @endvar
  @var        iof
  @vdesc      the HDF5 file handle
  @vtype      hid_t
  @vio        in
  @endvar
@@*/
int IOHDF5_DumpVar (cGH *GH, int index, int timelevel, ioHDF5Geo_t *geo, hid_t iof)
{
  pGH *pughGH;
  ioHDF5GH *h5GH;
  ioGH *ioUtilGH;
  int variable_type;
  dumpInfo_t info;
  int retval = 0;

  /* Get the handle for PUGH, IOUtil, and IOHDF5 extensions */
  pughGH = PUGH_pGH (GH);
  ioUtilGH = (ioGH *) GH->extensions [CCTK_GHExtensionHandle ("IO")];
  h5GH = (ioHDF5GH *) GH->extensions [CCTK_GHExtensionHandle ("IOHDF5")];

  variable_type = CCTK_VarTypeI (index);

  switch (variable_type) {
    case CCTK_VARIABLE_CHAR:
      info.iohdf5_type = IOHDF5_CHAR;
      info.element_size = sizeof (CCTK_CHAR);
#ifdef CCTK_MPI
      info.mpi_type = PUGH_MPI_CHAR;
#endif
      break;
    case CCTK_VARIABLE_INT:
      info.iohdf5_type = IOHDF5_INT;
      info.element_size = sizeof (CCTK_INT);
#ifdef CCTK_MPI
      info.mpi_type = PUGH_MPI_INT;
#endif
      break;
    case CCTK_VARIABLE_REAL: 
      info.iohdf5_type  = ioUtilGH->out_single ? IOHDF5_REAL4 : IOHDF5_REAL;
      info.element_size = ioUtilGH->out_single ?
                              sizeof (CCTK_REAL4) : sizeof (CCTK_REAL);
#ifdef CCTK_MPI
      info.mpi_type     = ioUtilGH->out_single ? PUGH_MPI_REAL4 : PUGH_MPI_REAL;
#endif
      break;
    case CCTK_VARIABLE_COMPLEX:
      info.iohdf5_type  = h5GH->IOHDF5_COMPLEX;
      info.element_size = sizeof (CCTK_COMPLEX);
#ifdef CCTK_MPI
      info.mpi_type     = pughGH->PUGH_mpi_complex;
#endif
      break;
    default:
      CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING,
                  "Unsupported variable type %d", variable_type);
      return(-1);
  }

  switch (CCTK_GroupTypeFromVarI (index)) {
    case CCTK_SCALAR:
      IOHDF5_DumpGS (GH, index, timelevel, iof, info.iohdf5_type);
      break;
    case CCTK_ARRAY:
    case CCTK_GF:
      retval = IOHDF5_DumpGA (GH, index, timelevel, geo, iof, &info);
      break;
    default:
      CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING,
                  "Invalid group type %d in IOHDF5_DumpVar",
                  CCTK_GroupTypeFromVarI (index));
      retval = -1;
  }
  return(retval);
}


/*@@
  @routine    IOHDF5_DumpGS
  @date       May 21 1999
  @author     Thomas Radke
  @desc       Dumps a SCALAR type variable into a HDF5 file
  @enddesc 
  @calledby   IOHDF5_DumpVar
  @var        GH
  @vdesc      Pointer to CCTK grid hierarchy
  @vtype      cGH
  @vio        in
  @endvar
  @var        index
  @vdesc      global index of the variable to be dumped
  @vtype      int
  @vio        in
  @endvar
  @var        timelevel
  @vdesc      the timelevel to store
  @vtype      int
  @vio        in
  @endvar
  @var        iof
  @vdesc      the HDF5 file to dump to
  @vtype      hid_t
  @vio        in
  @endvar
  @var        iohdf5_type
  @vdesc      the HDF5 datatype
  @vtype      int
  @vio        in
  @endvar
  @history 
  @endhistory 

@@*/
static void IOHDF5_DumpGS (cGH *GH, int index, int timelevel, hid_t iof,
                           int iohdf5_type)
{
  CCTK_INT global_shape [1] = {0};
  ioGH *ioUtilGH;
  ioHDF5GH *h5GH;
  hid_t dataset, dataspace;
  char *name, *datasetname;

  /* Get the handles for IOHDF5 and IOUtil extensions */
  ioUtilGH = (ioGH *) GH->extensions [CCTK_GHExtensionHandle ("IO")];
  h5GH = (ioHDF5GH *) GH->extensions [CCTK_GHExtensionHandle ("IOHDF5")];

  if (CCTK_MyProc (GH) != ioUtilGH->ioproc)
    return;

  CACTUS_IOHDF5_ERROR (dataspace = H5Screate (H5S_SCALAR));
  name = CCTK_FullName (index);
  datasetname = (char *) malloc (strlen (name) + 30);
  sprintf (datasetname, "%s@%d@%d", name, GH->cctk_iteration, timelevel);
  free (name);

  /* if restart from recovery delete an already existing dataset
     so that it can be created as anew */
  if (h5GH->checkForExistingObjects [index]) {
    CACTUS_IOHDF5_ERROR (H5Eset_auto (NULL, NULL));
    H5Gunlink (iof, datasetname);
    CACTUS_IOHDF5_ERROR(H5Eset_auto(h5GH->printErrorFn, h5GH->printErrorFnArg));
  }
  CACTUS_IOHDF5_ERROR (dataset = H5Dcreate (iof, datasetname,
                       iohdf5_type, dataspace, H5P_DEFAULT));
  free (datasetname);
  CACTUS_IOHDF5_ERROR (H5Dwrite (dataset, iohdf5_type, H5S_ALL, H5S_ALL,
                       H5P_DEFAULT,
                       CCTK_VarDataPtrI (GH, timelevel, index)));
    
  IOHDF5_AddCommonAttributes (GH, index, timelevel, global_shape, dataset);

  CACTUS_IOHDF5_ERROR (H5Dclose (dataset));
  CACTUS_IOHDF5_ERROR (H5Sclose (dataspace));
}


/*@@
  @routine    IOHDF5_DumpGA
  @date       May 21 1999
  @author     Paul Walker
  @desc       Dumps a GA type variable into a HDF5 file
  @enddesc 
  @calledby   IOHDF5_DumpVar
  @var        GH
  @vdesc      Pointer to CCTK grid hierarchy
  @vtype      cGH
  @vio        in
  @endvar
  @var        index
  @vdesc      global index of the variable to be dumped
  @vtype      int
  @vio        in
  @endvar
  @var        timelevel
  @vdesc      the timelevel to store
  @vtype      int
  @vio        in
  @endvar
  @var        iof
  @vdesc      the HDF5 file to dump to
  @vtype      hid_t
  @vio        in
  @endvar
  @var        info
  @vdesc      info structure describing
                - the HDF5 datatype to store
                - the size of an element of variable
                - the MPI datatype to communicate
  @vtype      dumpInfo_t
  @vio        in
  @endvar

  @history 
  @endhistory 

@@*/
static int IOHDF5_DumpGA (cGH *GH, int index, int timelevel, ioHDF5Geo_t *slab_geo,  
			   hid_t iof, dumpInfo_t *info)
{
  DECLARE_CCTK_PARAMETERS
  ioGH *ioUtilGH;
  pGH *pughGH;
  int myproc;
  int nprocs;
  int sdim;
  CCTK_INT *geom; /* Lower bounds, sizes and global shape of data */
  void *outme;    /* The data pointer to dump ... */
  int free_outme; /* and whether it needs freeing */
  int retval=0;
#ifdef CCTK_MPI
  int i, j;
  MPI_Status ms;
#endif

  /* Get the handles for PUGH and IOUtil extensions */
  pughGH   = PUGH_pGH (GH);
  ioUtilGH = (ioGH *) GH->extensions [CCTK_GHExtensionHandle ("IO")];

  sdim= slab_geo->sdim;

  /* allocate the geometry buffer */
  geom = (CCTK_INT *) malloc (3*sdim * sizeof (CCTK_INT));
  
  /* Get the pointer to the data we want to output (includes downsampling) */
  if (IOHDF5_getDumpData_ND(GH, index, timelevel, slab_geo, 
			    &outme, &free_outme, geom,
			    info->element_size)<0)
  {
    return(-1);
  }

  myproc = CCTK_MyProc (GH);
  nprocs = CCTK_nProcs (GH);

#ifdef CCTK_MPI
#ifdef HAVE_PARALLEL
  if (ioUtilGH->unchunked) {
    IOHDF5_collectiveDump (GH, index, timelevel, slab_geo,
			   outme, geom,
                           info->iohdf5_type, iof);
    if (free_outme)
      free (outme);
    free (geom);
    return;
  }
#endif
#endif

  /* Dump data held on IO processor */
  if (myproc == ioUtilGH->ioproc) {
    IOHDF5_procDump (GH, index, timelevel, slab_geo, 
		     outme, geom, myproc,
                     info->iohdf5_type,iof);
  }
#ifdef CCTK_MPI
  if (myproc == ioUtilGH->ioproc) {

    /* Dump data from all other processors */
    for (i = myproc+1; i < myproc+ioUtilGH->ioproc_every && i < nprocs; i++) {

      /* Allocate temp to the right size (based on GH->ub [i] et..) */
      void *tmpd;
      int incoming;

      /* receive geometry */
      CACTUS_MPI_ERROR (MPI_Recv (geom, 3*sdim, PUGH_MPI_INT, i, 2*i+IOTAGBASE+1,
                        pughGH->PUGH_COMM_WORLD, &ms));

      incoming = geom [sdim];
      for (j = 1; j < sdim; j++)
        incoming *= geom [sdim + j];

      /* receive data if incoming>0 */
      if (incoming>0) {

	tmpd = malloc (incoming * info->element_size);
	CACTUS_MPI_ERROR (MPI_Recv (tmpd, incoming, info->mpi_type, i,
				    2*i+IOTAGBASE, pughGH->PUGH_COMM_WORLD, &ms));
	
	IOHDF5_procDump (GH, index, timelevel, slab_geo,
			 tmpd, geom, i,
			 info->iohdf5_type, iof);
	free (tmpd);
      }
      
    } /* End loop over processors */

  } else {

    /* Send the geometry and data from the non-IO processors
     * to the IO processors
     */

    int outgoing = geom [sdim];

      for (i = 1; i < sdim; i++)
	outgoing *= geom [sdim + i];

      /* send geometry */
      CACTUS_MPI_ERROR (MPI_Send (geom, 3*sdim, PUGH_MPI_INT, ioUtilGH->ioproc,
				  2*myproc+IOTAGBASE+1, pughGH->PUGH_COMM_WORLD));
    if (outgoing>0) {

      /* send data if outgoing>0 */
	CACTUS_MPI_ERROR (MPI_Send (outme, outgoing, info->mpi_type, ioUtilGH->ioproc,
				    2*myproc+IOTAGBASE, pughGH->PUGH_COMM_WORLD));
#ifdef IOHDF5_DEBUG
	printf ("Processor %d sent %d data points\n", myproc, outgoing);
#endif
    }
  }

  CCTK_Barrier (GH);
#endif

  if (free_outme)
    free (outme);

  free (geom);
  return(retval);
}


/***************************** local functions ****************************/

/*@@
  @routine    IOHDF5_AddCommonAttributes
  @date       May 21 1999
  @author     Thomas Radke
  @desc 
  Add "Common" attributes, these are:
  <ul>
    <li> the variable's groupname
    <li> the grouptype
    <li> number of timelevels
    <li> the current date
         Note that the datestamp should be turned of if you are byte 
         comparing two output files.
    <li> simulation time
    <li> origin
    <li> bounding box
    <li> gridspacings (both downsampled and evolution)
    <li> global grid size
  </ul>
  @enddesc 
  @calledby   IOHDF5_DumpGS, IOHDF5_collectiveDump, IOHDF5_procDump
  @history 
  @endhistory 
  @var        GH
  @vdesc      Pointer to CCTK grid hierarchy
  @vtype      cGH
  @vio        in
  @endvar
  @var        index
  @vdesc      global index of the variable to be dumped
  @vtype      int
  @vio        in
  @endvar
  @var        timelevel
  @vdesc      the timelevel to store
  @vtype      int
  @vio        in
  @endvar
  @var        global_shape
  @vdesc      the global shape of the dataset
  @vtype      CCTK_INT [dim]
  @vio        in
  @endvar
  @var        dataset
  @vdesc      the dataset handle where the attributes should be attached to
  @vtype      hid_t
  @vio        in
  @endvar
@@*/
static void IOHDF5_AddCommonAttributes (cGH *GH, int index, int timelevel,
                                        CCTK_INT *global_shape, hid_t dataset)
{
  DECLARE_CCTK_PARAMETERS
  int vdim;
  CCTK_INT intAttr;
  CCTK_REAL dummy;
  CCTK_REAL realAttr [6];
  char *groupname;
  ioGH *ioUtilGH;
  ioHDF5GH *h5GH;

  /* Get the handles for IOUtil and IOHDF5 extensions */
  ioUtilGH = (ioGH *) GH->extensions [CCTK_GHExtensionHandle ("IO")];
  h5GH = (ioHDF5GH *) GH->extensions [CCTK_GHExtensionHandle ("IOHDF5")];

  /* get the dimension of the variable */
  vdim = CCTK_GroupDimI (CCTK_GroupIndexFromVarI (index));

  groupname = CCTK_GroupNameFromVarI (index);
  WRITE_ATTRIBUTE ("groupname", groupname,
                   dataset, h5GH->scalarDataspace, 0, h5GH->IOHDF5_STRING);
  free (groupname);
  intAttr = CCTK_GroupTypeFromVarI (index);
  WRITE_ATTRIBUTE ("grouptype", &intAttr, dataset, h5GH->scalarDataspace, 0,
                   IOHDF5_INT);
  intAttr = CCTK_NumTimeLevelsFromVarI (index);
  WRITE_ATTRIBUTE ("ntimelevels", &intAttr, dataset, h5GH->scalarDataspace, 0,
                   IOHDF5_INT);

  if (char_time_date && out3D_datestamp) 
    WRITE_ATTRIBUTE ("date", char_time_date, dataset, h5GH->scalarDataspace, 0,
                     h5GH->IOHDF5_STRING);

  WRITE_ATTRIBUTE ("time", &GH->cctk_time, dataset, h5GH->scalarDataspace, 0,
                   IOHDF5_REAL);

  /* NOTE: the attributes "origin", "min_ext", "max_ext", and "delta"
           are always stored as 3-dimensional points */
  CCTK_CoordRange (GH, &realAttr [0], &dummy, -1, "x", "cart3d");
  CCTK_CoordRange (GH, &realAttr [1], &dummy, -1, "y", "cart3d");
  CCTK_CoordRange (GH, &realAttr [2], &dummy, -1, "z", "cart3d");
  WRITE_ATTRIBUTE ("origin", realAttr, dataset, h5GH->arrayDataspace, 3,
                   IOHDF5_REAL);

  CCTK_CoordRange (GH, &realAttr [0], &realAttr [3], -1, "x", "cart3d");
  CCTK_CoordRange (GH, &realAttr [1], &realAttr [4], -1, "y", "cart3d");
  CCTK_CoordRange (GH, &realAttr [2], &realAttr [5], -1, "z", "cart3d");
  WRITE_ATTRIBUTE ("min_ext", realAttr, dataset, h5GH->arrayDataspace, 3,
                   IOHDF5_REAL);
  WRITE_ATTRIBUTE ("max_ext", realAttr + 3, dataset, h5GH->arrayDataspace, 3,
                   IOHDF5_REAL);
 

  realAttr [0] = GH->cctk_delta_space [0] * ioUtilGH->downsample [0]; 
  realAttr [1] = GH->cctk_delta_space [1] * ioUtilGH->downsample [1]; 
  realAttr [2] = GH->cctk_delta_space [2] * ioUtilGH->downsample [2];
  WRITE_ATTRIBUTE ("delta", realAttr, dataset, h5GH->arrayDataspace, 3,
                   IOHDF5_REAL);

  if (ioUtilGH->downsample [0] > 1 ||
      ioUtilGH->downsample [1] > 1 ||
      ioUtilGH->downsample [2] > 1)
    WRITE_ATTRIBUTE ("evolution_delta", GH->cctk_delta_space, dataset,
                     h5GH->arrayDataspace, 3, IOHDF5_REAL);

  WRITE_ATTRIBUTE ("global_size", global_shape, dataset, h5GH->arrayDataspace,
                   vdim, IOHDF5_INT);

}

static void IOHDF5_AddCommonAttributes_ND (cGH *GH, int index, ioHDF5Geo_t *geo,
					   CCTK_INT *global_shape, hid_t dataset)
{
  DECLARE_CCTK_PARAMETERS
  CCTK_INT intAttr;
  CCTK_REAL dummy;
  CCTK_REAL realAttr [6];
  char *groupname;

  ioGH *ioUtilGH;
  ioHDF5GH *h5GH;

  int idim;

  /* Get the handles for IOUtil and IOHDF5 extensions */
  ioUtilGH = (ioGH *) GH->extensions [CCTK_GHExtensionHandle ("IO")];
  h5GH = (ioHDF5GH *) GH->extensions [CCTK_GHExtensionHandle ("IOHDF5")];

  /* VARIABLE INFO */
  groupname = CCTK_GroupNameFromVarI (index);
  WRITE_ATTRIBUTE ("groupname", groupname,
                   dataset, h5GH->scalarDataspace, 0, h5GH->IOHDF5_STRING);
  free (groupname);
  intAttr = CCTK_GroupTypeFromVarI (index);
  WRITE_ATTRIBUTE ("grouptype", &intAttr, dataset, h5GH->scalarDataspace, 0,
                   IOHDF5_INT);

  /* TIMLEVELS (grid) */
  intAttr = CCTK_NumTimeLevelsFromVarI (index);
  WRITE_ATTRIBUTE ("ntimelevels", &intAttr, dataset, h5GH->scalarDataspace, 0,
                   IOHDF5_INT);

  /* DATE  */
  if (char_time_date && out3D_datestamp) 
    WRITE_ATTRIBUTE ("date", char_time_date, dataset, h5GH->scalarDataspace, 0,
                     h5GH->IOHDF5_STRING);

  /* TIME (grid) */
  WRITE_ATTRIBUTE ("time", &GH->cctk_time, dataset, h5GH->scalarDataspace, 0,
                   IOHDF5_REAL);

  /* ORIGIN (grid) */
  CCTK_CoordRange (GH, &realAttr [0], &dummy, -1, "x", "cart3d");
  CCTK_CoordRange (GH, &realAttr [1], &dummy, -1, "y", "cart3d");
  CCTK_CoordRange (GH, &realAttr [2], &dummy, -1, "z", "cart3d");
  WRITE_ATTRIBUTE ("origin", realAttr, dataset, h5GH->arrayDataspace, geo->vdim,
                   IOHDF5_REAL);

  /* ORIGIN (slab) -  slab may be offset by (slab_start) grid points */
  for (idim=0;idim<geo->sdim;idim++)
    realAttr[idim]=geo->slab_start[geo->direction[idim]]*
                   GH->cctk_delta_space[geo->direction[idim]];
  WRITE_ATTRIBUTE ("origin_slab", realAttr, dataset, h5GH->arrayDataspace, geo->sdim,
                   IOHDF5_REAL);
    
  /* MIN/MAX extent (grid) */
  CCTK_CoordRange (GH, &realAttr [0], &realAttr [3], -1, "x", "cart3d");
  CCTK_CoordRange (GH, &realAttr [1], &realAttr [4], -1, "y", "cart3d");
  CCTK_CoordRange (GH, &realAttr [2], &realAttr [5], -1, "z", "cart3d");
  WRITE_ATTRIBUTE ("min_ext", realAttr, dataset, h5GH->arrayDataspace, geo->vdim,
                   IOHDF5_REAL);
  WRITE_ATTRIBUTE ("max_ext", realAttr + 3, dataset, h5GH->arrayDataspace, geo->vdim,
                   IOHDF5_REAL);

  /* MIN/MAX extent (slab) -  slab may be smaller than grid*/
  for (idim=0;idim<geo->sdim;idim++) {
    realAttr[idim  ]=(geo->slab_start[idim]*GH->cctk_delta_space[idim]);
    realAttr[idim+3]=(geo->slab_start[idim]+geo->actlen[idim]-1)*
      (GH->cctk_delta_space[idim]*geo->downs[idim]);
  }
  WRITE_ATTRIBUTE ("min_ext_slab", realAttr, dataset, h5GH->arrayDataspace, geo->sdim,
                   IOHDF5_REAL);
  WRITE_ATTRIBUTE ("max_ext_slab", realAttr + 3, dataset, h5GH->arrayDataspace, geo->sdim,
                   IOHDF5_REAL);

  /* DELTA SPACING (grid) */
  realAttr [0] = GH->cctk_delta_space [0];
  realAttr [1] = GH->cctk_delta_space [1];
  realAttr [2] = GH->cctk_delta_space [2];
  WRITE_ATTRIBUTE ("delta", realAttr, dataset, h5GH->arrayDataspace, geo->vdim,
                   IOHDF5_REAL);

  /* DELTA SPACING (slab) */
  for (idim=0;idim<geo->sdim;idim++)
    realAttr[idim]=GH->cctk_delta_space[geo->direction[idim]]*
                   geo->downs[geo->direction[idim]];
  WRITE_ATTRIBUTE ("delta_slab", realAttr, dataset, h5GH->arrayDataspace, geo->sdim,
                   IOHDF5_REAL);

  /* FIXME: WHAT IS THIS ? */
  if (ioUtilGH->downsample [0] > 1 ||
      ioUtilGH->downsample [1] > 1 ||
      ioUtilGH->downsample [2] > 1)
    WRITE_ATTRIBUTE ("evolution_delta", GH->cctk_delta_space, dataset,
                     h5GH->arrayDataspace, 3, IOHDF5_REAL);

  /* GLOBAL SIZE (slab) */
  WRITE_ATTRIBUTE ("global_size_grid", GH->cctk_gsh, dataset, h5GH->arrayDataspace,
                   geo->vdim, IOHDF5_INT);

  /* GLOBAL SIZE (slab) */
  WRITE_ATTRIBUTE ("global_size", global_shape, dataset, h5GH->arrayDataspace,
                   geo->sdim, IOHDF5_INT);

  /* GLOBAL SIZE (slab) */
  WRITE_ATTRIBUTE ("global_size_slab", geo->actlen, dataset, h5GH->arrayDataspace,
                   geo->sdim, IOHDF5_INT);
}


static int IOHDF5_getDumpData_ND (cGH *GH, int index, int timelevel, ioHDF5Geo_t *slab_geo, 
				   void **outme, int *free_outme, 
				   CCTK_INT *geom, int element_size)
{
  ioGH *ioUtilGH;
  pGExtras *extras=NULL;
  
  CCTK_REAL4 *single_ptr;

  int myproc, do_downsample;
  int idim, sdim, vdim;                      /*iteration, slab dim, variable dim */
  int *hsizes,*hsizes_global,*hsizes_offset; /* geometry information */
  int sdir[3];                               /* slab direction FIXME */
  int locdowns[3];                           /* Holds the downsampling in the order 
						specified in slab_geo->direction[] */
  int slabstart[3];                          /* global start index for slab to extract */
  int ip;

  void *data = CCTK_VarDataPtrI (GH, timelevel, index);
  

  myproc = CCTK_MyProc (GH);

  /* get GH extension for IOUtil */
  ioUtilGH = (ioGH *) GH->extensions [CCTK_GHExtensionHandle ("IO")];

  /* get the pointer to PUGH specific structures */
  extras = ((pGA***) PUGH_pGH(GH)->variables) [index][timelevel]->extras;

  /* Shortcuts */
  vdim = slab_geo->vdim;
  sdim = slab_geo->sdim;
  
  if (vdim>3) {
    CCTK_WARN(1,"Cannot treat GFs with dim>3, check with cactus support\n");
    return(-1);
  }

  for (idim=0;idim<sdim;idim++)
    locdowns[idim]=slab_geo->downs[slab_geo->direction[idim]];

  /* Set downsample flag */
  do_downsample = 0;
  for (idim = 0; idim < sdim; idim++)
    do_downsample |= locdowns[idim] > 1;

  /* Simple case if no downsampling and sdim==vdim */
  if ((! do_downsample)&& (sdim==vdim)) 
  {
    if (ioUtilGH->out_single) {
      single_ptr = (CCTK_REAL4 *) malloc (extras->npoints*sizeof (CCTK_REAL4));
      for (ip = 0; ip < extras->npoints; ip++)
        single_ptr [ip] = (CCTK_REAL4) ((CCTK_REAL *) data) [ip];

      *outme      = single_ptr;
      *free_outme = 1;
    } 
    else 
    {
      *outme      = data;
      *free_outme = 0;
    }

    for (idim = 0; idim < sdim; idim++) {
      geom [idim + 0*sdim] = extras->lb [myproc][idim];
      geom [idim + 1*sdim] = extras->lnsize [idim];        
      geom [idim + 2*sdim] = extras->nsize [idim];
    }
  }
  else  /* Call Hyperslab else */ {  

    /* allocate array to hold size information of slab */
    hsizes = (int*)malloc(sdim*sizeof(int));
    hsizes_global= (int*)malloc(sdim*sizeof(int));
    hsizes_offset= (int*)malloc(sdim*sizeof(int));

    for (idim=0;idim<sdim;idim++) {
      hsizes[idim]=0;
      hsizes_global[idim]=0;
      hsizes_offset[idim]=0;
    }

    /* TEMPORARY FIX DUE TO INCONSISTENT DIR SPECIFICATION */
    /* direction vector of 1d line in 3d volume */
    if (sdim==1) {
      sdir[0] = (slab_geo->direction[0]==0) ? 1:0;
      sdir[1] = (slab_geo->direction[0]==1) ? 1:0;
      sdir[2] = (slab_geo->direction[0]==2) ? 1:0;
    }
    /* norm vector of 2d surface in 3d volume */
    else if (sdim==2)
      {
	sdir[0] = ((slab_geo->direction[0]==1)&&(slab_geo->direction[1]==2)) ? 1:0;
	sdir[1] = ((slab_geo->direction[0]==0)&&(slab_geo->direction[1]==2)) ? 1:0;
	sdir[2] = ((slab_geo->direction[0]==0)&&(slab_geo->direction[1]==1)) ? 1:0;
      }
    /* spanning directions for 3d */
    else if (sdim==3)
      {
	sdir[0] = slab_geo->direction[0];
	sdir[1] = slab_geo->direction[1];
	sdir[2] = slab_geo->direction[2];
      }
    
    /* Get the start indeces for the slab from slab_start. */
    /* Slab_start from the geo structure is not a true start index, it is currently
       used as the intersection point of three surfaces , that's why two entries
       have to be set to zero in the case of a surface (one entry for line) FIXME */
    for (idim=0;idim<vdim;idim++)
      slabstart[idim] = slab_geo->slab_start[idim];
    for (idim=0;idim<sdim;idim++)
      slabstart[slab_geo->direction[idim]]=0;
          

    if (Hyperslab_GetLocalHyperslab (GH, index, timelevel, sdim, slab_geo->slab_start,
				     sdir, slab_geo->length, locdowns, outme,
				     hsizes, hsizes_global, hsizes_offset)< 0) {
      char *fullname = CCTK_FullName (index);
      CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING,
		  "Failed to extract hyperslab for variable '%s'", fullname);
      free (fullname); 
      *free_outme = 0;
      return(-1);
      
    }
    
    *free_outme = 1;

    for (idim = 0; idim < sdim; idim++) {
      geom [idim + 0*sdim]  = hsizes_offset[idim];
      geom [idim + 1*sdim]  = hsizes[idim];        
      geom [idim + 2*sdim]  = hsizes_global[idim];
      slab_geo->actlen[idim]= hsizes_global[idim];
    }
  }
#ifdef IOHDF5_DEBUG
  printf ("***** %dD \n",sdim);
  printf ("Lower bound: %d", (int) geom [0]);
  for (idim = 1; idim < sdim; idim++)
    printf (" %d", (int) geom [idim]);
  printf ("\n");
  printf ("Chunk size : %d", (int) geom [1*sdim + 0]);
  for (idim = 1; idim < sdim; idim++)
    printf (" %d", (int) geom [1*sdim + idim]);
  printf ("\n");
  printf ("Global size: %d", (int) geom [2*sdim + 0]);
  for (idim = 1; idim < sdim; idim++) 
    printf (" %d", (int) geom [2*sdim + idim]);
  printf ("\n");
  printf ("Downsampling: ");
  for (idim = 0; idim < sdim; idim++)
    printf (" %d", locdowns[idim]);
  printf ("\n"); 
  printf ("Direction: "); 
  for (idim = 0; idim < sdim; idim++)
    printf (" %d",slab_geo->direction[idim]);
  printf ("\n");  
  printf ("Slab Start: %d", (int) slabstart [0]);
  for (idim = 1; idim < vdim; idim++)
    printf (" %d", (int) slabstart[idim]);
  printf ("\n"); 
  printf ("Slab intersect cntr: %d", (int) slab_geo->slab_start[0]);
  for (idim = 1; idim < vdim; idim++)
    printf (" %d", (int) slab_geo->slab_start[idim]);
  printf ("\n");  
  printf ("SlabDelta: "); 
  for (idim = 0; idim < sdim; idim++)
    printf (" %f",GH->cctk_delta_space[slab_geo->direction[idim]]*locdowns[idim]);
  printf ("\n");
#endif

  return(1);
}


/*@@
  @routine IOHDF5_procDump
  @author  Thomas Radke
  @date    May 21 1999
  @desc    All IO processors dump the data of their group into the file,
           either as one chunk per processor or unchunked
           (depending on parameter "unchunked").
  @enddesc
  @calledby   IOHDF5_DumpGA
  @history 
  @endhistory 
  @var        GH
  @vdesc      Pointer to CCTK grid hierarchy
  @vtype      cGH
  @vio        in
  @endvar
  @var        index
  @vdesc      global index of the variable to be dumped
  @vtype      int
  @vio        in
  @endvar
  @var        timelevel
  @vdesc      the timelevel to store
  @vtype      int
  @vio        in
  @endvar
  @var        outme
  @vdesc      pointer to the chunk to dump
  @vtype      void *
  @vio        in
  @endvar
  @var        geom
  @vdesc      bounds and sizes of the output
  @vtype      CCTK_INT [3*dim]
  @vio        in
  @vcomment   geom [0*dim..1*dim-1]: lower bounds
              geom [1*dim..2*dim-1]: (local) data sizes
              geom [2*dim..3*dim-1]: global space
  @endvar 
  @var        proc
  @vdesc      the processor which's chunk is to be dumped
  @vtype      int
  @vio        in
  @endvar
  @var        hdf5io_type
  @vdesc      the HDF5 datatype identifier
  @vtype      int
  @vio        in
  @endvar
  @var        iof
  @vdesc      the HDF5 file handle
  @vtype      hid_t
  @vio        in
  @endvar
@@*/
static void IOHDF5_procDump (cGH *GH, int index, int timelevel, ioHDF5Geo_t *slab_geo,
			     void *outme, CCTK_INT *geom, int proc, 
			     int iohdf5_type, hid_t iof)
{
  DECLARE_CCTK_PARAMETERS
  int i, vdim,sdim, idim;
  ioGH *ioUtilGH;
  ioHDF5GH *h5GH;
  hid_t group, dataset, memspace, filespace;
  char *name, datasetname [128], chunkname [138];
  hssize_t *chunk_origin;
  hsize_t *chunk_dims, *file_dims;
  int locpoints=0;


  ioUtilGH = (ioGH *) GH->extensions [CCTK_GHExtensionHandle ("IO")];
  h5GH = (ioHDF5GH *) GH->extensions [CCTK_GHExtensionHandle ("IOHDF5")];

  vdim = slab_geo->vdim;
  sdim = slab_geo->sdim;

  /* Check if our proc has any points to output */
  locpoints = geom [sdim];
  for (i = 1; i < sdim; i++)
    locpoints *= geom [sdim + i];
  if (locpoints<=0) return;

  chunk_origin = (hssize_t *) malloc (sdim * sizeof (hssize_t));
  chunk_dims   = (hsize_t *)  malloc (2*sdim * sizeof (hsize_t));
  file_dims    = chunk_dims + sdim;

  /* HDF5 needs it in reverse order */
  for (i = 0; i < sdim; i++) {
    chunk_origin [i] = geom [1*sdim - 1 - i];
    chunk_dims   [i] = geom [2*sdim - 1 - i];
    file_dims    [i] = geom [3*sdim - 1 - i];
  }

  /* build the unique dataset name: eg. WAVETOY::phi@v3s01@3@1  
     volume (v) 3-dim ; slab (s) spanned by 0,1 (->2-dim) 
     at(@) iteration at(@) timelevel
  */
  name = CCTK_FullName (index);
  sprintf (datasetname, "%s@v%ds",name,vdim);
  for (idim=0;idim<sdim;idim++)
   sprintf (datasetname, "%s%d",datasetname, slab_geo->direction[idim]);
  sprintf (datasetname, "%s@%d@%d",datasetname, GH->cctk_iteration, timelevel);
  free (name);

  printf("Datasetname: >%s< \n",datasetname);

  /* create the memspace according to chunk dims */
  CACTUS_IOHDF5_ERROR (memspace = H5Screate_simple (sdim, chunk_dims, NULL));
 
  if (ioUtilGH->unchunked) {

    /* create the (global) filespace and set the hyperslab for the chunk */
    CACTUS_IOHDF5_ERROR (filespace = H5Screate_simple (sdim, file_dims, NULL));
    CACTUS_IOHDF5_ERROR (H5Sselect_hyperslab (filespace, H5S_SELECT_SET,
                         chunk_origin, NULL, chunk_dims, NULL));
   
    /* the IO processor creates the dataset and adds the common attributes
       when writing its own data, otherwise the dataset is reopened */
    if (CCTK_MyProc (GH) == proc) {
      /* if restart from recovery delete an already existing dataset
         so that it can be created as anew */
      if (h5GH->checkForExistingObjects [index]) {
        CACTUS_IOHDF5_ERROR (H5Eset_auto (NULL, NULL));
        H5Gunlink (iof, datasetname);
        H5Eset_auto (h5GH->printErrorFn, h5GH->printErrorFnArg);
      }
      CACTUS_IOHDF5_ERROR (dataset = H5Dcreate (iof, datasetname,
						iohdf5_type, filespace, H5P_DEFAULT));
      IOHDF5_AddCommonAttributes_ND (GH, index, slab_geo, 
				     &geom [2*sdim], dataset);
    } else
      CACTUS_IOHDF5_ERROR (dataset = H5Dopen (iof, datasetname));

    /* write the data */
    CACTUS_IOHDF5_ERROR (H5Dwrite (dataset, iohdf5_type, memspace,
                         filespace, H5P_DEFAULT, outme));

    /* and close the file dataspace */
    CACTUS_IOHDF5_ERROR (H5Sclose (filespace));

  } else {

    /* the IO processor creates the chunk group and adds common attributes */
    if (CCTK_MyProc (GH) == proc) {
      /* if restart from recovery delete an already existing group
         so that it can be created as anew */
      if (h5GH->checkForExistingObjects [index]) {
        CACTUS_IOHDF5_ERROR (H5Eset_auto (NULL, NULL));
        H5Gunlink (iof, datasetname);
        H5Eset_auto (h5GH->printErrorFn, h5GH->printErrorFnArg);
      }
      CACTUS_IOHDF5_ERROR (group = H5Gcreate (iof, datasetname, 0));
      IOHDF5_AddCommonAttributes_ND (GH, index, slab_geo, 
				     &geom [2*sdim], group);
      CACTUS_IOHDF5_ERROR (H5Gclose (group));
    }

    /* now the chunk dataset for each processor is created within the group */
    sprintf (chunkname, "%s/chunk%d", datasetname, proc - CCTK_MyProc (GH));
    printf("chunkname: >%s< \n",chunkname);
    /* create the chunk dataset and dump the chunk data */
    CACTUS_IOHDF5_ERROR (dataset = H5Dcreate (iof, chunkname,
                         iohdf5_type, memspace, H5P_DEFAULT));
    CACTUS_IOHDF5_ERROR (H5Dwrite (dataset, iohdf5_type, H5S_ALL,
                         H5S_ALL, H5P_DEFAULT, outme));

    /* add the "origin" attribute for the chunk */
    WRITE_ATTRIBUTE ("chunk_origin", geom, dataset, h5GH->arrayDataspace, sdim,
                     IOHDF5_INT);
  }

  /* close the dataset and the memspace */
  CACTUS_IOHDF5_ERROR (H5Dclose (dataset));
  CACTUS_IOHDF5_ERROR (H5Sclose (memspace));

  free (chunk_origin);
  free (chunk_dims);
}


#ifdef CCTK_MPI
#ifdef HAVE_PARALLEL
/*@@
  @routine IOHDF5_collectiveDump
  @author  Thomas Radke
  @date    May 21 1999
  @desc    All processors dump their data into an unchunked file.
  @enddesc
  @calledby   IOHDF5_DumpGA
  @history
  @endhistory
  @var        GH
  @vdesc      Pointer to CCTK grid hierarchy
  @vtype      cGH
  @vio        in
  @endvar
  @var        index
  @vdesc      global index of the variable to be dumped
  @vtype      int
  @vio        in
  @endvar
  @var        timelevel
  @vdesc      the timelevel to store
  @vtype      int
  @vio        in
  @endvar
  @var        outme
  @vdesc      pointer to the chunk to dump
  @vtype      void *
  @vio        in
  @endvar
  @var        geom
  @vdesc      bounds and sizes of the output
  @vtype      CCTK_INT [3*dim]
  @vio        in
  @vcomment   geom [0*dim..1*dim-1]: lower bounds
              geom [1*dim..2*dim-1]: (local) data sizes
              geom [2*dim..3*dim-1]: global space
  @endvar 
  @var        hdf5io_type
  @vdesc      the HDF5 datatype identifier
  @vtype      int
  @vio        in
  @endvar
  @var        iof
  @vdesc      the HDF5 file handle
  @vtype      hid_t
  @vio        in
  @endvar
@@*/
static void IOHDF5_collectiveDump (cGH *GH, int index, int timelevel, ioHDF5Geo_t *slab_geo,  
				   void *outme, CCTK_INT *geom, 
				   int hdf5io_type, hid_t iof)
{
  DECLARE_CCTK_PARAMETERS
  int i, dim;
  ioHDF5GH *h5GH;
  hid_t dataset, memspace, filespace;
  char *name, datasetname [128];
  hssize_t *chunk_origin;
  hsize_t *chunk_dims, *file_dims;

  h5GH = (ioHDF5GH *) GH->extensions [CCTK_GHExtensionHandle ("IOHDF5")];

  /* get the dimension of the variable */
  sdim = slab_geo->sdim;

  chunk_origin = (hssize_t *) malloc (sdim * sizeof (hssize_t));
  chunk_dims   = (hsize_t *)  malloc (2*sdim * sizeof (hsize_t));
  file_dims    = chunk_dims + sdim;

  /* HDF5 needs it in reverse order */
  for (i = 0; i < sdim; i++) {
    chunk_origin [i] = geom [1*sdim - 1 - i];
    chunk_dims   [i] = geom [2*sdim - 1 - i];
    file_dims    [i] = geom [3*sdim - 1 - i];
  }

  /* build the unique dataset name */
  name = CCTK_FullName (index);
  sprintf (datasetname, "%s@%d@%d", name, GH->cctk_iteration, timelevel);
  free (name);

  /* create the memspace according to chunk dims */
  CACTUS_IOHDF5_ERROR (memspace = H5Screate_simple (sdim, chunk_dims, NULL));

  /* create the (global) filespace and set the hyperslab for the chunk */
  CACTUS_IOHDF5_ERROR (filespace = H5Screate_simple (sdim, file_dims, NULL));
  CACTUS_IOHDF5_ERROR (H5Sselect_hyperslab (filespace, H5S_SELECT_SET,
                       chunk_origin, NULL, chunk_dims, NULL));

  /* the IO processor creates the dataset and adds the common attributes
     when writing its own data, otherwise the dataset is reopened */
  /* if restart from recovery delete an already existing group
     so that it can be created as anew */
  if (h5GH->checkForExistingObjects [index]) {
    CACTUS_IOHDF5_ERROR (H5Eset_auto (NULL, NULL));
    H5Gunlink (iof, datasetname);
    H5Eset_auto (h5GH->printErrorFn, h5GH->printErrorFnArg);
  }
  CACTUS_IOHDF5_ERROR (dataset = H5Dcreate (iof, datasetname,
                       hdf5io_type, filespace, H5P_DEFAULT));
  if (CCTK_MyProc (GH) == 0)
    IOHDF5_AddCommonAttributes_ND (GH, index, slab_geo, 
				   &geom [2*sdim], dataset);

  /* write the data */
  CACTUS_IOHDF5_ERROR (H5Dwrite (dataset, hdf5io_type, memspace,
                       filespace, H5P_DEFAULT, outme));

  /* and close the file dataspace */
  CACTUS_IOHDF5_ERROR (H5Sclose (filespace));

  /* close the dataset and the memspace */
  CACTUS_IOHDF5_ERROR (H5Dclose (dataset));
  CACTUS_IOHDF5_ERROR (H5Sclose (memspace));

  free (chunk_origin);
  free (chunk_dims);
}
#endif /* HAVE_PARALLEL */
#endif /* MPI */