/*@@ @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 #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 H5_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, int check_exisiting_objects, hid_t file); #endif /* H5_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) { 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 DECLARE_CCTK_PARAMETERS /* 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 H5_HAVE_PARALLEL if (ioUtilGH->unchunked) { IOHDF5Util_collectiveDump (GH, vindex, timelevel, request, outme, geom, info->iohdf5_type, check_exisiting_objects,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)); } else { tmpd = NULL; } IOHDF5Util_procDump (GH, vindex, timelevel, request, tmpd, geom, i, info->iohdf5_type, check_exisiting_objects, file); if (tmpd) { 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) { int i, sdim; int myproc; ioGH *ioUtilGH; ioHDF5UtilGH *myGH; hid_t group, dataset, memspace, filespace, plist; char *fullname, *datasetname, *chunkname; hssize_t *chunk_origin; hsize_t *chunk_dims, *file_dims; hsize_t buffersize; DECLARE_CCTK_PARAMETERS ioUtilGH = (ioGH *) CCTK_GHExtension (GH, "IO"); myGH = (ioHDF5UtilGH *) CCTK_GHExtension (GH, "IOHDF5Util"); sdim = request->sdim; myproc = CCTK_MyProc (GH); 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); /* 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 (outme) { /* create the memspace according to chunk dims */ IOHDF5_ERROR (memspace = H5Screate_simple (sdim, chunk_dims, NULL)); } else { memspace = -1; } 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 (plist = H5Pcreate (H5P_DATASET_CREATE)); /* enable compression for chunked dataset if compression was requested */ if (compression_level) { IOHDF5_ERROR (H5Pset_chunk (plist, sdim, chunk_dims)); IOHDF5_ERROR (H5Pset_deflate (plist, compression_level)); } IOHDF5_ERROR (dataset = H5Dcreate (file, datasetname, iohdf5_type, filespace, plist)); IOHDF5_ERROR (H5Pclose (plist)); IOHDF5Util_DumpCommonAttributes (GH, vindex, timelevel, &geom[2*sdim], request, dataset); } else { IOHDF5_ERROR (dataset = H5Dopen (file, datasetname)); } if (memspace >= 0) { /* increase the buffer size if the default isn't sufficient */ IOHDF5_ERROR (plist = H5Pcreate (H5P_DATASET_XFER)); buffersize = H5Dget_storage_size (dataset); if (buffersize > H5Pget_buffer (plist, NULL, NULL)) { IOHDF5_ERROR (H5Pset_buffer (plist, buffersize, NULL, NULL)); } /* write the data */ IOHDF5_ERROR (H5Dwrite (dataset, iohdf5_type, memspace, filespace, plist, outme)); /* close the transfer property list */ IOHDF5_ERROR (H5Pclose (plist)); } /* close the file dataspace */ 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)); } if (memspace >= 0) { /* 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); IOHDF5_ERROR (plist = H5Pcreate (H5P_DATASET_CREATE)); /* enable compression for chunked dataset if compression was requested */ if (compression_level) { IOHDF5_ERROR (H5Pset_chunk (plist, sdim, chunk_dims)); IOHDF5_ERROR (H5Pset_deflate (plist, compression_level)); } /* create the chunk dataset and dump the chunk data */ IOHDF5_ERROR (dataset = H5Dcreate (file, chunkname, iohdf5_type, memspace, plist)); IOHDF5_ERROR (H5Pclose (plist)); 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); } else { dataset = -1; } } /* close the dataset and the memspace */ if (dataset >= 0) { IOHDF5_ERROR (H5Dclose (dataset)); } if (memspace >= 0) { IOHDF5_ERROR (H5Sclose (memspace)); } /* free allocated resources */ free (datasetname); free (fullname); free (chunk_origin); free (chunk_dims); } #ifdef CCTK_MPI #ifdef H5_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, const CCTK_INT *geom, int hdf5io_type, int check_exisiting_objects, hid_t file) { int i, dim; ioHDF5UtilGH *myGH; hid_t dataset, memspace, filespace, plist; char *fullname, *datasetname; hssize_t *chunk_origin; hsize_t *chunk_dims, *file_dims; hsize_t buffersize; DECLARE_CCTK_PARAMETERS myGH = (ioHDF5UtilGH *) CCTK_GHExtension (GH, "IOHDF5Util"); /* get the dimension of the variable */ chunk_origin = (hssize_t *) malloc (request->sdim * sizeof (hssize_t)); chunk_dims = (hsize_t *) malloc (2*request->sdim * sizeof (hsize_t)); file_dims = chunk_dims + request->sdim; /* HDF5 needs it in reverse order */ for (i = 0; i < request->sdim; i++) { chunk_origin[i] = geom[1*request->sdim - 1 - i]; chunk_dims [i] = geom[2*request->sdim - 1 - i]; file_dims [i] = geom[3*request->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 (request->sdim, chunk_dims, NULL)); /* create the (global) filespace and set the hyperslab for the chunk */ IOHDF5_ERROR (filespace = H5Screate_simple (request->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)); } /* enable compression for chunked dataset if compression was requested */ IOHDF5_ERROR (plist = H5Pcreate (H5P_DATASET_CREATE)); if (compression_level) { IOHDF5_ERROR (H5Pset_chunk (plist, request->sdim, chunk_dims)); IOHDF5_ERROR (H5Pset_deflate (plist, compression_level)); } IOHDF5_ERROR (dataset = H5Dcreate (file, datasetname, hdf5io_type, filespace, plist)); IOHDF5_ERROR (H5Pclose (plist)); if (CCTK_MyProc (GH) == 0) { IOHDF5Util_DumpCommonAttributes (GH, vindex, timelevel, &geom[2*request->sdim], request, dataset); } /* increase the buffer size if the default isn't sufficient */ IOHDF5_ERROR (plist = H5Pcreate (H5P_DATASET_XFER)); buffersize = H5Dget_storage_size (dataset); if (buffersize > H5Pget_buffer (plist, NULL, NULL)) { IOHDF5_ERROR (H5Pset_buffer (plist, buffersize, NULL, NULL)); } /* write the data */ IOHDF5_ERROR (H5Dwrite (dataset, hdf5io_type, memspace, filespace, plist, outme)); /* close resources */ IOHDF5_ERROR (H5Pclose (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 /* H5_HAVE_PARALLEL */ #endif /* MPI */