Refactoring for EQP training (#36)

* moved NS advection operators to hyperreduction integrator.

* added header to sketch file

* main workflow: split CollectSamples for TrainEQP.

* separated TrainEQP from TrainROM.

* SteadyNSSolver: comp_tensors are allocated within AllocateROMElements.

* name change from ROMElement to ROMProjElems.

* name change from ProjElems to LinElems. tensor will be factored out.

* name change back from LinElems to ProjElems. tensor is handled together with linear elemets.

* name change from ProjElems to LinElems....

* component tensor is factored out from linear ROM.

* ROM operator filename is always specified at main workflow level.

* interface for component-wise EQP training.

bin/main.cpp

@@ -31,6 +31,7 @@ int main(int argc, char *argv[])
       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);
       else
       {

include/advdiff_solver.hpp

@@ -46,7 +46,7 @@ friend class ParameterizedProblem;
    void BuildDomainOperators() override;
 
    // Component-wise assembly
-   void BuildCompROMElement(Array<FiniteElementSpace *> &fes_comp) override;
+   void BuildCompROMLinElems(Array<FiniteElementSpace *> &fes_comp) override;
 
    bool Solve() override;
 

include/hyperreduction_integ.hpp

@@ -158,6 +158,71 @@ class VectorConvectionTrilinearFormIntegrator : virtual public HyperReductionInt
                                      const Vector &x, DenseMatrix &jac) override;
 };
 
+/*
+   < \nabla v, uu > domain integrator
+*/
+class IncompressibleInviscidFluxNLFIntegrator :
+   public HyperReductionIntegrator
+{
+private:
+   int dim;
+   Coefficient *Q{};
+   DenseMatrix dshape, dshapex, EF, uu, ELV, elmat_comp;
+   Vector shape;
+
+public:
+   IncompressibleInviscidFluxNLFIntegrator(Coefficient &q)
+      : HyperReductionIntegrator(true), Q(&q) { }
+
+   IncompressibleInviscidFluxNLFIntegrator() = default;
+
+   static const IntegrationRule &GetRule(const FiniteElement &fe,
+                                         ElementTransformation &T);
+
+   void AssembleElementVector(const FiniteElement &el,
+                              ElementTransformation &trans,
+                              const Vector &elfun,
+                              Vector &elvect) override;
+
+   void AssembleElementGrad(const FiniteElement &el,
+                           ElementTransformation &trans,
+                           const Vector &elfun,
+                           DenseMatrix &elmat) override;
+};
+
+/*
+   Interior face, interface integrator
+   < [v], {uu \dot n} + \Lambda * [u] >
+   \Lambda = max( | u- \dot n |, | u+ \dot n | )
+*/
+class DGLaxFriedrichsFluxIntegrator : public HyperReductionIntegrator
+{
+private:
+   int dim, ndofs1, ndofs2, nvdofs;
+   double w, un1, un2, un;
+   Coefficient *Q{};
+
+   Vector nor, flux, shape1, shape2, u1, u2;
+   DenseMatrix udof1, udof2, elv1, elv2;
+   DenseMatrix elmat_comp11, elmat_comp12, elmat_comp21, elmat_comp22;
+public:
+   DGLaxFriedrichsFluxIntegrator(Coefficient &q)
+      : HyperReductionIntegrator(true), Q(&q) {}
+
+   void AssembleFaceVector(const FiniteElement &el1,
+                           const FiniteElement &el2,
+                           FaceElementTransformations &Tr,
+                           const Vector &elfun, Vector &elvect) override;
+
+   /// @brief Assemble the local action of the gradient of the
+   /// NonlinearFormIntegrator resulting from a face integral term.
+    void AssembleFaceGrad(const FiniteElement &el1,
+                           const FiniteElement &el2,
+                           FaceElementTransformations &Tr,
+                           const Vector &elfun, DenseMatrix &elmat) override;
+
+};
+
 } // namespace mfem
 
 #endif

include/interfaceinteg.hpp

@@ -11,7 +11,7 @@
 namespace mfem
 {
 
-   class InterfaceNonlinearFormIntegrator : virtual public HyperReductionIntegrator
+class InterfaceNonlinearFormIntegrator : virtual public HyperReductionIntegrator
 {
 protected:
   InterfaceNonlinearFormIntegrator(const bool precomputable_ = false, const IntegrationRule *ir = NULL)

include/linelast_solver.hpp

@@ -89,9 +89,9 @@ class LinElastSolver : public MultiBlockSolver
    virtual void SetupDomainBCOperators();
 
    // Component-wise assembly
-   virtual void BuildCompROMElement(Array<FiniteElementSpace *> &fes_comp);
-   virtual void BuildBdrROMElement(Array<FiniteElementSpace *> &fes_comp);
-   virtual void BuildInterfaceROMElement(Array<FiniteElementSpace *> &fes_comp);
+   virtual void BuildCompROMLinElems(Array<FiniteElementSpace *> &fes_comp);
+   virtual void BuildBdrROMLinElems(Array<FiniteElementSpace *> &fes_comp);
+   virtual void BuildItfaceROMLinElems(Array<FiniteElementSpace *> &fes_comp);
 
    virtual void ProjectOperatorOnReducedBasis();
 

include/main_workflow.hpp

@@ -18,10 +18,13 @@ SampleGenerator* InitSampleGenerator(MPI_Comm comm);
 std::vector<std::string> GetGlobalBasisTagList(const TopologyHandlerMode &topol_mode, const TrainMode &train_mode, bool separate_variable_basis);
 
 void GenerateSamples(MPI_Comm comm);
+void CollectSamples(SampleGenerator *sample_generator);
 void BuildROM(MPI_Comm comm);
 void TrainROM(MPI_Comm comm);
-// supremizer-enrichment, hypre-reduction optimization, etc..
+// supremizer-enrichment etc..
 void AuxiliaryTrainROM(MPI_Comm comm, SampleGenerator *sample_generator);
+// 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 std::string &basis_tag, const std::string &default_filename, std::vector<std::string> &file_list);
 // return relative error if comparing solution.

include/multiblock_solver.hpp

@@ -183,36 +183,36 @@ friend class ParameterizedProblem;
    virtual void AssembleInterfaceMatrices() = 0;
 
    // Global ROM operator Loading.
-   virtual void SaveROMOperator(const std::string input_prefix="")
-   { rom_handler->SaveOperator(input_prefix); }
-   virtual void LoadROMOperatorFromFile(const std::string input_prefix="")
-   { rom_handler->LoadOperatorFromFile(input_prefix); }
+   virtual void SaveROMOperator(const std::string filename)
+   { rom_handler->SaveOperator(filename); }
+   virtual void LoadROMOperatorFromFile(const std::string filename)
+   { rom_handler->LoadOperatorFromFile(filename); }
 
    // Component-wise assembly
    void GetComponentFESpaces(Array<FiniteElementSpace *> &comp_fes);
-   void AllocateROMElements();
+   virtual void AllocateROMLinElems();
 
-   void BuildROMElements();
-   virtual void BuildCompROMElement(Array<FiniteElementSpace *> &fes_comp) = 0;
-   virtual void BuildBdrROMElement(Array<FiniteElementSpace *> &fes_comp) = 0;
+   void BuildROMLinElems();
+   virtual void BuildCompROMLinElems(Array<FiniteElementSpace *> &fes_comp) = 0;
+   virtual void BuildBdrROMLinElems(Array<FiniteElementSpace *> &fes_comp) = 0;
    // TODO(kevin): part of this can be transferred to InterfaceForm.
-   virtual void BuildInterfaceROMElement(Array<FiniteElementSpace *> &fes_comp) = 0;
+   virtual void BuildItfaceROMLinElems(Array<FiniteElementSpace *> &fes_comp) = 0;
 
-   void SaveROMElements(const std::string &filename);
+   void SaveROMLinElems(const std::string &filename);
    // Save ROM Elements in a hdf5-format file specified with file_id.
    // TODO: add more arguments to support more general data structures of ROM Elements.
-   virtual void SaveCompBdrROMElement(hid_t &file_id);
-   void SaveBdrROMElement(hid_t &comp_grp_id, const int &comp_idx);
-   void SaveInterfaceROMElement(hid_t &file_id);
+   virtual void SaveCompBdrROMLinElems(hid_t &file_id);
+   void SaveBdrROMLinElems(hid_t &comp_grp_id, const int &comp_idx);
+   void SaveItfaceROMLinElems(hid_t &file_id);
 
-   void LoadROMElements(const std::string &filename);
+   void LoadROMLinElems(const std::string &filename);
    // Load ROM Elements in a hdf5-format file specified with file_id.
    // TODO: add more arguments to support more general data structures of ROM Elements.
-   virtual void LoadCompBdrROMElement(hid_t &file_id);
-   void LoadBdrROMElement(hid_t &comp_grp_id, const int &comp_idx);
-   void LoadInterfaceROMElement(hid_t &file_id);
+   virtual void LoadCompBdrROMLinElems(hid_t &file_id);
+   void LoadBdrROMLinElems(hid_t &comp_grp_id, const int &comp_idx);
+   void LoadItfaceROMLinElems(hid_t &file_id);
 
-   void AssembleROM();
+   void AssembleROMMat();
 
    virtual bool Solve() = 0;
 
@@ -225,12 +225,31 @@ friend class ParameterizedProblem;
    void LoadSolution(const std::string &filename);
    void CopySolution(BlockVector *input_sol);
 
-   virtual void TrainEQP(SampleGenerator *sample_generator)
-   { mfem_error("Abstract method MultiBlockSolver::TrainEQP!\n"); }
-   virtual void SaveEQP()
-   { mfem_error("Abstract method MultiBlockSolver::TrainEQP!\n"); }
-   virtual void LoadEQP()
-   { mfem_error("Abstract method MultiBlockSolver::TrainEQP!\n"); }
+   virtual void AllocateROMTensorElems()
+   { mfem_error("Abstract method MultiBlockSolver::AllocateROMTensorElems!\n"); }
+   virtual void BuildROMTensorElems()
+   { mfem_error("Abstract method MultiBlockSolver::BuildROMTensorElems!\n"); }
+   virtual void SaveROMTensorElems(const std::string &filename)
+   { mfem_error("Abstract method MultiBlockSolver::SaveROMTensorElems!\n"); }
+   virtual void LoadROMTensorElems(const std::string &filename)
+   { mfem_error("Abstract method MultiBlockSolver::LoadROMTensorElems!\n"); }
+   virtual void AssembleROMTensorOper()
+   { mfem_error("Abstract method MultiBlockSolver::AssembleROMTensorOper!\n"); }
+
+   virtual void AllocateROMEQPElems()
+   { mfem_error("Abstract method MultiBlockSolver::AllocateROMEQPElems!\n"); }
+   virtual void TrainEQPElems(SampleGenerator *sample_generator)
+   { mfem_error("Abstract method MultiBlockSolver::TrainEQPElems!\n"); }
+   virtual void SaveEQPElems(const std::string &filename)
+   { mfem_error("Abstract method MultiBlockSolver::SaveEQP!\n"); }
+   virtual void LoadEQPElems(const std::string &filename)
+   { mfem_error("Abstract method MultiBlockSolver::LoadEQP!\n"); }
+   virtual void AssembleROMEQPOper()
+   { mfem_error("Abstract method MultiBlockSolver::AssembleROMEQPOper!\n"); }
+
+   void AllocateROMNlinElems();
+   void LoadROMNlinElems(const std::string &input_prefix);
+   void AssembleROMNlinOper();
 
    void InitROMHandler();
    void GetBasisTags(std::vector<std::string> &basis_tags);

include/nonlinear_integ.hpp

@@ -10,31 +10,6 @@
 namespace mfem
 {
 
-class IncompressibleInviscidFluxNLFIntegrator :
-   public VectorConvectionNLFIntegrator
-{
-private:
-   int dim;
-   Coefficient *Q{};
-   DenseMatrix dshape, dshapex, EF, uu, ELV, elmat_comp;
-   Vector shape;
-
-public:
-   IncompressibleInviscidFluxNLFIntegrator(Coefficient &q): Q(&q) { }
-
-   IncompressibleInviscidFluxNLFIntegrator() = default;
-
-   virtual void AssembleElementVector(const FiniteElement &el,
-                                      ElementTransformation &trans,
-                                      const Vector &elfun,
-                                      Vector &elvect);
-
-   virtual void AssembleElementGrad(const FiniteElement &el,
-                                    ElementTransformation &trans,
-                                    const Vector &elfun,
-                                    DenseMatrix &elmat);
-};
-
 /*
    TrilinearForm domain integrator that computes Temam's domain integral:
    For the trial velocity u and the test velocity v,
@@ -70,33 +45,6 @@ class TemamTrilinearFormIntegrator :
                                     DenseMatrix &elmat);
 };
 
-class DGLaxFriedrichsFluxIntegrator : public NonlinearFormIntegrator
-{
-private:
-   int dim, ndofs1, ndofs2, nvdofs;
-   double w, un1, un2;
-   Coefficient *Q{};
-
-   Vector nor, nh, shape1, shape2, u1, u2, un;
-   DenseMatrix udof1, udof2, elv1, elv2, uu;
-   DenseMatrix elmat_comp11, elmat_comp12, elmat_comp21, elmat_comp22;
-public:
-   DGLaxFriedrichsFluxIntegrator(Coefficient &q) : Q(&q) {}
-
-   virtual void AssembleFaceVector(const FiniteElement &el1,
-                                   const FiniteElement &el2,
-                                   FaceElementTransformations &Tr,
-                                   const Vector &elfun, Vector &elvect);
-
-   /// @brief Assemble the local action of the gradient of the
-   /// NonlinearFormIntegrator resulting from a face integral term.
-   virtual void AssembleFaceGrad(const FiniteElement &el1,
-                                 const FiniteElement &el2,
-                                 FaceElementTransformations &Tr,
-                                 const Vector &elfun, DenseMatrix &elmat);
-
-};
-
 /*
    TrilinearForm face integrator that computes Temam's flux:
    For the trial velocity u^1, u^2 and the test velocity v^1, v^2 on the interior face with the normal vector n^1,

include/poisson_solver.hpp

@@ -88,9 +88,9 @@ friend class ParameterizedProblem;
    virtual void AssembleInterfaceMatrices();
 
    // Component-wise assembly
-   virtual void BuildCompROMElement(Array<FiniteElementSpace *> &fes_comp);
-   virtual void BuildBdrROMElement(Array<FiniteElementSpace *> &fes_comp);
-   virtual void BuildInterfaceROMElement(Array<FiniteElementSpace *> &fes_comp);
+   virtual void BuildCompROMLinElems(Array<FiniteElementSpace *> &fes_comp);
+   virtual void BuildBdrROMLinElems(Array<FiniteElementSpace *> &fes_comp);
+   virtual void BuildItfaceROMLinElems(Array<FiniteElementSpace *> &fes_comp);
 
    virtual bool Solve();
 

include/rom_handler.hpp

@@ -186,8 +186,8 @@ class ROMHandlerBase
    virtual void Solve(BlockVector* U) = 0;
    virtual void NonlinearSolve(Operator &oper, BlockVector* U, Solver *prec=NULL) = 0;   
 
-   virtual void SaveOperator(const std::string input_prefix="") = 0;
-   virtual void LoadOperatorFromFile(const std::string input_prefix="") = 0;
+   virtual void SaveOperator(const std::string filename) = 0;
+   virtual void LoadOperatorFromFile(const std::string filename) = 0;
    virtual void SetRomMat(BlockMatrix *input_mat) = 0;
    virtual void SaveRomSystem(const std::string &input_prefix, const std::string type="mm") = 0;
 

include/sample_generator.hpp

@@ -96,11 +96,18 @@ class SampleGenerator
 
    void ReportStatus(const int &sample_idx);
 
-   // Perform SVD over snapshot for basis_tag. Calculate the coverage for ref_num_basis (optional).
-   void FormReducedBasis(const std::string &basis_prefix,
+   /*
+      Collect snapshot matrices from the file list to the specified basis tag.
+   */
+   void CollectSnapshots(const std::string &basis_prefix,
                          const std::string &basis_tag,
-                         const std::vector<std::string> &file_list,
-                         const int &ref_num_basis = -1);
+                         const std::vector<std::string> &file_list);
+   /*
+      Perform SVD over snapshot for basis_tag.
+      Calculate the energy fraction for num_basis.
+      CollectSnapshots must be executed before this.
+   */
+   void FormReducedBasis(const std::string &basis_prefix);
 
 private:
    const int GetDimFromSnapshots(const std::string &filename);

include/steady_ns_solver.hpp

@@ -162,26 +162,26 @@ friend class SteadyNSOperator;
    void SaveROMOperator(const std::string input_prefix="") override;
    void LoadROMOperatorFromFile(const std::string input_prefix="") override;
 
-   // Component-wise assembly
-   void BuildCompROMElement(Array<FiniteElementSpace *> &fes_comp) override;
-   // virtual void BuildBdrROMElement(Array<FiniteElementSpace *> &fes_comp);
-   // virtual void BuildInterfaceROMElement(Array<FiniteElementSpace *> &fes_comp);
-   void SaveCompBdrROMElement(hid_t &file_id) override;
-   void LoadCompBdrROMElement(hid_t &file_id) override;
-
    bool Solve() override;
 
    void ProjectOperatorOnReducedBasis() override;
 
    void SolveROM() override;
 
-   void TrainEQP(SampleGenerator *sample_generator) override;
-   void SaveEQP() override;
-   void LoadEQP() override;
+   void AllocateROMTensorElems() override;
+   void BuildROMTensorElems() override;
+   void SaveROMTensorElems(const std::string &filename) override;
+   void LoadROMTensorElems(const std::string &filename) override;
+   void AssembleROMTensorOper() override;
+
+   void AllocateROMEQPElems() override;
+   void TrainEQPElems(SampleGenerator *sample_generator) override;
+   void SaveEQPElems(const std::string &filename) override;
+   void LoadEQPElems(const std::string &filename) override;
+   void AssembleROMEQPOper() override;
 
 private:
    DenseTensor* GetReducedTensor(DenseMatrix *basis, FiniteElementSpace *fespace);
-   void SetupEQPOperators();
 };
 
 #endif

include/stokes_solver.hpp

@@ -183,9 +183,9 @@ friend class ParameterizedProblem;
    virtual void SetupPressureMassMatrix();
 
    // Component-wise assembly
-   virtual void BuildCompROMElement(Array<FiniteElementSpace *> &fes_comp);
-   virtual void BuildBdrROMElement(Array<FiniteElementSpace *> &fes_comp);
-   virtual void BuildInterfaceROMElement(Array<FiniteElementSpace *> &fes_comp);
+   virtual void BuildCompROMLinElems(Array<FiniteElementSpace *> &fes_comp);
+   virtual void BuildBdrROMLinElems(Array<FiniteElementSpace *> &fes_comp);
+   virtual void BuildItfaceROMLinElems(Array<FiniteElementSpace *> &fes_comp);
 
    virtual bool Solve();
    virtual void Solve_obsolete();

sketches/ns_mms.cpp

@@ -8,6 +8,7 @@
 #include <algorithm>
 #include "linalg_utils.hpp"
 #include "nonlinear_integ.hpp"
+#include "hyperreduction_integ.hpp"
 #include "dg_mixed_bilin.hpp"
 #include "dg_bilinear.hpp"
 #include "dg_linear.hpp"