/*@@ @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 #include #ifdef SGI #include #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:
  • the variable's groupname
  • the grouptype
  • number of timelevels
  • the current date Note that the datestamp should be turned of if you are byte comparing two output files.
  • simulation time
  • origin
  • bounding box
  • gridspacings (both downsampled and evolution)
  • global grid size
@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;idimsdim;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;idimsdim;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;idimsdim;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;idimdowns[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;idimdirection[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;idimslab_start[idim]; for (idim=0;idimdirection[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;idimdirection[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 */