Minor fix

bin/main.cpp

@@ -28,11 +28,13 @@ int main(int argc, char *argv[])
    { if (rank == 0) RunExample(); }
    else 
    {
-      if (mode == "sample_generation") GenerateSamples(MPI_COMM_WORLD);
-      else if (mode == "build_rom")    BuildROM(MPI_COMM_WORLD);
-      else if (mode == "train_rom")    TrainROM(MPI_COMM_WORLD);
-      else if (mode == "train_eqp")    TrainEQP(MPI_COMM_WORLD);
-      else if (mode == "single_run")   double dump = SingleRun(MPI_COMM_WORLD, output_file);
+      if (mode == "sample_generation")    GenerateSamples(MPI_COMM_WORLD);
+      else if (mode == "build_rom")       BuildROM(MPI_COMM_WORLD);
+      else if (mode == "train_rom")       TrainROM(MPI_COMM_WORLD);
+      else if (mode == "aux_train_rom")   AuxiliaryTrainROM(MPI_COMM_WORLD);
+      else if (mode == "train_eqp")       TrainEQP(MPI_COMM_WORLD);
+      else if (mode == "single_run")      SingleRun(MPI_COMM_WORLD, output_file);
+      else if (mode == "print_eqp")       PrintEQPCoords(MPI_COMM_WORLD);
       else
       {
          if (rank == 0) printf("Unknown mode %s!\n", mode.c_str());

include/etc.hpp

@@ -6,6 +6,7 @@
 #define ETC_HPP
 
 #include "mfem.hpp"
+#include "input_parser.hpp"
 
 using namespace mfem;
 using namespace std;
@@ -38,4 +39,113 @@ std::string string_format( const std::string& format, Args ... args )
     return std::string( buf.get(), buf.get() + size - 1 ); // We don't want the '\0' inside
 }
 
+class TimeProfiler
+{
+private:
+    std::string title = "";
+
+    Array<StopWatch *> timers;
+
+    Array<int> calls;
+    Array<bool> starts;
+    std::vector<std::string> names;
+    std::unordered_map<std::string, int> indices;
+
+    const MPI_Comm comm;
+    int rank;
+
+public:
+    TimeProfiler(const std::string &title_, MPI_Comm comm_=MPI_COMM_WORLD)
+        : title(title_), comm(comm_)
+    {
+        MPI_Comm_rank(MPI_COMM_WORLD, &rank);
+        calls.SetSize(0);
+        names.clear();
+        indices.clear();
+    }
+
+    virtual ~TimeProfiler()
+    {
+        std::string option = "time_profile/" + title;
+        if (config.GetOption<bool>(option, false))
+            Print(title);
+        DeletePointers(timers);
+    }
+
+    void SetTitle(const std::string &title_)
+    { title = title_; }
+
+    void Start(const std::string &name)
+    {
+        if (!indices.count(name))
+        {
+            indices[name] = timers.Size();
+            timers.Append(new StopWatch);
+            names.push_back(name);
+            calls.Append(0);
+            starts.Append(false);
+        }
+
+        assert(indices.count(name));
+
+        int idx = indices[name];
+        if (starts[idx])
+        {
+            printf("TimeProfiler: %s is already started!\n", name.c_str());
+            assert(!starts[idx]);
+        }
+
+        timers[idx]->Start();
+        starts[idx] = true;
+    }
+
+    void Stop(const std::string &name)
+    {
+        assert(indices.count(name));
+
+        int idx = indices[name];
+        if (!starts[idx])
+        {
+            printf("TimeProfiler: %s is not started!\n", name.c_str());
+            assert(starts[idx]);
+        }
+
+        timers[idx]->Stop();
+        calls[idx] += 1;
+        starts[idx] = false;
+    }
+
+    void Print(const std::string &title_, const bool compute_sum=false)
+    {
+        int nfunc = timers.Size();
+        for (int k = 0; k < nfunc; k++)
+            assert(!starts[k]);
+
+        Array<double> times(nfunc);
+        for (int k = 0; k < nfunc; k++)
+            times[k] = timers[k]->RealTime();
+        double total = times.Sum();
+
+        MPI_Reduce(MPI_IN_PLACE, times.GetData(), nfunc, MPI_DOUBLE, MPI_SUM, 0, comm);
+        MPI_Reduce(MPI_IN_PLACE, calls.GetData(), nfunc, MPI_INT, MPI_SUM, 0, comm);
+
+        if (rank == 0)
+        {
+            printf("%s", (title_ + "\n").c_str());
+
+            std::string line = std::string(100, '=');
+            line += "\n";
+            printf("%s", line.c_str());
+            printf("%20s\t%20s\t%20s\t%20s\n", "Function", "Total time", "Calls", "Time per call");
+            for (int k = 0; k < nfunc; k++)
+            {
+                printf("%20s\t%20.5e\t%20d\t%20.5e\n", names[k].c_str(), times[k], calls[k], times[k] / calls[k]);
+            }
+            if (compute_sum)
+                printf("%20s\t%20.5e\t%20d\t%20.5e\n", "Total time", total, 1, total);
+            printf("%s", line.c_str());
+        }
+    }
+};
+
 #endif

include/interface_form.hpp

@@ -16,6 +16,8 @@ namespace mfem
 class InterfaceForm
 {
 protected:
+   mutable TimeProfiler timer;
+
    int numSub = -1;
    int skip_zeros = 1;
 

include/interfaceinteg.hpp

@@ -7,6 +7,7 @@
 
 #include "mfem.hpp"
 #include "hyperreduction_integ.hpp"
+#include "etc.hpp"
 
 namespace mfem
 {
@@ -298,9 +299,13 @@ class DGLaxFriedrichsFluxIntegrator : public InterfaceNonlinearFormIntegrator
    DenseMatrix udof1, udof2, elv1, elv2;
    DenseMatrix elmat_comp11, elmat_comp12, elmat_comp21, elmat_comp22;
 
+   mutable TimeProfiler timer;
+
 public:
    DGLaxFriedrichsFluxIntegrator(Coefficient &q, VectorCoefficient *ud = NULL, const IntegrationRule *ir = NULL)
-      : InterfaceNonlinearFormIntegrator(ir), Q(&q), UD(ud) {}
+      : InterfaceNonlinearFormIntegrator(ir), Q(&q), UD(ud), timer("DGLaxFriedrichsFluxIntegrator") {}
+
+   virtual ~DGLaxFriedrichsFluxIntegrator() {};
 
    void AssembleFaceVector(const FiniteElement &el1,
                            const FiniteElement &el2,

include/main_workflow.hpp

@@ -24,12 +24,15 @@ void CollectSamplesByBasis(SampleGenerator *sample_generator, const std::string
 void BuildROM(MPI_Comm comm);
 void TrainROM(MPI_Comm comm);
 // supremizer-enrichment etc..
-void AuxiliaryTrainROM(MPI_Comm comm, SampleGenerator *sample_generator);
+void AuxiliaryTrainROM(MPI_Comm comm);
 // EQP training, could include hypre-reduction optimization.
 void TrainEQP(MPI_Comm comm);
 // Input parsing routine to list out all snapshot files for training a basis.
 void FindSnapshotFilesForBasis(const BasisTag &basis_tag, const std::string &default_filename, std::vector<std::string> &file_list);
 // return relative error if comparing solution.
 double SingleRun(MPI_Comm comm, const std::string output_file = "");
 
+// Auxiliary function to print out EQP point coordinates
+void PrintEQPCoords(MPI_Comm comm);
+
 #endif

include/multiblock_solver.hpp

@@ -260,6 +260,8 @@ friend class ParameterizedProblem;
    void ComputeRelativeError(Array<GridFunction *> fom_sols, Array<GridFunction *> rom_sols, Vector &error);
    void CompareSolution(BlockVector &test_U, Vector &error);
 
+   virtual void SaveEQPCoords(const std::string &filename) {}
+
 protected:
    virtual void AssembleROMMat(BlockMatrix &romMat);
 };

include/rom_interfaceform.hpp

@@ -10,6 +10,7 @@
 #include "hdf5_utils.hpp"
 #include "hyperreduction_integ.hpp"
 #include "linalg/NNLS.h"
+#include "etc.hpp"
 
 namespace mfem
 {

include/rom_nonlinearform.hpp

@@ -16,9 +16,11 @@ namespace mfem
 
 class ROMNonlinearForm : public NonlinearForm
 {
+// private:
+//    static const int Nt = 6;
+//    mutable StopWatch *jac_timers[Nt];
 private:
-   static const int Nt = 6;
-   mutable StopWatch *jac_timers[Nt];
+   mutable TimeProfiler timer;
 
 protected:
    /// ROM basis for projection.
@@ -231,6 +233,8 @@ class ROMNonlinearForm : public NonlinearForm
        The state @a x must be a true-dof vector. */
    virtual Operator &GetGradient(const Vector &x) const;
 
+   void SaveDomainEQPCoords(const int k, hid_t file_id, const std::string &dsetname);
+
 private:
    void PrecomputeDomainEQPSample(const IntegrationRule &ir, const DenseMatrix &basis, EQPSample &eqp_sample);
    void PrecomputeFaceEQPSample(const IntegrationRule &ir, const DenseMatrix &basis,

include/steady_ns_solver.hpp

@@ -8,6 +8,7 @@
 #include "stokes_solver.hpp"
 #include "rom_nonlinearform.hpp"
 #include "rom_interfaceform.hpp"
+#include "etc.hpp"
 
 // By convention we only use mfem namespace as default, not CAROM.
 using namespace mfem;
@@ -17,6 +18,9 @@ using namespace mfem;
 // Ultimately, we should implement InterfaceForm to pass, not the MultiBlockSolver itself.
 class SteadyNSOperator : public Operator
 {
+private:
+   mutable TimeProfiler timer;
+
 protected:
    bool direct_solve;
 
@@ -99,6 +103,9 @@ class SteadyNSTensorROM : public SteadyNSROM
 
 class SteadyNSEQPROM : public SteadyNSROM
 {
+private:
+   mutable TimeProfiler timer;
+
 protected:
    Array<ROMNonlinearForm *> hs; // not owned by SteadyNSEQPROM.
    ROMInterfaceForm *itf = NULL; // not owned by SteadyNSEQPROM.
@@ -189,6 +196,8 @@ friend class SteadyNSOperator;
    virtual void LoadROMNlinElems(const std::string &input_prefix) override;
    virtual void AssembleROMNlinOper() override;
 
+   void SaveEQPCoords(const std::string &filename) override;
+
 private:
    DenseTensor* GetReducedTensor(DenseMatrix *basis, FiniteElementSpace *fespace);
 

include/unsteady_ns_solver.hpp

@@ -16,6 +16,9 @@ class UnsteadyNSSolver : public SteadyNSSolver
 friend class ParameterizedProblem;
 friend class SteadyNSOperator;
 
+private:
+   mutable TimeProfiler timer;
+
 protected:
 
    // number of timesteps
@@ -67,6 +70,9 @@ friend class SteadyNSOperator;
    Array<int> rom_u_offsets;
    BlockMatrix *rom_mass = NULL;
 
+private:
+   double times[10];
+
 public:
    UnsteadyNSSolver();