diff options
Diffstat (limited to 'Carpet/CarpetIOHDF5/src/iohdf5.cc')
-rw-r--r-- | Carpet/CarpetIOHDF5/src/iohdf5.cc | 1522 |
1 files changed, 738 insertions, 784 deletions
diff --git a/Carpet/CarpetIOHDF5/src/iohdf5.cc b/Carpet/CarpetIOHDF5/src/iohdf5.cc index aead94071..90bc2fa76 100644 --- a/Carpet/CarpetIOHDF5/src/iohdf5.cc +++ b/Carpet/CarpetIOHDF5/src/iohdf5.cc @@ -17,7 +17,7 @@ #include "cctk_Parameters.h" extern "C" { - static const char* rcsid = "$Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetIOHDF5/src/iohdf5.cc,v 1.37 2004/07/06 15:33:01 schnetter Exp $"; + static const char* rcsid = "$Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetIOHDF5/src/iohdf5.cc,v 1.38 2004/07/07 11:01:05 tradke Exp $"; CCTK_FILEVERSION(Carpet_CarpetIOHDF5_iohdf5_cc); } @@ -38,386 +38,387 @@ extern "C" { namespace CarpetIOHDF5 { - using namespace std; - using namespace Carpet; +using namespace std; +using namespace Carpet; - // Variable definitions - int GHExtension; - int IOMethod; - vector<bool> do_truncate; // [var] - vector<vector<vector<int> > > last_output; // [ml][rl][var] +// Variable definitions +vector<bool> do_truncate; // [var] +vector<vector<vector<int> > > last_output; // [ml][rl][var] - int CarpetIOHDF5Startup () - { - DECLARE_CCTK_PARAMETERS; - - int ierr; - - CCTK_RegisterBanner ("AMR 3D HDF5 I/O provided by CarpetIOHDF5"); - - GHExtension = CCTK_RegisterGHExtension ("CarpetIOHDF5"); - CCTK_RegisterGHExtensionSetupGH (GHExtension, SetupGH); +void CarpetIOHDF5Startup (void) +{ + DECLARE_CCTK_PARAMETERS - IOMethod = CCTK_RegisterIOMethod ("IOHDF5"); - CCTK_RegisterIOMethodOutputGH (IOMethod, OutputGH); - CCTK_RegisterIOMethodOutputVarAs (IOMethod, OutputVarAs); - CCTK_RegisterIOMethodTimeToOutput (IOMethod, TimeToOutput); - CCTK_RegisterIOMethodTriggerOutput (IOMethod, TriggerOutput); + CCTK_RegisterBanner ("AMR 3D HDF5 I/O provided by CarpetIOHDF5"); - /* initial I/O parameter check */ - int numvars = CCTK_NumVars (); - vector<bool> flags(numvars); + int GHExtension = CCTK_RegisterGHExtension ("CarpetIOHDF5"); + CCTK_RegisterGHExtensionSetupGH (GHExtension, SetupGH); - if (CCTK_TraverseString (out3D_vars, SetFlag, &flags,CCTK_GROUP_OR_VAR) < 0) - { - CCTK_VWarn (strict_io_parameter_check ? 0 : 1, - __LINE__, __FILE__, CCTK_THORNSTRING, - "error while parsing parameter 'IOHDF5::out3D_vars'"); - } + int IOMethod = CCTK_RegisterIOMethod ("IOHDF5"); + CCTK_RegisterIOMethodOutputGH (IOMethod, OutputGH); + CCTK_RegisterIOMethodOutputVarAs (IOMethod, OutputVarAs); + CCTK_RegisterIOMethodTimeToOutput (IOMethod, TimeToOutput); + CCTK_RegisterIOMethodTriggerOutput (IOMethod, TriggerOutput); + /* initial I/O parameter check */ + int numvars = CCTK_NumVars (); + vector<bool> flags(numvars); - // Christian's Recovery routine - if ( !(CCTK_Equals(recover,"no")) ) { - ierr = IOUtil_RegisterRecover ("CarpetIOHDF5 recovery", CarpetIOHDF5_Recover); - assert (! ierr); - } else { - // Erik's Recovery routine - ierr = IOUtil_RegisterRecover ("CarpetIOHDF5", Recover); - assert (! ierr); - } - + if (CCTK_TraverseString (out3D_vars, SetFlag, &flags,CCTK_GROUP_OR_VAR) < 0) + { + CCTK_VWarn (strict_io_parameter_check ? 0 : 1, + __LINE__, __FILE__, CCTK_THORNSTRING, + "error while parsing parameter 'IOHDF5::out3D_vars'"); + } - return 0; +#if 0 + // Christian's Recovery routine + if ( !(CCTK_Equals(recover,"no")) ) { + ierr = IOUtil_RegisterRecover ("CarpetIOHDF5 recovery", Recover); + assert (! ierr); + } else { + // Erik's Recovery routine + ierr = IOUtil_RegisterRecover ("CarpetIOHDF5", Recover); + assert (! ierr); + } +#else + if (IOUtil_RegisterRecover ("CarpetIOHDF5 recovery", Recover) < 0) + { + CCTK_WARN (1, "Failed to register IOFlexIO recovery routine"); } +#endif +} - void CarpetIOHDF5Init (CCTK_ARGUMENTS) - { - DECLARE_CCTK_ARGUMENTS; +int CarpetIOHDF5Init (const cGH* const cctkGH) +{ + DECLARE_CCTK_ARGUMENTS; - *this_iteration = -1; - *next_output_iteration = 0; - *next_output_time = cctk_time; - } + *this_iteration = -1; + *next_output_iteration = 0; + *next_output_time = cctk_time; + return (0); +} - void* SetupGH (tFleshConfig* const fc, - const int convLevel, cGH* const cctkGH) - { - DECLARE_CCTK_PARAMETERS; - CarpetIOHDF5GH* myGH; +void* SetupGH (tFleshConfig* const fc, + const int convLevel, cGH* const cctkGH) +{ + DECLARE_CCTK_PARAMETERS; - const void *dummy; - dummy = &fc; - dummy = &convLevel; - dummy = &cctkGH; - dummy = &dummy; + CarpetIOHDF5GH* myGH; - // Truncate all files if this is not a restart - do_truncate.resize(CCTK_NumVars(), true); + const void *dummy; + dummy = &fc; + dummy = &convLevel; + dummy = &cctkGH; + dummy = &dummy; - // No iterations have yet been output - last_output.resize(mglevels); - for (int ml=0; ml<mglevels; ++ml) { - last_output.at(ml).resize(maxreflevels); - for (int rl=0; rl<maxreflevels; ++rl) { - last_output.at(ml).at(rl).resize(CCTK_NumVars(), INT_MIN); - } + // Truncate all files if this is not a restart + do_truncate.resize(CCTK_NumVars(), true); + + // No iterations have yet been output + last_output.resize(mglevels); + for (int ml=0; ml<mglevels; ++ml) { + last_output.at(ml).resize(maxreflevels); + for (int rl=0; rl<maxreflevels; ++rl) { + last_output.at(ml).at(rl).resize(CCTK_NumVars(), INT_MIN); } + } - // We register only once, ergo we get only one handle. We store - // that statically, so there is no need to pass anything to - // Cactus. + // We register only once, ergo we get only one handle. We store + // that statically, so there is no need to pass anything to + // Cactus. - /* allocate a new GH extension structure */ + /* allocate a new GH extension structure */ - CCTK_INT numvars = CCTK_NumVars (); + CCTK_INT numvars = CCTK_NumVars (); - myGH = (CarpetIOHDF5GH*) malloc (sizeof (CarpetIOHDF5GH)); - myGH->out_last = (int *) malloc (numvars * sizeof (int)); - myGH->requests = (ioRequest **) calloc (numvars, sizeof (ioRequest *)); - myGH->cp_filename_list = (char **) calloc (abs (checkpoint_keep), sizeof (char *)); - myGH->cp_filename_index = 0; - myGH->out_vars = strdup (""); - myGH->out_every_default = out_every - 1; + myGH = (CarpetIOHDF5GH*) malloc (sizeof (CarpetIOHDF5GH)); + myGH->out_last = (int *) malloc (numvars * sizeof (int)); + myGH->requests = (ioRequest **) calloc (numvars, sizeof (ioRequest *)); + myGH->cp_filename_list = (char **) calloc (abs (checkpoint_keep), sizeof (char *)); + myGH->cp_filename_index = 0; + myGH->out_vars = strdup (""); + myGH->out_every_default = out_every - 1; - for (int i = 0; i < numvars; i++) - { - myGH->out_last [i] = -1; - } + for (int i = 0; i < numvars; i++) + { + myGH->out_last [i] = -1; + } - myGH->open_output_files = NULL; + myGH->open_output_files = NULL; - // Now set hdf5 complex datatypes - // Stolen from Thomas Radke - HDF5_ERROR (myGH->HDF5_COMPLEX = - H5Tcreate (H5T_COMPOUND, sizeof (CCTK_COMPLEX))); - HDF5_ERROR (H5Tinsert (myGH->HDF5_COMPLEX, "real", - offsetof (CCTK_COMPLEX, Re), HDF5_REAL)); - HDF5_ERROR (H5Tinsert (myGH->HDF5_COMPLEX, "imag", - offsetof (CCTK_COMPLEX, Im), HDF5_REAL)); + // Now set hdf5 complex datatypes + // Stolen from Thomas Radke + HDF5_ERROR (myGH->HDF5_COMPLEX = + H5Tcreate (H5T_COMPOUND, sizeof (CCTK_COMPLEX))); + HDF5_ERROR (H5Tinsert (myGH->HDF5_COMPLEX, "real", + offsetof (CCTK_COMPLEX, Re), HDF5_REAL)); + HDF5_ERROR (H5Tinsert (myGH->HDF5_COMPLEX, "imag", + offsetof (CCTK_COMPLEX, Im), HDF5_REAL)); #ifdef CCTK_REAL4 - HDF5_ERROR (myGH->HDF5_COMPLEX8 = - H5Tcreate (H5T_COMPOUND, sizeof (CCTK_COMPLEX8))); - HDF5_ERROR (H5Tinsert (myGH->HDF5_COMPLEX8, "real", - offsetof (CCTK_COMPLEX8, Re), HDF5_REAL4)); - HDF5_ERROR (H5Tinsert (myGH->HDF5_COMPLEX8, "imag", - offsetof (CCTK_COMPLEX8, Im), HDF5_REAL4)); + HDF5_ERROR (myGH->HDF5_COMPLEX8 = + H5Tcreate (H5T_COMPOUND, sizeof (CCTK_COMPLEX8))); + HDF5_ERROR (H5Tinsert (myGH->HDF5_COMPLEX8, "real", + offsetof (CCTK_COMPLEX8, Re), HDF5_REAL4)); + HDF5_ERROR (H5Tinsert (myGH->HDF5_COMPLEX8, "imag", + offsetof (CCTK_COMPLEX8, Im), HDF5_REAL4)); #endif #ifdef CCTK_REAL8 - HDF5_ERROR (myGH->HDF5_COMPLEX16 = - H5Tcreate (H5T_COMPOUND, sizeof (CCTK_COMPLEX16))); - HDF5_ERROR (H5Tinsert (myGH->HDF5_COMPLEX16, "real", - offsetof (CCTK_COMPLEX16, Re), HDF5_REAL8)); - HDF5_ERROR (H5Tinsert (myGH->HDF5_COMPLEX16, "imag", - offsetof (CCTK_COMPLEX16, Im), HDF5_REAL8)); + HDF5_ERROR (myGH->HDF5_COMPLEX16 = + H5Tcreate (H5T_COMPOUND, sizeof (CCTK_COMPLEX16))); + HDF5_ERROR (H5Tinsert (myGH->HDF5_COMPLEX16, "real", + offsetof (CCTK_COMPLEX16, Re), HDF5_REAL8)); + HDF5_ERROR (H5Tinsert (myGH->HDF5_COMPLEX16, "imag", + offsetof (CCTK_COMPLEX16, Im), HDF5_REAL8)); #endif #ifdef CCTK_REAL16 - HDF5_ERROR (myGH->HDF5_COMPLEX32 = - H5Tcreate (H5T_COMPOUND, sizeof (CCTK_COMPLEX32))); - HDF5_ERROR (H5Tinsert (myGH->HDF5_COMPLEX32, "real", - offsetof (CCTK_COMPLEX32, Re), HDF5_REAL16)); - HDF5_ERROR (H5Tinsert (myGH->HDF5_COMPLEX32, "imag", - offsetof (CCTK_COMPLEX32, Im), HDF5_REAL16)); + HDF5_ERROR (myGH->HDF5_COMPLEX32 = + H5Tcreate (H5T_COMPOUND, sizeof (CCTK_COMPLEX32))); + HDF5_ERROR (H5Tinsert (myGH->HDF5_COMPLEX32, "real", + offsetof (CCTK_COMPLEX32, Re), HDF5_REAL16)); + HDF5_ERROR (H5Tinsert (myGH->HDF5_COMPLEX32, "imag", + offsetof (CCTK_COMPLEX32, Im), HDF5_REAL16)); #endif - return (myGH); - } + return (myGH); +} - int OutputGH (const cGH* const cctkGH) { - for (int vindex=0; vindex<CCTK_NumVars(); ++vindex) { - if (TimeToOutput(cctkGH, vindex)) { - TriggerOutput(cctkGH, vindex); - } +int OutputGH (const cGH* const cctkGH) { + for (int vindex=0; vindex<CCTK_NumVars(); ++vindex) { + if (TimeToOutput(cctkGH, vindex)) { + TriggerOutput(cctkGH, vindex); } - return 0; } + return 0; +} - int OutputVarAs (const cGH* const cctkGH, const char* const varname, - const char* const alias) { - DECLARE_CCTK_ARGUMENTS; - DECLARE_CCTK_PARAMETERS; - - herr_t herr; - - const int n = CCTK_VarIndex(varname); - assert (n>=0 && n<CCTK_NumVars()); - const int group = CCTK_GroupIndexFromVarI (n); - assert (group>=0 && group<(int)Carpet::arrdata.size()); - const int n0 = CCTK_FirstVarIndexI(group); - assert (n0>=0 && n0<CCTK_NumVars()); - const int var = n - n0; - assert (var>=0 && var<CCTK_NumVars()); +int OutputVarAs (const cGH* const cctkGH, const char* const varname, + const char* const alias) { + DECLARE_CCTK_ARGUMENTS; + DECLARE_CCTK_PARAMETERS; - // Check for storage - if (! CCTK_QueryGroupStorageI(cctkGH, group)) { - CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING, - "Cannot output variable \"%s\" because it has no storage", - varname); - return 0; - } + herr_t herr; - const int grouptype = CCTK_GroupTypeI(group); - switch (grouptype) { - case CCTK_SCALAR: - case CCTK_ARRAY: - assert (do_global_mode); - break; - case CCTK_GF: - /* do nothing */ - break; - default: - assert (0); - } + const int n = CCTK_VarIndex(varname); + assert (n>=0 && n<CCTK_NumVars()); + const int group = CCTK_GroupIndexFromVarI (n); + assert (group>=0 && group<(int)Carpet::arrdata.size()); + const int n0 = CCTK_FirstVarIndexI(group); + assert (n0>=0 && n0<CCTK_NumVars()); + const int var = n - n0; + assert (var>=0 && var<CCTK_NumVars()); - /* get the default I/O request for this variable */ - ioRequest* request = IOUtil_DefaultIORequest (cctkGH, n, 1); - - // Get grid hierarchy extentsion from IOUtil - const ioGH * const iogh = (const ioGH *)CCTK_GHExtension (cctkGH, "IO"); - assert (iogh); + // Check for storage + if (! CCTK_QueryGroupStorageI(cctkGH, group)) { + CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING, + "Cannot output variable \"%s\" because it has no storage", + varname); + return 0; + } - // Create the output directory - const char* myoutdir = out3D_dir; - if (CCTK_EQUALS(myoutdir, "")) { - myoutdir = out_dir; - } - if (CCTK_MyProc(cctkGH)==0) { - CCTK_CreateDirectory (0755, myoutdir); - } + const int grouptype = CCTK_GroupTypeI(group); + switch (grouptype) { + case CCTK_SCALAR: + case CCTK_ARRAY: + assert (do_global_mode); + break; + case CCTK_GF: + /* do nothing */ + break; + default: + assert (0); + } - // Invent a file name - ostringstream filenamebuf; - filenamebuf << myoutdir << "/" << alias << out3D_extension; - string filenamestr = filenamebuf.str(); - const char * const filename = filenamestr.c_str(); + /* get the default I/O request for this variable */ + ioRequest* request = IOUtil_DefaultIORequest (cctkGH, n, 1); - hid_t writer = -1; + // Get grid hierarchy extentsion from IOUtil + const ioGH * const iogh = (const ioGH *)CCTK_GHExtension (cctkGH, "IO"); + assert (iogh); - // Write the file only on the root processor - if (CCTK_MyProc(cctkGH)==0) { + // Create the output directory + const char* myoutdir = out3D_dir; + if (CCTK_EQUALS(myoutdir, "")) { + myoutdir = out_dir; + } + if (CCTK_MyProc(cctkGH)==0) { + CCTK_CreateDirectory (0755, myoutdir); + } - // If this is the first time, then create and truncate the file - if (do_truncate.at(n)) { - struct stat fileinfo; - if (! iogh->recovered || stat(filename, &fileinfo)!=0) { - writer = H5Fcreate (filename, H5F_ACC_TRUNC, H5P_DEFAULT, H5P_DEFAULT); - assert (writer>=0); - herr = H5Fclose (writer); - assert (!herr); - writer = -1; - } + // Invent a file name + ostringstream filenamebuf; + filenamebuf << myoutdir << "/" << alias << out3D_extension; + string filenamestr = filenamebuf.str(); + const char * const filename = filenamestr.c_str(); + + hid_t writer = -1; + + // Write the file only on the root processor + if (CCTK_MyProc(cctkGH)==0) { + + // If this is the first time, then create and truncate the file + if (do_truncate.at(n)) { + struct stat fileinfo; + if (! iogh->recovered || stat(filename, &fileinfo)!=0) { + writer = H5Fcreate (filename, H5F_ACC_TRUNC, H5P_DEFAULT, H5P_DEFAULT); + assert (writer>=0); + herr = H5Fclose (writer); + assert (!herr); + writer = -1; } - - // Open the file - writer = H5Fopen (filename, H5F_ACC_RDWR, H5P_DEFAULT); - assert (writer>=0); - } - if(verbose) { - CCTK_VInfo (CCTK_THORNSTRING, - "Writing variable %s on refinement level %d",varname,reflevel); - } + // Open the file + writer = H5Fopen (filename, H5F_ACC_RDWR, H5P_DEFAULT); + assert (writer>=0); - WriteVar(cctkGH,writer,request,0); + } - // Close the file - if (CCTK_MyProc(cctkGH)==0) { - herr = H5Fclose (writer); - assert (!herr); - writer = -1; - } + if(verbose) { + CCTK_VInfo (CCTK_THORNSTRING, + "Writing variable %s on refinement level %d",varname,reflevel); + } - // Don't truncate again - do_truncate.at(n) = false; + WriteVar(cctkGH,writer,request,0); - return 0; + // Close the file + if (CCTK_MyProc(cctkGH)==0) { + herr = H5Fclose (writer); + assert (!herr); + writer = -1; } - int WriteVar (const cGH* const cctkGH, const hid_t writer, const ioRequest* request, - const int called_from_checkpoint) { + // Don't truncate again + do_truncate.at(n) = false; - DECLARE_CCTK_ARGUMENTS; - DECLARE_CCTK_PARAMETERS; + return 0; +} - herr_t herr=0; +int WriteVar (const cGH* const cctkGH, const hid_t writer, const ioRequest* request, + const int called_from_checkpoint) { - void * h5data=NULL; + DECLARE_CCTK_ARGUMENTS; + DECLARE_CCTK_PARAMETERS; - const int n = request->vindex; - assert (n>=0 && n<CCTK_NumVars()); - const char * varname = CCTK_FullName(n); - const int group = CCTK_GroupIndexFromVarI (n); - assert (group>=0 && group<(int)Carpet::arrdata.size()); - const int n0 = CCTK_FirstVarIndexI(group); - assert (n0>=0 && n0<CCTK_NumVars()); - const int var = n - n0; - assert (var>=0 && var<CCTK_NumVars()); - int tl = 0; + herr_t herr=0; - const int gpdim = CCTK_GroupDimI(group); + void * h5data=NULL; - // Check for storage - if (! CCTK_QueryGroupStorageI(cctkGH, group)) { - CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING, - "Cannot output variable \"%s\" because it has no storage", - varname); - return 0; - } + const int n = request->vindex; + assert (n>=0 && n<CCTK_NumVars()); + const char * varname = CCTK_FullName(n); + const int group = CCTK_GroupIndexFromVarI (n); + assert (group>=0 && group<(int)Carpet::arrdata.size()); + const int n0 = CCTK_FirstVarIndexI(group); + assert (n0>=0 && n0<CCTK_NumVars()); + const int var = n - n0; + assert (var>=0 && var<CCTK_NumVars()); + int tl = 0; - const int grouptype = CCTK_GroupTypeI(group); - switch (grouptype) { - case CCTK_SCALAR: - case CCTK_ARRAY: - assert (do_global_mode); - break; - case CCTK_GF: - /* do nothing */ - break; - default: - assert (0); - } - const int rl = grouptype==CCTK_GF ? reflevel : 0; + const int gpdim = CCTK_GroupDimI(group); + + // Check for storage + if (! CCTK_QueryGroupStorageI(cctkGH, group)) { + CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING, + "Cannot output variable \"%s\" because it has no storage", + varname); + return 0; + } - cGroup cgdata; - int ierr = CCTK_GroupData(group,&cgdata); - assert(ierr==0); + const int grouptype = CCTK_GroupTypeI(group); + switch (grouptype) { + case CCTK_SCALAR: + case CCTK_ARRAY: + assert (do_global_mode); + break; + case CCTK_GF: + /* do nothing */ + break; + default: + assert (0); + } + const int rl = grouptype==CCTK_GF ? reflevel : 0; + + cGroup cgdata; + int ierr = CCTK_GroupData(group,&cgdata); + assert(ierr==0); - // Select memory (source) and file (destination) datatypes - int cctkDataType = CCTK_VarTypeI(n); - const hid_t memdatatype = h5DataType(cctkGH, cctkDataType); - assert(memdatatype >= 0); - if (out_single_precision && ! called_from_checkpoint) + // Select memory (source) and file (destination) datatypes + int cctkDataType = CCTK_VarTypeI(n); + const hid_t memdatatype = h5DataType(cctkGH, cctkDataType); + assert(memdatatype >= 0); + if (out_single_precision && ! called_from_checkpoint) + { + if (cctkDataType == CCTK_VARIABLE_REAL) { - if (cctkDataType == CCTK_VARIABLE_REAL) - { - cctkDataType = CCTK_VARIABLE_REAL4; - } - else if (cctkDataType == CCTK_VARIABLE_COMPLEX) - { - cctkDataType = CCTK_VARIABLE_COMPLEX8; - } + cctkDataType = CCTK_VARIABLE_REAL4; + } + else if (cctkDataType == CCTK_VARIABLE_COMPLEX) + { + cctkDataType = CCTK_VARIABLE_COMPLEX8; + } #ifdef CCTK_INT2 - else if (cctkDataType == CCTK_VARIABLE_INT) - { - cctkDataType = CCTK_VARIABLE_INT2; - } -#endif + else if (cctkDataType == CCTK_VARIABLE_INT) + { + cctkDataType = CCTK_VARIABLE_INT2; } - const hid_t filedatatype = h5DataType(cctkGH, cctkDataType); - assert(filedatatype >= 0); +#endif + } + const hid_t filedatatype = h5DataType(cctkGH, cctkDataType); + assert(filedatatype >= 0); - // let's get the correct Carpet time level (which is the (-1) * Cactus timelevel): - tl = - request->timelevel; + // let's get the correct Carpet time level (which is the (-1) * Cactus timelevel): + tl = - request->timelevel; - // Traverse all components - BEGIN_MAP_LOOP(cctkGH, grouptype) { - BEGIN_COMPONENT_LOOP(cctkGH, grouptype) { + // Traverse all components + BEGIN_MAP_LOOP(cctkGH, grouptype) { + BEGIN_COMPONENT_LOOP(cctkGH, grouptype) { - const ggf<dim>* ff = 0; + const ggf<dim>* ff = 0; - ff = (ggf<dim>*)arrdata.at(group).at(Carpet::map).data.at(var); + ff = (ggf<dim>*)arrdata.at(group).at(Carpet::map).data.at(var); - const gdata<dim>* const data = (*ff) (tl, rl, component, mglevel); + const gdata<dim>* const data = (*ff) (tl, rl, component, mglevel); - // Make temporary copy on processor 0 - const ibbox ext = data->extent(); + // Make temporary copy on processor 0 + const ibbox ext = data->extent(); // vect<int,dim> lo = ext.lower(); // vect<int,dim> hi = ext.upper(); // vect<int,dim> str = ext.stride(); - gdata<dim>* const tmp = data->make_typed (n); - tmp->allocate (ext, 0); + gdata<dim>* const tmp = data->make_typed (n); + tmp->allocate (ext, 0); - if ( !((cgdata.disttype == CCTK_DISTRIB_CONSTANT) && - (arrdata.at(group).at(Carpet::map).hh->processors.at(reflevel).at(component)!=0))) { + if ( !((cgdata.disttype == CCTK_DISTRIB_CONSTANT) && + (arrdata.at(group).at(Carpet::map).hh->processors.at(reflevel).at(component)!=0))) { - if (cgdata.disttype == CCTK_DISTRIB_CONSTANT) { - assert(grouptype == CCTK_ARRAY || grouptype == CCTK_SCALAR); - if(component!=0) goto skip; - h5data = CCTK_VarDataPtrI(cctkGH,tl,n); - } else { - for (comm_state<dim> state; !state.done(); state.step()) { - tmp->copy_from (state, data, ext); - } + if (cgdata.disttype == CCTK_DISTRIB_CONSTANT) { + assert(grouptype == CCTK_ARRAY || grouptype == CCTK_SCALAR); + if(component!=0) goto skip; + h5data = CCTK_VarDataPtrI(cctkGH,tl,n); + } else { + for (comm_state<dim> state; !state.done(); state.step()) { + tmp->copy_from (state, data, ext); } - // Write data - if (CCTK_MyProc(cctkGH)==0) { - - int ldim=0; - if ( grouptype==CCTK_SCALAR ) { - ldim = 1; - } else { - ldim = gpdim; - } + } + // Write data + if (CCTK_MyProc(cctkGH)==0) { + + int ldim=0; + if ( grouptype==CCTK_SCALAR ) { + ldim = 1; + } else { + ldim = gpdim; + } // hsize_t shape[ldim]; @@ -575,629 +576,582 @@ namespace CarpetIOHDF5 { } - herr = H5Sclose (dataspace); - assert (!herr); - - } // if on root processor - } // if ! CCTK_DISTRIB_BLAH - - skip: - // Delete temporary copy - delete tmp; - - } END_COMPONENT_LOOP; - } END_MAP_LOOP; + herr = H5Sclose (dataspace); + assert (!herr); - return 0; - } + } // if on root processor + } // if ! CCTK_DISTRIB_BLAH + skip: + // Delete temporary copy + delete tmp; + } END_COMPONENT_LOOP; + } END_MAP_LOOP; - int TimeToOutput (const cGH* const cctkGH, const int vindex) - { - DECLARE_CCTK_ARGUMENTS; - DECLARE_CCTK_PARAMETERS; - - assert (vindex>=0 && vindex<CCTK_NumVars()); + return 0; +} - const int grouptype = CCTK_GroupTypeFromVarI(vindex); - switch (grouptype) { - case CCTK_SCALAR: - case CCTK_ARRAY: - if (! do_global_mode) return 0; - break; - case CCTK_GF: - // do nothing - break; - default: - assert (0); - } +int TimeToOutput (const cGH* const cctkGH, const int vindex) +{ + DECLARE_CCTK_ARGUMENTS; + DECLARE_CCTK_PARAMETERS; + assert (vindex>=0 && vindex<CCTK_NumVars()); - // check whether to output at this iteration - bool output_this_iteration; - const char* myoutcriterion = out3D_criterion; - if (CCTK_EQUALS(myoutcriterion, "default")) { - myoutcriterion = out_criterion; - } + const int grouptype = CCTK_GroupTypeFromVarI(vindex); + switch (grouptype) { + case CCTK_SCALAR: + case CCTK_ARRAY: + if (! do_global_mode) return 0; + break; + case CCTK_GF: + // do nothing + break; + default: + assert (0); + } - if (CCTK_EQUALS (myoutcriterion, "never")) { - // Never output - output_this_iteration = false; - } else if (CCTK_EQUALS (myoutcriterion, "iteration")) { + // check whether to output at this iteration + bool output_this_iteration; - int myoutevery = out3D_every; - if (myoutevery == -2) { - myoutevery = out_every; - } - if (myoutevery <= 0) { - // output is disabled - output_this_iteration = false; - } else if (cctk_iteration == *this_iteration) { - // we already decided to output this iteration - output_this_iteration = true; - } else if (cctk_iteration >= *next_output_iteration) { - // it is time for the next output - output_this_iteration = true; - *next_output_iteration = cctk_iteration + myoutevery; - *this_iteration = cctk_iteration; - } else { - // we want no output at this iteration - output_this_iteration = false; - } + const char* myoutcriterion = out3D_criterion; + if (CCTK_EQUALS(myoutcriterion, "default")) { + myoutcriterion = out_criterion; + } - } else if (CCTK_EQUALS (myoutcriterion, "divisor")) { + if (CCTK_EQUALS (myoutcriterion, "never")) { - int myoutevery = out3D_every; - if (myoutevery == -2) { - myoutevery = out_every; - } - if (myoutevery <= 0) { - // output is disabled - output_this_iteration = false; - } else if ((cctk_iteration % myoutevery) == 0) { - // we already decided to output this iteration - output_this_iteration = true; - } else { - // we want no output at this iteration - output_this_iteration = false; - } + // Never output + output_this_iteration = false; - } else if (CCTK_EQUALS (myoutcriterion, "time")) { - - CCTK_REAL myoutdt = out3D_dt; - if (myoutdt == -2) { - myoutdt = out_dt; - } - if (myoutdt < 0) { - // output is disabled - output_this_iteration = false; - } else if (myoutdt == 0) { - // output all iterations - output_this_iteration = true; - } else if (cctk_iteration == *this_iteration) { - // we already decided to output this iteration - output_this_iteration = true; - } else if (cctk_time / cctk_delta_time - >= *next_output_time / cctk_delta_time - 1.0e-12) { - // it is time for the next output - output_this_iteration = true; - *next_output_time = cctk_time + myoutdt; - *this_iteration = cctk_iteration; - } else { - // we want no output at this iteration - output_this_iteration = false; - } + } else if (CCTK_EQUALS (myoutcriterion, "iteration")) { + int myoutevery = out3D_every; + if (myoutevery == -2) { + myoutevery = out_every; + } + if (myoutevery <= 0) { + // output is disabled + output_this_iteration = false; + } else if (cctk_iteration == *this_iteration) { + // we already decided to output this iteration + output_this_iteration = true; + } else if (cctk_iteration >= *next_output_iteration) { + // it is time for the next output + output_this_iteration = true; + *next_output_iteration = cctk_iteration + myoutevery; + *this_iteration = cctk_iteration; } else { + // we want no output at this iteration + output_this_iteration = false; + } - assert (0); - - } // select output criterion - - if (! output_this_iteration) return 0; - - + } else if (CCTK_EQUALS (myoutcriterion, "divisor")) { - if (! CheckForVariable(cctkGH, out3D_vars, vindex)) { - // This variable should not be output - return 0; + int myoutevery = out3D_every; + if (myoutevery == -2) { + myoutevery = out_every; + } + if (myoutevery <= 0) { + // output is disabled + output_this_iteration = false; + } else if ((cctk_iteration % myoutevery) == 0) { + // we already decided to output this iteration + output_this_iteration = true; + } else { + // we want no output at this iteration + output_this_iteration = false; } + } else if (CCTK_EQUALS (myoutcriterion, "time")) { - - if (last_output.at(mglevel).at(reflevel).at(vindex) == cctk_iteration) { - // Has already been output during this iteration - char* varname = CCTK_FullName(vindex); - CCTK_VWarn (5, __LINE__, __FILE__, CCTK_THORNSTRING, - "Skipping output for variable \"%s\", because this variable " - "has already been output during the current iteration -- " - "probably via a trigger during the analysis stage", - varname); - free (varname); - return 0; + CCTK_REAL myoutdt = out3D_dt; + if (myoutdt == -2) { + myoutdt = out_dt; + } + if (myoutdt < 0) { + // output is disabled + output_this_iteration = false; + } else if (myoutdt == 0) { + // output all iterations + output_this_iteration = true; + } else if (cctk_iteration == *this_iteration) { + // we already decided to output this iteration + output_this_iteration = true; + } else if (cctk_time / cctk_delta_time + >= *next_output_time / cctk_delta_time - 1.0e-12) { + // it is time for the next output + output_this_iteration = true; + *next_output_time = cctk_time + myoutdt; + *this_iteration = cctk_iteration; + } else { + // we want no output at this iteration + output_this_iteration = false; } - assert (last_output.at(mglevel).at(reflevel).at(vindex) < cctk_iteration); + } else { - // Should be output during this iteration - return 1; - } + assert (0); + } // select output criterion + if (! output_this_iteration) return 0; - int TriggerOutput (const cGH* const cctkGH, const int vindex) { - assert (vindex>=0 && vindex<CCTK_NumVars()); - char* varname = CCTK_FullName(vindex); - const int retval = OutputVarAs (cctkGH, varname, CCTK_VarName(vindex)); - free (varname); - last_output.at(mglevel).at(reflevel).at(vindex) = cctkGH->cctk_iteration; - - return retval; + if (! CheckForVariable(out3D_vars, vindex)) { + // This variable should not be output + return 0; } - - int InputGH (const cGH* const cctkGH) { - DECLARE_CCTK_PARAMETERS; - int retval = 0; - for (int vindex=0; vindex<CCTK_NumVars(); ++vindex) { - if (CheckForVariable(cctkGH, in3D_vars, vindex)) { - char* varname = CCTK_FullName(vindex); - retval = InputVarAs (cctkGH, varname, CCTK_VarName(vindex)); - free (varname); - if (retval != 0) return retval; - } - } - return retval; + if (last_output.at(mglevel).at(reflevel).at(vindex) == cctk_iteration) { + // Has already been output during this iteration + char* varname = CCTK_FullName(vindex); + CCTK_VWarn (5, __LINE__, __FILE__, CCTK_THORNSTRING, + "Skipping output for variable \"%s\", because this variable " + "has already been output during the current iteration -- " + "probably via a trigger during the analysis stage", + varname); + free (varname); + return 0; } + assert (last_output.at(mglevel).at(reflevel).at(vindex) < cctk_iteration); + // Should be output during this iteration + return 1; +} - int ReadVar (const cGH* const cctkGH, const hid_t reader, const char* const varname, - const hid_t dataset, vector<ibset> ®ions_read, - const int called_from_recovery) { - DECLARE_CCTK_PARAMETERS; - const void *dummy; - dummy = &reader; - dummy = &dummy; +int TriggerOutput (const cGH* const cctkGH, const int vindex) { + assert (vindex>=0 && vindex<CCTK_NumVars()); - const int n = CCTK_VarIndex(varname); - assert (n>=0 && n<CCTK_NumVars()); - const int group = CCTK_GroupIndexFromVarI (n); - assert (group>=0 && group<(int)Carpet::arrdata.size()); - const int n0 = CCTK_FirstVarIndexI(group); - assert (n0>=0 && n0<CCTK_NumVars()); - const int var = n - n0; - assert (var>=0 && var<CCTK_NumVars()); - int tl = 0; + char* varname = CCTK_FullName(vindex); + const int retval = OutputVarAs (cctkGH, varname, CCTK_VarName(vindex)); + free (varname); - herr_t herr = 1; - bool did_read_something = false; + last_output.at(mglevel).at(reflevel).at(vindex) = cctkGH->cctk_iteration; - // Stuff needed for Recovery - int recovery_tl = -1; - int recovery_mglevel = -1; - int recovery_rl = -1; - int recovery_comp = -1; + return retval; +} - void * h5data = NULL; - // Check for storage - if (! CCTK_QueryGroupStorageI(cctkGH, group)) { - CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING, - "Cannot input variable \"%s\" because it has no storage", - varname); - return 0; - } - const int grouptype = CCTK_GroupTypeI(group); - const int rl = grouptype==CCTK_GF ? reflevel : 0; - const int gpdim = CCTK_GroupDimI(group); +int ReadVar (const cGH* const cctkGH, const char* const varname, + const hid_t dataset, vector<ibset> ®ions_read, + const int called_from_recovery) +{ + DECLARE_CCTK_PARAMETERS; - int amr_level; - int amr_origin[dim]; - int amr_dims[dim]; - - // We need to initialize amr_dims to 0 - for(int i=0;i<dim;i++) amr_dims[i] = 0; + const int n = CCTK_VarIndex(varname); + assert (n>=0 && n<CCTK_NumVars()); + const int group = CCTK_GroupIndexFromVarI (n); + assert (group>=0 && group<(int)Carpet::arrdata.size()); + const int n0 = CCTK_FirstVarIndexI(group); + assert (n0>=0 && n0<CCTK_NumVars()); + const int var = n - n0; + assert (var>=0 && var<CCTK_NumVars()); + int tl = 0; - if (CCTK_MyProc(cctkGH)==0) { + bool did_read_something = false; - // get dataset dimensions - const hid_t dataspace = H5Dget_space(dataset); - assert (dataspace>=0); - int rank = (int) H5Sget_simple_extent_ndims(dataspace); - vector<hsize_t> shape(rank); - int rank2 = H5Sget_simple_extent_dims (dataspace, &shape[0], NULL); - herr = H5Sclose(dataspace); - assert(!herr); - assert (rank2 == rank); + // Stuff needed for Recovery - if(grouptype == CCTK_SCALAR) { - assert (gpdim+1 == rank); - } else { - assert (gpdim == rank); - } + void *h5data = NULL; - int datalength=1; + // Check for storage + if (! CCTK_QueryGroupStorageI(cctkGH, group)) + { + CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING, + "Cannot input variable \"%s\" because it has no storage", + varname); + return 0; + } - for(int i=0;i<rank;i++) { - datalength=datalength*shape[i]; - amr_dims[i]=shape[rank-i-1]; - } + const int grouptype = CCTK_GroupTypeI(group); + if ((grouptype == CCTK_SCALAR || grouptype == CCTK_ARRAY) && reflevel) + { + return 0; + } - if(grouptype == CCTK_ARRAY) { + const int gpdim = CCTK_GroupDimI(group); - for(int i=rank;i<dim;i++) { - amr_dims[i]=1; - } - } - const int cctkDataType = CCTK_VarTypeI(n); - const hid_t datatype = h5DataType(cctkGH,cctkDataType); + int intbuffer[2 + 2*dim]; + int &group_timelevel = intbuffer[0], + &amr_level = intbuffer[1]; + int *amr_origin = &intbuffer[2], + *amr_dims = &intbuffer[2+dim]; - //cout << "datalength: " << datalength << " rank: " << rank << "\n"; - //cout << shape[0] << " " << shape[1] << " " << shape[2] << "\n"; + if (CCTK_MyProc(cctkGH)==0) + { + // get dataset dimensions + hid_t dataspace; + HDF5_ERROR (dataspace = H5Dget_space (dataset)); + int rank = (int) H5Sget_simple_extent_ndims (dataspace); + assert (0 < rank && rank <= dim); + assert ((grouptype == CCTK_SCALAR ? gpdim+1 : gpdim) == rank); + vector<hsize_t> shape(rank); + HDF5_ERROR (H5Sget_simple_extent_dims (dataspace, &shape[0], NULL)); + HDF5_ERROR (H5Sclose (dataspace)); + + for (int i = 0; i < dim; i++) + { + amr_dims[i] = 1; + amr_origin[i] = 0; + } + int datalength = 1; + for (int i = 0; i < rank; i++) + { + datalength *= shape[i]; + amr_dims[i] = shape[rank-i-1]; + } - // to do: read in an allocate with correct datatype + const int cctkDataType = CCTK_VarTypeI(n); + const hid_t datatype = h5DataType(cctkGH,cctkDataType); - h5data = (void*) malloc( CCTK_VarTypeSize(cctkDataType) *datalength); + //cout << "datalength: " << datalength << " rank: " << rank << "\n"; + //cout << shape[0] << " " << shape[1] << " " << shape[2] << "\n"; - herr = H5Dread(dataset,datatype,H5S_ALL, H5S_ALL, H5P_DEFAULT,(void*)h5data); - assert(!herr); + // to do: read in an allocate with correct datatype + h5data = malloc (CCTK_VarTypeSize (cctkDataType) * datalength); + HDF5_ERROR (H5Dread (dataset, datatype, H5S_ALL, H5S_ALL, H5P_DEFAULT, + h5data)); - herr = ReadAttribute(dataset,"level",amr_level); - assert(herr>=0); - ReadAttribute(dataset,"iorigin",amr_origin,dim); - assert(herr>=0); + ReadAttribute (dataset, "level", amr_level); + ReadAttribute (dataset, "iorigin", amr_origin, rank); - if(called_from_recovery) { - recovery_rl = amr_level; - ReadAttribute(dataset,"carpet_component",recovery_comp); - ReadAttribute(dataset,"group_timelevel",recovery_tl); - ReadAttribute(dataset,"carpet_mglevel",recovery_mglevel); - } + if(called_from_recovery) + { + ReadAttribute(dataset,"group_timelevel", group_timelevel); + } + } // MyProc == 0 - } // MyProc == 0 + MPI_Bcast (intbuffer, sizeof (intbuffer) / sizeof (*intbuffer), MPI_INT, 0, dist::comm); + if (verbose) cout << "amr_level: " << amr_level << " reflevel: " << + reflevel << endl; + if (amr_level == reflevel) + { + // Traverse all components on all levels + BEGIN_MAP_LOOP(cctkGH, grouptype) + { + BEGIN_COMPONENT_LOOP(cctkGH, grouptype) + { + did_read_something = true; - if(called_from_recovery) { - MPI_Bcast (&recovery_rl, 1, MPI_INT, 0, dist::comm); - MPI_Bcast (&recovery_tl, 1, MPI_INT, 0, dist::comm); - MPI_Bcast (&recovery_mglevel, 1, MPI_INT, 0, dist::comm); - MPI_Bcast (&recovery_comp, 1, MPI_INT, 0, dist::comm); - } + ggf<dim>* ff = 0; - MPI_Bcast (&amr_level, 1, MPI_INT, 0, dist::comm); - MPI_Bcast (amr_origin, dim, MPI_INT, 0, dist::comm); - MPI_Bcast (amr_dims, dim, MPI_INT, 0, dist::comm); + assert (var < (int)arrdata.at(group).at(Carpet::map).data.size()); + ff = (ggf<dim>*)arrdata.at(group).at(Carpet::map).data.at(var); - if ((grouptype == CCTK_SCALAR || grouptype == CCTK_ARRAY) && reflevel != 0) { - return 0; - } + if(called_from_recovery) tl = group_timelevel; - if (verbose) cout << "amr_level: " << amr_level << " reflevel: " << reflevel << endl; + gdata<dim>* const data = (*ff) (tl, reflevel, component, mglevel); - if (amr_level == rl) { + // Create temporary data storage on processor 0 + vect<int,dim> str = vect<int,dim>(maxreflevelfact/reflevelfact); - // Traverse all components on all levels - BEGIN_MAP_LOOP(cctkGH, grouptype) { - BEGIN_COMPONENT_LOOP(cctkGH, grouptype) { + if(grouptype == CCTK_SCALAR || grouptype == CCTK_ARRAY) + str = vect<int,dim> (1); - did_read_something = true; + vect<int,dim> lb = vect<int,dim>::ref(amr_origin) * str; + vect<int,dim> ub + = lb + (vect<int,dim>::ref(amr_dims) - 1) * str; - ggf<dim>* ff = 0; + gdata<dim>* const tmp = data->make_typed (n); - assert (var < (int)arrdata.at(group).at(Carpet::map).data.size()); - ff = (ggf<dim>*)arrdata.at(group).at(Carpet::map).data.at(var); - if(called_from_recovery) tl = recovery_tl; + cGroup cgdata; + int ierr = CCTK_GroupData(group,&cgdata); + assert(ierr==0); + //cout << "lb_before: " << lb << endl; + //cout << "ub_before: " << ub << endl; + if (cgdata.disttype == CCTK_DISTRIB_CONSTANT) { + if (verbose) cout << "CCTK_DISTRIB_CONSTANT: " << varname << endl; + assert(grouptype == CCTK_ARRAY || grouptype == CCTK_SCALAR); + if (grouptype == CCTK_SCALAR) { + lb[0] = arrdata.at(group).at(Carpet::map).hh->processors.at(reflevel).at(component); + ub[0] = arrdata.at(group).at(Carpet::map).hh->processors.at(reflevel).at(component); + for(int i=1;i<dim;i++) { + lb[i]=0; + ub[i]=0; + } + } else { + const int newlb = lb[gpdim-1] + + (ub[gpdim-1]-lb[gpdim-1]+1)* + (arrdata.at(group).at(Carpet::map).hh->processors.at(reflevel).at(component)); + const int newub = ub[gpdim-1] + + (ub[gpdim-1]-lb[gpdim-1]+1)* + (arrdata.at(group).at(Carpet::map).hh->processors.at(reflevel).at(component)); + lb[gpdim-1] = newlb; + ub[gpdim-1] = newub; + } + if (verbose) cout << "lb: " << lb << endl; + if (verbose) cout << "ub: " << ub << endl; + } + const bbox<int,dim> ext(lb,ub,str); + if (verbose) cout << "ext: " << ext << endl; - gdata<dim>* const data = (*ff) (tl, rl, component, mglevel); + if (CCTK_MyProc(cctkGH)==0) { + tmp->allocate (ext, 0, h5data); + } else { + tmp->allocate (ext, 0); + } - // Create temporary data storage on processor 0 - vect<int,dim> str - = vect<int,dim>(maxreflevelfact/reflevelfact); + // Initialise with what is found in the file -- this does + // not guarantee that everything is initialised. + const bbox<int,dim> overlap = tmp->extent() & data->extent(); + regions_read.at(Carpet::map) |= overlap; + + if (verbose) { + cout << "working on component: " << component << endl; + cout << "tmp->extent " << tmp->extent() << endl; + cout << "data->extent " << data->extent() << endl; + cout << "overlap " << overlap << endl; + cout << "-----------------------------------------------------" << endl; + } + MPI_Barrier(MPI_COMM_WORLD); - if(grouptype == CCTK_SCALAR || grouptype == CCTK_ARRAY) - str = vect<int,dim> (1); + // Copy into grid function + for (comm_state<dim> state; !state.done(); state.step()) { + data->copy_from (state, tmp, overlap); + } - vect<int,dim> lb = vect<int,dim>::ref(amr_origin) * str; - vect<int,dim> ub - = lb + (vect<int,dim>::ref(amr_dims) - 1) * str; - gdata<dim>* const tmp = data->make_typed (n); + // Delete temporary copy + delete tmp; - cGroup cgdata; - int ierr = CCTK_GroupData(group,&cgdata); - assert(ierr==0); - //cout << "lb_before: " << lb << endl; - //cout << "ub_before: " << ub << endl; - if (cgdata.disttype == CCTK_DISTRIB_CONSTANT) { - if (verbose) cout << "CCTK_DISTRIB_CONSTANT: " << varname << endl; - assert(grouptype == CCTK_ARRAY || grouptype == CCTK_SCALAR); - if (grouptype == CCTK_SCALAR) { - lb[0] = arrdata.at(group).at(Carpet::map).hh->processors.at(rl).at(component); - ub[0] = arrdata.at(group).at(Carpet::map).hh->processors.at(rl).at(component); - for(int i=1;i<dim;i++) { - lb[i]=0; - ub[i]=0; - } - } else { - const int newlb = lb[gpdim-1] + - (ub[gpdim-1]-lb[gpdim-1]+1)* - (arrdata.at(group).at(Carpet::map).hh->processors.at(rl).at(component)); - const int newub = ub[gpdim-1] + - (ub[gpdim-1]-lb[gpdim-1]+1)* - (arrdata.at(group).at(Carpet::map).hh->processors.at(rl).at(component)); - lb[gpdim-1] = newlb; - ub[gpdim-1] = newub; - } - if (verbose) cout << "lb: " << lb << endl; - if (verbose) cout << "ub: " << ub << endl; - } - const bbox<int,dim> ext(lb,ub,str); + } END_COMPONENT_LOOP; - if (verbose) cout << "ext: " << ext << endl; + if (called_from_recovery) { + arrdata.at(group).at(Carpet::map).tt->set_time(reflevel,mglevel, + (CCTK_REAL) ((cctkGH->cctk_time - cctk_initial_time) + / (delta_time * mglevelfact)) ); + } - if (CCTK_MyProc(cctkGH)==0) { - tmp->allocate (ext, 0, h5data); - } else { - tmp->allocate (ext, 0); - } - // Initialise with what is found in the file -- this does - // not guarantee that everything is initialised. - const bbox<int,dim> overlap = tmp->extent() & data->extent(); - regions_read.at(Carpet::map) |= overlap; + } END_MAP_LOOP; - if (verbose) { - cout << "working on component: " << component << endl; - cout << "tmp->extent " << tmp->extent() << endl; - cout << "data->extent " << data->extent() << endl; - cout << "overlap " << overlap << endl; - cout << "-----------------------------------------------------" << endl; - } - MPI_Barrier(MPI_COMM_WORLD); + } // if amr_level == reflevel - // Copy into grid function - for (comm_state<dim> state; !state.done(); state.step()) { - data->copy_from (state, tmp, overlap); - } + free (h5data); + return did_read_something; +} - // Delete temporary copy - delete tmp; - } END_COMPONENT_LOOP; +static int InputVarAs (const cGH* const cctkGH, const char* const varname, + const char* const alias) +{ + DECLARE_CCTK_PARAMETERS; - if (called_from_recovery) { - arrdata.at(group).at(Carpet::map).tt->set_time(reflevel,mglevel, - (CCTK_REAL) ((cctkGH->cctk_time - cctk_initial_time) - / (delta_time * mglevelfact)) ); - } + const int n = CCTK_VarIndex(varname); + assert (n>=0 && n<CCTK_NumVars()); + const int group = CCTK_GroupIndexFromVarI (n); + assert (group>=0 && group<(int)Carpet::arrdata.size()); + const int n0 = CCTK_FirstVarIndexI(group); + assert (n0>=0 && n0<CCTK_NumVars()); + const int var = n - n0; + assert (var>=0 && var<CCTK_NumVars()); + herr_t herr = 1; + int want_dataset = 0; + bool did_read_something = false; + int ndatasets = 0; + hid_t dataset = 0; - } END_MAP_LOOP; - - } // if amr_level == rl - if (CCTK_MyProc(cctkGH)==0) { - free (h5data); - } - - return did_read_something; + char datasetname[1024]; + // Check for storage + if (! CCTK_QueryGroupStorageI(cctkGH, group)) { + CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING, + "Cannot input variable \"%s\" because it has no storage", + varname); + return 0; } + const int grouptype = CCTK_GroupTypeI(group); + const int rl = grouptype==CCTK_GF ? reflevel : 0; + //cout << "want level " << rl << " reflevel " << reflevel << endl; + // Find the input directory + const char* const myindir = in3D_dir; - int InputVarAs (const cGH* const cctkGH, const char* const varname, - const char* const alias) { - DECLARE_CCTK_PARAMETERS; + // Invent a file name + ostringstream filenamebuf; + filenamebuf << myindir << "/" << alias << in3D_extension; + string filenamestr = filenamebuf.str(); + const char * const filename = filenamestr.c_str(); - const int n = CCTK_VarIndex(varname); - assert (n>=0 && n<CCTK_NumVars()); - const int group = CCTK_GroupIndexFromVarI (n); - assert (group>=0 && group<(int)Carpet::arrdata.size()); - const int n0 = CCTK_FirstVarIndexI(group); - assert (n0>=0 && n0<CCTK_NumVars()); - const int var = n - n0; - assert (var>=0 && var<CCTK_NumVars()); + hid_t reader = -1; - herr_t herr = 1; - int want_dataset = 0; - bool did_read_something = false; - int ndatasets = 0; - hid_t dataset = 0; + // Read the file only on the root processor + if (CCTK_MyProc(cctkGH)==0) { - char datasetname[1024]; - - // Check for storage - if (! CCTK_QueryGroupStorageI(cctkGH, group)) { - CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING, - "Cannot input variable \"%s\" because it has no storage", - varname); - return 0; + // Open the file + if (verbose) CCTK_VInfo (CCTK_THORNSTRING, "Opening file \"%s\"", filename); + reader = H5Fopen (filename, H5F_ACC_RDONLY, H5P_DEFAULT); + if (reader<0) { + CCTK_VWarn (0, __LINE__, __FILE__, CCTK_THORNSTRING, + "Could not open file \"%s\" for reading", filename); } + assert (reader>=0); + // get the number of datasets in the file + ndatasets=GetnDatasets(reader); + } - const int grouptype = CCTK_GroupTypeI(group); - const int rl = grouptype==CCTK_GF ? reflevel : 0; - //cout << "want level " << rl << " reflevel " << reflevel << endl; + vector<ibset> regions_read(Carpet::maps); - // Find the input directory - const char* const myindir = in3D_dir; + // Broadcast number of datasets + MPI_Bcast (&ndatasets, 1, MPI_INT, 0, dist::comm); + assert (ndatasets>=0); - // Invent a file name - ostringstream filenamebuf; - filenamebuf << myindir << "/" << alias << in3D_extension; - string filenamestr = filenamebuf.str(); - const char * const filename = filenamestr.c_str(); - hid_t reader = -1; + for (int datasetid=0; datasetid<ndatasets; ++datasetid) { + if (verbose) CCTK_VInfo (CCTK_THORNSTRING, "Handling dataset #%d", datasetid); - // Read the file only on the root processor + + // Read data if (CCTK_MyProc(cctkGH)==0) { + GetDatasetName(reader,datasetid,datasetname); + // cout << datasetname << "\n"; - // Open the file - if (verbose) CCTK_VInfo (CCTK_THORNSTRING, "Opening file \"%s\"", filename); - reader = H5Fopen (filename, H5F_ACC_RDONLY, H5P_DEFAULT); - if (reader<0) { - CCTK_VWarn (0, __LINE__, __FILE__, CCTK_THORNSTRING, - "Could not open file \"%s\" for reading", filename); - } - assert (reader>=0); - // get the number of datasets in the file - ndatasets=GetnDatasets(reader); + dataset = H5Dopen (reader, datasetname); + assert(dataset); } - vector<ibset> regions_read(Carpet::maps); - - // Broadcast number of datasets - MPI_Bcast (&ndatasets, 1, MPI_INT, 0, dist::comm); - assert (ndatasets>=0); - - - for (int datasetid=0; datasetid<ndatasets; ++datasetid) { - if (verbose) CCTK_VInfo (CCTK_THORNSTRING, "Handling dataset #%d", datasetid); - - - // Read data - if (CCTK_MyProc(cctkGH)==0) { - GetDatasetName(reader,datasetid,datasetname); - // cout << datasetname << "\n"; - - dataset = H5Dopen (reader, datasetname); - assert(dataset); - } - #if 0 - int amr_level; - int amr_origin[dim]; - int amr_dims[dim]; + int amr_level; + int amr_origin[dim]; + int amr_dims[dim]; #endif - if (CCTK_MyProc(cctkGH)==0) { - - // Read data - char * name; - ReadAttribute (dataset, "name", name); - // cout << "dataset name is " << name << endl; - if (verbose) { - if (name) { - CCTK_VInfo (CCTK_THORNSTRING, "Dataset name is \"%s\"", name); - } - } - want_dataset = name && CCTK_EQUALS(name, varname); - free (name); - } // myproc == 0 - - MPI_Bcast (&want_dataset, 1, MPI_INT, 0, dist::comm); - - - if(want_dataset) { - - did_read_something = ReadVar(cctkGH,reader,varname,dataset,regions_read,0); - - } // want_dataset - - } // loop over datasets - - // Close the file if (CCTK_MyProc(cctkGH)==0) { - if (verbose) CCTK_VInfo (CCTK_THORNSTRING, "Closing file"); - herr = H5Fclose(reader); - // cout << "blah! " << reader << "\n"; - // cout << "closing file " << herr << "\n"; - assert(!herr); - reader=-1; - } - // Was everything initialised? - if (did_read_something) { - for (int m=0; m<Carpet::maps; ++m) { - dh<dim>& thedd = *arrdata.at(group).at(m).dd; - ibset all_exterior; - for (size_t c=0; c<thedd.boxes.at(rl).size(); ++c) { - all_exterior |= thedd.boxes.at(rl).at(c).at(mglevel).exterior; - } - if (regions_read.at(m) != all_exterior) { - cout << "read: " << regions_read.at(m) << endl - << "want: " << all_exterior << endl; - CCTK_VWarn (0, __LINE__, __FILE__, CCTK_THORNSTRING, - "Variable \"%s\" could not be initialised from file -- the file may be missing data", - varname); - } + // Read data + char * name; + ReadAttribute (dataset, "name", name); + // cout << "dataset name is " << name << endl; + if (verbose && name) { + CCTK_VInfo (CCTK_THORNSTRING, "Dataset name is \"%s\"", name); + } + want_dataset = name && CCTK_EQUALS(name, varname); + free (name); + } // myproc == 0 + + MPI_Bcast (&want_dataset, 1, MPI_INT, 0, dist::comm); + + if(want_dataset) { + did_read_something = ReadVar(cctkGH,varname,dataset,regions_read,0); + } // want_dataset + + } // loop over datasets + + // Close the file + if (CCTK_MyProc(cctkGH)==0) { + if (verbose) CCTK_VInfo (CCTK_THORNSTRING, "Closing file"); + herr = H5Fclose(reader); + // cout << "blah! " << reader << "\n"; + // cout << "closing file " << herr << "\n"; + assert(!herr); + reader=-1; + } + + // Was everything initialised? + if (did_read_something) { + for (int m=0; m<Carpet::maps; ++m) { + dh<dim>& thedd = *arrdata.at(group).at(m).dd; + ibset all_exterior; + for (size_t c=0; c<thedd.boxes.at(rl).size(); ++c) { + all_exterior |= thedd.boxes.at(rl).at(c).at(mglevel).exterior; + } + if (regions_read.at(m) != all_exterior) { + cout << "read: " << regions_read.at(m) << endl + << "want: " << all_exterior << endl; + CCTK_VWarn (0, __LINE__, __FILE__, CCTK_THORNSTRING, + "Variable \"%s\" could not be initialised from file -- the file may be missing data", + varname); } - } // if did_read_something - // CCTK_VWarn (0, __LINE__, __FILE__, CCTK_THORNSTRING,"stop!"); + } + } // if did_read_something + // CCTK_VWarn (0, __LINE__, __FILE__, CCTK_THORNSTRING,"stop!"); - return did_read_something ? 0 : -1; + return did_read_something ? 0 : -1; - } +} - /** returns the number of recovered variables */ - int Recover (cGH* const cctkGH, const char *basefilename, - const int called_from) - { - assert (cctkGH); - assert (basefilename); - assert (called_from == CP_INITIAL_DATA - || called_from == CP_EVOLUTION_DATA - || called_from == CP_RECOVER_PARAMETERS - || called_from == CP_RECOVER_DATA - || called_from == FILEREADER_DATA); - - // the other modes are not supported yet - assert (called_from == FILEREADER_DATA); - - const ioGH * const iogh = (const ioGH *) CCTK_GHExtension (cctkGH, "IO"); - assert (iogh); - - int num_vars_read = 0; - assert (iogh->do_inVars); - for (int n=0; n<CCTK_NumVars(); ++n) { - if (iogh->do_inVars[n]) { - char * const fullname = CCTK_FullName(n); - assert (fullname); - const int ierr = InputVarAs (cctkGH, fullname, basefilename); - if (! ierr) { - ++ num_vars_read; - } - free (fullname); - } +int CarpetIOHDF5ReadData (const cGH* const cctkGH) +{ + DECLARE_CCTK_PARAMETERS; + int retval = 0; + for (int vindex=0; vindex<CCTK_NumVars(); ++vindex) { + if (CheckForVariable(in3D_vars, vindex)) { + char* varname = CCTK_FullName(vindex); + retval = InputVarAs (cctkGH, varname, CCTK_VarName(vindex)); + free (varname); + if (retval != 0) return retval; } - - return num_vars_read; } + return retval; +} - void CarpetIOHDF5ReadData (CCTK_ARGUMENTS) - { - InputGH (cctkGH); +#if 0 +/** returns the number of recovered variables */ +int Recover (cGH* const cctkGH, const char *basefilename, + const int called_from) +{ + assert (cctkGH); + assert (basefilename); + assert (called_from == CP_INITIAL_DATA + || called_from == CP_EVOLUTION_DATA + || called_from == CP_RECOVER_PARAMETERS + || called_from == CP_RECOVER_DATA + || called_from == FILEREADER_DATA); + + // the other modes are not supported yet + assert (called_from == FILEREADER_DATA); + + const ioGH * const iogh = (const ioGH *) CCTK_GHExtension (cctkGH, "IO"); + assert (iogh); + + int num_vars_read = 0; + assert (iogh->do_inVars); + for (int n=0; n<CCTK_NumVars(); ++n) { + if (iogh->do_inVars[n]) { + char * const fullname = CCTK_FullName(n); + assert (fullname); + const int ierr = InputVarAs (cctkGH, fullname, basefilename); + if (! ierr) { + ++ num_vars_read; + } + free (fullname); + } } - + return num_vars_read; +} +#endif } // namespace CarpetIOHDF5 |