draft for SSV

Apps/Common/SIMSolver.h

@@ -19,6 +19,7 @@
 #include "TimeStep.h"
 #include "HDF5Restart.h"
 #include "HDF5Writer.h"
+#include "HDF5WriterSSV.h"
 #include "tinyxml.h"
 
 
@@ -47,17 +48,23 @@ template<class T1> class SIMSolverStat : public SIMadmin
   //! \param[in] hdf5file The file to save to
   //! \param[in] modelAdm Process administrator to use
   //! \param[in] saveInterval The stride in the output file
+  //! \param[in] ssv Enable SSV writer
   void handleDataOutput(const std::string& hdf5file,
                         const ProcessAdm& modelAdm,
-                        int saveInterval = 1)
+                        int saveInterval = 1,
+                        bool vis = true,
+                        bool ssv = false)
   {
     if (IFEM::getOptions().discretization == ASM::Spectral && !hdf5file.empty())
       IFEM::cout <<"\n  ** HDF5 output is available for spline/lagrangian discretization"
                  <<" only. Deactivating...\n"<< std::endl;
-    else
+    else if (vis || ssv)
     {
       exporter = new DataExporter(true,saveInterval);
-      exporter->registerWriter(new HDF5Writer(hdf5file,modelAdm));
+      if (vis)
+        exporter->registerWriter(new HDF5Writer(hdf5file,modelAdm));
+      if (ssv)
+        exporter->registerWriter(new HDF5WriterSSV(hdf5file,modelAdm));
       S1.registerFields(*exporter);
       IFEM::registerCallback(*exporter);
     }
@@ -169,15 +176,18 @@ template<class T1> class SIMSolver : public SIMSolverStat<T1>
   //! \param[in] modelAdm Process administrator to use
   //! \param[in] saveInterval The stride in the output file
   //! \param[in] restartInterval The stride in the restart file
+  //! \param[in] ssv Enable SSV writer
   void handleDataOutput(const std::string& hdf5file,
                         const ProcessAdm& modelAdm,
                         int saveInterval = 1,
-                        int restartInterval = 0)
+                        int restartInterval = 0,
+                        bool vis = true,
+                        bool ssv = false)
   {
     if (restartInterval > 0)
       restartAdm = new HDF5Restart(hdf5file+"_restart",modelAdm,restartInterval);
 
-    this->SIMSolverStat<T1>::handleDataOutput(hdf5file, modelAdm, saveInterval);
+    this->SIMSolverStat<T1>::handleDataOutput(hdf5file, modelAdm, saveInterval, vis, ssv);
   }
 
   //! \brief Serialize internal state for restarting purposes.

src/SIM/SIMoptions.C

@@ -52,6 +52,8 @@ SIMoptions::SIMoptions ()
   nViz[0] = nViz[1] = nViz[2] = 2;
 
   printPid = 0;
+  saveHDF5viz = true;
+  saveHDF5ssv = false;
 }
 
 
@@ -174,6 +176,9 @@ bool SIMoptions::parseOutputTag (const TiXmlElement* elem)
     }
     else // use the default output file name
       hdf5 = "(default)";
+
+    utl::getAttribute(elem,"ssv",saveHDF5ssv);
+    utl::getAttribute(elem,"viz",saveHDF5viz);
   }
 
   else if (!strcasecmp(elem->Value(),"primarySolOnly"))

src/SIM/SIMoptions.h

@@ -91,6 +91,8 @@ class SIMoptions
   bool pSolOnly; //!< If \e true, don't save secondary solution variables
   bool saveNorms;//!< If \e true, save element norms
   bool saveLog;  //!< If \e true, export the log
+  bool saveHDF5viz; //!< If \e true, save HDF5 visualization file
+  bool saveHDF5ssv; //!< If \e true, save HDF5 SSV file
 
   std::string hdf5; //!< Prefix for HDF5-file
   std::string vtf;  //!< Prefix for VTF-file

src/Utility/DataExporter.C

@@ -121,6 +121,13 @@ bool DataExporter::dumpTimeLevel (const TimeStep* tp, bool geometryUpdated)
 
   for (DataWriter* writer : m_writers) {
     writer->openFile(m_level);
+    for (const DataEntry& it : m_entry)
+      if (it.second.field == SIM)
+        if (!writer->prepare(m_level,it,geometryUpdated,it.second.prefix)) {
+          std::cerr << "** DataExporter::dumpTimeLevel: Error preparing output for field " << it.first << std::endl;
+          return false;
+        }
+
     for (const DataEntry& it : m_entry) {
       if (!it.second.data)
         return false;

src/Utility/DataExporter.h

@@ -189,6 +189,16 @@ class DataWriter
   //! \param[in] level Level we just wrote to the file
   virtual void closeFile(int level) = 0;
 
+  //! \brief Prepare writer for writing a time level.
+  //! \param[in] level The time level to write the data at
+  //! \param[in] entry The DataEntry describing the vector
+  //! \param[in] geometryUpdated Whether or not geometries should be written
+  //! \param[in] prefix Field name prefix
+  //! \details Used by SVV where we have to preallocate datasets
+  virtual bool prepare(int level, const DataEntry& entry,
+                       bool geometryUpdated, const std::string& prefix)
+  { return true; }
+
   //! \brief Writes a vector to file.
   //! \param[in] level The time level to write the vector at
   //! \param[in] entry The DataEntry describing the vector

src/Utility/HDF5Base.C

@@ -35,27 +35,40 @@ HDF5Base::HDF5Base (const std::string& name, const ProcessAdm& adm)
 }
 
 
-bool HDF5Base::openFile (unsigned int flags)
+bool HDF5Base::openFile (unsigned int flags, bool latest)
 {
 #ifdef HAS_HDF5
   if (m_file != -1)
     return true;
 
+  hid_t acc_tpl = H5P_DEFAULT;
+  bool acc_tpl_alloc = false;
+  if (latest
+#ifdef HAVE_MPI
+      && m_adm.getProcId() == 0 || !m_adm.dd.isPartitioned()
+ #endif
+      ) {
+    acc_tpl = H5Pcreate(H5P_FILE_ACCESS);
+    acc_tpl_alloc = true;
+  }
+
 #ifdef HAVE_MPI
-  hid_t acc_tpl;
   if (m_adm.dd.isPartitioned()) {
     if (m_adm.getProcId() != 0)
       return true;
-    acc_tpl = H5P_DEFAULT;
   } else {
     MPI_Info info = MPI_INFO_NULL;
-    acc_tpl = H5Pcreate(H5P_FILE_ACCESS);
+    if (!acc_tpl_alloc) {
+        acc_tpl = H5Pcreate(H5P_FILE_ACCESS);
+        acc_tpl_alloc = true;
+    }
     H5Pset_fapl_mpio(acc_tpl, *m_adm.getCommunicator(), info);
   }
-#else
-  hid_t acc_tpl = H5P_DEFAULT;
 #endif
 
+  if (latest)
+    H5Pset_libver_bounds(acc_tpl, H5F_LIBVER_LATEST, H5F_LIBVER_LATEST);
+
   if (flags == H5F_ACC_TRUNC)
     m_file = H5Fcreate(m_hdf5_name.c_str(),flags,H5P_DEFAULT,acc_tpl);
   else
@@ -64,13 +77,14 @@ bool HDF5Base::openFile (unsigned int flags)
   if (m_file < 0)
   {
     std::cerr <<" *** HDF5Base: Failed to open "<< m_hdf5_name << std::endl;
+    if (acc_tpl_alloc)
+      H5Pclose(acc_tpl);
     return false;
   }
 
-#ifdef HAVE_MPI
-  if (!m_adm.dd.isPartitioned())
+  if (acc_tpl_alloc)
     H5Pclose(acc_tpl);
-#endif
+
   return true;
 #else
   return false;

src/Utility/HDF5Base.h

@@ -39,7 +39,8 @@ class HDF5Base
 protected:
   //! \brief Opens the HDF5 file.
   //! \param[in] flag Mode to open file using
-  bool openFile(unsigned int flag);
+  //! \param[in] latest If true, open using latest version
+  bool openFile(unsigned int flag, bool latest = false);
   //! \brief Closes the HDF5 file.
   void closeFile();
 

src/Utility/HDF5Writer.C

@@ -168,12 +168,18 @@ void HDF5Writer::writeArray(hid_t group, const std::string& name, int patch,
     space = H5Screate_simple(1,&siz,nullptr);
     std::stringstream str;
     str << patch;
-    set = H5Dcreate2(group1,str.str().c_str(),
-                     type,space,H5P_DEFAULT,H5P_DEFAULT,H5P_DEFAULT);
+    if (checkGroupExistence(group1,str.str().c_str()))
+      set = H5Dopen2(group1,str.str().c_str(),H5P_DEFAULT);
+    else
+      set = H5Dcreate2(group1,str.str().c_str(),
+                       type,space,H5P_DEFAULT,H5P_DEFAULT,H5P_DEFAULT);
   } else {
     space = H5Screate_simple(1,&siz,nullptr);
-    set = H5Dcreate2(group,name.c_str(),
-                     type,space,H5P_DEFAULT,H5P_DEFAULT,H5P_DEFAULT);
+    if (checkGroupExistence(group,name.c_str()))
+      set = H5Dopen2(group,name.c_str(),H5P_DEFAULT);
+    else
+      set = H5Dcreate2(group,name.c_str(),
+                       type,space,H5P_DEFAULT,H5P_DEFAULT,H5P_DEFAULT);
   }
   if (len > 0) {
     hid_t file_space = H5Dget_space(set);
@@ -185,6 +191,7 @@ void HDF5Writer::writeArray(hid_t group, const std::string& name, int patch,
     H5Sclose(mem_space);
     H5Sclose(file_space);
   }
+  H5Dflush(set);
   H5Dclose(set);
   H5Sclose(space);
   if (group1 != -1)
@@ -241,6 +248,14 @@ void HDF5Writer::writeBasis (int level, const DataEntry& entry,
 
 void HDF5Writer::writeSIM (int level, const DataEntry& entry,
                            bool geometryUpdated, const std::string& prefix)
+{
+  this->writeSIMInt(level, entry, geometryUpdated, prefix, false);
+}
+
+
+void HDF5Writer::writeSIMInt (int level, const DataEntry& entry,
+                              bool geometryUpdated,
+                              const std::string& prefix, bool noData)
 {
   if (!entry.second.enabled || !entry.second.data || entry.second.data2.empty())
     return;
@@ -336,7 +351,7 @@ void HDF5Writer::writeSIM (int level, const DataEntry& entry,
         (!(results & DataExporter::REDUNDANT) || sim->getGlobalProcessID() == 0)) // we own the patch
     {
       ASMbase* pch = sim->getPatch(loc);
-      if (results & DataExporter::PRIMARY && !sol->empty()) {
+      if (results & DataExporter::PRIMARY && !sol->empty() && !noData) {
         Vector psol;
         size_t ndof1 = sim->extractPatchSolution(*sol,psol,pch,entry.second.ncmps,
                                                  usedescription ? 1 : 0);
@@ -357,7 +372,7 @@ void HDF5Writer::writeSIM (int level, const DataEntry& entry,
                      i+1, ndof1, data, H5T_NATIVE_DOUBLE);
       }
 
-      if (results & DataExporter::SECONDARY && !sol->empty()) {
+      if (results & DataExporter::SECONDARY && !sol->empty() && !noData) {
         Matrix field;
         SIM::SolutionMode mode = prob->getMode();
         const_cast<SIMbase*>(sim)->setMode(SIM::RECOVERY);
@@ -369,7 +384,7 @@ void HDF5Writer::writeSIM (int level, const DataEntry& entry,
                      i+1, field.cols(), field.getRow(j+1).ptr(), H5T_NATIVE_DOUBLE);
       }
 
-      if (proj)
+      if (proj && !noData)
         for (size_t p = 0; p < proj->size(); ++p) {
           if (proj->at(p).empty())
             continue;
@@ -394,7 +409,7 @@ void HDF5Writer::writeSIM (int level, const DataEntry& entry,
                        i+1, field.cols(), field.getRow(j+1).ptr(), H5T_NATIVE_DOUBLE);
         }
 
-      if (results & DataExporter::NORMS && eNorm) {
+      if (results & DataExporter::NORMS && eNorm && !noData) {
         Matrix patchEnorm;
         sim->extractPatchElmRes(*eNorm,patchEnorm,loc-1);
         for (size_t j = 1, l = 1; l < eNorm->rows(); j++)
@@ -405,7 +420,7 @@ void HDF5Writer::writeSIM (int level, const DataEntry& entry,
                          i+1, patchEnorm.cols(), patchEnorm.getRow(l++).ptr(), H5T_NATIVE_DOUBLE);
       }
 
-      if (results & DataExporter::EIGENMODES) {
+      if (results & DataExporter::EIGENMODES && !noData) {
         size_t iMode = 0;
         const std::vector<Mode>* modes = static_cast<const std::vector<Mode>*>(entry.second.data2.front());
         for (const Mode& mode : *modes)
@@ -440,7 +455,7 @@ void HDF5Writer::writeSIM (int level, const DataEntry& entry,
     else // must write empty dummy records for the other patches
     {
       double dummy=0.0;
-      if (results & DataExporter::PRIMARY) {
+      if (results & DataExporter::PRIMARY && !noData) {
         if (usedescription)
           writeArray(group.front(), entry.second.description,
                      i+1, 0, &dummy, H5T_NATIVE_DOUBLE);
@@ -453,18 +468,18 @@ void HDF5Writer::writeSIM (int level, const DataEntry& entry,
                      i+1, 0, &dummy, H5T_NATIVE_DOUBLE);
       }
 
-      if (results & DataExporter::SECONDARY)
+      if (results & DataExporter::SECONDARY && !noData)
         for (size_t j = 0; j < prob->getNoFields(2); j++)
           writeArray(group.front(), prefix+prob->getField2Name(j),
                      i+1, 0, &dummy,H5T_NATIVE_DOUBLE);
 
-      if (proj)
+      if (proj && !noData)
         for (size_t p = 0; p < proj->size(); p++)
           for (size_t j = 0; j < prob->getNoFields(2); j++)
             writeArray(group.front(), m_prefix[p]+" "+prob->getField2Name(j),
                        i+1, 0, &dummy,H5T_NATIVE_DOUBLE);
 
-      if (results & DataExporter::NORMS && eNorm)
+      if (results & DataExporter::NORMS && eNorm && !noData)
         for (size_t j = 1; j <= norm->getNoFields(0); j++)
           for (size_t k = 1; k <= norm->getNoFields(j); k++)
             if (norm->hasElementContributions(j,k))

src/Utility/HDF5Writer.h

@@ -99,6 +99,19 @@ class HDF5Writer : public DataWriter, public HDF5Base
 
 #ifdef HAS_HDF5
 protected:
+  //! \brief Writes data from a SIM to file.
+  //! \param[in] level The time level to write the data at
+  //! \param[in] entry The DataEntry describing the data to write
+  //! \param[in] geometryUpdated Whether or not geometries should be written
+  //! \param[in] prefix Field name prefix
+  //! \param[in] noData True to not write data (only open groups)
+  //!
+  //! \details If prefix is non-empty and we are asked to dump secondary
+  //! solutions, we assume they come from different projections
+  //! \sa SIMbase::project
+  void writeSIMInt(int level, const DataEntry& entry,
+                   bool geometryUpdated, const std::string& prefix, bool noData);
+
   //! \brief Internal helper function writing a data array to file.
   //! \param[in] group The HDF5 group to write data into
   //! \param[in] name The name of the array
@@ -120,7 +133,6 @@ class HDF5Writer : public DataWriter, public HDF5Base
                   int basis, int level, bool redundant = false,
                   bool l2g = false);
 
-private:
   unsigned int m_flag; //!< The file flags to open HDF5 file with
 #endif
 };