path: root/src/DumpVar.c

/*@@
   @file      DumpVar.c
   @date      Fri May 21 1999
   @author    Thomas Radke
   @desc 
              Routines for writing variables into HDF5 data or checkpoint files.
              These routines are used by other HDF5 IO methods.
   @enddesc
   @version   $Id$
 @@*/


#include <stdlib.h>

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


/* the rcs ID and its dummy function to use it */
static const char *rcsid = "$Header$";
CCTK_FILEVERSION(BetaThorns_IOHDF5Util_DumpVar_c)

/* #define DEBUG_ME 1 */

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


/* 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;


/* prototypes of routines defined in this source file */
static int IOHDF5Util_DumpGS (const cGH *GH,
                              int vindex,
                              int timelevel,
                              int iohdf5_type,
                              const ioHDF5Geo_t *request,  
                              int check_exisiting_objects,
                              hid_t file);
static int IOHDF5Util_DumpGA (const cGH *GH,
                              int vindex,
                              int timelevel,
                              const ioHDF5Geo_t *request,  
                              const dumpInfo *info,
                              int check_exisiting_objects,
                              hid_t file);
static int IOHDF5Util_getDumpData (const cGH *GH,
                                   int vindex,
                                   int timelevel,
                                   const ioHDF5Geo_t *request,  
                                   void **outme,
                                   int *free_outme, 
                                   CCTK_INT *geom);
static void IOHDF5Util_procDump (const cGH *GH,
                                 int vindex,
                                 int timelevel,
                                 const ioHDF5Geo_t *request,  
                                 const void *outme,
                                 const CCTK_INT *geom,
                                 int proc, 
                                 int iohdf5_type,
                                 int check_exisiting_objects,
                                 hid_t file);
#ifdef CCTK_MPI
#ifdef HAVE_PARALLEL
static void IOHDF5Util_collectiveDump (const cGH *GH,
                                       int vindex,
                                       int timelevel,
                                       const ioHDF5Geo_t *request,  
                                       const void *outme,
                                       const CCTK_INT *geom, 
                                       int hdf5io_type,
                                       hid_t file);
#endif /* HAVE_PARALLEL */
#endif /* MPI */



/*@@
   @routine    IOHDF5Util_DumpVar
   @date       16 Apr 1999
   @author     Thomas Radke
   @desc
               Generic dump routine, just calls the type-specific dump routine.
   @enddesc

   @calls      IOHDF5Util_DumpGS
               IOHDF5Util_DumpGA

   @var        GH
   @vdesc      Pointer to CCTK grid hierarchy
   @vtype      const cGH *
   @vio        in
   @endvar
   @var        vindex
   @vdesc      index of the variable to be dumped
   @vtype      int
   @vio        in
   @endvar
   @var        timelevel
   @vdesc      the timelevel to dump
   @vtype      int
   @vio        in
   @endvar
   @var        request
   @vdesc      pointer to the IO request structure
   @vtype      const ioHDF5Geo_t *
   @vio        in
   @endvar
   @var        file
   @vdesc      the HDF5 file handle
   @vtype      hid_t
   @vio        in
   @endvar
   @var        check_exisiting_objects
   @vdesc      flag indicating if we should check for already existing objects
               before these are created anew
               This might happen for existing files reopened after recovery.
   @vtype      int
   @vio        in
   @endvar

   @returntype int
   @returndesc
               -1 if variable has unknown type
               or return code of either
                 @seeroutine IOHDF5Util_DumpGS or
                 @seeroutine IOHDF5Util_DumpGA
   @endreturndesc
@@*/
int IOHDF5Util_DumpVar (const cGH *GH,
                        int vindex,
                        int timelevel,
                        const ioHDF5Geo_t *request,
                        hid_t file,
                        int check_exisiting_objects)
{
  int vtype, gtype;
  pGH *pughGH;
  ioHDF5UtilGH *myGH;
  ioGH *ioUtilGH;
  dumpInfo info;
  int retval;


  /* Get the handle for PUGH, IOUtil, and IOHDF5Util extensions */
  pughGH = PUGH_pGH (GH);
  ioUtilGH = (ioGH *) CCTK_GHExtension (GH, "IO");
  myGH = (ioHDF5UtilGH *) CCTK_GHExtension (GH, "IOHDF5Util");

  vtype = CCTK_VarTypeI (vindex);
  /* downgrade the output precision if requested */
  if (ioUtilGH->out_single)
  {
    /* Do we only want to downgrade generic CCTK types ? */
#ifdef CCTK_REAL4
    if (vtype == CCTK_VARIABLE_REAL)
    {
      vtype = CCTK_VARIABLE_REAL4;
    }
    else if (vtype == CCTK_VARIABLE_COMPLEX)
    {
      vtype = CCTK_VARIABLE_COMPLEX8;
    }
#endif
#ifdef CCTK_INT4
    if (vtype == CCTK_VARIABLE_INT)
    {
      vtype = CCTK_VARIABLE_INT4;
    }
#endif
  }
  info.element_size = CCTK_VarTypeSize (vtype);
#ifdef CCTK_MPI
  info.mpi_type = PUGH_MPIDataType (pughGH, vtype);
#endif
  info.iohdf5_type = IOHDF5Util_DataType (myGH, vtype);

  if (info.element_size <= 0 ||
#ifdef CCTK_MPI
      info.mpi_type == MPI_DATATYPE_NULL ||
#endif
      info.iohdf5_type < 0)
  {
    CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING,
                "Unsupported variable type %d", vtype);
    return (-1);
  }

  /* now branch to the appropriate dump routine */
  gtype = CCTK_GroupTypeFromVarI (vindex);
  if (gtype == CCTK_SCALAR)
  {
    retval = IOHDF5Util_DumpGS (GH, vindex, timelevel, info.iohdf5_type,
                                request, check_exisiting_objects, file);
  }
  else if (gtype == CCTK_ARRAY || gtype == CCTK_GF)
  {
    retval = IOHDF5Util_DumpGA (GH, vindex, timelevel, request, &info,
                                check_exisiting_objects, file);
  }
  else
  {
    CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING,
                "Invalid group type %d", gtype);
    retval = -1;
  }

  return (retval);
}


/*@@
   @routine    IOHDF5Util_DumpGS
   @date       May 21 1999
   @author     Thomas Radke
   @desc
               Dumps a CCTK_SCALAR type variable into a HDF5 file.
   @enddesc 
   @var        GH
   @vdesc      Pointer to CCTK grid hierarchy
   @vtype      const cGH *
   @vio        in
   @endvar
   @var        vindex
   @vdesc      index of the variable to be dumped
   @vtype      int
   @vio        in
   @endvar
   @var        timelevel
   @vdesc      the timelevel to dump
   @vtype      int
   @vio        in
   @endvar
   @var        iohdf5_type
   @vdesc      the HDF5 datatype
   @vtype      int
   @vio        in
   @endvar
   @var        request
   @vdesc      pointer to the IO request structure
   @vtype      const ioHDF5Geo_t *
   @vio        in
   @endvar
   @var        check_exisiting_objects
   @vdesc      flag indicating if we should check for already existing objects
               before these are created anew
               This might happen for existing files reopened after recovery.
   @vtype      int
   @vio        in
   @endvar
   @var        file
   @vdesc      the HDF5 file to dump to
   @vtype      hid_t
   @vio        in
   @endvar

   @returntype int
   @returndesc
               always 0
   @endreturndesc
@@*/
static int IOHDF5Util_DumpGS (const cGH *GH,
                              int vindex,
                              int timelevel,
                              int iohdf5_type,
                              const ioHDF5Geo_t *request,  
                              int check_exisiting_objects,
                              hid_t file)
{
  ioGH *ioUtilGH;
  ioHDF5UtilGH *myGH;
  hid_t dataset;
  char *fullname;
  char *datasetname;
  static CCTK_INT global_shape[1] = {0};


  /* Get the handles for IOHDF5Util and IOUtil extensions */
  ioUtilGH = (ioGH *) CCTK_GHExtension (GH, "IO");
  myGH = (ioHDF5UtilGH *) CCTK_GHExtension (GH, "IOHDF5Util");

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

  fullname = CCTK_FullName (vindex);
  datasetname = (char *) malloc (strlen (fullname) + 80);
  sprintf (datasetname, "%s timelevel %d at iteration %d",
           fullname, timelevel, GH->cctk_iteration);

  /* if restart from recovery delete an already existing dataset
     so that it can be created as anew */
  if (check_exisiting_objects)
  {
    IOHDF5_ERROR (H5Eset_auto (NULL, NULL));
    H5Gunlink (file, datasetname);
    IOHDF5_ERROR (H5Eset_auto (myGH->print_error_fn, myGH->print_error_fn_arg));
  }

  IOHDF5_ERROR (dataset = H5Dcreate (file, datasetname, iohdf5_type,
                                     myGH->scalar_dataspace, H5P_DEFAULT));
  IOHDF5_ERROR (H5Dwrite (dataset, iohdf5_type, H5S_ALL, H5S_ALL, H5P_DEFAULT,
                          CCTK_VarDataPtrI (GH, timelevel, vindex)));
    
  IOHDF5Util_DumpCommonAttributes (GH, vindex, timelevel, global_shape,
                                   request, dataset);

  IOHDF5_ERROR (H5Dclose (dataset));

  free (datasetname);
  free (fullname);

  return (0);
}


/*@@
   @routine    IOHDF5Util_DumpGA
   @date       May 21 1999
   @author     Paul Walker
   @desc
               Dumps a CCTK_GF or CCTK_ARRAY type variable into a HDF5 file.
   @enddesc 

   @calls      IOHDF5Util_getDumpData
               IOHDF5Util_procDump
               IOHDF5Util_collectiveDump

   @var        GH
   @vdesc      Pointer to CCTK grid hierarchy
   @vtype      const cGH *
   @vio        in
   @endvar
   @var        vindex
   @vdesc      index of the variable to be dumped
   @vtype      int
   @vio        in
   @endvar
   @var        timelevel
   @vdesc      the timelevel to dump
   @vtype      int
   @vio        in
   @endvar
   @var        request
   @vdesc      pointer to the IO request structure
   @vtype      const ioHDF5Geo_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      const dumpInfo *
   @vio        in
   @endvar
   @var        check_exisiting_objects
   @vdesc      flag indicating if we should check for already existing objects
               before these are created anew
               This might happen for existing files reopened after recovery.
   @vtype      int
   @vio        in
   @endvar
   @var        file
   @vdesc      the HDF5 file to dump to
   @vtype      hid_t
   @vio        in
   @endvar

   @returntype int
   @returndesc
               0 for success
               or returncode of @seeroutine IOHDF5Util_getDumpData
   @endreturndesc
@@*/
static int IOHDF5Util_DumpGA (const cGH *GH,
                              int vindex,
                              int timelevel,
                              const ioHDF5Geo_t *request,  
                              const dumpInfo *info,
                              int check_exisiting_objects,
                              hid_t file)
{
  DECLARE_CCTK_PARAMETERS
  ioGH *ioUtilGH;
  pGH *pughGH;
  int myproc;
  int nprocs;
  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;
#ifdef CCTK_MPI
  void *tmpd;
  int incoming;
  int outgoing;
  int i, j;
  MPI_Status ms;
#endif


  /* Get the handles for PUGH and IOUtil extensions */
  pughGH   = PUGH_pGH (GH);
  ioUtilGH = (ioGH *) CCTK_GHExtension (GH, "IO");

  /* allocate the geometry buffer */
  geom = (CCTK_INT *) malloc (3 * request->sdim * sizeof (CCTK_INT));
  
  /* Get the pointer to the data we want to output (includes downsampling) */
  retval = IOHDF5Util_getDumpData (GH, vindex, timelevel, request, &outme,
                                   &free_outme, geom);
  if (retval < 0)
  {
    return (retval);
  }

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

#ifdef CCTK_MPI
#ifdef HAVE_PARALLEL
  if (ioUtilGH->unchunked)
  {
    IOHDF5Util_collectiveDump (GH, vindex, timelevel, request, outme, geom,
                               info->iohdf5_type, file);
    if (free_outme)
    {
      free (outme);
    }
    free (geom);
    return (0);
  }
#endif
#endif

  /* Dump data held on IO processor */
  if (myproc == ioUtilGH->ioproc)
  {
    IOHDF5Util_procDump (GH, vindex, timelevel, request, outme, geom,
                         myproc, info->iohdf5_type, check_exisiting_objects,
                         file);
  }
#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++)
    {
      /* receive geometry */
      CACTUS_MPI_ERROR (MPI_Recv (geom, 3*request->sdim, PUGH_MPI_INT, i,
                                  2*i + IOTAGBASE + 1, pughGH->PUGH_COMM_WORLD,
                                  &ms));

      incoming = geom[request->sdim];
      for (j = 1; j < request->sdim; j++)
      {
        incoming *= geom[request->sdim + j];
      }

      /* receive data */
      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));

        IOHDF5Util_procDump (GH, vindex, timelevel, request, tmpd, geom, i,
                             info->iohdf5_type, check_exisiting_objects, file);
        free (tmpd);
      }

    } /* End loop over processors */

  }
  else
  {

    /* Send the geometry and data from the non-IO processors
     * to the IO processors
     */
    outgoing = geom[request->sdim];
    for (i = 1; i < request->sdim; i++)
    {
      outgoing *= geom[request->sdim + i];
    }

    /* send geometry */
    CACTUS_MPI_ERROR (MPI_Send (geom, 3*request->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 DEBUG_ME
      printf ("Processor %d sent %d data points\n", myproc, outgoing);
#endif
    }
  }

  /* wait for every processor to catch up */
  CCTK_Barrier (GH);
#endif

  /* free allocated resources */
  if (free_outme)
  {
    free (outme);
  }
  free (geom);

  return (retval);
}


/**************************************************************************/
/*                             local functions                            */
/**************************************************************************/
static int IOHDF5Util_getDumpData (const cGH *GH,
                                   int vindex,
                                   int timelevel,
                                   const ioHDF5Geo_t *request, 
                                   void **outme,
                                   int *free_outme, 
                                   CCTK_INT *geom)
{
  int i, retval;
  int cctk_output_type;
  ioGH *ioUtilGH;
  int *hsizes, *hsizes_global, *hsizes_offset; /* geometry information */
  char *fullname;
  

  /* get GH extension for IOUtil */
  ioUtilGH = (ioGH *) CCTK_GHExtension (GH, "IO");

  /* determine the datatype for the hyperslab */
  cctk_output_type = CCTK_VarTypeI (vindex);
  if (ioUtilGH->out_single)
  {
#ifdef CCTK_INT4
    if (cctk_output_type == CCTK_VARIABLE_INT)
    {
      cctk_output_type = CCTK_VARIABLE_INT4;
    }
#endif
#ifdef CCTK_REAL4
    if (cctk_output_type == CCTK_VARIABLE_REAL)
    {
      cctk_output_type = CCTK_VARIABLE_REAL4;
    }
    else if (cctk_output_type == CCTK_VARIABLE_COMPLEX)
    {
      cctk_output_type = CCTK_VARIABLE_COMPLEX8;
    }
#endif
  }

  /* FIXME */
  /* allocate array to hold size information of request */
  hsizes        = (int *) malloc (3 * request->sdim * sizeof (int));
  hsizes_global = hsizes + 1*request->sdim;
  hsizes_offset = hsizes + 2*request->sdim;

#if 0
  FIXME: what is this for ???
  /* Get the start indeces for the request from origin. */
  /* 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 (i = 0; i < vdim; i++)
  {
    slabstart[i] = request->origin[i];
  }
  for (i = 0; i < request->sdim; i++)
  {
    slabstart[request->direction[i]] = 0;
  }     
#endif

  retval = NewHyperslab_GetLocalHyperslab (GH, vindex, timelevel, request->sdim,
                                      cctk_output_type, NULL, request->origin,
                                      request->direction, request->length,
                                      request->downsample, outme, free_outme,
                                      hsizes, hsizes_global, hsizes_offset);
  if (retval == 0)
  {
    for (i = 0; i < request->sdim; i++)
    {
      geom[i + 0*request->sdim] = hsizes_offset[i];
      geom[i + 1*request->sdim] = hsizes[i];        
      geom[i + 2*request->sdim] = hsizes_global[i];
      request->actlen[i]  = hsizes_global[i];
    }
  }
  else
  {
    fullname = CCTK_FullName (vindex);
    CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING,
                "Failed to extract hyperslab for variable '%s'", fullname);
    free (fullname);
  }

  free (hsizes);

  return (retval);
}


/*@@
   @routine    IOHDF5Util_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

   @calls      IOHDF5Util_DumpCommonAttributes

   @var        GH
   @vdesc      Pointer to CCTK grid hierarchy
   @vtype      const cGH *
   @vio        in
   @endvar
   @var        vindex
   @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        request
   @vdesc      pointer to the IO request structure
   @vtype      const ioHDF5Geo_t *
   @vio        in
   @endvar
   @var        outme
   @vdesc      pointer to the chunk to dump
   @vtype      const void *
   @vio        in
   @endvar
   @var        geom
   @vdesc      bounds and sizes of the output
   @vtype      const 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        check_exisiting_objects
   @vdesc      flag indicating if we should check for already existing objects
               before these are created anew
               This might happen for existing files reopened after recovery.
   @vtype      int
   @vio        in
   @endvar
   @var        file
   @vdesc      the HDF5 file handle
   @vtype      hid_t
   @vio        in
   @endvar
@@*/
static void IOHDF5Util_procDump (const cGH *GH,
                                 int vindex,
                                 int timelevel,
                                 const ioHDF5Geo_t *request,
                                 const void *outme,
                                 const CCTK_INT *geom,
                                 int proc, 
                                 int iohdf5_type,
                                 int check_exisiting_objects,
                                 hid_t file)
{
  DECLARE_CCTK_PARAMETERS
  int i, sdim;
  int myproc;
  ioGH *ioUtilGH;
  ioHDF5UtilGH *myGH;
  hid_t group, dataset, memspace, filespace, xfer_plist;
  char *fullname, *datasetname, *chunkname;
  hssize_t *chunk_origin;
  hsize_t *chunk_dims, *file_dims;
  hsize_t buffersize;
  int locpoints;


  ioUtilGH = (ioGH *) CCTK_GHExtension (GH, "IO");
  myGH = (ioHDF5UtilGH *) CCTK_GHExtension (GH, "IOHDF5Util");

  sdim = request->sdim;
  myproc = CCTK_MyProc (GH);

  /* Check if there are 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 from the variable's full name,
     the timelevel and the current iteration number */
  fullname = CCTK_FullName (vindex);
  datasetname = (char *) malloc (strlen (fullname) + 80);
  sprintf (datasetname, "%s timelevel %d at iteration %d",
           fullname, timelevel, GH->cctk_iteration);

  /* create the memspace according to chunk dims */
  IOHDF5_ERROR (memspace = H5Screate_simple (sdim, chunk_dims, NULL));
 
  /* if restart from recovery delete an already existing dataset
      so that it can be created as anew */
  if (proc == myproc && check_exisiting_objects)
  {
    IOHDF5_ERROR (H5Eset_auto (NULL, NULL));
    H5Gunlink (file, datasetname);
    IOHDF5_ERROR (H5Eset_auto (myGH->print_error_fn, myGH->print_error_fn_arg));
  }

  if (ioUtilGH->unchunked)
  {
    /* create the (global) filespace and set the hyperslab for the chunk */
    IOHDF5_ERROR (filespace = H5Screate_simple (sdim, file_dims, NULL));
    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 (proc == myproc)
    {
      IOHDF5_ERROR (dataset = H5Dcreate (file, datasetname,
                                         iohdf5_type, filespace, H5P_DEFAULT));
      IOHDF5Util_DumpCommonAttributes (GH, vindex, timelevel, &geom[2*sdim],
                                       request, dataset);
    }
    else
    {
      IOHDF5_ERROR (dataset = H5Dopen (file, datasetname));
    }

    /* increase the buffer size if the default isn't sufficient */
    IOHDF5_ERROR (xfer_plist = H5Pcreate (H5P_DATASET_XFER));
    buffersize = H5Dget_storage_size (dataset);
    if (buffersize > H5Pget_buffer (xfer_plist, NULL, NULL))
    {
      IOHDF5_ERROR (H5Pset_buffer (xfer_plist, buffersize, NULL, NULL));
    }

    /* write the data */
    IOHDF5_ERROR (H5Dwrite (dataset, iohdf5_type, memspace, filespace,
                            xfer_plist, outme));

    /* and close the transfer property list and the file dataspace */
    IOHDF5_ERROR (H5Pclose (xfer_plist));
    IOHDF5_ERROR (H5Sclose (filespace));
  }
  else
  {
    /* the IO processor creates the chunk group and adds common attributes */
    if (proc == myproc)
    {
      IOHDF5_ERROR (group = H5Gcreate (file, datasetname, 0));
      IOHDF5Util_DumpCommonAttributes (GH, vindex, timelevel, &geom[2*sdim],
                                       request, group);
      IOHDF5_ERROR (H5Gclose (group));
    }

    /* now the chunk dataset for each processor is created within the group */
    chunkname = (char *) malloc (strlen (datasetname) + 20);
    sprintf (chunkname, "%s/chunk%d", datasetname, proc - myproc);
    /* create the chunk dataset and dump the chunk data */
    IOHDF5_ERROR (dataset = H5Dcreate (file, chunkname,
                                       iohdf5_type, memspace, H5P_DEFAULT));
    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, myGH->array_dataspace, sdim,
                     IOHDF5_INT);

    free (chunkname);
  }

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

  /* free allocated resources */
  free (datasetname);
  free (fullname);
  free (chunk_origin);
  free (chunk_dims);
}


#ifdef CCTK_MPI
#ifdef HAVE_PARALLEL
/*@@
   @routine    IOHDF5Util_collectiveDump
   @author     Thomas Radke
   @date       May 21 1999
   @desc
               All processors dump their data into an unchunked file.
   @enddesc

   @calls      IOHDF5Util_DumpCommonAttributes

   @var        GH
   @vdesc      Pointer to CCTK grid hierarchy
   @vtype      const cGH *
   @vio        in
   @endvar
   @var        vindex
   @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        request
   @vdesc      pointer to the IO request structure
   @vtype      const ioHDF5Geo_t *
   @vio        in
   @endvar
   @var        outme
   @vdesc      pointer to the chunk to dump
   @vtype      const void *
   @vio        in
   @endvar
   @var        geom
   @vdesc      bounds and sizes of the output
   @vtype      const 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        check_exisiting_objects
   @vdesc      flag indicating if we should check for already existing objects
               before these are created anew
               This might happen for existing files reopened after recovery.
   @vtype      int
   @vio        in
   @endvar
   @var        file
   @vdesc      the HDF5 file handle
   @vtype      hid_t
   @vio        in
   @endvar
@@*/
static void IOHDF5Util_collectiveDump (const cGH *GH,
                                       int vindex,
                                       int timelevel,
                                       const ioHDF5Geo_t *request,  
                                       const void *outme,
                                       CCTK_INT *geom, 
                                       int hdf5io_type,
                                       hid_t file)
{
  DECLARE_CCTK_PARAMETERS
  int i, dim;
  ioHDF5GH *myGH;
  hid_t dataset, memspace, filespace, xfer_plist;
  char *fullname, *datasetname;
  hssize_t *chunk_origin;
  hsize_t *chunk_dims, *file_dims;
  hsize_t buffersize;


  myGH = (ioHDF5GH *) CCTK_GHExtension (GH, "IOHDF5Util");

  /* get the dimension of the variable */
  sdim = request->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 */
  fullname = CCTK_FullName (vindex);
  datasetname = (char *) malloc (strlen (fullname) + 80);
  sprintf (datasetname, "%s timelevel %d at iteration %d",
           fullname, timelevel, GH->cctk_iteration);

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

  /* create the (global) filespace and set the hyperslab for the chunk */
  IOHDF5_ERROR (filespace = H5Screate_simple (sdim, file_dims, NULL));
  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 (check_exisiting_objects)
  {
    IOHDF5_ERROR (H5Eset_auto (NULL, NULL));
    H5Gunlink (file, datasetname);
    IOHDF5_ERROR (H5Eset_auto (myGH->print_error_fn, myGH->print_error_fn_arg));
  }
  IOHDF5_ERROR (dataset = H5Dcreate (file, datasetname, hdf5io_type, filespace,
                                     H5P_DEFAULT));
  if (CCTK_MyProc (GH) == 0)
  {
    IOHDF5Util_DumpCommonAttributes (GH, vindex, timelevel, &geom[2*sdim],
                                     request, dataset);
  }

  /* increase the buffer size if the default isn't sufficient */
  IOHDF5_ERROR (xfer_plist = H5Pcreate (H5P_DATASET_XFER));
  buffersize = H5Dget_storage_size (dataset);
  if (buffersize > H5Pget_buffer (xfer_plist, NULL, NULL))
  {
    IOHDF5_ERROR (H5Pset_buffer (xfer_plist, buffersize, NULL, NULL));
  }

  /* write the data */
  IOHDF5_ERROR (H5Dwrite (dataset, hdf5io_type, memspace,
                       filespace, xfer_plist, outme));

  /* close resources */
  IOHDF5_ERROR (H5Pclose (xfer_plist));
  IOHDF5_ERROR (H5Sclose (filespace));
  IOHDF5_ERROR (H5Dclose (dataset));
  IOHDF5_ERROR (H5Sclose (memspace));

  /* free allocated resources */
  free (datasetname);
  free (fullname);
  free (chunk_origin);
  free (chunk_dims);
}
#endif /* HAVE_PARALLEL */
#endif /* MPI */