diff options
Diffstat (limited to 'CarpetDev/CarpetIOF5/src/iof5.cc.old')
| -rw-r--r-- | CarpetDev/CarpetIOF5/src/iof5.cc.old | 643 |
1 files changed, 0 insertions, 643 deletions
diff --git a/CarpetDev/CarpetIOF5/src/iof5.cc.old b/CarpetDev/CarpetIOF5/src/iof5.cc.old deleted file mode 100644 index fe780030b..000000000 --- a/CarpetDev/CarpetIOF5/src/iof5.cc.old +++ /dev/null @@ -1,643 +0,0 @@ -#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 |