path: root/CarpetDev/CarpetIOF5/src/iof5.cc.old

#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