EQP ROM workflow for unsteady NS

include/multiblock_solver.hpp

@@ -259,6 +259,9 @@ friend class ParameterizedProblem;
    void ComputeSubdomainErrorAndNorm(GridFunction *fom_sol, GridFunction *rom_sol, double &error, double &norm);
    void ComputeRelativeError(Array<GridFunction *> fom_sols, Array<GridFunction *> rom_sols, Vector &error);
    void CompareSolution(BlockVector &test_U, Vector &error);
+
+protected:
+   virtual void AssembleROMMat(BlockMatrix &romMat);
 };
 
 #endif

include/rom_element_collection.hpp

@@ -50,6 +50,8 @@ class ROMLinearElement : public ROMElementCollection
    Array<Array<MatrixBlocks *> *> bdr;
    Array<MatrixBlocks *> port;   // reference ports.
 
+   Array<MatrixBlocks *> mass;     // Size(num_components);
+
 public:
    ROMLinearElement(TopologyHandler *topol_handler_,
                     const Array<FiniteElementSpace *> &fes_,

include/rom_handler.hpp

@@ -31,6 +31,12 @@ enum NonlinearHandling
    NUM_NLNHNDL
 };
 
+enum class ROMOrderBy
+{
+   VARIABLE,
+   DOMAIN
+};
+
 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="");
 
@@ -90,16 +96,11 @@ class ROMHandlerBase
    /*
       offset for the global domain ROM blocks.
       For i-th subdomain and j-th variable,
-         index = i * num_var + j
+         ROMOrderBy::DOMAIN:   index = i * num_var + j
+         ROMOrderBy::VARIABLE:     index = j * numSub + i
    */
    Array<int> num_basis;
    Array<int> rom_block_offsets;
-   /*
-      the global domain ROM blocks with index order switched.
-      For i-th subdomain and j-th variable,
-         index = j * num_var + i
-   */
-   Array<int> rom_varblock_offsets;
 
    CAROM::Options* rom_options;
    CAROM::BasisGenerator *basis_generator;
@@ -109,6 +110,11 @@ class ROMHandlerBase
    bool update_right_SV = false;
    bool incremental = false;
 
+   ROMOrderBy ordering = ROMOrderBy::DOMAIN;
+
+   mfem::BlockVector *reduced_rhs = NULL;
+   mfem::BlockVector *reduced_sol = NULL;
+
    void ParseInputs();
 public:
    ROMHandlerBase(TopologyHandler *input_topol, const Array<int> &input_var_offsets,
@@ -129,7 +135,6 @@ class ROMHandlerBase
    const std::string GetBasisPrefix() { return basis_prefix; }
    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;
    const bool GetNonlinearMode() { return nonlinear_mode; }
    void SetNonlinearMode(const bool nl_mode)
@@ -139,6 +144,12 @@ class ROMHandlerBase
       nonlinear_mode = nl_mode;
    }
    const NonlinearHandling GetNonlinearHandling() { return nlin_handle; }
+   const ROMOrderBy GetOrdering() { return ordering; }
+   const int GetBlockIndex(const int m, const int v=-1);
+   void GetDomainAndVariableIndex(const int &rom_block_index, int &m, int &v);
+
+   mfem::BlockVector* GetReducedSolution() { return reduced_sol; }
+   mfem::BlockVector* GetReducedRHS() { return reduced_rhs; }
 
    /* parse inputs for supremizer. only for Stokes/SteadyNS Solver. */
    void ParseSupremizerInput(Array<int> &num_ref_supreme, Array<int> &num_supreme);
@@ -157,9 +168,9 @@ class ROMHandlerBase
    virtual void ProjectToDomainBasis(const Array<int> &idx_i, const Array<int> &idx_j, const Array2D<Operator*> &mats, Array2D<SparseMatrix *> &rom_mats) = 0;
 
    virtual void ProjectComponentToRefBasis(const int &c, const Array2D<Operator*> &mats, Array2D<SparseMatrix *> &rom_mats) = 0;
-   virtual void ProjectComponentToDomainBasis(const int &c, const Array2D<Operator*> &mats, Array2D<SparseMatrix *> &rom_mats) = 0;
+   virtual void ProjectComponentToDomainBasis(const int &m, const Array2D<Operator*> &mats, Array2D<SparseMatrix *> &rom_mats) = 0;
    virtual void ProjectInterfaceToRefBasis(const int &c1, const int &c2, const Array2D<Operator*> &mats, Array2D<SparseMatrix *> &rom_mats) = 0;
-   virtual void ProjectInterfaceToDomainBasis(const int &c1, const int &c2, const Array2D<Operator*> &mats, Array2D<SparseMatrix *> &rom_mats) = 0;
+   virtual void ProjectInterfaceToDomainBasis(const int &m1, const int &m2, const Array2D<Operator*> &mats, Array2D<SparseMatrix *> &rom_mats) = 0;
    virtual void ProjectVariableToDomainBasis(const int &vi, const int &vj, const Array2D<Operator*> &mats, Array2D<SparseMatrix *> &rom_mats) = 0;
    virtual void ProjectGlobalToDomainBasis(const Array2D<Operator*> &mats, Array2D<SparseMatrix *> &rom_mats) = 0;
    virtual void ProjectOperatorOnReducedBasis(const Array2D<Operator*> &mats) = 0;
@@ -173,12 +184,13 @@ class ROMHandlerBase
    virtual void LiftUpFromDomainBasis(const int &i, const Vector &rom_vec, Vector &vec) = 0;
    virtual void LiftUpGlobal(const BlockVector &rom_vec, BlockVector &vec) = 0;
 
+   virtual void Solve(BlockVector &rhs, BlockVector &sol) = 0;
    virtual void Solve(BlockVector* U) = 0;
    virtual void NonlinearSolve(Operator &oper, BlockVector* U, Solver *prec=NULL) = 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 SetRomMat(BlockMatrix *input_mat, const bool init_direct_solver=true) = 0;
    virtual void SaveRomSystem(const std::string &input_prefix, const std::string type="mm") = 0;
 
    virtual void SaveBasisVisualization(const Array<FiniteElementSpace *> &fes, const std::vector<std::string> &var_names) = 0;
@@ -217,9 +229,6 @@ class MFEMROMHandler : public ROMHandlerBase
    HYPRE_BigInt sys_row_starts[2];
    HypreParMatrix *romMat_hypre = NULL;
    MUMPSSolver *mumps = NULL;
-
-   mfem::BlockVector *reduced_rhs = NULL;
-   mfem::BlockVector *reduced_sol = NULL;
 
 public:
    MFEMROMHandler(TopologyHandler *input_topol, const Array<int> &input_var_offsets,
@@ -241,9 +250,9 @@ class MFEMROMHandler : public ROMHandlerBase
    virtual void ProjectToDomainBasis(const Array<int> &idx_i, const Array<int> &idx_j, const Array2D<Operator*> &mats, Array2D<SparseMatrix *> &rom_mats);
 
    virtual void ProjectComponentToRefBasis(const int &c, const Array2D<Operator*> &mats, Array2D<SparseMatrix *> &rom_mats);
-   virtual void ProjectComponentToDomainBasis(const int &c, const Array2D<Operator*> &mats, Array2D<SparseMatrix *> &rom_mats);
+   virtual void ProjectComponentToDomainBasis(const int &m, const Array2D<Operator*> &mats, Array2D<SparseMatrix *> &rom_mats);
    virtual void ProjectInterfaceToRefBasis(const int &c1, const int &c2, const Array2D<Operator*> &mats, Array2D<SparseMatrix *> &rom_mats);
-   virtual void ProjectInterfaceToDomainBasis(const int &c1, const int &c2, const Array2D<Operator*> &mats, Array2D<SparseMatrix *> &rom_mats);
+   virtual void ProjectInterfaceToDomainBasis(const int &m1, const int &m2, const Array2D<Operator*> &mats, Array2D<SparseMatrix *> &rom_mats);
    virtual void ProjectVariableToDomainBasis(const int &vi, const int &vj, const Array2D<Operator*> &mats, Array2D<SparseMatrix *> &rom_mats);
    virtual void ProjectGlobalToDomainBasis(const Array2D<Operator*> &mats, Array2D<SparseMatrix *> &rom_mats);
    virtual void ProjectOperatorOnReducedBasis(const Array2D<Operator*> &mats);
@@ -258,12 +267,13 @@ class MFEMROMHandler : public ROMHandlerBase
    virtual void LiftUpFromDomainBasis(const int &i, const Vector &rom_vec, Vector &vec);
    virtual void LiftUpGlobal(const BlockVector &rom_vec, BlockVector &vec);
 
-   virtual void Solve(BlockVector* U);
-   virtual void NonlinearSolve(Operator &oper, BlockVector* U, Solver *prec=NULL) override;
+   void Solve(BlockVector &rhs, BlockVector &sol) override;
+   void Solve(BlockVector* U) override;
+   void NonlinearSolve(Operator &oper, BlockVector* U, Solver *prec=NULL) override;
 
    virtual void SaveOperator(const std::string input_prefix="");
    virtual void LoadOperatorFromFile(const std::string input_prefix="");
-   virtual void SetRomMat(BlockMatrix *input_mat);
+   virtual void SetRomMat(BlockMatrix *input_mat, const bool init_direct_solver=true);
    virtual void SaveRomSystem(const std::string &input_prefix, const std::string type="mm");
 
    virtual void SaveBasisVisualization(const Array<FiniteElementSpace *> &fes, const std::vector<std::string> &var_names);

include/steady_ns_solver.hpp

@@ -60,7 +60,10 @@ class SteadyNSROM : public Operator
    const int num_var = 2;
    int numSub = -1;
 
+   /* block offsets of velocity */
    Array<int> block_offsets;
+   /* block index of velocity on FOM variable */
+   Array<int> midxs;
    Array<Array<int> *> block_idxs;
    SparseMatrix *linearOp = NULL;
 
@@ -71,7 +74,7 @@ class SteadyNSROM : public Operator
    HYPRE_BigInt sys_glob_size;
    mutable HYPRE_BigInt sys_row_starts[2];
 public:
-   SteadyNSROM(SparseMatrix *linearOp_, const int numSub_, const Array<int> &block_offsets_, const bool direct_solve_=true);
+   SteadyNSROM(const int numSub_, ROMHandlerBase *rom_handler, const bool direct_solve_=true);
 
    virtual ~SteadyNSROM();
 
@@ -85,8 +88,8 @@ class SteadyNSTensorROM : public SteadyNSROM
    Array<DenseTensor *> hs; // not owned by SteadyNSTensorROM.
 
 public:
-   SteadyNSTensorROM(SparseMatrix *linearOp_, Array<DenseTensor *> &hs_, const Array<int> &block_offsets_, const bool direct_solve_=true)
-      : SteadyNSROM(linearOp_, hs_.Size(), block_offsets_, direct_solve_), hs(hs_) {}
+   SteadyNSTensorROM(ROMHandlerBase *rom_handler, Array<DenseTensor *> &hs_, const bool direct_solve_=true)
+      : SteadyNSROM(hs_.Size(), rom_handler, direct_solve_), hs(hs_) {}
 
    virtual ~SteadyNSTensorROM() {}
 
@@ -101,13 +104,12 @@ class SteadyNSEQPROM : public SteadyNSROM
    ROMInterfaceForm *itf = NULL; // not owned by SteadyNSEQPROM.
 
    Array<int> u_offsets;
-   Array<int> u_idxs;
 
    mutable Vector x_u, y_u;
 
 public:
-   SteadyNSEQPROM(SparseMatrix *linearOp_, Array<ROMNonlinearForm *> &hs_, ROMInterfaceForm *itf_,
-                  const Array<int> &block_offsets_, const bool direct_solve_=true);
+   SteadyNSEQPROM(ROMHandlerBase *rom_handler, Array<ROMNonlinearForm *> &hs_,
+                  ROMInterfaceForm *itf_, const bool direct_solve_=true);
 
    virtual ~SteadyNSEQPROM() {}
 

include/unsteady_ns_solver.hpp

@@ -59,6 +59,14 @@ friend class SteadyNSOperator;
 
    BlockOperator *Hop = NULL;
 
+   /* function coefficients for initial condition */
+   VectorCoefficient *u_ic = NULL;
+   VectorCoefficient *p_ic = NULL;
+
+   /* mass matrix operator for ROM */
+   Array<int> rom_u_offsets;
+   BlockMatrix *rom_mass = NULL;
+
 public:
    UnsteadyNSSolver();
 
@@ -81,25 +89,23 @@ friend class SteadyNSOperator;
 
    void SetParameterizedProblem(ParameterizedProblem *problem) override;
 
+   BlockVector* PrepareSnapshots(std::vector<BasisTag> &basis_tags) override;
+
    void ProjectOperatorOnReducedBasis() override
    { mfem_error("UnsteadyNSSolver::ProjectOperatorOnReducedBasis is not implemented yet!\n"); }
 
-   void SolveROM() override
-   { mfem_error("UnsteadyNSSolver::SolveROM is not implemented yet!\n"); }
+   void BuildCompROMLinElems() override;
+
+   void SolveROM() override;
+
+   void InitROMHandler() override;
 
    void BuildROMTensorElems() override
    { mfem_error("UnsteadyNSSolver::BuildROMTensorElems is not implemented yet!\n"); }
-   void TrainROMEQPElems(SampleGenerator *sample_generator) override
-   { mfem_error("UnsteadyNSSolver::TrainROMEQPElems is not implemented yet!\n"); }
-   void SaveROMNlinElems(const std::string &input_prefix) override
-   { mfem_error("UnsteadyNSSolver::SaveROMNlinElems is not implemented yet!\n"); }
-   void LoadROMNlinElems(const std::string &input_prefix) override
-   { mfem_error("UnsteadyNSSolver::LoadROMNlinElems is not implemented yet!\n"); }
-   void AssembleROMNlinOper() override
-   { mfem_error("UnsteadyNSSolver::AssembleROMNlinOper is not implemented yet!\n"); }
 
 private:
    void InitializeTimeIntegration();
+   void SetupInitialCondition(int &initial_step, double &time);
    void Step(double &time, int step);
 
    void SanityCheck(const int step)
@@ -113,6 +119,8 @@ friend class SteadyNSOperator;
    double ComputeCFL(const double dt);
    void SetTime(const double time);
 
+   void AssembleROMMat(BlockMatrix &romMat) override;
+
 };
 
 #endif

src/main_workflow.cpp

@@ -288,12 +288,14 @@ void AuxiliaryTrainROM(MPI_Comm comm, SampleGenerator *sample_generator)
    bool separate_variable_basis = config.GetOption<bool>("model_reduction/separate_variable_basis", false);
 
    /* Supremizer enrichment */
-   if ((separate_variable_basis) && ((solver_type == "stokes") || (solver_type == "steady-ns")))
+   if ((separate_variable_basis) &&
+       ((solver_type == "stokes") || (solver_type == "steady-ns") || (solver_type == "unsteady-ns")))
    {
       ParameterizedProblem *problem = InitParameterizedProblem();
       StokesSolver *solver = NULL;
       if (solver_type == "stokes")           { solver = new StokesSolver; }
       else if (solver_type == "steady-ns")   { solver = new SteadyNSSolver; }
+      else if (solver_type == "unsteady-ns")   { solver = new UnsteadyNSSolver; }
 
       if (!solver->UseRom()) mfem_error("ROM must be enabled for supremizer enrichment!\n");
 

src/multiblock_solver.cpp

@@ -272,31 +272,33 @@ void MultiBlockSolver::BuildROMLinElems()
 
 void MultiBlockSolver::AssembleROMMat()
 {
-   assert(topol_mode == TopologyHandlerMode::COMPONENT);
-   assert(rom_elems);
-
    const Array<int> *rom_block_offsets = rom_handler->GetBlockOffsets();
    BlockMatrix *romMat = new BlockMatrix(*rom_block_offsets);
    romMat->owns_blocks = true;
 
+   AssembleROMMat(*romMat);  
+
+   romMat->Finalize();
+   rom_handler->SetRomMat(romMat);
+}
+
+void MultiBlockSolver::AssembleROMMat(BlockMatrix &romMat)
+{
+   assert(topol_mode == TopologyHandlerMode::COMPONENT);
+   assert(rom_elems);
+
    // component domain matrix.
    for (int m = 0; m < numSub; m++)
    {
       int c_type = topol_handler->GetMeshType(m);
-      int num_block = 1;
-      int midx0 = m;
-      if (separate_variable_basis)
-      {
-         num_block *= num_var;
-         midx0 *= num_var;
-      }
+      int num_block = (separate_variable_basis) ? num_var : 1;
 
       Array<int> midx(num_block);
-      for (int k = 0; k < num_block; k++)
-         midx[k] = midx0 + k;
+      for (int v = 0; v < num_block; v++)
+         midx[v] = rom_handler->GetBlockIndex(m, v);
 
       MatrixBlocks *comp_mat = rom_elems->comp[c_type];
-      AddToBlockMatrix(midx, midx, *comp_mat, *romMat);
+      AddToBlockMatrix(midx, midx, *comp_mat, romMat);
 
       // boundary matrices of each component.
       Array<int> *bdr_c2g = topol_handler->GetBdrAttrComponentToGlobalMap(m);
@@ -312,7 +314,7 @@ void MultiBlockSolver::AssembleROMMat()
             continue;
 
          MatrixBlocks *bdr_mat = (*(rom_elems->bdr[c_type]))[b];
-         AddToBlockMatrix(midx, midx, *bdr_mat, *romMat);
+         AddToBlockMatrix(midx, midx, *bdr_mat, romMat);
       }  // for (int b = 0; b < bdr_c2g->Size(); b++)
    }  // for (int m = 0; m < numSub; m++)
 
@@ -325,38 +327,19 @@ void MultiBlockSolver::AssembleROMMat()
 
       const int m1 = pInfo->Mesh1;
       const int m2 = pInfo->Mesh2;
-      const int c1 = topol_handler->GetMeshType(m1);
-      const int c2 = topol_handler->GetMeshType(m2);
 
-      int num_block = 2;
-      int c1idx0 = c1, c2idx0 = c2;
-      int m1idx0 = m1, m2idx0 = m2;
-      if (separate_variable_basis)
-      {
-         num_block *= num_var;
-         c1idx0 *= num_var;
-         c2idx0 *= num_var;
-         m1idx0 *= num_var;
-         m2idx0 *= num_var;
-      }
+      int num_block = (separate_variable_basis) ? num_var : 1;
 
-      Array<int> midx(num_block), num_basis(num_block);
-      for (int k = 0, cidx = c1idx0, m = m1idx0; k < num_block/2; k++, cidx++, m++)
-      {
-         num_basis[k] = rom_handler->GetRefNumBasis(cidx);
-         midx[k] = m;
-      }
-      for (int k = num_block/2, cidx = c2idx0, m = m2idx0; k < num_block; k++, cidx++, m++)
-      {
-         num_basis[k] = rom_handler->GetRefNumBasis(cidx);
-         midx[k] = m;
-      }
+      Array<int> midx(0);
 
-      AddToBlockMatrix(midx, midx, *port_mat, *romMat);
-   }
+      for (int v = 0; v < num_block; v++)
+         midx.Append(rom_handler->GetBlockIndex(m1, v));
 
-   romMat->Finalize();
-   rom_handler->SetRomMat(romMat);
+      for (int v = 0; v < num_block; v++)
+         midx.Append(rom_handler->GetBlockIndex(m2, v));
+
+      AddToBlockMatrix(midx, midx, *port_mat, romMat);
+   }
 }
 
 void MultiBlockSolver::InitVisualization(const std::string& output_path)