diff options
Diffstat (limited to 'Carpet')
-rw-r--r-- | Carpet/CarpetIOHDF5/param.ccl | 3 | ||||
-rw-r--r-- | Carpet/CarpetIOHDF5/schedule.ccl | 5 | ||||
-rw-r--r-- | Carpet/CarpetIOHDF5/src/iohdf5.cc | 495 | ||||
-rw-r--r-- | Carpet/CarpetIOHDF5/src/iohdf5.hh | 39 | ||||
-rw-r--r-- | Carpet/CarpetIOHDF5/src/iohdf5chckpt_recover.cc | 310 |
5 files changed, 633 insertions, 219 deletions
diff --git a/Carpet/CarpetIOHDF5/param.ccl b/Carpet/CarpetIOHDF5/param.ccl index 3d9d55b28..62912c491 100644 --- a/Carpet/CarpetIOHDF5/param.ccl +++ b/Carpet/CarpetIOHDF5/param.ccl @@ -1,5 +1,5 @@ # Parameter definitions for thorn CarpetIOHDF5 -# $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetIOHDF5/param.ccl,v 1.2 2004/03/10 21:23:49 cott Exp $ +# $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetIOHDF5/param.ccl,v 1.3 2004/03/12 00:13:25 cott Exp $ @@ -16,6 +16,7 @@ USES BOOLEAN recover_and_remove USES BOOLEAN checkpoint_on_terminate USES KEYWORD recover USES STRING recover_file +USES STRING recover_dir diff --git a/Carpet/CarpetIOHDF5/schedule.ccl b/Carpet/CarpetIOHDF5/schedule.ccl index 8ab0101d5..6ed265c02 100644 --- a/Carpet/CarpetIOHDF5/schedule.ccl +++ b/Carpet/CarpetIOHDF5/schedule.ccl @@ -1,5 +1,5 @@ # Schedule definitions for thorn CarpetIOHDF5 -# $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetIOHDF5/schedule.ccl,v 1.4 2004/03/11 09:33:23 cott Exp $ +# $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetIOHDF5/schedule.ccl,v 1.5 2004/03/12 00:13:25 cott Exp $ schedule CarpetIOHDF5Startup at STARTUP after IOUtil_Startup { @@ -15,7 +15,7 @@ schedule CarpetIOHDF5ReadData at INITIAL schedule CarpetIOHDF5_EvolutionCheckpoint at CHECKPOINT { LANG: C - OPTIONS: global + OPTIONS: meta } "Do checkpointing" if (! CCTK_Equals (recover, "no") && *recover_file) @@ -23,5 +23,6 @@ if (! CCTK_Equals (recover, "no") && *recover_file) schedule CarpetIOHDF5_RecoverParameters at RECOVER_PARAMETERS { LANG:C + OPTIONS: meta } "Parameter recovery routine" } diff --git a/Carpet/CarpetIOHDF5/src/iohdf5.cc b/Carpet/CarpetIOHDF5/src/iohdf5.cc index 543bfd843..a11519714 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.13 2004/03/11 11:50:51 cott Exp $"; + static const char* rcsid = "$Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetIOHDF5/src/iohdf5.cc,v 1.14 2004/03/12 00:13:25 cott Exp $"; CCTK_FILEVERSION(Carpet_CarpetIOHDF5_iohdf5_cc); } @@ -280,6 +280,8 @@ namespace CarpetIOHDF5 { const int var = n - n0; assert (var>=0 && var<CCTK_NumVars()); int tl = 0; + + const int gpdim = CCTK_GroupDimI(group); // Check for storage if (! CCTK_QueryGroupStorageI(cctkGH, group)) { @@ -303,6 +305,10 @@ namespace CarpetIOHDF5 { } const int rl = grouptype==CCTK_GF ? reflevel : 0; + cGroup cgdata; + int ierr = CCTK_GroupData(group,&cgdata); + assert(ierr==0); + // let's get the correct Carpet time level (which is the (-1) * Cactus timelevel): if (request->timelevel==0) tl = 0; @@ -327,181 +333,203 @@ namespace CarpetIOHDF5 { gdata<dim>* const tmp = data->make_typed (n); tmp->allocate (ext, 0); - for (comm_state<dim> state; !state.done(); state.step()) { - tmp->copy_from (state, data, ext); - } - - // Write data - if (CCTK_MyProc(cctkGH)==0) { - - hsize_t shape[dim]; - for (int d=0; d<dim; ++d) { - shape[dim-1-d] = (ext.shape() / ext.stride())[d]; - } - const hid_t dataspace = H5Screate_simple (dim, shape, NULL); - assert (dataspace>=0); + + if ( !((cgdata.disttype == CCTK_DISTRIB_CONSTANT) && + (arrdata[group][Carpet::map].hh->processors[reflevel][component]!=0))) { + + if (cgdata.disttype == CCTK_DISTRIB_CONSTANT) { + assert(grouptype == CCTK_ARRAY || grouptype == CCTK_SCALAR); + if(component!=0) continue; + } 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; + if ( grouptype==CCTK_SCALAR ) { + ldim = 1; + } else { + ldim = gpdim; + } + + hsize_t shape[ldim]; + for (int d=0; d<ldim; ++d) { + shape[ldim-1-d] = (ext.shape() / ext.stride())[d]; + } + const hid_t dataspace = H5Screate_simple (ldim, shape, NULL); + assert (dataspace>=0); + +// hsize_t shape[dim]; +// for (int d=0; d<dim; ++d) { +// shape[dim-1-d] = (ext.shape() / ext.stride())[d]; +// } +// const hid_t dataspace = H5Screate_simple (dim, shape, NULL); +// assert (dataspace>=0); +// // Select datatype #if 0 - assert (true - || (CCTK_VarTypeI(n) == CCTK_VARIABLE_REAL8 - && sizeof(CCTK_REAL8) == sizeof(double)) - || (CCTK_VarTypeI(n) == CCTK_VARIABLE_REAL - && sizeof(CCTK_REAL) == sizeof(double))); - // TODO: Set datatype correctly + assert (true + || (CCTK_VarTypeI(n) == CCTK_VARIABLE_REAL8 + && sizeof(CCTK_REAL8) == sizeof(double)) + || (CCTK_VarTypeI(n) == CCTK_VARIABLE_REAL + && sizeof(CCTK_REAL) == sizeof(double))); + // TODO: Set datatype correctly #endif - const hid_t datatype = h5DataType(CCTK_VarTypeI(n)); - - ostringstream datasetnamebuf; - datasetnamebuf << varname - << " it=" << cctk_iteration - << " tl=" << tl - << " ml=" << mglevel - << " rl=" << rl - << " m=" << Carpet::map - << " c=" << component; - string datasetnamestr = datasetnamebuf.str(); - const char * const datasetname = datasetnamestr.c_str(); - const hid_t dataset = H5Dcreate (writer, datasetname, datatype, dataspace, H5P_DEFAULT); - assert (dataset>=0); - - const void * const data = (void*)tmp->storage(); - herr = H5Dwrite (dataset, datatype, H5S_ALL, H5S_ALL, H5P_DEFAULT, data); - assert (!herr); - - // Write FlexIO attributes - WriteAttribute (dataset, "level", rl); - { - CCTK_REAL origin[dim], delta[dim]; - CCTK_REAL min_ext[dim], max_ext[dim]; - for (int d=0; d<dim; ++d) { - origin[d] = CCTK_ORIGIN_SPACE(d) + cctk_lbnd[d] * delta[d]; - delta[d] = CCTK_DELTA_SPACE(d); - min_ext[d] = origin[d]; - max_ext[d] = origin[d] + cctk_lsh[d] * delta[d]; - } - WriteAttribute (dataset, "origin", origin, dim); - WriteAttribute (dataset, "delta", delta, dim); - WriteAttribute (dataset, "min_ext", min_ext, dim); - WriteAttribute (dataset, "max_ext", max_ext, dim); - } - WriteAttribute (dataset, "time", cctk_time); - WriteAttribute (dataset, "timestep", cctk_iteration); - WriteAttribute (dataset, "level_timestep", cctk_iteration / reflevelfact); - WriteAttribute (dataset, "persistence", maxreflevelfact / reflevelfact); - { - int time_refinement; - int spatial_refinement[dim]; - int grid_placement_refinement[dim]; - time_refinement = reflevelfact; - for (int d=0; d<dim; ++d) { - spatial_refinement[d] = reflevelfact; - grid_placement_refinement[d] = reflevelfact; - } - WriteAttribute (dataset, "time_refinement", time_refinement); - WriteAttribute (dataset, "spatial_refinement", spatial_refinement, dim); - WriteAttribute (dataset, "grid_placement_refinement", grid_placement_refinement, dim); - } - { - int iorigin[dim]; - for (int d=0; d<dim; ++d) { - iorigin[d] = (ext.lower() / ext.stride())[d]; - } - WriteAttribute (dataset, "iorigin", iorigin, dim); - } - - // Write some additional attributes - - // Legacy arguments - { - char * fullname = CCTK_FullName(n); - assert (fullname); - WriteAttribute (dataset, "name", fullname); - free (fullname); - } - - // Group arguments - WriteAttribute (dataset, "group_version", 1); - { - char * fullname = CCTK_FullName(n); - assert (fullname); - WriteAttribute (dataset, "group_fullname", fullname); - free (fullname); - } - WriteAttribute (dataset, "group_varname", CCTK_VarName(n)); - { - char * groupname = CCTK_GroupName(group); - assert (groupname); - WriteAttribute (dataset, "group_groupname", groupname); - free (groupname); - } - switch (grouptype) { - case CCTK_GF: - WriteAttribute (dataset, "group_grouptype", "CCTK_GF"); - break; - case CCTK_ARRAY: - WriteAttribute (dataset, "group_grouptype", "CCTK_ARRAY"); - break; - case CCTK_SCALAR: - WriteAttribute (dataset, "group_grouptype", "CCTK_SCALAR"); - break; - default: - assert (0); - } - WriteAttribute (dataset, "group_dim", CCTK_GroupDimI(group)); - WriteAttribute (dataset, "group_timelevel", tl); - WriteAttribute (dataset, "group_numtimelevels", CCTK_NumTimeLevelsI(group)); - - // Cactus arguments - WriteAttribute (dataset, "cctk_version", 1); - WriteAttribute (dataset, "cctk_dim", cctk_dim); - WriteAttribute (dataset, "cctk_iteration", cctk_iteration); -// TODO: disable temporarily -// WriteAttribute (dataset, "cctk_nmaps", cctk_nmaps); -// WriteAttribute (dataset, "cctk_map", cctk_map); - WriteAttribute (dataset, "cctk_gsh", cctk_gsh, dim); - WriteAttribute (dataset, "cctk_lsh", cctk_lsh, dim); - WriteAttribute (dataset, "cctk_lbnd", cctk_lbnd, dim); - WriteAttribute (dataset, "cctk_delta_time", cctk_delta_time); - WriteAttribute (dataset, "cctk_delta_space", cctk_delta_space, dim); - WriteAttribute (dataset, "cctk_origin_space", cctk_origin_space, dim); - WriteAttribute (dataset, "cctk_bbox", cctk_bbox, 2*dim); - WriteAttribute (dataset, "cctk_levfac", cctk_levfac, dim); - WriteAttribute (dataset, "cctk_levoff", cctk_levoff, dim); - WriteAttribute (dataset, "cctk_levoffdenom", cctk_levoffdenom, dim); - WriteAttribute (dataset, "cctk_timefac", cctk_timefac); - WriteAttribute (dataset, "cctk_convlevel", cctk_convlevel); - WriteAttribute (dataset, "cctk_convfac", cctk_convfac); - WriteAttribute (dataset, "cctk_nghostzones", cctk_nghostzones, dim); - WriteAttribute (dataset, "cctk_time", cctk_time); + const hid_t datatype = h5DataType(CCTK_VarTypeI(n)); - // Carpet arguments - WriteAttribute (dataset, "carpet_version", 1); - WriteAttribute (dataset, "carpet_dim", dim); - WriteAttribute (dataset, "carpet_basemglevel", basemglevel); - WriteAttribute (dataset, "carpet_mglevel", mglevel); - WriteAttribute (dataset, "carpet_mglevels", mglevels); - WriteAttribute (dataset, "carpet_mgface", mgfact); - WriteAttribute (dataset, "carpet_reflevel", reflevel); - WriteAttribute (dataset, "carpet_reflevels", reflevels); - WriteAttribute (dataset, "carpet_reffact", reffact); - WriteAttribute (dataset, "carpet_map", Carpet::map); - WriteAttribute (dataset, "carpet_maps", maps); - WriteAttribute (dataset, "carpet_component", component); - WriteAttribute (dataset, "carpet_components", vhh.at(Carpet::map)->components(reflevel)); + ostringstream datasetnamebuf; + datasetnamebuf << varname + << " it=" << cctk_iteration + << " tl=" << tl + << " ml=" << mglevel + << " rl=" << rl + << " m=" << Carpet::map + << " c=" << component; + string datasetnamestr = datasetnamebuf.str(); + const char * const datasetname = datasetnamestr.c_str(); + const hid_t dataset = H5Dcreate (writer, datasetname, datatype, dataspace, H5P_DEFAULT); + assert (dataset>=0); - herr = H5Dclose (dataset); - assert (!herr); - - herr = H5Sclose (dataspace); - assert (!herr); + const void * const data = (void*)tmp->storage(); + herr = H5Dwrite (dataset, datatype, H5S_ALL, H5S_ALL, H5P_DEFAULT, data); + assert (!herr); - } // if on root processor - - // Delete temporary copy - delete tmp; - + // Write FlexIO attributes + WriteAttribute (dataset, "level", rl); + { + CCTK_REAL origin[dim], delta[dim]; + CCTK_REAL min_ext[dim], max_ext[dim]; + for (int d=0; d<dim; ++d) { + origin[d] = CCTK_ORIGIN_SPACE(d) + cctk_lbnd[d] * delta[d]; + delta[d] = CCTK_DELTA_SPACE(d); + min_ext[d] = origin[d]; + max_ext[d] = origin[d] + cctk_lsh[d] * delta[d]; + } + WriteAttribute (dataset, "origin", origin, dim); + WriteAttribute (dataset, "delta", delta, dim); + WriteAttribute (dataset, "min_ext", min_ext, dim); + WriteAttribute (dataset, "max_ext", max_ext, dim); + } + WriteAttribute (dataset, "time", cctk_time); + WriteAttribute (dataset, "timestep", cctk_iteration); + WriteAttribute (dataset, "level_timestep", cctk_iteration / reflevelfact); + WriteAttribute (dataset, "persistence", maxreflevelfact / reflevelfact); + { + int time_refinement; + int spatial_refinement[dim]; + int grid_placement_refinement[dim]; + time_refinement = reflevelfact; + for (int d=0; d<dim; ++d) { + spatial_refinement[d] = reflevelfact; + grid_placement_refinement[d] = reflevelfact; + } + WriteAttribute (dataset, "time_refinement", time_refinement); + WriteAttribute (dataset, "spatial_refinement", spatial_refinement, dim); + WriteAttribute (dataset, "grid_placement_refinement", grid_placement_refinement, dim); + } + { + int iorigin[dim]; + for (int d=0; d<dim; ++d) { + iorigin[d] = (ext.lower() / ext.stride())[d]; + } + WriteAttribute (dataset, "iorigin", iorigin, dim); + } + + // Write some additional attributes + + // Legacy arguments + { + char * fullname = CCTK_FullName(n); + assert (fullname); + WriteAttribute (dataset, "name", fullname); + free (fullname); + } + + // Group arguments + WriteAttribute (dataset, "group_version", 1); + { + char * fullname = CCTK_FullName(n); + assert (fullname); + WriteAttribute (dataset, "group_fullname", fullname); + free (fullname); + } + WriteAttribute (dataset, "group_varname", CCTK_VarName(n)); + { + char * groupname = CCTK_GroupName(group); + assert (groupname); + WriteAttribute (dataset, "group_groupname", groupname); + free (groupname); + } + switch (grouptype) { + case CCTK_GF: + WriteAttribute (dataset, "group_grouptype", "CCTK_GF"); + break; + case CCTK_ARRAY: + WriteAttribute (dataset, "group_grouptype", "CCTK_ARRAY"); + break; + case CCTK_SCALAR: + WriteAttribute (dataset, "group_grouptype", "CCTK_SCALAR"); + break; + default: + assert (0); + } + WriteAttribute (dataset, "group_dim", CCTK_GroupDimI(group)); + WriteAttribute (dataset, "group_timelevel", tl); + WriteAttribute (dataset, "group_numtimelevels", CCTK_NumTimeLevelsI(group)); + + // Cactus arguments + WriteAttribute (dataset, "cctk_version", 1); + WriteAttribute (dataset, "cctk_dim", cctk_dim); + WriteAttribute (dataset, "cctk_iteration", cctk_iteration); + // TODO: disable temporarily + // WriteAttribute (dataset, "cctk_nmaps", cctk_nmaps); + // WriteAttribute (dataset, "cctk_map", cctk_map); + WriteAttribute (dataset, "cctk_gsh", cctk_gsh, dim); + WriteAttribute (dataset, "cctk_lsh", cctk_lsh, dim); + WriteAttribute (dataset, "cctk_lbnd", cctk_lbnd, dim); + WriteAttribute (dataset, "cctk_delta_time", cctk_delta_time); + WriteAttribute (dataset, "cctk_delta_space", cctk_delta_space, dim); + WriteAttribute (dataset, "cctk_origin_space", cctk_origin_space, dim); + WriteAttribute (dataset, "cctk_bbox", cctk_bbox, 2*dim); + WriteAttribute (dataset, "cctk_levfac", cctk_levfac, dim); + WriteAttribute (dataset, "cctk_levoff", cctk_levoff, dim); + WriteAttribute (dataset, "cctk_levoffdenom", cctk_levoffdenom, dim); + WriteAttribute (dataset, "cctk_timefac", cctk_timefac); + WriteAttribute (dataset, "cctk_convlevel", cctk_convlevel); + WriteAttribute (dataset, "cctk_convfac", cctk_convfac); + WriteAttribute (dataset, "cctk_nghostzones", cctk_nghostzones, dim); + WriteAttribute (dataset, "cctk_time", cctk_time); + + // Carpet arguments + WriteAttribute (dataset, "carpet_version", 1); + WriteAttribute (dataset, "carpet_dim", dim); + WriteAttribute (dataset, "carpet_basemglevel", basemglevel); + WriteAttribute (dataset, "carpet_mglevel", mglevel); + WriteAttribute (dataset, "carpet_mglevels", mglevels); + WriteAttribute (dataset, "carpet_mgface", mgfact); + WriteAttribute (dataset, "carpet_reflevel", reflevel); + WriteAttribute (dataset, "carpet_reflevels", reflevels); + WriteAttribute (dataset, "carpet_reffact", reffact); + WriteAttribute (dataset, "carpet_map", Carpet::map); + WriteAttribute (dataset, "carpet_maps", maps); + WriteAttribute (dataset, "carpet_component", component); + WriteAttribute (dataset, "carpet_components", vhh.at(Carpet::map)->components(reflevel)); + + herr = H5Dclose (dataset); + assert (!herr); + + herr = H5Sclose (dataspace); + assert (!herr); + + } // if on root processor + + // Delete temporary copy + delete tmp; + } //if( !((CCTK_DISTRIB_BLAH) } END_COMPONENT_LOOP; } END_MAP_LOOP; @@ -595,7 +623,7 @@ namespace CarpetIOHDF5 { } -#if 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) { @@ -610,11 +638,17 @@ namespace CarpetIOHDF5 { assert (n0>=0 && n0<CCTK_NumVars()); const int var = n - n0; assert (var>=0 && var<CCTK_NumVars()); - const int tl = 0; + int tl = 0; herr_t herr = 1; bool did_read_something = false; + // Stuff needed for Recovery + int recovery_tl = -1; + int recovery_mglevel = -1; + int recovery_rl = -1; + int recovery_comp = -1; + CCTK_REAL *h5data; // Check for storage @@ -636,24 +670,34 @@ namespace CarpetIOHDF5 { 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; + if (CCTK_MyProc(cctkGH)==0) { // get dataset dimensions const hid_t dataspace = H5Dget_space(dataset); assert (dataspace>=0); - hsize_t rank=H5Sget_simple_extent_ndims(dataspace); + hsize_t rank = H5Sget_simple_extent_ndims(dataspace); hsize_t shape[rank]; int rank2 = H5Sget_simple_extent_dims (dataspace, shape, NULL); herr = H5Sclose(dataspace); assert(!herr); assert (rank2 == rank); - assert (gpdim == rank); - + + if(grouptype == CCTK_SCALAR) { + assert (gpdim+1 == rank); + } else { + assert (gpdim == rank); + } + int datalength=1; for(int i=0;i<rank;i++) { datalength=datalength*shape[i]; + amr_dims[i]=shape[i]; + // cout << "amr_dims[" << i << "]: "<<amr_dims[i] << endl; } const hid_t datatype = H5T_NATIVE_DOUBLE; @@ -664,6 +708,7 @@ namespace CarpetIOHDF5 { herr = H5Dread(dataset,datatype,H5S_ALL, H5S_ALL, H5P_DEFAULT,(void*)h5data); assert(!herr); + // cout << h5data[100] << endl; //cout << datasetname << endl; //cout << name << endl; @@ -673,19 +718,26 @@ namespace CarpetIOHDF5 { ReadAttribute(dataset,"iorigin",amr_origin,dim); assert(herr>=0); // cout << "amr_origin[0] " << amr_origin[0] << "\n"; - - herr = H5Dclose(dataset); - assert(!herr); - for (int d=0; d<gpdim; ++d) { - amr_dims[d] = shape[gpdim-1-d]; + + 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); } - } // MyProc == 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); + + 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); + } + + 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); if (amr_level == rl) { // cout << "I want this" << endl; @@ -700,34 +752,71 @@ namespace CarpetIOHDF5 { 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; + gdata<dim>* const data = (*ff) (tl, rl, component, mglevel); // Create temporary data storage on processor 0 - const vect<int,dim> str + vect<int,dim> str = vect<int,dim>(maxreflevelfact/reflevelfact); - const vect<int,dim> lb = vect<int,dim>::ref(amr_origin) * str; - const vect<int,dim> ub + vect<int,dim> lb = vect<int,dim>::ref(amr_origin) * str; + vect<int,dim> ub = lb + (vect<int,dim>::ref(amr_dims) - 1) * str; - const bbox<int,dim> ext(lb,ub,str); - + + // cout << "lb: " << lb << endl; + // cout << "ub: " << ub << endl; + // cout << "str: " << str << endl; gdata<dim>* const tmp = data->make_typed (n); + //cout << "overlap " << overlap << endl; + //cout << "tmp->extent " << tmp->extent() << endl; + //cout << "data->extent " << data->extent() << endl; + + cGroup cgdata; + int ierr = CCTK_GroupData(group,&cgdata); + assert(ierr==0); + + // cout << "trying to handle DISTRIB=const data " << endl; + + // handle distrib=constant vars + + if(grouptype == CCTK_SCALAR || grouptype == CCTK_ARRAY) + str = vect<int,dim> (1); + + if (cgdata.disttype == CCTK_DISTRIB_CONSTANT) { + assert(grouptype == CCTK_ARRAY || grouptype == CCTK_SCALAR); + if (grouptype == CCTK_SCALAR) { + lb[0] = arrdata[group][Carpet::map].hh->processors.at(rl).at(component); + ub[0] = arrdata[group][Carpet::map].hh->processors.at(rl).at(component); + for(int i=1;i<dim;i++) { + lb[i]=0; + ub[i]=0; + } + } else { + lb[dim-1] = lb[dim-1] + (ub[dim-1]-lb[dim-1]+1)*(arrdata[group][Carpet::map].hh->processors.at(rl).at(component)); + ub[dim-1] = ub[dim-1] + (ub[dim-1]-lb[dim-1]+1)*(arrdata[group][Carpet::map].hh->processors.at(rl).at(component)); + // cout << "lb: " << lb << endl; + // cout << "ub: " << ub << endl; + } + } + const bbox<int,dim> ext(lb,ub,str); + + // cout << ext << endl; + 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; - - //cout << "overlap " << overlap << endl; - //cout << "tmp->extent " << tmp->extent() << endl; - //cout << "data->extent " << data->extent() << endl; - + + + // Copy into grid function for (comm_state<dim> state; !state.done(); state.step()) { data->copy_from (state, tmp, overlap); @@ -735,7 +824,13 @@ namespace CarpetIOHDF5 { // Delete temporary copy delete tmp; - + + // set tt (ask Erik why...) + if (called_from_recovery) { + arrdata[group][Carpet::map].tt->set_time(reflevel,mglevel, + (CCTK_REAL) cctkGH->cctk_iteration/maxreflevelfact); + } + } END_COMPONENT_LOOP; } END_MAP_LOOP; @@ -749,7 +844,7 @@ namespace CarpetIOHDF5 { } -#endif + int InputVarAs (const cGH* const cctkGH, const char* const varname, const char* const alias) { @@ -821,8 +916,8 @@ namespace CarpetIOHDF5 { MPI_Bcast (&ndatasets, 1, MPI_INT, 0, dist::comm); assert (ndatasets>=0); - - for (int datasetid=0; datasetid<ndatasets; ++datasetid) { + // added +1 to ndatasets since the number seems to be 1 to small + for (int datasetid=0; datasetid<ndatasets+1; ++datasetid) { if (verbose) CCTK_VInfo (CCTK_THORNSTRING, "Handling dataset #%d", datasetid); diff --git a/Carpet/CarpetIOHDF5/src/iohdf5.hh b/Carpet/CarpetIOHDF5/src/iohdf5.hh index fd4c6db8c..a0a5ca44b 100644 --- a/Carpet/CarpetIOHDF5/src/iohdf5.hh +++ b/Carpet/CarpetIOHDF5/src/iohdf5.hh @@ -1,4 +1,4 @@ -// $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetIOHDF5/src/iohdf5.hh,v 1.7 2004/03/11 10:00:16 cott Exp $ +// $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetIOHDF5/src/iohdf5.hh,v 1.8 2004/03/12 00:13:25 cott Exp $ #ifndef CARPETIOHDF5_HH #define CARPETIOHDF5_HH @@ -13,6 +13,43 @@ #include "iohdf5.h" #include "CactusBase/IOUtil/src/ioutil_Utils.h" +// Some MPI Datatypes we need for Recovery +// Originally written by Thomas Radke. + +#ifdef CCTK_INT4 +#define CARPET_MPI_INT4 (sizeof (CCTK_INT4) == sizeof (int) ? MPI_INT : \ + sizeof (CCTK_INT4) == sizeof (short) ? MPI_SHORT : \ + MPI_DATATYPE_NULL) +#endif + +#define CARPET_MPI_CHAR MPI_CHAR + +/* floating point types are architecture-independent, + ie. a float has always 4 bytes, and a double has 8 bytes + + PUGH_MPI_REAL is used for communicating reals of the generic CCTK_REAL type + PUGH_MPI_REALn is used to explicitely communicate n-byte reals */ +#ifdef CCTK_REAL4 +#define CARPET_MPI_REAL4 MPI_FLOAT +#endif +#ifdef CCTK_REAL8 +#define CARPET_MPI_REAL8 MPI_DOUBLE +#endif +#ifdef CCTK_REAL16 +#define CARPET_MPI_REAL16 (sizeof (CCTK_REAL16) == sizeof (long double) ? \ + MPI_LONG_DOUBLE : MPI_DATATYPE_NULL) +#endif + + +#ifdef CCTK_REAL_PRECISION_16 +#define CARPET_MPI_REAL CARPET_MPI_REAL16 +#elif CCTK_REAL_PRECISION_8 +#define CARPET_MPI_REAL CARPET_MPI_REAL8 +#elif CCTK_REAL_PRECISION_4 +#define CARPET_MPI_REAL CARPET_MPI_REAL4 +#endif + + /*** Define the different datatypes used for HDF5 I/O NOTE: the complex datatype SHOULD be [is] defined dynamically at runtime in Startup.c diff --git a/Carpet/CarpetIOHDF5/src/iohdf5chckpt_recover.cc b/Carpet/CarpetIOHDF5/src/iohdf5chckpt_recover.cc index e692c7701..4f95e0c51 100644 --- a/Carpet/CarpetIOHDF5/src/iohdf5chckpt_recover.cc +++ b/Carpet/CarpetIOHDF5/src/iohdf5chckpt_recover.cc @@ -51,9 +51,9 @@ namespace CarpetIOHDF5 { int Checkpoint (const cGH* const cctkGH, int called_from); int DumpParametersGHExtentions (const cGH *cctkGH, int all, hid_t writer); - // int RecoverParameters (IObase* reader); - //int RecoverGHextensions (cGH* cgh, IObase* reader); - //int RecoverVariables (cGH* cgh, IObase* reader, AmrGridReader* amrreader); + int RecoverParameters (hid_t reader); + int RecoverGHextensions (cGH* cctkGH, hid_t reader); + int RecoverVariables (cGH* cctkGH, hid_t reader); void CarpetIOHDF5_EvolutionCheckpoint( const cGH* const cgh){ @@ -92,8 +92,10 @@ namespace CarpetIOHDF5 { int result,myproc; CarpetIOHDF5GH *myGH; - char filename[1024]; + + static hid_t reader; //this thing absolutely needs to be static!!! + myGH = NULL; result = 0; @@ -101,29 +103,304 @@ namespace CarpetIOHDF5 { if (called_from == CP_RECOVER_PARAMETERS) { - CCTK_WARN (-1,"Sorry, this feature is not implemented yet."); - } + // Okay, let's see what we can do about the parameters + + // Invent a file name + ostringstream filenamebuf; + + if(CCTK_nProcs(cctkGH) == 1) + filenamebuf << recover_dir << "/" << basefilename << ".h5"; + else + filenamebuf << recover_dir << "/" << basefilename << ".file_0.h5"; + + string filenamestr = filenamebuf.str(); + const char * const filename = filenamestr.c_str(); + + if (myproc == 0) { + // First, open the file + if (verbose) + CCTK_VInfo(CCTK_THORNSTRING, "Opening Checkpoint file %s for recovery",filename); + reader = H5Fopen (filename, H5F_ACC_RDONLY, H5P_DEFAULT); + if (reader<0) { + CCTK_VWarn (0, __LINE__, __FILE__, CCTK_THORNSTRING, + "Could not recover from \"%s\"", filename); + } + assert (reader>=0); + } // myproc == 0 + } else { /* This is the case for CP_RECOVER_DATA. CCTK_RECOVER_PARAMETERS must have been called before and set up the file info structure. */ if (myproc == 0) { - CCTK_WARN (-1,"Sorry, this feature is not implemented yet."); + assert(reader>=0); } } + if (called_from == CP_RECOVER_PARAMETERS) + { + return (RecoverParameters (reader)); + } if (called_from == CP_RECOVER_DATA) { + CCTK_INT4 numberofmgtimes; + + if (myproc == 0) { + + /* we need all the times on the individual levels */ + // these are a bit messy to extract + + hid_t group = H5Gopen (reader, PARAMETERS_GLOBAL_ATTRIBUTES_GROUP); + assert(group>=0); + hid_t dataset = H5Dopen (group, ALL_PARAMETERS); + assert(dataset >= 0); + hid_t attr = H5Aopen_name (dataset, "numberofmgtimes"); + assert(attr >= 0); + hid_t atype = H5Aget_type (attr); + if(H5Tequal(atype, H5T_NATIVE_INT)) { + herr_t herr = H5Aread(attr, atype, &numberofmgtimes); + assert(!herr); + herr = H5Aclose(attr); + assert(numberofmgtimes==mglevels); + char buffer[100]; + for(int lcv=0;lcv<numberofmgtimes;lcv++) { + sprintf(buffer,"mgleveltimes %d",lcv); + attr = H5Aopen_name(dataset, buffer); + assert (attr>=0); + atype = H5Aget_type (attr); + assert (atype>=0); + herr = H5Aread (attr, atype, &leveltimes[lcv][0]); + assert(!herr); + herr = H5Aclose(attr); + assert(!herr); + } + herr = H5Dclose(dataset); + assert(!herr); + herr = H5Gclose(group); + assert(!herr); + } else { + CCTK_WARN(0,"BAD BAD BAD! Can't read leveltimes!!"); + } + } // myproc == 0 + + int mpierr = MPI_Bcast (&numberofmgtimes, 1, CARPET_MPI_INT4, 0,MPI_COMM_WORLD); + assert(!mpierr); + + + for(int i=0;i<numberofmgtimes;i++) { + mpierr = MPI_Bcast (&leveltimes[i][0], numberofmgtimes, CARPET_MPI_REAL, 0, MPI_COMM_WORLD); + assert(!mpierr); + } + + cout << "leveltimes: " << leveltimes << endl; + + + BEGIN_MGLEVEL_LOOP(cctkGH) { + BEGIN_REFLEVEL_LOOP(cctkGH) { + + // set tt (ask Erik why...) + // arrdata[0][0].tt->set_time(reflevel,mglevel,(CCTK_REAL) cctkGH->cctk_iteration/maxreflevelfact); + + // Now let us recover the GHextentions + result += RecoverGHextensions(cctkGH,reader); + + result += RecoverVariables (cctkGH,reader); + + } END_REFLEVEL_LOOP; + + } END_MGLEVEL_LOOP; + + CCTK_VInfo (CCTK_THORNSTRING, "Restarting simulation at iteration %d (physical time %g)", cctkGH->cctk_iteration, (double) cctkGH->cctk_time); - } + } // called_from == CP_RECOVER_DATA - CCTK_WARN (-1,"STOPSTOPSTOP2"); + if (myproc == 0) + H5Fclose(reader); + + + // CCTK_WARN (-1,"Sorry, this feature is not implemented yet."); return (result); } // CarpetIOHDF5_Recover + int RecoverVariables (cGH* cctkGH, hid_t reader) { + + int retval = 0; + int myproc = CCTK_MyProc (cctkGH); + int currdataset,ndatasets,namelength; + char datasetname[1024]; + char * name; + + hid_t dataset; + herr_t herr; + + CCTK_VInfo(CCTK_THORNSTRING,"Starting to recover data!!!"); + + if(myproc==0) { + ndatasets=GetnDatasets(reader); + assert (ndatasets>=0); + } + + // Broadcast number of datasets + MPI_Bcast (&ndatasets, 1, MPI_INT, 0, dist::comm); + assert (ndatasets>=0); + + for (currdataset=0;currdataset < ndatasets+1;currdataset++) { + + if (myproc==0) { + GetDatasetName(reader,currdataset,datasetname); + dataset = H5Dopen (reader, datasetname); + assert(dataset); + // Read data + ReadAttribute (dataset, "name", name); + namelength=strlen(name)+1; + } + MPI_Bcast (&namelength, 1, MPI_INT, 0, dist::comm); + + + if(myproc!=0) { + name = (char *) malloc (sizeof(char)*namelength); + } + + MPI_Bcast (name,namelength, MPI_CHAR, 0, dist::comm); + + cout << name << endl; + vector<ibset> regions_read(Carpet::maps); + + ReadVar(cctkGH,reader,name,dataset,regions_read,0); + + if(myproc==0) { + herr = H5Dclose(dataset); + assert(!herr); + } + } + return retval; + } + + + + + int RecoverGHextensions (cGH *cctkGH, hid_t reader) + { + const int myproc = CCTK_MyProc(cctkGH); + CCTK_INT4 int4Buffer[2]; + CCTK_REAL realBuffer; + CCTK_REAL realBuffer2; + CCTK_INT4 intbuffer; + + int mpierr = 0; + + if (myproc==0) + { + + // First open group and dataset + hid_t group = H5Gopen (reader, PARAMETERS_GLOBAL_ATTRIBUTES_GROUP); + assert(group>=0); + hid_t dataset = H5Dopen (group, ALL_PARAMETERS); + assert(dataset >= 0); + + ReadAttribute(dataset,"GH$iteration",int4Buffer[0]); + ReadAttribute(dataset,"main loop index",int4Buffer[1]); + ReadAttribute(dataset,"carpet_global_time",realBuffer); + + herr_t herr = H5Dclose(dataset); + assert(!herr); + herr = H5Gclose(group); + assert(!herr); + + } + /* Broadcast the GH extensions to all processors */ + /* NOTE: We have to use MPI_COMM_WORLD here + because PUGH_COMM_WORLD is not yet set up at parameter recovery time. + We also assume that PUGH_MPI_INT4 is a compile-time defined datatype. */ + + mpierr = MPI_Bcast (int4Buffer, 2, CARPET_MPI_INT4, 0,MPI_COMM_WORLD); + assert(!mpierr); + mpierr = MPI_Bcast (int4Buffer, 2, CARPET_MPI_INT4, 0,MPI_COMM_WORLD); + assert(!mpierr); + mpierr = MPI_Bcast (&realBuffer, 1, CARPET_MPI_REAL,0,MPI_COMM_WORLD); + assert(!mpierr); + + global_time = (CCTK_REAL) realBuffer; + cctkGH->cctk_iteration = (int) int4Buffer[0]; + CCTK_SetMainLoopIndex ((int) int4Buffer[1]); + + + return (0); + + } // RecoverGHExtensions + + int RecoverParameters(hid_t reader){ + + DECLARE_CCTK_PARAMETERS; + + int myproc, retval; + char *parameters; + CCTK_INT4 parameterSize; + + hid_t group,dataset; + herr_t herr; + + int mpierr; + + myproc = CCTK_MyProc (NULL); + + if (myproc == 0){ + CCTK_VInfo (CCTK_THORNSTRING, "Recovering parameters from checkpoint "); + + group = H5Gopen (reader, PARAMETERS_GLOBAL_ATTRIBUTES_GROUP); + assert(group >= 0); + dataset = H5Dopen (group, ALL_PARAMETERS); + assert(dataset>= 0); + + parameterSize = H5Dget_storage_size (dataset); + parameters = (char *) malloc (parameterSize); + herr = H5Dread (dataset, H5T_NATIVE_CHAR, H5S_ALL, H5S_ALL, + H5P_DEFAULT, parameters); + assert(!herr); + herr = H5Dclose(dataset); + assert(!herr); + herr = H5Gclose(group); + assert(!herr); + + if(verbose) + CCTK_VInfo (CCTK_THORNSTRING, "\n%s\n",parameters); + + CCTK_VInfo(CCTK_THORNSTRING, "Successfully recovered parameters!"); + } // myproc == 0 + + /* Broadcast the parameter buffer size to all processors */ + /* NOTE: We have to use MPI_COMM_WORLD here + because CARPET_COMM_WORLD is not yet set up at parameter recovery time. + We also assume that CARPET_MPI_INT4 is a compile-time defined datatype. */ + mpierr = MPI_Bcast (¶meterSize, 1, CARPET_MPI_INT4, 0, + MPI_COMM_WORLD); + assert(!mpierr); + + if (parameterSize > 0) { + if (myproc) { + parameters = (char*) malloc (parameterSize + 1); + } + + mpierr = MPI_Bcast (parameters, parameterSize + 1, CARPET_MPI_CHAR, + 0,MPI_COMM_WORLD); + assert(!mpierr); + + IOUtil_SetAllParameters (parameters); + + free (parameters); + } + + /* return positive value for success otherwise negative */ + retval = (parameterSize > 0 ? 1 : -1); + + return (retval); + + } // RecoverParameters + + int Checkpoint (const cGH* const cctkGH, int called_from) { @@ -392,19 +669,22 @@ namespace CarpetIOHDF5 { dtmp = cctkGH->cctk_time; WriteAttribute(dataset,"GH$time", dtmp); + dtmp = global_time; + WriteAttribute(dataset,"carpet_global_time", dtmp); + version = CCTK_FullVersion(); WriteAttribute(dataset,"Cactus version", version); /* finally, we need all the times on the individual refinement levels */ - const int numberoftimes=leveltimes[0].size(); - WriteAttribute(dataset,"numberoftimes",numberoftimes); - for(int i=0;i < numberoftimes;i++) { + + const int numberofmgtimes=mglevels; + WriteAttribute(dataset,"numberofmgtimes",numberofmgtimes); + for(int i=0;i < numberofmgtimes;i++) { char buffer[100]; - sprintf(buffer,"leveltime%d",i); - WriteAttribute(dataset,buffer,leveltimes[0][i]); + sprintf(buffer,"mgleveltimes %d",i); + WriteAttribute(dataset,buffer,(double *) &leveltimes[i][0], reflevels); } - herr = H5Dclose (dataset); assert(!herr); herr = H5Sclose (dataspace); |