hdf5 writing working so far...

py/src/test_py_unsio.py

@@ -201,8 +201,10 @@ def process(args):
   
   ok,insnap=readSnap(args.simname,args.component,args.float)
   ok,svsnap=saveSnap(insnap,args.component,"gadget2",args.float)
+  ok,hdf5snap=saveSnap(insnap,args.component,"gadget3",args.float)
   ok,nemosnap=saveSnap(insnap,args.component,"nemo",args.float)
   compare(insnap,svsnap)
+  compare(insnap,hdf5snap)
   compare(insnap,nemosnap)
 
 # -----------------------------------------------------

src/snapshotgadget.cc

@@ -1748,13 +1748,8 @@ int CSnapshotGadgetOut<T>::setData(std::string name,std::string array,  const in
   bool ok=true;
   int status=0;
   
-  switch(CunsOut2<T>::s_mapStringValues[name]) {
-  case uns::Gas   :
-  case uns::Halo  :
-  case uns::Disk  :
-  case uns::Bulge :
-  case uns::Stars :
-  case uns::Bndry :
+  switch(CunsOut2<T>::s_mapStringValues[array]) {
+  case uns::Id   :
     status = setId(name, n, data, _addr);
     break;
   default: ok=false;

src/snapshotgadgeth5.cc

@@ -17,6 +17,7 @@
 #include "unsdebug.h"
 #include "uns.h"
 #include <sstream>
+#include <limits>
 
 namespace  uns {
 
@@ -540,6 +541,390 @@ bool CSnapshotGadgetH5In<T>::loadCommonDataset(std::string tag, std::vector<U> &
   }
   return ok;
 }
+// ----------------------------------------------------------------------------
+//
+//                     CSnapshotGadgetH5Out CLASS implementation
+//
+// ----------------------------------------------------------------------------
+
+// ----------------------------------------------------------------------------
+// WRITING constructor
+template <class T>
+CSnapshotGadgetH5Out<T>::CSnapshotGadgetH5Out(const std::string _n, const std::string _t, const bool _v):CSnapshotInterfaceOut<T>(_n, _t, _v)
+{
+  myH5 = NULL;
+  try {
+    /*
+     * Turn off the auto-printing when failure occurs so that we can
+     * handle the errors appropriately
+     */
+    Exception::dontPrint();
+    myH5 = new GH5<T>(this->simname,H5F_ACC_TRUNC,this->verbose);
+    this->interface_type = "Gadget3";
+    this->file_structure = "component"; // "component" like file
+    if (this->verbose)
+      std::cerr << "CSnapshotGadgetH5Out::CSnapshotGadgetH5Out simname = " << this->simname <<"\n";
+    bzero(&header,sizeof(t_h5_header));
+
+    // setup header
+    header.MassTable.resize(6);
+    header.NumFilesPerSnapshot=1;
+    header.NumPart_ThisFile.resize(6);
+    header.NumPart_Total.resize(6);
+    header.NumPart_Total_HighWord.resize(6);
+
+   // detect precision
+   if ((std::numeric_limits<T>::max()) == (std::numeric_limits<double>::max())) {
+      header.Flag_DoublePrecision=1;
+    };
+
+  } catch (...) {
+    if (this->verbose) {
+      std::cerr << "WARNING: impossible to create Gadget3 (HDF5) output file...\n";
+    }
+    throw(-1);
+    std::exit(1);
+  }
+}
+// ----------------------------------------------------------------------------
+// destructor
+template <class T>
+CSnapshotGadgetH5Out<T>::~CSnapshotGadgetH5Out()
+{
+}
+// ----------------------------------------------------------------------------
+// putHeader
+template <class T>
+int CSnapshotGadgetH5Out<T>::setHeader(void * _header)
+{
+  memcpy(&header,_header,sizeof(t_h5_header));
+  return 1;
+}
+// ============================================================================
+// setNbody:
+template <class T>
+int CSnapshotGadgetH5Out<T>::setNbody(const int _nbody)
+{
+  return 1;
+}
+// ----------------------------------------------------------------------------
+// setData
+template <class T>
+int CSnapshotGadgetH5Out<T>::setData(std::string name,T  data)
+{
+  int status=0;
+
+  switch(CunsOut2<T>::s_mapStringValues[name]) {
+  case uns::Time :
+    status = 1;
+    header.Time = data;
+    //status=myH5->setAttribute("Time",&header.Time,1);
+    break;
+  default: status=0;
+  }
+
+  if (this->verbose) {
+    if (status) {
+      std::cerr << "CSnapshotGadgetH5Out::setData name["<<name<<"]=" << CunsOut2<T>::s_mapStringValues[name] << "\n";
+    } else {
+      std::cerr << "** WARNING ** SnapshotGadgetH5Out::setData Value ["<<name<<"] does not exist.....\n";
+    }
+  }
+  return status;
+}
+// ----------------------------------------------------------------------------
+// setData
+// setData("gas","pos",n,gas_array,true)
+template <class T>
+int CSnapshotGadgetH5Out<T>::setData(std::string name,std::string array,  const int n ,T * data,const bool _addr)
+{
+  int status=0;
+  try {
+    switch(CunsOut2<T>::s_mapStringValues[array]) {
+    case uns::Pos  :
+      status = saveCommonDataset(name,"Coordinates",n,data,3);
+      break;
+    case uns::Vel  :
+      status = saveCommonDataset(name,"Velocities",n,data,3);
+      break;
+    case uns::Mass :
+      status = saveCommonDataset(name,"Masses",n,data,1);
+      break;
+    case uns::Pot  :
+      status = saveCommonDataset(name,"Potential",n,data,1);
+      break;
+    case uns::Acc  :
+      status = saveCommonDataset(name,"Acceleration",n,data,3);
+      break;
+    case uns::Hsml :
+      status = saveCommonDataset(name,"SmoothingLength",n,data,1);
+      break;
+    case uns::Rho  :
+      status = saveCommonDataset(name,"Density",n,data,1);
+      break;
+    case uns::U  :
+      if (name=="gas") {
+        status = saveCommonDataset(name,"InternalEnergy",n,data,1);
+      }
+      break;
+    case uns::Age  :
+      if (name=="stars") {
+        status = saveCommonDataset(name,"StellarFormationTime",n,data,1);
+      }
+      break;
+    case uns::GasMetal  :
+      if (name=="gas")
+        status = saveCommonDataset(name,"Metallicity",n,data,1);
+      break;
+    case uns::StarsMetal  :
+      if (name=="stars")
+      status = saveCommonDataset(name,"Metallicity",n,data,1);
+      break;
+    case uns::Metal  :
+      if (name=="gas" || name=="stars") {
+        status = saveCommonDataset(name,"Metallicity",n,data,1);
+      }
+      break;
+    default: status=0;
+    }
+  } catch (...) {
+    std::cerr << "Error during CSnapshotGadgetH5Out<T>::setData(..)\n";
+  }
+  if (this->verbose) {
+    if (status) {
+      std::cerr << "CSnapshotGadgetH5Out::setData name["<<name<<"]=" << CunsOut2<T>::s_mapStringValues[name] << "\n";
+    } else {
+      if (name != "EXTRA") {
+        std::cerr << "** WARNING ** CSnapshotGadgetH5Out::setData  Value ["<<name<<"]=<"<< array<<"] does not exist.....\n";
+      } else {
+        std::cerr << "CSnapshotGadgetH5Out::setData EXTRA tags["<<array<<"]\n";
+      }
+    }
+  }
+  return status;
+}
+// ----------------------------------------------------------------------------
+// setData
+template <class T>
+int CSnapshotGadgetH5Out<T>::setData(std::string name, const int n ,T * data,const bool _addr)
+{
+  int status=0;
+  switch(CunsOut2<T>::s_mapStringValues[name]) {
+  case uns::Rho :
+    status = setData("gas",name,n,data,_addr);
+    break;
+  case uns::Hsml:
+    status = setData("gas",name,n,data,_addr);
+    break;
+  case uns::U:
+    status = setData("gas",name,n,data,_addr);
+    break;
+  case uns::Temp  :
+    status = setData("gas",name,n,data,_addr);
+    break;
+  case uns::Age  :
+    status = setData("stars",name,n,data,_addr);
+    break;
+  case uns::GasMetal  :
+    status = setData("gas","metal",n,data,_addr);
+    break;
+  case uns::StarsMetal  :
+    status = setData("stars","metal",n,data,_addr);
+    break;
+  default: status=0;
+  }
+  if (this->verbose) {
+    if (status) { std::cerr << "CSnapshotGadgetH5Out::setData name["<<name<<"]=" << CunsOut2<T>::s_mapStringValues[name] << "\n";
+    } else {
+      std::cerr << "** WARNING ** CSnapshotGadgetH5Out::setData Value ["<<name<<"] does not exist.....\n";
+    }
+  }
+  return status;
+}
+// ----------------------------------------------------------------------------
+// setData
+// setData("all","id",n,gas_array,true)
+template <class T>
+int CSnapshotGadgetH5Out<T>::setData(std::string name,std::string array,  const int n ,int * data,const bool _addr)
+{
+  int status=0;
+  switch(CunsOut2<T>::s_mapStringValues[array]) {
+  case uns::Id   :
+    status = saveCommonDataset(name,"ParticleIDs",n,data,1);
+    break;
+  default: status=false;
+  }
+
+  if (this->verbose) {
+    if (status) {
+      std::cerr << "CSnapshotGadgetH5Out::setData name["<<name<<"]=" << CunsOut2<T>::s_mapStringValues[array] << "\n";
+    } else {
+      std::cerr << "** WARNING ** CSnapshotGadgetH5Out::setData Value ["<<name<<"] does not exist.....\n";
+    }
+  }
+  return status;
+}
+// ----------------------------------------------------------------------------
+// setData
+template <class T>
+int CSnapshotGadgetH5Out<T>::setData(std::string name, const int n ,int * data,const bool _addr)
+{
+  int status=0;
+  return status;
+}
+// ----------------------------------------------------------------------------
+// setData
+// setData("gas",n,pos,vel,mass,true)
+template <class T>
+int CSnapshotGadgetH5Out<T>::setData(std::string name, const int n ,T * data, T * data1, T * data2, const bool _addr)
+{
+  int status=0;
+  return status;
+}
+// ============================================================================
+// save:
+template <class T>
+int CSnapshotGadgetH5Out<T>::save()
+{
+  int fail=0;
+  // write headers
+  myH5->setAttribute("MassTable",&header.MassTable[0], header.MassTable.size());
+  myH5->setAttribute("Time",&header.Time,1);
+  myH5->setAttribute("Redshift",&header.Redshift,1);
+  myH5->setAttribute("BoxSize",&header.BoxSize,1);
+  myH5->setAttribute("Omega0",&header.Omega0,1);
+  myH5->setAttribute("OmegaLambda",&header.OmegaLambda,1);
+  myH5->setAttribute("HubbleParam",&header.HubbleParam,1);
+  myH5->setAttribute("Flag_Cooling",&header.Flag_Cooling,1);
+  myH5->setAttribute("Flag_DoublePrecision",&header.Flag_DoublePrecision,1);
+  myH5->setAttribute("Flag_IC_Info",&header.Flag_IC_Info,1);
+  myH5->setAttribute("Flag_Metals",&header.Flag_Metals,1);
+  myH5->setAttribute("Flag_Sfr",&header.Flag_Sfr,1);
+  myH5->setAttribute("Flag_StellarAge",&header.Flag_StellarAge,1);
+  myH5->setAttribute("NumFilesPerSnapshot",&header.NumFilesPerSnapshot,1);
+  myH5->setAttribute("NumPart_ThisFile",&header.NumPart_ThisFile[0],header.NumPart_ThisFile.size());
+  myH5->setAttribute("NumPart_Total",&header.NumPart_Total[0],header.NumPart_Total.size());
+  myH5->setAttribute("NumPart_Total_HighWord",&header.NumPart_Total_HighWord[0],header.NumPart_Total_HighWord.size());
+
+  return fail;
+}
+// ----------------------------------------------------------------------------
+// checkMasses
+// check out if masses are differents. If not save mass value in attribute MassTable
+// return ok if NOT same masses
+template <class T>
+template <class U>
+bool CSnapshotGadgetH5Out<T>::checkMasses(const int n, U *data, const int comp_id)
+{
+  bool same_mass=true;
+  U massref=data[0];
+  // check if same masses
+  for ( int i=1; i<n;i++) {
+    if (massref != data[i]) { // not same mass
+      same_mass=false;
+      break;
+    }
+  }
+  if (same_mass ) { // we must same unique mass value in Attribute
+    header.MassTable[comp_id]=massref;
+  } else {
+    header.MassTable[comp_id]=0.0;
+  }
+  // save attribute MassTable
+  //myH5->setAttribute("MassTable",&header.MassTable[0], header.MassTable.size());
+  return !same_mass;
+}
+
+// ----------------------------------------------------------------------------
+// saveCommonDataset
+template <class T>
+template <class U>
+bool CSnapshotGadgetH5Out<T>::saveCommonDataset(std::string comp,std::string tag,
+                                                const int n ,U * data, const unsigned int second_dim)
+{
+  bool ok=false;
+
+  std::map<std::string,int> compo_indx;
+  compo_indx["gas"  ]=0;
+  compo_indx["halo" ]=1;
+  compo_indx["dm"   ]=1;
+  compo_indx["disk" ]=2;
+  compo_indx["bulge"]=3;
+  compo_indx["stars"]=4;
+  compo_indx["bndry"]=5;
+
+  try {
+
+    std::map<std::string,int>::iterator it;
+    it = compo_indx.find(comp);
+    if (it != compo_indx.end()) { // key exist
+      bool save=true;
+      if (tag=="Masses") { // special process for masses
+        save=checkMasses(n,data,(*it).second);
+      }
+      if (save) {
+        std::ostringstream myid;
+        myid << (*it).second;
+        std::string dataset="/PartType"+myid.str()+"/"+tag;
+        if (this->verbose)  std::cerr << dataset << "\n";
+        try {
+          Exception::dontPrint();
+          ok = myH5->setDataset(dataset,data,n,second_dim);
+          header.NumPart_ThisFile[(*it).second]=n;
+          header.NumPart_Total[(*it).second]=n;
+          header.NumPart_Total_HighWord[(*it).second]=n;
+//          bool ok2 = myH5->setAttribute("NumPart_ThisFile",
+//                                        &header.NumPart_ThisFile[0],header.NumPart_ThisFile.size());
+//          bool ok3 = myH5->setAttribute("NumPart_Total",
+//                                        &header.NumPart_Total[0],header.NumPart_Total.size());
+//          bool ok4 = myH5->setAttribute("NumPart_Total_HighWord",
+//                                        &header.NumPart_Total_HighWord[0],header.NumPart_Total_HighWord.size());
+
+        }
+        // catch failure caused by the H5File operations
+        catch( FileIException error )
+        {
+           error.printError();
+           throw -1;
+        }
+        // catch failure caused by the DataSet operations
+        catch( DataSetIException error )
+        {
+           error.printError();
+           throw -1;
+        }
+        // catch failure caused by the DataSpace operations
+        catch( DataSpaceIException error )
+        {
+           error.printError();
+           throw -1;
+        }
+        // catch failure caused by the DataSpace operations
+        catch( DataTypeIException error )
+        {
+           error.printError();
+           throw -1;
+        }
+        catch( GroupIException error )
+        {
+          error.printError();
+          return -1;
+        }
+        catch (...) {
+          if (this->verbose) {
+            std::cerr << "WARNING !!! : error while saving Dataset["
+                      << dataset << "]\n";
+          }
+          throw -1;
+        }
+      }
+    }
+  }
+  catch (int e) {
+    std::cerr << "WARNING !!! : enable to saveCommonDataset....\n";
+  }
+  return ok;
+}
 
 // ============================================================================
 // Implementation of class GH5
@@ -548,14 +933,20 @@ bool CSnapshotGadgetH5In<T>::loadCommonDataset(std::string tag, std::vector<U> &
 // ============================================================================
 // Constructor
 template <class T>
-GH5<T>::GH5(const std::string _f_name, int mode, const bool verb)
+GH5<T>::GH5(const std::string _f_name,unsigned int mode, const bool verb)
 {
   verbose = verb;
   f_name=_f_name;
   myfile = NULL;
   myfile = new H5File(f_name,mode);
 
-  readHeaderAttributes();
+  if (mode==H5F_ACC_RDONLY) { // reading mode only
+    readHeaderAttributes();
+  } else
+    if (mode==H5F_ACC_TRUNC) { // read and write mode
+      // we create header group
+      header_group=myfile->createGroup("/Header");
+    }
 }
 // ============================================================================
 // Destructor
@@ -610,7 +1001,7 @@ template <class U> std::vector<U> GH5<T>::getAttribute(std::string attr_name)
 {
   if (verbose) {
     std::cerr << "= = = = = = = = = = = = = = = = = =\n";
-    std::cerr << "Attribute ["<<attr_name<< "]\n";
+    std::cerr << "Read Attribute ["<<attr_name<< "]\n";
   }
   Group     grp  = myfile->openGroup("Header"  );
   Attribute attr = grp.openAttribute(attr_name);
@@ -708,7 +1099,7 @@ template <class U> std::vector<U> GH5<T>::getDataset(std::string dset_name, U du
    default :
     std::cerr << "We should not be here.....\n";
     assert(0);
-  }  
+  }
 
   // read vector of data
   dataset.read(&vret[0],mem_type,H5S_ALL,H5S_ALL,H5P_DEFAULT);
@@ -719,6 +1110,125 @@ template <class U> std::vector<U> GH5<T>::getDataset(std::string dset_name, U du
   return vret;
 
 }
+// ============================================================================
+// setDataset
+// Write the corresponding dataset
+//
+template <class T>
+template <class U>
+bool GH5<T>::setDataset(std::string dset_name, U * data, const unsigned int n,
+                        const unsigned int second_dim)
+{
+  bool ok=true;
+  assert(second_dim==1 || second_dim==3);
+  if (verbose) {
+    std::cerr << "= = = = = = = = = = = = = = = = = =\n";
+    std::cerr << "Set Dataset ["<<dset_name<< "]\n";
+  }
+
+  // first create group
+  std::size_t o1=dset_name.find("/");
+  std::size_t o2=dset_name.find("/",1);
+  if (o1==std::string::npos || o1==std::string::npos ) {
+    std::cerr << "GH5<T>::setDataset no '/' in datasetname....\n";
+    throw -1;
+  }
+  std::string groupname=dset_name.substr(o1,o2-o1); // find groupname
+
+  if (histo_group[groupname]==false) { // check if not created yet
+    Group  group=myfile->createGroup(groupname);
+    histo_group[groupname]=true;
+  }
+  // set the number of dimensions in the dataspace
+  int rank=1;
+
+  hsize_t dims_out[2];
+  dims_out[0] = n;
+  if (second_dim>1) {
+    rank=2;
+    dims_out[1]=second_dim;
+  }
+  if (verbose) {
+    std::cerr << "rank " << rank << "\n";
+  }
+  // set dataspace of the dataset
+  DataSpace dataspace(rank,dims_out);
+  // guess U type
+  DataType mem_type=guessType(U (1));
+  // Create dataset
+  DataSet dataset = myfile->createDataSet(dset_name,mem_type,dataspace);
+  // write datasey
+  dataset.write(data,mem_type);
+
+  if (verbose) {
+    std::cerr << "rank " << rank << ", dimensions " ;
+  }
+  return ok;
+}
+// ============================================================================
+// setAttribute
+// Write the corresponding attribute
+//
+template <class T>
+template <class U>
+bool GH5<T>::setAttribute(std::string attr_name, U * attr_value, const int n)
+{
+  if (verbose) {
+    std::cerr << "= = = = = = = = = = = = = = = = = =\n";
+    std::cerr << "set Attribute ["<<attr_name<< "]\n";
+  }
+  // guess type
+  DataType mem_type=guessType(U (1));
+  // set dataspace for the attribute
+  hsize_t dims_out[1];
+  dims_out[0] = n;
+  DataSpace dataspace(1,dims_out);
+  // create attribute
+  Attribute attr=header_group.createAttribute(attr_name,mem_type,dataspace);
+  // write attribute
+  attr.write(mem_type,attr_value);
+
+  return true;
+}
+// ============================================================================
+// return HDF5 native type corresponding to U parameter
+//
+template <class T>
+template <class U>
+DataType GH5<T>::guessType(U value)
+{
+  DataType ret;
+  if ((std::numeric_limits<U>::max()) == (std::numeric_limits<int>::max())) {
+    if (verbose) std::cerr << "U is an INT\n";
+    ret = PredType::NATIVE_INT;
+  } else
+  if ((std::numeric_limits<U>::max()) == (std::numeric_limits<float>::max())) {
+    if (verbose) std::cerr << "U is an FLOAT\n";
+    ret = PredType::NATIVE_FLOAT;
+  } else
+  if ((std::numeric_limits<U>::max()) == (std::numeric_limits<double>::max())) {
+    if (verbose) std::cerr << "U is an double\n";
+    ret = PredType::NATIVE_DOUBLE;
+  } else
+  if ((std::numeric_limits<U>::max()) == (std::numeric_limits<long int>::max())) {
+    if (verbose) std::cerr << "U is an LONG INT\n";
+    ret = PredType::NATIVE_LONG;
+  } else
+  if ((std::numeric_limits<U>::max()) == (std::numeric_limits<long long int>::max())) {
+    if (verbose) std::cerr << "U is an LONG LONG INT\n";
+    ret = PredType::NATIVE_LLONG;
+  } else
+    if ((std::numeric_limits<U>::max()) == (std::numeric_limits<long double>::max())) {
+      if (verbose) std::cerr << "U is an LONG LONG INT\n";
+      ret = PredType::NATIVE_LDOUBLE;
+    } else {
+      std::cerr << "GH5<T>::guessType, unknown type !!!!\n";
+      throw(-1);
+      std::exit(1);
+    }
+  return ret;
+}
+
 // !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 // Templates instantiation MUST be declared **AFTER** templates declaration
 // !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
@@ -729,6 +1239,12 @@ template class CSnapshotGadgetH5In<float>;
 //extern template class CSnapshotGadgetH5In<double>;
 template class CSnapshotGadgetH5In<double>;
 
+//extern template class CSnapshotGadgetH5Out<float>;
+template class CSnapshotGadgetH5Out<float>;
+
+//extern template class CSnapshotGadgetH5Out<double>;
+template class CSnapshotGadgetH5Out<double>;
+
 //extern template class GH5<float>;
 template class GH5<float>;