/*@@ @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 #include "cctk.h" #include "util_Table.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_BYTE: retval = FLEXIO_BYTE; 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 */ retval = 0; if (myproc == ioUtilGH->ioproc) { if (IOisValid (file)) { /* first dump the data then add the attributes */ FLEXIO_ERROR (IOwrite (file, IOFlexIO_DataType (request->hdatatype), dim, &dim, buffer)); AddCommonAttributes (GH, request, file); } else { retval = -1; } } 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_LocalMappingByIndex 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, table, 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 */ table = -1; if (request->with_ghostzones) { table = Util_TableCreateFromString ("with_ghostzones = 1"); if (table < 0) { CCTK_WARN (1, "Failed to hyperslab parameter create table from string"); } } mapping = Hyperslab_LocalMappingByIndex (GH, request->vindex, request->hdim, request->direction, request->origin, request->extent, request->downsample, table, NULL, request->hsize_chunk, request->hsize, request->hoffset); if (table >= 0) { Util_TableDestroy (table); } 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: - the GV variable and group name, type, timelevel info - the current iteration number and simulation time - bbox information (only if the GV has a coordinate system associated with it) - chunk information @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 hdim, vdim, coord_system_handle; CCTK_REAL *dtmp; CCTK_INT4 *itmp; CCTK_INT *coord_handles; char *name; const ioGH *ioUtilGH; DECLARE_CCTK_PARAMETERS ioUtilGH = CCTK_GHExtension (GH, "IO"); /* 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", CHAR, strlen (name) + 1, name)); free (name); name = CCTK_GroupNameFromVarI (request->vindex); FLEXIO_ERROR (IOwriteAttribute (file, "groupname", CHAR, strlen (name) + 1, name)); coord_system_handle = -1; if (CCTK_IsFunctionAliased ("Coord_GroupSystem")) { coord_system_handle = Coord_GroupSystem (GH, 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)); /* write bbox attributes if we have coordinate system info */ coord_handles = malloc (request->vdim * sizeof (CCTK_INT)); if (coord_system_handle >= 0 && Util_TableGetIntArray (coord_system_handle, request->vdim, coord_handles, "COORDINATES") >= 0) { dtmp = calloc (3 * request->vdim, sizeof (CCTK_REAL)); for (vdim = 0; vdim < request->vdim; vdim++) { for (hdim = 0; hdim < request->hdim; hdim++) { if (request->direction[hdim*request->hdim + vdim]) { Util_TableGetReal (coord_handles[vdim], &dtmp[hdim + 0*request->vdim], "COMPMIN"); if (Util_TableGetReal (coord_handles[vdim], &dtmp[hdim+2*request->vdim], "DELTA") >= 0) { dtmp[hdim+2*request->vdim] *= request->downsample[hdim]; dtmp[hdim+1*request->vdim] = dtmp[hdim+0*request->vdim]; dtmp[hdim+1*request->vdim] += ((request->extent[hdim] + request->downsample[hdim]-1) / request->downsample[hdim] - 1) * dtmp[hdim+2*request->vdim]; } } } } FLEXIO_ERROR (IOwriteAttribute (file, "origin", FLEXIO_REAL, request->hdim, dtmp + 0*request->vdim)); FLEXIO_ERROR (IOwriteAttribute (file, "min_ext", FLEXIO_REAL, request->hdim, dtmp + 0*request->vdim)); FLEXIO_ERROR (IOwriteAttribute (file, "max_ext", FLEXIO_REAL, request->hdim, dtmp + 1*request->vdim)); FLEXIO_ERROR (IOwriteAttribute (file, "delta", FLEXIO_REAL, request->hdim, dtmp + 2*request->vdim)); #if 0 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)); } #endif free (dtmp); } free (coord_handles); itmp[0] = 0; for (hdim = 0; hdim < request->hdim; hdim++) { itmp[hdim] = request->hsize[hdim]; } FLEXIO_ERROR (IOwriteAttribute (file, "global_size", FLEXIO_INT4, hdim ? hdim : 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", 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