/*@@ @file DumpVar.c @author Paul Walker, Gabrielle Allen, Tom Goodale, Thomas Radke @desc Routines for writing variables into IOFlexIO data or checkpoint files. @enddesc @version $Id$ @@*/ #include #include #include "cctk.h" #include "cctk_Parameters.h" #include "CactusPUGH/PUGH/src/include/pugh.h" #include "CactusBase/IOUtil/src/ioGH.h" #include "ioFlexGH.h" /* the rcs ID and its dummy funtion to use it */ static const char *rcsid = "$Id$"; CCTK_FILEVERSION(CactusPUGHIO_IOFlexIO_DumpVar_c) /******************************************************************** ******************** Macro Definitions ************************ ********************************************************************/ /* #define IOFLEXIO_DEBUG 1 */ /* tag base for MPI messages (this may break on more than 2000 processors) */ #define MPITAGBASE 20000 /******************************************************************** ******************** Internal Routines ************************ ********************************************************************/ static int WriteGS (const cGH *GH, const ioRequest *request, IOFile file); static int WriteGA (const cGH *GH, const ioRequest *request, IOFile file); static void EachProcDump (const cGH *GH, const ioRequest *request, void *hdata, IOFile file); static void AddChunkAttributes (const cGH *GH, const ioRequest *request, IOFile file); static void AddCommonAttributes (const cGH *GH, const ioRequest *request, IOFile file); #ifdef CCTK_MPI static void ProcDump (const cGH *GH, const ioRequest *request, void *hdata, IOFile file); #endif /*@@ @routine IOFlexIO_DumpVar @date 16 Apr 1999 @author Thomas Radke @desc Generic dump routine, just calls the appropriate dump routine @enddesc @calls @var GH @vdesc pointer to CCTK grid hierarchy @vtype const cGH * @vio in @endvar @var request @vdesc reference to the I/O request description @vtype const ioRequest * @vio in @endvar @var file @vdesc the IEEEIO file to dump to @vtype IOFile @vio in @endvar @returntype int @returndesc return code of either @seeroutine WriteGS or @seeroutine WriteGA or -1 if variable has unknown type @endreturndesc @@*/ int IOFlexIO_DumpVar (const cGH *GH, const ioRequest *request, IOFile file) { int gtype, retval; /* branch to the appropriate dump routine */ gtype = CCTK_GroupTypeFromVarI (request->vindex); if (gtype == CCTK_SCALAR) { retval = WriteGS (GH, request, file); } else if (gtype == CCTK_ARRAY || gtype == CCTK_GF) { retval = WriteGA (GH, request, file); } else { CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING, "Invalid group type %d", gtype); retval = -1; } return (retval); } /*@@ @routine IOFlexIO_DataType @author Thomas Radke @date Mon 11 June 2001 @desc Returns the FlexIO datatype for a given CCTK datatype @enddesc @var cctk_type @vdesc CCTK datatype @vtype int @vio in @endvar @returntype int @returndesc the appropriate FlexIO datatype for success, or -1 otherwise @endreturndesc @@*/ int IOFlexIO_DataType (int cctk_type) { int retval; switch (cctk_type) { case CCTK_VARIABLE_CHAR: retval = FLEXIO_CHAR; break; case CCTK_VARIABLE_INT: retval = FLEXIO_INT; break; case CCTK_VARIABLE_REAL: retval = FLEXIO_REAL; break; #ifdef CCTK_INT1 case CCTK_VARIABLE_INT1: retval = FLEXIO_INT1; break; #endif #ifdef CCTK_INT2 case CCTK_VARIABLE_INT2: retval = FLEXIO_INT2; break; #endif #ifdef CCTK_INT4 case CCTK_VARIABLE_INT4: retval = FLEXIO_INT4; break; #endif #ifdef CCTK_INT8 case CCTK_VARIABLE_INT8: retval = FLEXIO_INT8; break; #endif #ifdef CCTK_REAL4 case CCTK_VARIABLE_REAL4: retval = FLEXIO_REAL4; break; #endif #ifdef CCTK_REAL8 case CCTK_VARIABLE_REAL8: retval = FLEXIO_REAL8; break; #endif #ifdef CCTK_REAL16 case CCTK_VARIABLE_REAL16: retval = FLEXIO_REAL16; break; #endif default: CCTK_VWarn (0, __LINE__, __FILE__, CCTK_THORNSTRING, "Unsupported CCTK variable datatype %d", cctk_type); retval = -1; break; } return (retval); } /******************************************************************** ******************** Internal Routines ************************ ********************************************************************/ /*@@ @routine WriteGS @date 16 Apr 1999 @author Thomas Radke @desc Writes a CCTK_SCALAR type variable into a IEEEIO file. Only the value from processor 0 is written here. If other processors have different values for this scalar variable a warning will be issued. @enddesc @calls AddCommonAttributes @var GH @vdesc Pointer to CCTK grid hierarchy @vtype const cGH * @vio in @endvar @var request @vdesc reference to the I/O request description @vtype const ioRequest * @vio in @endvar @var file @vdesc the IEEEIO file to dump to @vtype IOFile @vio in @endvar @returntype int @returndesc 0 for success, or -1 if file handle is invalid @endreturndesc @@*/ static int WriteGS (const cGH *GH, const ioRequest *request, IOFile file) { int i, myproc, nprocs, hdatatypesize, retval; char *buffer, *fullname; /* const */ int dim = 1; const ioGH *ioUtilGH; ioUtilGH = CCTK_GHExtension (GH, "IO"); myproc = CCTK_MyProc (GH); nprocs = CCTK_nProcs (GH); fullname = CCTK_FullName (request->vindex); hdatatypesize = CCTK_VarTypeSize (request->hdatatype); buffer = calloc (nprocs, hdatatypesize); memcpy (buffer + myproc*hdatatypesize, CCTK_VarDataPtrI (GH, request->timelevel, request->vindex), hdatatypesize); if (nprocs > 1) { i = CCTK_ReductionHandle ("sum"); if (i >= 0) { i = CCTK_ReduceArray (GH, -1, i, nprocs, request->hdatatype, buffer, 1, 1, request->hdatatype, nprocs, buffer); } if (i < 0) { CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING, "WriteGS: Cannot check whether values on different " "processors are the same for grid scalar '%s'", fullname); /* copy this processor's value to the start of buffer */ memcpy (buffer, buffer + myproc*hdatatypesize, hdatatypesize); } else { retval = 0; for (i = 1; i < nprocs; i++) { retval |= memcmp (buffer, buffer + i*hdatatypesize, hdatatypesize); } if (retval) { CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING, "WriteGS: value of grid scalar variable '%s' (timelevel %d)" " differs between processors, only value from processor 0 " "will be written", fullname, request->timelevel); } } } /* only I/O processors write data */ if (myproc != ioUtilGH->ioproc) { retval = 0; } else if (! IOisValid (file)) { retval = -1; } else { /* first dump the data then add the attributes */ FLEXIO_ERROR (IOwrite (file, IOFlexIO_DataType (request->hdatatype), dim, &dim, buffer)); AddCommonAttributes (GH, request, file); retval = 0; } free (buffer); free (fullname); return (retval); } /*@@ @routine WriteGA @date July 1998 @author Paul Walker @desc Dumps a grid array variable into a IEEEIO file. @enddesc @calls Hyperslab_DefineLocalMappingByIndex Hyperslab_FreeMapping Hyperslab_Get WriteData @var GH @vdesc pointer to CCTK grid hierarchy @vtype const cGH * @vio in @endvar @var request @vdesc reference to the I/O request description @vtype const ioRequest * @vio in @endvar @var file @vdesc the IEEEIO file to dump to @vtype IOFile @vio in @endvar @returntype int @returndesc 0 for success, or
-1 if hyperslab mapping couldn't be defined, or
-2 if hyperslab couldn't be extracted @endreturndesc @@*/ static int WriteGA (const cGH *GH, const ioRequest *request, IOFile file) { int i, mapping, hdatasize, retval; void *hdata; char *fullname; const ioGH *ioUtilGH; #ifdef CCTK_MPI int myproc, nprocs; int *hsize_chunk; void *tmpd; MPI_Comm comm; MPI_Status ms; MPI_Datatype mpitype; #endif DECLARE_CCTK_PARAMETERS /* define the hyperslab mapping */ mapping = Hyperslab_DefineLocalMappingByIndex (GH, request->vindex, request->hdim, request->direction, request->origin, request->extent, request->downsample, -1, NULL, request->hsize_chunk, request->hsize, request->hoffset); if (mapping < 0) { fullname = CCTK_FullName (request->vindex); CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING, "Failed to define hyperslab mapping for variable '%s'", fullname); free (fullname); return (-1); } /* calculate the size of the hyperslab */ request->hsize_chunk[request->hdim] = 1; for (i = 0; i < request->hdim; i++) { request->hsize_chunk[request->hdim] *= request->hsize_chunk[i]; } /* get the hyperslab */ hdatasize = CCTK_VarTypeSize (request->hdatatype); hdata = request->hsize_chunk[request->hdim] > 0 ? malloc (request->hsize_chunk[request->hdim] * hdatasize) : NULL; retval = Hyperslab_Get (GH, mapping, -1, request->vindex, request->timelevel, request->hdatatype, hdata); /* release the mapping structure */ Hyperslab_FreeMapping (mapping); if (retval) { fullname = CCTK_FullName (request->vindex); CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING, "Failed to extract hyperslab for variable '%s'", fullname); free (fullname); if (hdata) { free (hdata); } return (-2); } /* dump data held on I/O processor */ ioUtilGH = CCTK_GHExtension (GH, "IO"); if (ioUtilGH->ioproc_every == 1) { EachProcDump (GH, request, hdata, file); } #ifdef CCTK_MPI else { myproc = CCTK_MyProc (GH); nprocs = CCTK_nProcs (GH); comm = PUGH_pGH (GH)->PUGH_COMM_WORLD; mpitype = PUGH_MPIDataType (PUGH_pGH (GH), request->hdatatype); if (myproc == ioUtilGH->ioproc) { if (ioUtilGH->unchunked) { /* need to copy CCTK_INT[] to int[] */ hsize_chunk = malloc (request->hdim * sizeof (int)); for (i = 0; i < request->hdim; i++) { hsize_chunk[i] = request->hsize[i]; } IOreserveChunk (file, IOFlexIO_DataType (request->hdatatype), request->hdim, hsize_chunk); free (hsize_chunk); } /* write the data first then the attributes */ ProcDump (GH, request, hdata, file); /* delay adding attributes for unchunked files until all chunks were written (otherwise attributes get lost !) */ if (! ioUtilGH->unchunked) { AddCommonAttributes (GH, request, file); } /* dump data from all other processors */ for (i = myproc+1; i < myproc+ioUtilGH->ioproc_every && i < nprocs; i++) { /* receive geometry (this assumes the geometry arrays to be contiguous starting at hoffset */ CACTUS_MPI_ERROR (MPI_Recv (request->hoffset, 3*request->hdim + 1, PUGH_MPI_INT, i, 2*i + MPITAGBASE + 1, comm, &ms)); /* receive data */ tmpd = NULL; if (request->hsize_chunk[request->hdim] > 0) { tmpd = malloc (request->hsize_chunk[request->hdim] * hdatasize); CACTUS_MPI_ERROR (MPI_Recv (tmpd, request->hsize_chunk[request->hdim], mpitype, i, 2*i + MPITAGBASE, comm, &ms)); } /* write data */ ProcDump (GH, request, tmpd, file); if (tmpd) { free (tmpd); } } /* end loop over processors */ /* now add the attributes for unchunked files */ if (ioUtilGH->unchunked) { AddCommonAttributes (GH, request, file); } } else { /* send geometry (this assumes the geometry arrays to be contiguous starting at hoffset) */ CACTUS_MPI_ERROR (MPI_Send (request->hoffset, 3*request->hdim + 1, PUGH_MPI_INT, ioUtilGH->ioproc, 2*myproc + MPITAGBASE + 1, comm)); /* send data */ if (request->hsize_chunk[request->hdim] > 0) { CACTUS_MPI_ERROR (MPI_Send (hdata, request->hsize_chunk[request->hdim], mpitype, ioUtilGH->ioproc, 2*myproc + MPITAGBASE, comm)); #ifdef IOFLEXIO_DEBUG printf ("Processor %d sent %d data points\n", myproc, request->hsize_chunk[request->hdim]); #endif } } } /* wait for every processor to catch up */ CCTK_Barrier (GH); #endif /* MPI */ /* free allocated resources */ if (hdata) { free (hdata); } return (retval); } /*@@ @routine AddCommonAttributes @date July 1998 @author Paul Walker @desc Add "Common" attributes, these are the GF name, the current date, simulation time, origin, bounding box, and gridspacings (both downsampled and evolution). @enddesc @calls @history @hdate Wed Sep 2 10:15:22 1998 @hauthor Tom Goodale @hdesc Abstracted WRITE_ATTRIBUTE to merge in Szu-Wen's Panda changes. @hdate Apr 16 1999 @hauthor Thomas Radke @hdesc Added attributes groupname, grouptype, ntimelevels, and current timelevel to be stored @hdate May 02 1999 @hauthor Thomas Radke @hdesc Made chunked attribute nioprocs common so that it can be found by the recombiner even for unchunked datasets (eg. SCALARs) @hdate May 05 1999 @hauthor Thomas Radke @hdesc Added "unchunked" attribute to distinguish between chunked/unchunked output files @endhistory @@*/ static void AddCommonAttributes (const cGH *GH, const ioRequest *request, IOFile file) { int i; CCTK_REAL *dtmp; CCTK_INT4 *itmp; char *name; ioGH *ioUtilGH; char coord_system_name[20]; DECLARE_CCTK_PARAMETERS /* allocate at least one CCTK_INT4 if hdim is 0 */ itmp = malloc ((request->hdim+1) * sizeof (CCTK_INT4)); name = CCTK_FullName (request->vindex); FLEXIO_ERROR (IOwriteAttribute (file, "name", FLEXIO_CHAR, strlen (name) + 1, name)); free (name); name = CCTK_GroupNameFromVarI (request->vindex); FLEXIO_ERROR (IOwriteAttribute (file, "groupname", FLEXIO_CHAR, strlen (name) + 1, name)); free (name); itmp[0] = CCTK_GroupTypeFromVarI (request->vindex); FLEXIO_ERROR (IOwriteAttribute (file, "grouptype", FLEXIO_INT4, 1, itmp)); itmp[0] = CCTK_MaxTimeLevelsVI (request->vindex); FLEXIO_ERROR (IOwriteAttribute (file, "ntimelevels", FLEXIO_INT4, 1, itmp)); itmp[0] = request->timelevel; FLEXIO_ERROR (IOwriteAttribute (file, "timelevel", FLEXIO_INT4, 1, itmp)); FLEXIO_ERROR (IOwriteAttribute (file, "time", FLEXIO_REAL, 1,&GH->cctk_time)); /* attributes describing the underlying grid These are only stored for CCTK_GF variables if they are associated with coordinates. */ /* FIXME: This is hardcoded for cartesian coordinate systems. A better solution would be to be able to query the coordinate system which is associated with the variable. */ ioUtilGH = CCTK_GHExtension (GH, "IO"); sprintf (coord_system_name, "cart%dd", request->hdim); if (CCTK_GroupTypeFromVarI (request->vindex) == CCTK_GF && CCTK_CoordSystemHandle (coord_system_name) >= 0) { dtmp = malloc (3 * request->hdim * sizeof (CCTK_REAL)); for (i = 0; i < request->hdim; i++) { CCTK_CoordRange (GH, &dtmp[i], &dtmp[i + request->hdim], i + 1, NULL, coord_system_name); dtmp[i + 2*request->hdim] = GH->cctk_delta_space[i] * ioUtilGH->downsample[i]; } FLEXIO_ERROR (IOwriteAttribute (file, "origin", FLEXIO_REAL, request->hdim, dtmp)); FLEXIO_ERROR (IOwriteAttribute (file, "min_ext", FLEXIO_REAL, request->hdim, dtmp)); FLEXIO_ERROR (IOwriteAttribute (file, "max_ext", FLEXIO_REAL, request->hdim, dtmp + request->hdim)); FLEXIO_ERROR (IOwriteAttribute (file, "delta", FLEXIO_REAL, request->hdim, dtmp + 2*request->hdim)); if (ioUtilGH->downsample[0] > 1 || ioUtilGH->downsample[1] > 1 || ioUtilGH->downsample[2] > 1) { FLEXIO_ERROR (IOwriteAttribute (file, "evolution_delta", FLEXIO_REAL, request->hdim, GH->cctk_delta_space)); } free (dtmp); } if (request->hdim) { for (i = 0; i < request->hdim; i++) { itmp[i] = request->hsize[i]; } FLEXIO_ERROR (IOwriteAttribute (file, "global_size", FLEXIO_INT4, request->hdim, itmp)); } else { itmp[0] = 0; FLEXIO_ERROR (IOwriteAttribute (file, "global_size", FLEXIO_INT4, 1, itmp)); } itmp[0] = CCTK_nProcs (GH); FLEXIO_ERROR (IOwriteAttribute (file, "nprocs", FLEXIO_INT4, 1, itmp)); itmp[0] = ioUtilGH->ioproc_every; FLEXIO_ERROR (IOwriteAttribute (file, "ioproc_every", FLEXIO_INT4, 1, itmp)); itmp[0] = ioUtilGH->unchunked; FLEXIO_ERROR (IOwriteAttribute (file, "unchunked", FLEXIO_INT4, 1, itmp)); itmp[0] = GH->cctk_iteration; FLEXIO_ERROR (IOwriteAttribute (file, "iteration", FLEXIO_INT4, 1, itmp)); free (itmp); } /*@@ @routine AddChunkAttributes @author Paul Walker @date Feb 1997 @desc This routine adds chunk attributes to a data set. That is, is adds attributes to the chunks of data passed to the I/O processors. @enddesc @history @hauthor Gabrielle Allen @hdate Jul 10 1998 @hdesc Added the name of the grid function @hdate Wed Sep 2 10:08:31 1998 @hauthor Tom Goodale @hdesc Abstracted WRITE_ATTRIBUTE to merge in Szu-Wen's Panda calls. @endhistory @@*/ static void AddChunkAttributes (const cGH *GH, const ioRequest *request, IOFile file) { int i; char *fullname; CCTK_INT4 *itmp; /* there is nothing to do for a serial run */ if (CCTK_nProcs (GH) == 1) { return; } itmp = malloc (2 * request->hdim * sizeof (CCTK_INT4)); /* copy from CCTK_INT[] to CCTK_INT4[] */ for (i = 0; i < request->hdim; i++) { itmp[0*request->hdim + i] = request->hoffset[i]; itmp[1*request->hdim + i] = request->hsize[i]; } FLEXIO_ERROR (IOwriteAttribute (file, "chunk_origin", FLEXIO_INT4, request->hdim, &itmp[0])); FLEXIO_ERROR (IOwriteAttribute (file, "subchunk_lb", FLEXIO_INT4, request->hdim, &itmp[0])); FLEXIO_ERROR (IOwriteAttribute (file, "global_size", FLEXIO_INT4, request->hdim, &itmp[request->hdim])); itmp[0] = GH->cctk_iteration; FLEXIO_ERROR (IOwriteAttribute (file, "chunk_dataset", FLEXIO_INT4, 1, itmp)); fullname = CCTK_FullName (request->vindex); FLEXIO_ERROR (IOwriteAttribute (file, "name", FLEXIO_CHAR, strlen (fullname)+1, fullname)); free (fullname); free (itmp); } /*@@ @routine EachProcDump @author Paul Walker @date Feb 1997 @desc Dump data from each processor @enddesc @history @hauthor Gabrielle Allen @hdate Oct 5 1998 @hdesc Made into subroutine @endhistory @@*/ static void EachProcDump (const cGH *GH, const ioRequest *request, void *hdata, IOFile file) { int i, *chunk_dims; /* copy from CCTK_INT[] to int[] */ chunk_dims = malloc (request->hdim * sizeof (int)); for (i = 0; i < request->hdim; i++) { chunk_dims[i] = request->hsize_chunk[i]; } /* dump the data */ FLEXIO_ERROR (IOwrite (file, IOFlexIO_DataType (request->hdatatype), request->hdim, chunk_dims, hdata)); /* add attributes for global space */ AddCommonAttributes (GH, request, file); /* add chunk attributes */ AddChunkAttributes (GH, request, file); free (chunk_dims); } /*@@ @routine ProcDump @author Paul Walker @date Feb 1997 @desc Dump data @enddesc @history @hauthor Gabrielle Allen @hdate Oct 5 1998 @hdesc Made into subroutine @endhistory @@*/ #ifdef CCTK_MPI static void ProcDump (const cGH *GH, const ioRequest *request, void *hdata, IOFile file) { ioGH *ioUtilGH; int i, *chunk_dims, *chunk_origin; /* copy from CCTK_INT[] to int[] */ chunk_origin = malloc (2*request->hdim * sizeof (int)); chunk_dims = chunk_origin + request->hdim; for (i = 0; i < request->hdim; i++) { chunk_origin[i] = request->hoffset[i]; chunk_dims[i] = request->hsize_chunk[i]; } ioUtilGH = CCTK_GHExtension (GH, "IO"); if (ioUtilGH->unchunked) { /* unchunked data */ FLEXIO_ERROR (IOwriteChunk (file, chunk_dims, chunk_origin, hdata)); } else { /* chunked data */ FLEXIO_ERROR (IOwrite (file, IOFlexIO_DataType (request->hdatatype), request->hdim, chunk_dims, hdata)); /* write chunk attributes */ AddChunkAttributes (GH, request, file); } free (chunk_origin); } #endif