diff options
Diffstat (limited to 'CarpetDev/CarpetIOF5/src/iof5.cc.old')
| -rw-r--r-- | CarpetDev/CarpetIOF5/src/iof5.cc.old | 643 |
1 files changed, 643 insertions, 0 deletions
diff --git a/CarpetDev/CarpetIOF5/src/iof5.cc.old b/CarpetDev/CarpetIOF5/src/iof5.cc.old new file mode 100644 index 000000000..fe780030b --- /dev/null +++ b/CarpetDev/CarpetIOF5/src/iof5.cc.old @@ -0,0 +1,643 @@ +#include <cctk.h> +#include <cctk_Arguments.h> +#include <cctk_Parameters.h> + +#include <cassert> +#include <cstdlib> +#include <cstring> +#include <iomanip> +#include <iostream> +#include <limits> +#include <sstream> +#include <string> + +#include <hdf5.h> +#include <F5/F5F.h> +#include <F5/F5R.h> +#include <F5/F5iterate.h> +#include <F5/F5uniform.h> + +#include <bbox.hh> +#include <vect.hh> + +#include <carpet.hh> + + + +namespace CarpetIOF5 { + + using namespace std; + using namespace Carpet; + + + + // File mode for creating directories + int const mode = 0755; + + // A nan + CCTK_REAL const nan = numeric_limits<CCTK_REAL>::quiet_NaN(); + + // Indentation + inline string indent (int const n) + { + int const width = 3; + return string(width*n, ' '); + } + + + + typedef vect<hsize_t,dim> hvect; + typedef vect<float,dim> fvect; + + + + static inline + hvect v2h (ivect const& a) + { + return hvect(a); + } + + static inline + F5_vec3_point_t v2p (rvect const& a) + { + F5_vec3_point_t res; + res.x = a[0]; + res.y = a[1]; + res.z = a[2]; + return res; + } + + static inline + F5_vec3_float_t v2f (rvect const& a) + { + F5_vec3_float_t res; + res.x = a[0]; + res.y = a[1]; + res.z = a[2]; + return res; + } + + static inline + F5_vec3_double_t v2d (rvect const& a) + { + F5_vec3_double_t res; + res.x = a[0]; + res.y = a[1]; + res.z = a[2]; + return res; + } + + + + // Generate a good grid name (simulation name) + string + generate_gridname (cGH const *const cctkGH) + { +#if 0 + char const *const gridname = (char const*) UniqueSimulationID(cctkGH); + assert (gridname); + return gridname; +#endif + // Use the parameter file name, since the simulation ID is too + // long and doesn't look nice + char parfilename[10000]; + CCTK_ParameterFilename (sizeof parfilename, parfilename); + char *const p = strstr(parfilename, ".par"); + if (p) *p = '\0'; + char *const s = strrchr(parfilename, '/'); + char *const gridname = s ? s+1 : parfilename; + return gridname; + } + + + + // Generate a good file name ("alias") for a variable + string + generate_basename (cGH const *const cctkGH, + int const variable) + { + DECLARE_CCTK_PARAMETERS; + + assert (variable >= 0); + + ostringstream filename_buf; + + if (CCTK_EQUALS (file_content, "variable")) { + char *const varname = CCTK_FullName(variable); + assert (varname); + for (char *p = varname; *p; ++p) { + *p = tolower(*p); + } + filename_buf << varname; + free (varname); + } + else if (CCTK_EQUALS (file_content, "group")) { + char *const groupname = CCTK_GroupNameFromVarI(variable); + assert (groupname); + for (char *p = groupname; *p; ++p) { + *p = tolower(*p); + } + filename_buf << groupname; + free (groupname); + } + else if (CCTK_EQUALS (file_content, "thorn")) { + char const *const impname = CCTK_ImpFromVarI(variable); + char *const thornname = strdup(impname); + assert (thornname); + char *const colon = strchr(thornname, ':'); + assert (colon); + *colon = '\0'; + for (char *p = thornname; *p; ++p) { + *p = tolower(*p); + } + filename_buf << thornname; + free (thornname); + } + else if (CCTK_EQUALS (file_content, "everything")) { + filename_buf << out_filename; + } + else { + assert (0); + } + + if (out_timesteps_per_file > 0) { + int const iteration = (cctkGH->cctk_iteration + / out_timesteps_per_file * out_timesteps_per_file); + filename_buf << ".it" + << setw (iteration_digits) << setfill ('0') << iteration; + } + + return filename_buf.str(); + } + + + + // Create the final file name on a particular processor + string + create_filename (cGH const *const cctkGH, + string const basename, + bool const create_directories) + { + DECLARE_CCTK_PARAMETERS; + + int const proc = CCTK_MyProc(cctkGH); + int const nprocs = CCTK_nProcs(cctkGH); + + bool const use_IO_out_dir = strcmp(out_dir, "") == 0; + string path = use_IO_out_dir ? IO_out_dir : out_dir; + + if (create_subdirs) { + if (nprocs >= 10000) { + ostringstream buf; + buf << path + "/" + << basename + << ".p" << setw (max (0, processor_digits - 4)) << setfill ('0') + << proc / 10000 + << "nnnn/"; + path = buf.str(); + if (create_directories) { + check (CCTK_CreateDirectory (mode, path.c_str()) >= 0); + } + } + if (nprocs >= 100) { + ostringstream buf; + buf << path + "/" + << basename + << ".p" << setw (max (0, processor_digits - 2)) << setfill ('0') + << proc / 100 + << "nn/"; + path = buf.str(); + if (create_directories) { + check (CCTK_CreateDirectory (mode, path.c_str()) >= 0); + } + } + } + if (one_dir_per_file and nprocs >= 1) { + ostringstream buf; + buf << path + << basename + << ".p" << setw (processor_digits) << setfill ('0') + << proc + << "/"; + path = buf.str(); + if (create_directories) { + check (CCTK_CreateDirectory (mode, path.c_str()) >= 0); + } + } + ostringstream buf; + buf << path + "/" + << basename + << ".p" << setw (processor_digits) << setfill ('0') << proc + << out_extension; + return buf.str(); + } + + + + extern "C" + void F5_Output (CCTK_ARGUMENTS) + { + DECLARE_CCTK_ARGUMENTS; + DECLARE_CCTK_PARAMETERS; + + herr_t herr; + + + + assert (is_global_mode()); + CCTK_VInfo (CCTK_THORNSTRING, "F5_Output: iteration=%d", cctk_iteration); + + string const gridname = generate_gridname(cctkGH); + + + + // Open file + static bool first_time = true; + + string const basename = + generate_basename (cctkGH, CCTK_VarIndex("grid::r")); + string const name = + create_filename (cctkGH, basename, first_time); + + bool const truncate_file = first_time and IO_TruncateOutputFiles(cctkGH); + hid_t const file = + truncate_file ? + H5Fcreate (name.c_str(), H5F_ACC_TRUNC, H5P_DEFAULT, H5P_DEFAULT) : + H5Fopen (name.c_str(), H5F_ACC_RDWR, H5P_DEFAULT); + assert (file >= 0); + first_time = false; + + + + BEGIN_REFLEVEL_LOOP (cctkGH) { + DECLARE_CCTK_ARGUMENTS; + +#if 0 + // Write data +#warning "TODO: Save the coordinates in double precision" + F5Path *const path = + F5Fwrite_uniform_cartesian3D + (file, cctk_time, gridname, &v2p(lower), &v2f(delta), &v2h(lsh)[0], + "radius" /* group name */, H5T_NATIVE_DOUBLE, r, NULL, F5P_DEFAULT); +#endif + + // Define grid hierarchy + ivect const reffact = spacereffacts.AT(reflevel); + F5Path *const path = F5Rcreate_vertex_refinement3D + (file, cctk_time, gridname.c_str(), &v2h(reffact)[0], NULL); + + // Define iteration + F5Rset_timestep (path, cctk_iteration); + +#if 0 + F5Path *refpath = NULL; + static CCTK_REAL last_root_time = -1; + if (reflevel == 0) { + last_root_time = cctk_time; + } else { + refpath = F5Rcreate_relative_vertex_Qrefinement3D + (file, cctk_time, gridname, &v2h(reffact)[0], + last_root_time, &v2h(ivect(1))[0]); + } +#endif + + + + BEGIN_LOCAL_MAP_LOOP (cctkGH, CCTK_GF) { + DECLARE_CCTK_ARGUMENTS; + + // Determine level coordinates + ivect gsh; + rvect origin, delta; + rvect lower, upper; + for (int d=0; d<dim; ++d) { + gsh[d] = cctk_gsh[d]; + origin[d] = CCTK_ORIGIN_SPACE(d); + delta[d] = CCTK_DELTA_SPACE(d); + lower[d] = origin[d]; + upper[d] = lower[d] + (gsh[d]-1) * delta[d]; + } + + // Define level topology + if (reflevel == 0) { + // Choose (arbitrarily) the root level as default topology, + // for readers which don't understand AMR + F5Rlink_default_vertex_topology (path, &v2h(reffact)[0]); + } + + // Define level geometry + F5Fwrite_linear (path, FIBER_HDF5_POSITIONS_STRING, + dim, &v2h(gsh)[0], + F5T_COORD3_DOUBLE, &v2d(lower), &v2d(delta)); + F5Fset_range (path, &v2d(lower), &v2d(upper)); + + + + BEGIN_LOCAL_COMPONENT_LOOP (cctkGH, CCTK_GF) { + DECLARE_CCTK_ARGUMENTS; + + // Determine component coordinates + ivect lbnd, lsh, lghosts, ughosts; + rvect clower, cupper; +#if 0 + F5_refinement3D_point_t iorigin, idelta; +#endif + for (int d=0; d<dim; ++d) { + lbnd[d] = cctk_lbnd[d]; + lsh[d] = cctk_lsh[d]; + // F5 counts the total overlap, which is the sum of the + // ghost zones on this and the adjacent component + lghosts[d] = cctk_bbox[2*d ] ? 0 : 2*cctk_nghostzones[d]; + ughosts[d] = cctk_bbox[2*d+1] ? 0 : 2*cctk_nghostzones[d]; + clower[d] = origin[d] + cctk_lbnd[d] * delta[d]; + cupper[d] = clower[d] + (lsh[d]-1) * delta[d]; +#if 0 + iorigin.d[d].num = + cctk_levoff[d] + cctk_lbnd[d] * cctk_levoffdenom[d]; + iorigin.d[d].denom = cctk_levoffdenom[d] * reffact[d]; + idelta .d[d].num = 1; + idelta .d[d].denom = reffact[d]; +#endif + } + ostringstream buf2; + buf2 << "component=" << component; + string const cname = buf2.str(); + +#if 0 + if (reflevel > 0) { + F5Fwrite_linear_fraction (refpath, FIBER_HDF5_POSITIONS_STRING, + cctk_dim, &v2h(gsh)[0], &v2h(lsh)[0], + F5T_COORD3_INT_FRACTION32, + &iorigin, &idelta, &v2h(lbnd)[0], + cname.c_str()); + } +#endif + + // Define coordinates + // (This is redundant, since the level's overall bounding + // box was already defined above, but it provides the + // individual components' bounding boxes.) + F5Fwrite_linear_fraction (path, FIBER_HDF5_POSITIONS_STRING, + cctk_dim, &v2h(gsh)[0], &v2h(lsh)[0], + F5T_COORD3_DOUBLE, + &clower, &delta, &v2h(lbnd)[0], + cname.c_str()); + + // Write data + // scalar field + F5Fwrite_fraction (path, "radius" /* group name */, + cctk_dim, &v2h(gsh)[0], &v2h(lsh)[0], + H5T_NATIVE_DOUBLE, H5T_NATIVE_DOUBLE, + r, + &v2h(lbnd)[0], &v2h(lghosts)[0], &v2h(ughosts)[0], + cname.c_str(), H5P_DEFAULT); + + // vector field + void const* data[dim]; + data[0] = x; + data[1] = y; + data[2] = z; + F5FSwrite_fraction (path, "coordinates" /* group name */, + cctk_dim, &v2h(gsh)[0], &v2h(lsh)[0], + F5T_VEC3_DOUBLE, F5T_VEC3_DOUBLE, + data, + &v2h(lbnd)[0], &v2h(lghosts)[0], &v2h(ughosts)[0], + cname.c_str(), H5P_DEFAULT, 0); + + } END_LOCAL_COMPONENT_LOOP; + + } END_LOCAL_MAP_LOOP; + + // Close topology + F5close (path); + + } END_REFLEVEL_LOOP; + + // Close file + herr = H5Fclose (file); + assert (not herr); + } + + + + // Use a class for reading data, so that we have an easy way to pass + // variables between the various iterators + class iterator_t { + cGH *cctkGH; + double time; + char const *gridname; + char const *topologyname; + int index_depth; // 0=vertex, 1=cell + int topological_dimension; + char const *fieldname; + char const *fragmentname; + + public: + iterator_t (cGH *const cctkGH_) + : cctkGH(cctkGH_), + time(nan), + gridname(NULL), + topologyname(NULL), index_depth(-1), topological_dimension(-1), + fieldname(NULL), + fragmentname(NULL) + { + } + + + + private: + + void read_timeslice (F5Path *const path) + { + cout << indent(1) << "time=" << time << "\n"; + F5iterate_grids (path, NULL, grid_iterator, this, NULL, NULL); + } + + void read_grid (F5Path *const path) + { + cout << indent(2) << "grid=\"" << gridname << "\"\n"; + // F5iterate_vertex_fields (path, NULL, field_iterator, this, NULL, NULL); + F5iterate_topologies (path, NULL, topology_iterator, this); + } + + void read_topology (F5Path *const path) + { + herr_t herr; + + cout << indent(3) << "topology=\"" << topologyname << "\"" + << " (" << (index_depth==0 ? "vertex" : "cell") << ")\n"; + + // Ignore topologies that are only an alias for another topology + H5G_stat_t stat; + herr = H5Gget_objinfo (path->Grid_hid, topologyname, false, &stat); + assert (not herr); + if (stat.type == H5G_LINK) { + char linkval[100000]; + herr = H5Gget_linkval + (path->Grid_hid, topologyname, sizeof linkval, linkval); + assert (not herr); + cout << indent(4) << "alias for topology \"" << linkval << "\"\n" + << indent(4) << "ignoring this topology\n"; + return; + } + + F5iterate_topology_fields (path, NULL, field_iterator, this, NULL, NULL); + } + + void read_field (F5Path *const path) + { + cout << indent(4) << "field=\"" << fieldname << "\"\n"; + F5iterate_field_fragments (path, NULL, fragment_iterator, this); + } + + void read_fragment (F5Path *const path) + { + herr_t herr; + + cout << indent(5) << "fragment=\"" << fragmentname << "\"\n"; + + if (strcmp(fieldname, "radius") == 0) { + + hid_t const group = path->Field_hid; + read_variable (group, fragmentname, CCTK_VarIndex("grid::r")); + + } else if (strcmp(fieldname, "coordinates") == 0) { + + hid_t const group = + H5Gopen (path->Field_hid, fragmentname, H5P_DEFAULT); + assert (group); + read_variable (group, "x", CCTK_VarIndex("grid::x")); + read_variable (group, "y", CCTK_VarIndex("grid::y")); + read_variable (group, "z", CCTK_VarIndex("grid::z")); + herr = H5Gclose (group); + assert (not herr); + + } else { + // do nothing + } + } + + void read_variable (hid_t const group, char const *const name, + int const var) + { + assert (group >= 0); + assert (name); + assert (var >= 0); + + cout << indent(6) << "dataset=\"" << name << "\"\n"; + // hid_t const fragment = H5Dopen(group, name); + // assert (fragment); + // + // hid_t const space = H5Dget_space(fragment); + // assert (space); + // int const rank = H5Sget_simple_extent_dims(space, NULL, NULL); + // assert (rank>=0); + // + // assert (rank == cctk_dim); + // + // H5Sclose(space_id); + } + + + + public: + + void iterate (hid_t const object) + { + F5iterate_timeslices (object, NULL, timeslice_iterator, this); + } + + static + herr_t timeslice_iterator (F5Path *const path, double const time, + void *const userdata) + { + iterator_t* const iterator = (iterator_t*)userdata; + iterator->time = time; + iterator->read_timeslice (path); + return 0; + } + + static + herr_t grid_iterator (F5Path *const path, char const *const gridname, + void *const userdata) + { + iterator_t* const iterator = (iterator_t*)userdata; + iterator->gridname = gridname; + iterator->read_grid (path); + return 0; + } + + static + herr_t topology_iterator (F5Path *const path, + char const *const topologyname, + int const index_depth, + int const topological_dimension, + void *const userdata) + { + iterator_t* const iterator = (iterator_t*)userdata; + iterator->topologyname = topologyname; + iterator->index_depth = index_depth; + iterator->topological_dimension = topological_dimension; + iterator->read_topology (path); + return 0; + } + + static + herr_t field_iterator (F5Path *const path, char const *const fieldname, + void *const userdata) + { + iterator_t* const iterator = (iterator_t*)userdata; + iterator->fieldname = fieldname; + iterator->read_field (path); + return 0; + } + + static + herr_t fragment_iterator (F5Path *const path, + char const *const fragmentname, + void *const userdata) + { + iterator_t* const iterator = (iterator_t*)userdata; + iterator->fragmentname = fragmentname; + iterator->read_fragment (path); + return 0; + } + + }; + + + + extern "C" + void F5_Input (CCTK_ARGUMENTS) + { + DECLARE_CCTK_ARGUMENTS; + DECLARE_CCTK_PARAMETERS; + + herr_t herr; + + + + assert (is_global_mode()); + CCTK_VInfo (CCTK_THORNSTRING, "F5_Input: iteration=%d", cctk_iteration); + + + + // Open file + string const basename = + generate_basename (cctkGH, CCTK_VarIndex("grid::r")); + string const name = + create_filename (cctkGH, basename, false); + + hid_t const file = H5Fopen (name.c_str(), H5F_ACC_RDONLY, H5P_DEFAULT); + assert (file >= 0); + + // Iterate over all time slices + iterator_t iterator(cctkGH); + iterator.iterate (file); + + // Close file + herr = H5Fclose (file); + assert (not herr); + } + +} // end namespace CarpetIOF5 |