Separate variable ROM (#9)

* GetTopologyHandlerMode().

* SetTrainMode function

* sketches/static_members.

* GetGlobalBasisTagList

* ROMHandler - name change from num_basis_sets to num_rom_comp_blocks.

* ROMHandler: name change from num_basis to comp_num_basis.

* ROMHandler: name change from carom_spatialbasis to carom_comp_basis.

* array for component basis and domain basis.

* ROMHandler: name change from vdim to fom_vdim.

* ROMHandler initialization change. fom_num_vdofs still needs to be changed.

* ROMHandler::SetBlockSizes. we may need num_vdofs instead of offsets.

* GetBasisTag/GetBasisTagForComponent, LookUpFromDict.

* ROMHandler::LoadReducedBasis, MFEMROMHandler::LoadReducedBasis

* ROMHandler::GetBasis -> GetReferenceBasis.

* ROMHandler -> ROMHandlerBase. Only MFEMROMHandler is used in practice.

* name changes for reference basis.

* ROMHandler::ProjectToRefBasis.

* ROMHandler::ProjectToRef/DomainBasis.

* ROMHandler - upper level projection routines.

* ROMHandler::ProjectGlobalToDomainBasis

* main_workflow: refactored the workflow.

* ROMHandler::SaveOperator.

* ROMHandler::LoadOperator -> SetRomMat.

* MultiBlockSolver::ProjectOperatorOnReducedBasis.

* ROMHandler: projection routines for vector.

* bug fix over global rom, separate variable workflow.

* Orthonormalize(mat1, mat).

* ROMHandler::ParseSupremizerInput. minor fix for MultiBlockSolver::GetComponentFESpace.

* StokesSolver::EnrichSupremizer. not tested yet.

* AuxiliaryTrainROM for supremizer enrichment and hyper-reduction.

* Using MultiBlockSolver::LoadReducedBasis.

* supremizer enrichment verified over global operator.

* separate variable rom does not switch block indices any more for the sake of simplicity.

* MatrixBlocks for Array2D<SparseMatrix *> container.

* hdf5 I/O for MatrixBlocks.

* linalg_utils - AddToBlockMatrix.

* hdf5_utils::WriteDataset for MatrixBlocks bug fix.

* component wise rom workflow with separate variable basis. needs verification and bug fix.

* MUMPSSolver seems to not behave well with SYMMETRIC_INDEFINITE.

* component-wise separate variable workflow verified for stokes flow.

* linalg_utils::AddToBlockMatrix - changed the block initialization.

* steady ns solver with separate variable verified.

* Bug fix in StokesFlow solver for heterogeneous component meshes. steady-ns fully verified for both fom/rom.

* minor adjustment of threshold.

include/hdf5_utils.hpp

@@ -7,6 +7,7 @@
 
 #include "mfem.hpp"
 #include "hdf5.h"
+#include "linalg_utils.hpp"
 
 using namespace mfem;
 
@@ -15,6 +16,7 @@ namespace hdf5_utils
 
 inline hid_t GetType(int) { return (H5T_NATIVE_INT); }
 inline hid_t GetType(double) { return (H5T_NATIVE_DOUBLE); }
+inline hid_t GetType(bool) { return (H5T_NATIVE_HBOOL); }
 
 hid_t GetNativeType(hid_t type);
 
@@ -169,6 +171,9 @@ void WriteDataset(hid_t &source, std::string dataset, const DenseMatrix &value);
 void ReadDataset(hid_t &source, std::string dataset, DenseTensor &value);
 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);
+
 inline bool pathExists(hid_t id, const std::string& path)
 {
   return H5Lexists(id, path.c_str(), H5P_DEFAULT) > 0;

include/input_parser.hpp

@@ -143,6 +143,46 @@ class InputParser
    void SetOption(const std::string &keys, const T &value)
    { SetOptionInDict<T>(keys, value, dict_); }
 
+   template<class T>
+   YAML::Node LookUpFromDict(const std::string &lkey, const T &lval, YAML::Node input_dict)
+   {
+      YAML::Node tmp;
+      if (input_dict.size() == 0)
+         return tmp["none"];  // this returns a zombie node that is false.
+
+      for (int k = 0; k < input_dict.size(); k++)
+      {
+         YAML::Node lval_node = FindNodeFromDict(lkey, input_dict[k]);
+         if (lval_node && (lval_node.as<T>() == lval))
+         {
+            return input_dict[k];
+            break;
+         }
+      }
+
+      return tmp["none"];  // this returns a zombie node that is false.
+   }
+
+   template<class LT, class RT>
+   const RT LookUpFromDict(const std::string &lkey, const LT &lval, const std::string &rkey, const RT &rval_default, YAML::Node input_dict)
+   {
+      RT rval = rval_default;
+      if (input_dict.size() == 0)
+         return rval;
+
+      for (int k = 0; k < input_dict.size(); k++)
+      {
+         YAML::Node lval_node = FindNodeFromDict(lkey, input_dict[k]);
+         if (lval_node && (lval_node.as<LT>() == lval))
+         {
+            rval = GetOptionFromDict<RT>(rkey, rval_default, input_dict[k]);
+            break;
+         }
+      }
+
+      return rval;
+   }
+
 private:
    void OverwriteOption(const std::string &forced_input);
 

include/linalg_utils.hpp

@@ -7,6 +7,7 @@
 
 #include "mfem.hpp"
 #include "librom.h"
+#include "etc.hpp"
 
 using namespace mfem;
 
@@ -45,16 +46,95 @@ void PrintVector(const CAROM::Vector &vec,
 namespace mfem
 {
 
+struct MatrixBlocks
+{
+   int nrows;
+   int ncols;
+   Array2D<SparseMatrix *> blocks;
+
+   MatrixBlocks() : nrows(0), ncols(0), blocks(0, 0) {}
+   MatrixBlocks(const int nrows_, const int ncols_)
+      : nrows(nrows_), ncols(ncols_), blocks(nrows_, ncols_)
+   { blocks = NULL; }
+   MatrixBlocks(const Array2D<SparseMatrix *> &mats)
+   { *this = mats; }
+
+   ~MatrixBlocks() { DeletePointers(blocks); }
+
+   void SetSize(const int w, const int h)
+   {
+      for (int i = 0; i < nrows; i++)
+         for (int j = 0; j < ncols; j++)
+            if (blocks(i, j)) delete blocks(i, j);
+
+      nrows = w; ncols = h;
+      blocks.SetSize(w, h);
+      blocks = NULL;
+   }
+
+   SparseMatrix*& operator()(const int i, const int j)
+   {
+      assert((i >= 0) && (i < nrows));
+      assert((j >= 0) && (j < ncols));
+      return blocks(i, j);
+   }
+
+   SparseMatrix const * operator()(const int i, const int j) const
+   {
+      assert((i >= 0) && (i < nrows));
+      assert((j >= 0) && (j < ncols));
+      return blocks(i, j);
+   }
+
+   // MatrixBlocks has the ownership. no copy assignment.
+   MatrixBlocks& operator=(const Array2D<SparseMatrix *> &mats)
+   {
+printf("operator= const Array2D<SparseMatrix *>.\n");
+      nrows = mats.NumRows();
+      ncols = mats.NumCols();
+      blocks = mats;
+
+      return *this;
+   }
+
+   // Copy assignment.
+   MatrixBlocks& operator=(const MatrixBlocks &matblock)
+   {
+printf("operator= const MatrixBlocks.\n");
+      if (this == &matblock) return *this;
+      DeletePointers(blocks);
+
+      nrows = matblock.nrows;
+      ncols = matblock.ncols;
+      blocks.SetSize(nrows, ncols);
+      blocks = NULL;
+
+      for (int i = 0; i < nrows; i++)
+         for (int j = 0; j < ncols; j++)
+            if (matblock(i, j)) blocks(i, j) = new SparseMatrix(*matblock(i, j));
+
+      return *this;
+   }
+};
+
+void AddToBlockMatrix(const Array<int> &ridx, const Array<int> &cidx, const MatrixBlocks &mats, BlockMatrix &bmat);
+
 void GramSchmidt(DenseMatrix& mat);
 
+/* Orthonormalize mat over mat1 and itself. */
+void Orthonormalize(DenseMatrix& mat1, DenseMatrix& mat);
+
 // Compute Rt * A * P
 void RtAP(DenseMatrix& R,
         const Operator& A,
         DenseMatrix& P,
         DenseMatrix& RAP);
-SparseMatrix* RtAP(DenseMatrix& R,
-                   const Operator& A,
-                   DenseMatrix& P);
+DenseMatrix* DenseRtAP(DenseMatrix& R,
+                       const Operator& A,
+                       DenseMatrix& P);
+SparseMatrix* SparseRtAP(DenseMatrix& R,
+                        const Operator& A,
+                        DenseMatrix& P);
 
 template<typename T>
 void PrintMatrix(const T &mat,

include/main_workflow.hpp

@@ -14,11 +14,14 @@ double dbc4(const Vector &);
 void RunExample();
 
 MultiBlockSolver* InitSolver();
-
 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 BuildROM(MPI_Comm comm);
 void TrainROM(MPI_Comm comm);
+// supremizer-enrichment, hypre-reduction optimization, etc..
+void AuxiliaryTrainROM(MPI_Comm comm);
 // return relative error if comparing solution.
 double SingleRun(MPI_Comm comm, const std::string output_file = "");
 

include/multiblock_solver.hpp

@@ -84,16 +84,17 @@ friend class ParameterizedProblem;
    Array<GridFunction *> global_us_visual;
 
    // rom variables.
-   ROMHandler *rom_handler = NULL;
+   ROMHandlerBase *rom_handler = NULL;
    TrainMode train_mode = NUM_TRAINMODE;
    bool use_rom = false;
+   bool separate_variable_basis = false;
 
    // Used for bottom-up building, only with ComponentTopologyHandler.
    // For now, assumes ROM basis represents the entire vector solution.
-   Array<SparseMatrix *> comp_mats;     // Size(num_components);
+   Array<MatrixBlocks *> comp_mats;     // Size(num_components);
    // boundary condition is enforced via forcing term.
-   Array<Array<SparseMatrix *> *> bdr_mats;
-   Array<Array2D<SparseMatrix *> *> port_mats;   // reference ports.
+   Array<Array<MatrixBlocks *> *> bdr_mats;
+   Array<MatrixBlocks *> port_mats;   // reference ports.
    // DenseTensor objects from nonlinear operators
    // will be defined per each derived MultiBlockSolver.
 
@@ -113,7 +114,7 @@ friend class ParameterizedProblem;
    GridFunction* GetGridFunction(const int k) { return us[k]; }
    const int GetDiscretizationOrder() const { return order; }
    const bool UseRom() const { return use_rom; }
-   ROMHandler* GetROMHandler() const { return rom_handler; }
+   ROMHandlerBase* GetROMHandler() const { return rom_handler; }
    const TrainMode GetTrainMode() { return train_mode; }
    const bool IsVisualizationSaved() const { return save_visual; }
    const std::string GetSolutionFilePrefix() const { return sol_prefix; }
@@ -166,6 +167,8 @@ friend class ParameterizedProblem;
    virtual void AssembleInterfaceMatrixes() = 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); }
 
@@ -210,7 +213,7 @@ friend class ParameterizedProblem;
    void GetBasisTags(std::vector<std::string> &basis_tags);
 
    virtual void PrepareSnapshots(BlockVector* &U_snapshots, std::vector<std::string> &basis_tags);
-   void LoadReducedBasis() { rom_handler->LoadReducedBasis(); }
+   virtual void LoadReducedBasis() { rom_handler->LoadReducedBasis(); }
    virtual void ProjectOperatorOnReducedBasis() = 0;
    virtual void ProjectRHSOnReducedBasis();
    virtual void SolveROM();

include/poisson_solver.hpp

@@ -54,6 +54,13 @@ friend class ParameterizedProblem;
 
    virtual ~PoissonSolver();
 
+   static const std::vector<std::string> GetVariableNames()
+   {
+      std::vector<std::string> varnames(1);
+      varnames[0] = "solution";
+      return varnames;
+   }
+
    virtual void SetupBCVariables() override;
    virtual void AddBCFunction(std::function<double(const Vector &)> F, const int battr = -1);
    virtual void AddBCFunction(const double &F, const int battr = -1);