Make doxygen happy.

config/defaults.mk

@@ -24,7 +24,7 @@ EGREP_BIN = $(shell command -v egrep 2> /dev/null)
 CXX = g++
 MPICXX = mpicxx
 
-BASE_FLAGS  = -std=c++11
+BASE_FLAGS  = -std=c++17
 OPTIM_FLAGS = -O3 $(BASE_FLAGS)
 DEBUG_FLAGS = -g $(XCOMPILER)-Wall $(BASE_FLAGS)
 

examples/ex15p.cpp

@@ -222,6 +222,10 @@ int main(int argc, char *argv[])
    visit_dc.RegisterField("solution", &x);
    int vis_cycle = 0;
 
+   MFEMDataCollection mfem_dc("Example15-internal", pmesh, true);
+   mfem_dc.RegisterField("solution", &x);
+   mfem_dc.ResetMetadata();
+
    // 10. As in Example 6p, we set up an estimator that will be used to obtain
    //     element error indicators. The integrator needs to provide the method
    //     ComputeElementFlux. We supply an L2 space for the discontinuous flux
@@ -344,6 +348,10 @@ int main(int argc, char *argv[])
          }
          if (visit)
          {
+            mfem_dc.SetCycle(vis_cycle);
+            mfem_dc.SetTime(time);
+            mfem_dc.Save();
+
             visit_dc.SetCycle(vis_cycle++);
             visit_dc.SetTime(time);
             visit_dc.Save();

examples/ex9p.cpp

@@ -252,6 +252,7 @@ int main(int argc, char *argv[])
    double t_final = 10.0;
    double dt = 0.01;
    bool visualization = true;
+   bool internal = false;
    bool visit = false;
    bool paraview = false;
    bool adios2 = false;
@@ -300,6 +301,9 @@ int main(int argc, char *argv[])
    args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
                   "--no-visualization",
                   "Enable or disable GLVis visualization.");
+   args.AddOption(&internal, "-idf", "--mfem-datafiles", "-no-idf",
+                  "--no-mfem-datafiles",
+                  "Save data files in the internal formal.");
    args.AddOption(&visit, "-visit", "--visit-datafiles", "-no-visit",
                   "--no-visit-datafiles",
                   "Save data files for VisIt (visit.llnl.gov) visualization.");
@@ -505,6 +509,18 @@ int main(int argc, char *argv[])
       pd->Save();
    }
 
+   MFEMDataCollection *idc = NULL;
+   if (internal)
+   {
+      idc = new MFEMDataCollection("Example9P-internal", *pmesh, false);
+      idc->SetPrefixPath("internal");
+      idc->RegisterField("solution", u);
+      idc->ResetMetadata();
+      idc->SetCycle(0);
+      idc->SetTime(0.0);
+      idc->Save();
+   }
+
    // Optionally output a BP (binary pack) file using ADIOS2. This can be
    // visualized with the ParaView VTX reader.
 #ifdef MFEM_USE_ADIOS2
@@ -606,6 +622,13 @@ int main(int argc, char *argv[])
             pd->Save();
          }
 
+         if (internal)
+         {
+            idc->SetCycle(ti);
+            idc->SetTime(t);
+            idc->Save();
+         }
+
 #ifdef MFEM_USE_ADIOS2
          // transient solutions can be visualized with ParaView
          if (adios2)
@@ -639,6 +662,7 @@ int main(int argc, char *argv[])
    delete fes;
    delete pmesh;
    delete ode_solver;
+   delete idc;
    delete pd;
 #ifdef MFEM_USE_ADIOS2
    if (adios2)

fem/datacollection.cpp

@@ -1171,4 +1171,408 @@ const char *ParaViewDataCollection::GetDataTypeString() const
    }
 }
 
+std::istream& operator>>(std::istream& in,
+                         mfem::MFEMDataCollection::MetaInfo& i)
+{
+   std::string  line;
+   ///@TODO error management.
+   if (std::getline(in, line))
+   {
+      std::stringstream linestream(line);
+      linestream >> i.cycle;
+      linestream >> i.t;
+   }
+
+   return in;
+}
+
+//! Constructor to store stuff
+MFEMDataCollection::MFEMDataCollection(const std::string& collection_name,
+                                       Mesh& mesh_, bool adaptive_mesh_)
+   : DataCollection(collection_name, &mesh_), adaptive_mesh(adaptive_mesh_)
+{
+   own_data = false;
+   format = SERIAL_FORMAT;
+}
+
+#ifdef MFEM_USE_MPI
+//! Constructor to store stuff
+MFEMDataCollection::MFEMDataCollection(const std::string& collection_name,
+                                       ParMesh& mesh_, bool adaptive_mesh_)
+   : DataCollection(collection_name, &mesh_), adaptive_mesh(adaptive_mesh_)
+{
+   own_data = false;
+   format = PARALLEL_FORMAT;
+   appendRankToFileName = true;
+}
+#endif
+
+//! Constructor to load stuff
+MFEMDataCollection::MFEMDataCollection(const std::string& collection_name
+                                       , Format format_)
+   : DataCollection(collection_name)
+{
+   format = format_;
+   if (format == PARALLEL_FORMAT)
+   {
+      appendRankToFileName = true;
+      serial = false;
+   }
+   else
+   {
+      serial = true;
+   }
+   own_data = true;
+#ifdef MFEM_USE_MPI
+   m_comm = MPI_COMM_WORLD;
+   MPI_Comm_rank(m_comm, &myid);
+   MPI_Comm_size(m_comm, &num_procs);
+#endif
+}
+
+void MFEMDataCollection::ResetMetadata()
+{
+   // Clear file.
+   std::ofstream metadata_fs(GetMetafileName(), std::ios::trunc);
+   ///@TODO error management
+   metadata_fs
+         << (adaptive_mesh ? "amr_on" : "amr_off")
+         << std::endl;
+   for (auto [field_name, _] : GetFieldMap())
+   {
+      metadata_fs << field_name << " ";
+   }
+   metadata_fs << std::endl;
+   for (auto [field_name, _] : GetQFieldMap())
+   {
+      metadata_fs << field_name << " ";
+   }
+   metadata_fs << std::endl;
+
+   first_io = true;
+}
+
+
+void MFEMDataCollection::Save()
+{
+   //std::filesystem::path data_folder = CyclePathName(cycle);
+   std::string data_folder = CyclePathName(cycle);
+   // check if the directories are created
+   {
+      int err = create_directory(data_folder, mesh, myid);
+      if (err)
+      {
+         error = WRITE_ERROR;
+         MFEM_WARNING("Error creating directory: " << data_folder);
+         return; // do not even try to write the mesh
+      }
+   }
+
+   // Add metadata
+   if (myid == 0)
+   {
+      std::ofstream metadata_fs(GetMetafileName(),
+                                std::ios::out | std::ios::app);
+      metadata_fs << GetCycle() << " " << GetTime() << std::endl;
+   }
+
+   // If the mesh is not adaptive, then we just have to store the mesh once.
+   if (first_io || adaptive_mesh)
+   {
+      first_io = false;
+      //std::filesystem::path mesh_filename_raw;
+      std::string mesh_filename_raw;
+
+      if (!adaptive_mesh)
+      {
+         mesh_filename_raw = prefix_path + DataCollection::GetCollectionName();
+         mesh_filename_raw += "/mesh";
+      }
+      else
+      {
+         mesh_filename_raw =  data_folder + "/mesh";
+      }
+
+      auto mesh_filename = mesh_filename_raw;
+      if (appendRankToFileName)
+      {
+         mesh_filename += "." + to_padded_string(myid, pad_digits_rank);
+      }
+      mfem::ofgzstream mesh_fs(mesh_filename, compression);
+      mesh_fs.precision(precision);
+      if (format == PARALLEL_FORMAT)
+      {
+#ifdef MFEM_USE_MPI
+         ParMesh *pmesh = dynamic_cast<ParMesh*>(mesh);
+         if (pmesh && format == PARALLEL_FORMAT)
+         {
+            pmesh->ParPrint(mesh_fs);
+         }
+         else
+         {
+            MFEM_WARNING("Cannot save. Mesh is not a ParMesh." << data_folder);
+         }
+#else
+         MFEM_WARNING("Parallel output format not allowed in serial simulation." <<
+                      data_folder);
+         error = WRITE_ERROR;
+         return;
+#endif
+      }
+      else if (format == SERIAL_FORMAT)
+      {
+         mesh->Print(mesh_fs);
+      }
+      else
+      {
+         error = WRITE_ERROR;
+         MFEM_WARNING("Unknown output format: " << data_folder);
+         return;
+      }
+   }
+
+   for (FieldMapIterator it = field_map.begin(); it != field_map.end(); ++it)
+   {
+      auto field_filename = data_folder + "/" + it->first;
+      if (appendRankToFileName)
+      {
+         field_filename += "." + to_padded_string(myid, pad_digits_rank);
+      }
+      mfem::ofgzstream field_file(field_filename, compression);
+
+      field_file.precision(precision);
+      (it->second)->Save(field_file);
+      if (!field_file)
+      {
+         error = WRITE_ERROR;
+         MFEM_WARNING("Error writing field to file: " << it->first);
+      }
+      // Even if there is an error, try saving the other fields
+   }
+
+   for (QFieldMapIterator it = q_field_map.begin(); it != q_field_map.end();
+        ++it)
+   {
+      auto field_filename = data_folder + "/" + it->first;
+      if (appendRankToFileName)
+      {
+         field_filename += "." + to_padded_string(myid, pad_digits_rank);
+      }
+      mfem::ofgzstream field_file(field_filename, compression);
+
+      field_file.precision(precision);
+      (it->second)->Save(field_file);
+      if (!field_file)
+      {
+         error = WRITE_ERROR;
+         MFEM_WARNING("Error writing field to file: " << it->first);
+      }
+   }
+}
+
+
+void MFEMDataCollection::Load(int cycle_)
+{
+   this->cycle = cycle_;
+
+   std::string data_folder = CyclePathName(cycle);
+   // if(!std::filesystem::exists(data_folder))
+   // {
+   //    error = READ_ERROR;
+   //    MFEM_WARNING("Unable to find directory: " << data_folder);
+   //    return; // do not even try to read the data
+   // }
+
+   // Clear relevant data.
+   if (own_data && adaptive_mesh && !first_io)
+   {
+      delete this->mesh;
+      this->mesh = nullptr;
+   }
+   field_map.DeleteData(true);
+   q_field_map.DeleteData(true);
+
+   // Load mesh and take ownership
+   // std::filesystem::path mesh_filename;
+   std::string mesh_filename;
+   // If the mesh is not adaptive, then we just have to store the mesh once.
+   if ((first_io && own_data) || adaptive_mesh)
+   {
+      own_data = true;
+      if (!adaptive_mesh)
+      {
+         mesh_filename = prefix_path + DataCollection::GetCollectionName();
+         mesh_filename += "/mesh";
+      }
+      else
+      {
+         mesh_filename =  data_folder + "/mesh";
+      }
+
+      if (appendRankToFileName)
+      {
+         mesh_filename += "." + to_padded_string(myid, pad_digits_rank);
+      }
+
+      // if(!std::filesystem::exists(mesh_filename))
+      // {
+      //       error = READ_ERROR;
+      //       MFEM_WARNING("Unable to find mesh: " << mesh_filename);
+      //       return; // do not even try to read the data
+      // }
+
+      mfem::ifgzstream mesh_fs(mesh_filename);
+      if (adaptive_mesh || first_io)
+      {
+         if (format == PARALLEL_FORMAT)
+         {
+#ifdef MFEM_USE_MPI
+            mesh = new ParMesh(MPI_COMM_WORLD, mesh_fs);
+#else
+            error = READ_ERROR;
+            return;
+#endif
+         }
+         else if (format == SERIAL_FORMAT)
+         {
+            mesh = new Mesh(mesh_fs);
+         }
+         else
+         {
+            error = READ_ERROR;
+            return;
+         }
+      }
+   }
+   first_io = false;
+
+   // Load grids and take ownership
+   for (auto field_name : field_names_meta)
+   {
+      auto field_filename = data_folder + "/" + field_name;
+      if (appendRankToFileName)
+      {
+         field_filename += "." + to_padded_string(myid, pad_digits_rank);
+      }
+
+      if (format == PARALLEL_FORMAT)
+      {
+#ifdef MFEM_USE_MPI
+         mfem::ifgzstream ifs(field_filename);
+         this->RegisterField(field_name, new ParGridFunction(GetParMesh(), ifs));
+#else
+         error = READ_ERROR;
+         return;
+#endif
+      }
+      else if (format == SERIAL_FORMAT)
+      {
+         mfem::ifgzstream ifs(field_filename);
+         this->RegisterField(field_name, new GridFunction(mesh, ifs));
+      }
+   }
+
+   for (auto q_field_name : q_field_names_meta)
+   {
+      auto q_field_filename = data_folder + "/" + q_field_name;
+      if (appendRankToFileName)
+      {
+         q_field_filename += "." + to_padded_string(myid, pad_digits_rank);
+      }
+
+      mfem::ifgzstream ifs(q_field_filename);
+      this->RegisterQField(q_field_name, new QuadratureFunction(mesh, ifs));
+   }
+}
+
+
+std::optional<std::vector<MFEMDataCollection::MetaInfo>>
+                                                      MFEMDataCollection::ReloadMetaInfo()
+{
+   std::vector<MetaInfo> metavector;
+   std::ifstream in(GetMetafileName(), std::ios::in);
+
+   field_names_meta.clear();
+   q_field_names_meta.clear();
+
+   std::string line;
+   if (std::getline(in, line))
+   {
+      std::stringstream linestream(line);
+      std::string option_name;
+      ///@TODO refactor
+      while (linestream >> option_name)
+      {
+         if (option_name == "amr_on")
+         {
+            adaptive_mesh = true;
+         }
+         else if (option_name == "amr_off")
+         {
+            adaptive_mesh = false;
+         }
+         else
+         {
+            MFEM_WARNING("WARNING: Unknown option '"
+                         << option_name << "' in Metadata file "
+                         << GetMetafileName());
+         }
+      }
+   }
+   else
+   {
+      MFEM_WARNING("Corrputed Metadata. Option line missing.");
+      error = READ_ERROR;
+      return std::nullopt;
+   }
+   if (std::getline(in, line))
+   {
+      std::stringstream linestream(line);
+      std::string field_name;
+      while (linestream >> field_name)
+      {
+         field_names_meta.push_back(field_name);
+      }
+   }
+   else
+   {
+      MFEM_WARNING("Corrputed Metadata. Fields line missing.");
+      error = READ_ERROR;
+      return std::nullopt;
+   }
+   if (std::getline(in, line))
+   {
+      std::stringstream linestream(line);
+      std::string field_name;
+      while (linestream >> field_name)
+      {
+         q_field_names_meta.push_back(field_name);
+      }
+   }
+   else
+   {
+      MFEM_WARNING("Corrputed Metadata. QFields line missing.");
+      error = READ_ERROR;
+      return std::nullopt;
+   }
+
+   MetaInfo buffer;
+   while (in >> buffer)
+   {
+      metavector.push_back(buffer);
+   }
+   return metavector;
+}
+
+
+
+std::string MFEMDataCollection::CyclePathName(int cycle) const
+{
+   std::string out = "";
+   out = prefix_path + DataCollection::GetCollectionName();
+   out += "/Cycle" + to_padded_string(cycle,pad_digits_cycle);
+   return out;
+}
+
+
 }  // end namespace MFEM

fem/datacollection.hpp

@@ -20,6 +20,7 @@
 #include <string>
 #include <map>
 #include <fstream>
+#include <optional>
 
 namespace mfem
 {
@@ -556,5 +557,82 @@ public:
    virtual void Load(int cycle_ = 0) override;
 };
 
+
+
+/** Basically a helper to store a ".meta" file with triplets of the form (cycle,
+ * time, dt) along the fields for each cycle. */
+class MFEMDataCollection final : public DataCollection
+{
+public:
+   struct MetaInfo
+   {
+      int cycle;
+      double t;
+
+      friend std::istream& operator>>(std::istream& in,
+                                      mfem::MFEMDataCollection::MetaInfo& i);
+   };
+
+private:
+   std::string GetMetafileName()
+   {
+      return prefix_path + GetCollectionName() + ".meta";
+   }
+
+   // We require these two for loading.
+   std::vector<std::string> field_names_meta;
+   std::vector<std::string> q_field_names_meta;
+
+   bool adaptive_mesh = false;
+   bool first_io = true;
+
+public:
+   //! Constructor to store stuff
+   MFEMDataCollection(const std::string& collection_name,
+                      Mesh& mesh_, bool adaptive_mesh_);
+
+#ifdef MFEM_USE_MPI
+   //! Constructor to store stuff
+   MFEMDataCollection(const std::string& collection_name,
+                      ParMesh& mesh_, bool adaptive_mesh_);
+#endif
+
+#ifdef MFEM_USE_MPI
+   //! Constructor to load stuff
+   MFEMDataCollection(const std::string& collection_name
+                      , Format format_ = PARALLEL_FORMAT
+                     );
+#else
+   //! Constructor to load stuff
+   MFEMDataCollection(const std::string& collection_name
+                      , Format format_ = SERIAL_FORMAT
+                     );
+#endif
+   void ResetMetadata();
+
+#ifdef MFEM_USE_MPI
+   ParMesh* GetParMesh()
+   {
+      return dynamic_cast<ParMesh*>(GetMesh());
+   }
+#endif
+
+   void Save() override;
+
+   std::string CyclePathName(int cycle) const;
+
+   bool IsAMRData() const
+   {
+      return adaptive_mesh;
+   }
+
+   void Load(int cycle_) override;
+
+   std::optional<std::vector<MetaInfo>> ReloadMetaInfo();
+
+   /// Delete the mesh and fields if owned by the collection
+   ~MFEMDataCollection() = default;
+};
+
 }
 #endif

makefile

@@ -716,7 +716,7 @@ status info:
 	$(info MFEM_MPI_NP            = $(MFEM_MPI_NP))
 	@true
 
-ASTYLE_BIN = astyle
+ASTYLE_BIN ?= astyle
 ASTYLE = $(ASTYLE_BIN) --options=$(SRC)config/mfem.astylerc
 ASTYLE_VER = "Artistic Style Version 3.1"
 FORMAT_FILES = $(foreach dir,$(DIRS) $(EM_DIRS) config,$(dir)/*.?pp)

miniapps/tools/makefile

@@ -26,7 +26,7 @@ MFEM_LIB_FILE = mfem_is_not_built
 -include $(CONFIG_MK)
 
 SEQ_MINIAPPS = display-basis load-dc convert-dc get-values lor-transfer
-PAR_MINIAPPS =
+PAR_MINIAPPS = replay-simulation
 ifeq ($(MFEM_USE_MPI),NO)
    MINIAPPS = $(SEQ_MINIAPPS)
 else
@@ -76,7 +76,7 @@ RUN_MPI = $(MFEM_MPIEXEC) $(MFEM_MPIEXEC_NP) $(MFEM_MPI_NP)
 
 # Testing: Specific execution options
 # Do not test: display-basis, load-dc, convert-dc, get-values, lor-transfer
-NO_TEST_APPS = display-basis load-dc convert-dc get-values lor-transfer
+NO_TEST_APPS = display-basis load-dc convert-dc get-values lor-transfer replay-simulation
 $(foreach app,$(NO_TEST_APPS),$(app)-test-seq $(app)-test-par):
 	@true
 

miniapps/tools/replay-simulation.cpp

@@ -0,0 +1,71 @@
+#include "mfem.hpp"
+#include <memory>
+#include <iostream>
+#include <fstream>
+
+using namespace std;
+using namespace mfem;
+
+int main(int argc, char *argv[])
+{
+   MPI_Session mpi;
+   const int myid = mpi.WorldRank();
+
+   // 1. Parse command-line options
+   const char *simulation_name = "../../examples/internal/Example9P-internal";
+   const char *field_name = "solution";
+   double t0 = 0.0;
+   double T = std::numeric_limits<double>::infinity();
+   const char* vishost = "localhost";
+   int visport = 19916;
+   int precision = 8;
+
+   OptionsParser args(argc, argv);
+   args.AddOption(&simulation_name, "-sn", "--simulation-name",
+                  "Path with name of the MFEMDataCollection.");
+   args.AddOption(&field_name, "-fn", "--field-name",
+                  "Name of the field in the MFEMDataCollection.");
+   args.AddOption(&t0, "-t0", "--first-time-point",
+                  "Time point to begin the replay at.");
+   args.AddOption(&T, "-T", "--last-time-point",
+                  "Time point to end the replay at.");
+   args.Parse();
+   if (!args.Good())
+   {
+      if (myid == 0)
+      {
+         args.PrintUsage(cout);
+      }
+      return 1;
+   }
+   if (myid == 0)
+   {
+      args.PrintOptions(cout);
+   }
+
+   // Load data collection
+   auto dc = std::make_shared<MFEMDataCollection>(simulation_name);
+   auto metainfo = dc->ReloadMetaInfo();
+   if (!metainfo)
+   {
+      mfem_error("Failed to load meta info.");
+   }
+
+   socketstream vis;
+   vis.open(vishost, visport);
+   vis.precision(precision);
+
+   for (auto [cycle, t] : metainfo.value())
+   {
+      if (t < t0 || t > T) { continue; }
+
+      dc->Load(cycle);
+      auto pmesh = dc->GetParMesh();
+      auto u = dc->GetParField(field_name);
+      vis << "parallel " << pmesh->GetNRanks() << " " << pmesh->GetMyRank() << "\n";
+      vis << "solution\n" << *pmesh << *u << std::flush;
+      vis << "plot_caption 't=" << t << "'\n";
+   }
+
+   return 0;
+}