Eqp train (#48)

* SteadyNSSolver::InitROMHandler

* separated InitROMHandler from InitVariables

* cleaning up top-level routines. need some more thoughts on ROM elements..

* add random boundary on eqp test.

* LF flux does not work with constant rhs

* BasisTag struct

* hdf5_utils::Read/WriteAttributes for BasisTag

* Sample collection by port files

* interface test with periodic config.

* ROMInterfaceForm takes component basis and fes. subdomain basis and fes are view arrays.

* SteadyNSSolver::TrainEQPElems

* minor fix on face integraton rule.

* Read/WriteDataset for Array<SampleInfo>

* ChannelFlow minor fix

* ROMInterfaceForm::SaveEQPForIntegrator

* ROMInterfaceForm::LoadEQPForIntegrator

* Plug in rom interface form into SteadyNSEQPROM. Still needs to handle boundaries...

* SampleInfo consolidated el, be, face, itf into el.

* ROMNonlinearForm::GetEQPforIntegrator

* ROMNonlinearForm: decided to own bfnfi_markers, as opposed to the parent class.

* SteadyNSSolver::AssembleROMEQPOper. need verification..

* lid driven cavity problem

* randomizing rom nonlinear solve initial guess.

* SteadyNSSolver LF test without interface.

* interface eqp verified. interface to self is not working and needs a debugging.

* workflow with multiple port file verified.

* error flag for self interface

* adjusted threshold

include/hdf5_utils.hpp

@@ -8,6 +8,22 @@
 #include "mfem.hpp"
 #include "hdf5.h"
 #include "linalg_utils.hpp"
+#include "rom_handler.hpp"
+#include "hyperreduction_integ.hpp"
+
+namespace mfem
+{
+
+enum IntegratorType
+{
+   DOMAIN,
+   INTERIORFACE,
+   BDRFACE,
+   INTERFACE,
+   NUM_INTEG_TYPE
+};
+
+}
 
 using namespace mfem;
 
@@ -25,6 +41,9 @@ int GetDatasetSize(hid_t &source, std::string dataset, hsize_t* &dims);
 void ReadAttribute(hid_t &source, std::string attribute, std::string &value);
 void WriteAttribute(hid_t &source, std::string attribute, const std::string &value);
 
+void ReadAttribute(hid_t &source, std::string attribute, BasisTag &value);
+void WriteAttribute(hid_t &source, std::string attribute, const BasisTag &value);
+
 template <typename T>
 void ReadAttribute(hid_t &source, std::string attribute, T &value) {
    herr_t status;
@@ -68,8 +87,11 @@ void ReadDataset(hid_t &source, std::string dataset, Array<T> &value)
    assert(errf >= 0);
 
    value.SetSize(dims[0]);
-   hid_t dataType = hdf5_utils::GetType(value[0]);
-   H5Dread(dset_id, dataType, H5S_ALL, H5S_ALL, H5P_DEFAULT, value.Write());
+   if (dims[0] > 0)
+   {
+      hid_t dataType = hdf5_utils::GetType(value[0]);
+      H5Dread(dset_id, dataType, H5S_ALL, H5S_ALL, H5P_DEFAULT, value.Write());
+   }
 
    errf = H5Dclose(dset_id);
    assert(errf >= 0);
@@ -115,7 +137,8 @@ void WriteDataset(hid_t &source, std::string dataset, const Array<T> &value)
 {
    herr_t errf = 0;
 
-   hid_t dataType = hdf5_utils::GetType(value[0]);
+   T tmp;
+   hid_t dataType = hdf5_utils::GetType(tmp);
    hsize_t dims[1];
    dims[0] = value.Size();
 
@@ -125,8 +148,11 @@ void WriteDataset(hid_t &source, std::string dataset, const Array<T> &value)
    hid_t dset_id = H5Dcreate2(source, dataset.c_str(), dataType, dspace_id, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
    assert(dset_id >= 0);
 
-   errf = H5Dwrite(dset_id, dataType, H5S_ALL, H5S_ALL, H5P_DEFAULT, value.Read());
-   assert(errf >= 0);
+   if (dims[0] > 0)
+   {
+      errf = H5Dwrite(dset_id, dataType, H5S_ALL, H5S_ALL, H5P_DEFAULT, value.Read());
+      assert(errf >= 0);
+   }
 
    errf = H5Dclose(dset_id);
    assert(errf >= 0);
@@ -174,6 +200,9 @@ void WriteDataset(hid_t &source, std::string dataset, const DenseTensor &value);
 void ReadDataset(hid_t &source, std::string dataset, MatrixBlocks &value);
 void WriteDataset(hid_t &source, std::string dataset, const MatrixBlocks &value);
 
+void ReadDataset(hid_t &source, std::string dataset, const IntegratorType type, Array<SampleInfo> &value);
+void WriteDataset(hid_t &source, std::string dataset, const IntegratorType type, const Array<SampleInfo> &value);
+
 inline bool pathExists(hid_t id, const std::string& path)
 {
   return H5Lexists(id, path.c_str(), H5P_DEFAULT) > 0;

include/hyperreduction_integ.hpp

@@ -11,12 +11,13 @@ namespace mfem
 {
 
 struct SampleInfo {
-   // TODO(kevin): all integrators could use the same integer variable, reducing the memory size of this struct.
-   int el;         // element index (used for DomainIntegrator)
-   int face;       // face index (used for InteriorFaceIntegrator)
-   int be;         // boundary element index (used for BdrFaceIntegrator)
-   int itf;        // interface info index (used for InterfaceIntegrator)
-   // can add dofs for other hyper reductions.
+   /*
+      - For DomainIntegrator: element index
+      - For BdrFaceIntegrator: boundary element index
+      - For InteriorFaceIntegrator: face index
+      - For InterfaceIntegrator: interface info index
+   */
+   int el;
 
    int qp;         // quadrature point
    double qw;      // quadrature weight

include/main_workflow.hpp

@@ -15,18 +15,20 @@ void RunExample();
 
 MultiBlockSolver* InitSolver();
 SampleGenerator* InitSampleGenerator(MPI_Comm comm);
-std::vector<std::string> GetGlobalBasisTagList(const TopologyHandlerMode &topol_mode, bool separate_variable_basis);
+std::vector<BasisTag> GetGlobalBasisTagList(const TopologyHandlerMode &topol_mode, bool separate_variable_basis);
 
 void GenerateSamples(MPI_Comm comm);
 void CollectSamples(SampleGenerator *sample_generator);
+void CollectSamplesByPort(SampleGenerator *sample_generator, const std::string &basis_prefix);
+void CollectSamplesByBasis(SampleGenerator *sample_generator, const std::string &basis_prefix);
 void BuildROM(MPI_Comm comm);
 void TrainROM(MPI_Comm comm);
 // 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);
+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 = "");
 

include/multiblock_solver.hpp

@@ -228,36 +228,23 @@ friend class ParameterizedProblem;
    void LoadSolutionWithTime(const std::string &filename, int &step, double &time);
    void CopySolution(BlockVector *input_sol);
 
-   virtual void AllocateROMTensorElems()
-   { mfem_error("Abstract method MultiBlockSolver::AllocateROMTensorElems!\n"); }
+   virtual void AllocateROMNlinElems()
+   { mfem_error("Abstract method MultiBlockSolver::AllocateROMNlinElems!\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);
-
-   virtual BlockVector* PrepareSnapshots(std::vector<std::string> &basis_tags);
+   virtual void TrainROMEQPElems(SampleGenerator *sample_generator)
+   { mfem_error("Abstract method MultiBlockSolver::TrainROMEQPElems!\n"); }
+   virtual void SaveROMNlinElems(const std::string &input_prefix)
+   { mfem_error("Abstract method MultiBlockSolver::SaveROMNlinElems!\n"); }
+   virtual void LoadROMNlinElems(const std::string &input_prefix)
+   { mfem_error("Abstract method MultiBlockSolver::LoadROMNlinElems!\n"); }
+   virtual void AssembleROMNlinOper()
+   { mfem_error("Abstract method MultiBlockSolver::AssembleROMNlinOper!\n"); }
+
+   virtual void InitROMHandler();
+   void GetBasisTags(std::vector<BasisTag> &basis_tags);
+
+   virtual BlockVector* PrepareSnapshots(std::vector<BasisTag> &basis_tags);
    void SaveSnapshots(SampleGenerator *sample_generator);
 
    virtual void LoadReducedBasis() { rom_handler->LoadReducedBasis(); }

include/parameterized_problem.hpp

@@ -82,6 +82,13 @@ namespace backward_facing_step
    void pic(const Vector &x, double t, Vector &y);
 }
 
+namespace lid_driven_cavity
+{
+   extern double u0, L;
+
+   void ubdr(const Vector &x, double t, Vector &y);
+}
+
 namespace periodic_flow_past_array
 {
    extern Vector f;
@@ -339,12 +346,26 @@ class BackwardFacingStep : public FlowProblem
    virtual ~BackwardFacingStep() {};
 };
 
+class LidDrivenCavity : public FlowProblem
+{
+public:
+   LidDrivenCavity();
+   virtual ~LidDrivenCavity() {};
+};
+
 class PeriodicFlowPastArray : public FlowProblem
 {
 public:
    PeriodicFlowPastArray();
 };
 
+class ForceDrivenCorner : public PeriodicFlowPastArray
+{
+public:
+   ForceDrivenCorner();
+   virtual ~ForceDrivenCorner() {};
+};
+
 class LinElastProblem : public ParameterizedProblem
 {
 friend class LinElastSolver;

include/rom_element_collection.hpp

@@ -100,7 +100,7 @@ class ROMEQPElement : public ROMElementCollection
    Array<ROMNonlinearForm *> comp;     // Size(num_components);
    // boundary condition is enforced via forcing term.
    Array<Array<ROMNonlinearForm *> *> bdr;
-   Array<ROMInterfaceForm *> port;   // reference ports.
+   ROMInterfaceForm *port = NULL;   // reference ports.
 
 public:
    ROMEQPElement(TopologyHandler *topol_handler_,

include/rom_handler.hpp

@@ -30,8 +30,50 @@ enum NonlinearHandling
    NUM_NLNHNDL
 };
 
-const std::string GetBasisTagForComponent(const int &comp_idx, const TopologyHandler *topol_handler, const std::string var_name="");
-const std::string GetBasisTag(const int &subdomain_index, const TopologyHandler *topol_handler, const std::string var_name="");
+struct BasisTag
+{
+   /* component mesh name */
+   std::string comp = "";
+   /* variable name, if separate basis is used */
+   std::string var = "";
+
+   BasisTag() {}
+
+   BasisTag(const std::string &comp_, const std::string &var_="")
+      : comp(comp_), var(var_) {}
+
+   const std::string print() const
+   {
+      std::string tag = comp;
+      if (var != "")
+         tag += "_" + var;
+      return tag;
+   }
+
+   bool operator==(const BasisTag &tag) const
+   {
+      return ((comp == tag.comp) && (var == tag.var));
+   }
+
+   bool operator<(const BasisTag &tag) const
+   {
+      if (comp == tag.comp)
+         return (var < tag.var);
+
+      return (comp < tag.comp);
+   }
+
+   BasisTag& operator=(const BasisTag &tag)
+   {
+      comp = tag.comp;
+      var = tag.var;
+      return *this;
+   }
+};
+
+
+const BasisTag GetBasisTagForComponent(const int &comp_idx, const TopologyHandler *topol_handler, const std::string var_name="");
+const BasisTag GetBasisTag(const int &subdomain_index, const TopologyHandler *topol_handler, const std::string var_name="");
 
 class ROMHandlerBase
 {
@@ -81,7 +123,7 @@ class ROMHandlerBase
    Array<int> dim_ref_basis;
    Array<CAROM::Matrix*> carom_ref_basis;
    Array<int> rom_comp_block_offsets;
-   std::vector<std::string> basis_tags;
+   std::vector<BasisTag> basis_tags;
    bool basis_loaded;
    bool operator_loaded;
 
@@ -126,7 +168,7 @@ class ROMHandlerBase
    const bool SeparateVariable() { return separate_variable; }
    const std::string GetOperatorPrefix() { return operator_prefix; }
    const std::string GetBasisPrefix() { return basis_prefix; }
-   const std::string GetRefBasisTag(const int ref_idx) { return basis_tags[ref_idx]; }
+   const BasisTag GetRefBasisTag(const int ref_idx) { return basis_tags[ref_idx]; }
    const Array<int>* GetBlockOffsets() { return &rom_block_offsets; }
    const Array<int>* GetVarBlockOffsets() { return &rom_varblock_offsets; }
    virtual SparseMatrix* GetOperator() = 0;