minor changes before ns solver

include/multiblock_solver.hpp

@@ -92,6 +92,8 @@ friend class ParameterizedProblem;
    // boundary condition is enforced via forcing term.
    Array<Array<SparseMatrix *> *> bdr_mats;
    Array<Array2D<SparseMatrix *> *> port_mats;   // reference ports.
+   // DenseTensor objects from nonlinear operators
+   // will be defined per each derived MultiBlockSolver.
 
 public:
    MultiBlockSolver();
@@ -185,14 +187,14 @@ friend class ParameterizedProblem;
    void SaveROMElements(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.
-   void SaveCompBdrROMElement(hid_t &file_id);
+   virtual void SaveCompBdrROMElement(hid_t &file_id);
    void SaveBdrROMElement(hid_t &comp_grp_id, const int &comp_idx);
    void SaveInterfaceROMElement(hid_t &file_id);
 
    void LoadROMElements(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.
-   void LoadCompBdrROMElement(hid_t &file_id);
+   virtual void LoadCompBdrROMElement(hid_t &file_id);
    void LoadBdrROMElement(hid_t &comp_grp_id, const int &comp_idx);
    void LoadInterfaceROMElement(hid_t &file_id);
 

include/rom_handler.hpp

@@ -112,6 +112,8 @@ class ROMHandler
    int GetBasisIndexForSubdomain(const int &subdomain_index);
    void GetBasis(const int &basis_index, const CAROM::Matrix* &basis);
    virtual void GetBasisOnSubdomain(const int &subdomain_index, const CAROM::Matrix* &basis);
+   virtual Vector* GetBasisVector(const int &basis_index, const int &column_idx, const int size=-1, const int offset=0)
+   { mfem_error("ROMHandler::GetBasisVector is not supported! Use MFEMROMHandler.\n"); }
    virtual void SetBlockSizes();
    virtual void AllocROMMat();  // allocate matrixes for rom.
    // TODO: extension to nonlinear operators.
@@ -185,6 +187,7 @@ class MFEMROMHandler : public ROMHandler
    virtual void LoadReducedBasis();
    void GetBasis(const int &basis_index, DenseMatrix* &basis);
    void GetBasisOnSubdomain(const int &subdomain_index, DenseMatrix* &basis);
+   virtual Vector* GetBasisVector(const int &basis_index, const int &column_idx, const int size=-1, const int offset=0) override;
    // virtual void AllocROMMat();  // allocate matrixes for rom.
    // TODO: extension to nonlinear operators.
    virtual void ProjectOperatorOnReducedBasis(const Array2D<Operator*> &mats);

include/stokes_solver.hpp

@@ -175,6 +175,8 @@ friend class ParameterizedProblem;
    // For bilinear case.
    // system-specific.
    virtual void AssembleInterfaceMatrixes();
+   virtual void SetupMUMPSSolver();
+   virtual void SetupPressureMassMatrix();
 
    // Component-wise assembly
    virtual void BuildCompROMElement(Array<FiniteElementSpace *> &fes_comp);

src/rom_handler.cpp

@@ -547,6 +547,18 @@ void MFEMROMHandler::GetBasisOnSubdomain(const int &subdomain_index, DenseMatrix
    GetBasis(idx, basis);
 }
 
+const Vector* MFEMROMHandler::GetBasisVector(
+   const int &basis_index, const int &column_idx, const int size, const int offset)
+{
+   assert(num_basis_sets > 0);
+   assert((basis_index >= 0) && (basis_index < num_basis_sets));
+   assert((column_idx >= 0) && (column_idx < num_basis[basis_index]));
+
+   // This vector does not assume ownership of the basis data.
+   const int vec_size = (size > 0) ? size : spatialbasis[basis_index]->NumRows();
+   return new Vector(spatialbasis[basis_index]->GetColumn(column_idx) + offset, vec_size);
+}
+
 void MFEMROMHandler::ProjectOperatorOnReducedBasis(const Array2D<Operator*> &mats)
 {
    assert(mats.NumRows() == numSub);

src/stokes_solver.cpp

@@ -367,6 +367,12 @@ void StokesSolver::Assemble()
    AssembleRHS();
 
    AssembleOperator();
+
+   if (direct_solve)
+      SetupMUMPSSolver();
+   else
+      // pressure mass matrix for preconditioner.
+      SetupPressureMassMatrix();
 }
 
 void StokesSolver::AssembleRHS()
@@ -457,34 +463,43 @@ void StokesSolver::AssembleOperator()
    systemOp->SetBlock(0,0, M);
    systemOp->SetBlock(0,1, Bt);
    systemOp->SetBlock(1,0, B);
+}
 
-   if (direct_solve)
-   {
-      systemOp_mono = systemOp->CreateMonolithic();
+void StokesSolver::SetupMUMPSSolver()
+{
+   assert(systemOp);
+   delete systemOp_mono;
+   delete systemOp_hypre;
+   delete mumps;
 
-      // TODO: need to change when the actual parallelization is implemented.
-      sys_glob_size = systemOp_mono->NumRows();
-      sys_row_starts[0] = 0;
-      sys_row_starts[1] = systemOp_mono->NumRows();
-      systemOp_hypre = new HypreParMatrix(MPI_COMM_WORLD, sys_glob_size, sys_row_starts, systemOp_mono);
-
-      mumps = new MUMPSSolver();
-      mumps->SetMatrixSymType(MUMPSSolver::MatType::SYMMETRIC_POSITIVE_DEFINITE);
-      mumps->SetOperator(*systemOp_hypre);
-   }
-   else
+   systemOp_mono = systemOp->CreateMonolithic();
+
+   // TODO: need to change when the actual parallelization is implemented.
+   sys_glob_size = systemOp_mono->NumRows();
+   sys_row_starts[0] = 0;
+   sys_row_starts[1] = systemOp_mono->NumRows();
+   systemOp_hypre = new HypreParMatrix(MPI_COMM_WORLD, sys_glob_size, sys_row_starts, systemOp_mono);
+
+   mumps = new MUMPSSolver();
+   mumps->SetMatrixSymType(MUMPSSolver::MatType::SYMMETRIC_POSITIVE_DEFINITE);
+   mumps->SetOperator(*systemOp_hypre);
+}
+
+void StokesSolver::SetupPressureMassMatrix()
+{
+   assert(pms.Size() == numSub);
+   delete pmMat;
+
+   pmMat = new BlockMatrix(p_offsets);
+   for (int m = 0; m < numSub; m++)
    {
-      // pressure mass matrix for preconditioner.
-      pmMat = new BlockMatrix(p_offsets);
-      for (int m = 0; m < numSub; m++)
-      {
-         pms[m]->Assemble();
-         pms[m]->Finalize();
+      assert(pms[m]);
+      pms[m]->Assemble();
+      pms[m]->Finalize();
 
-         pmMat->SetBlock(m, m, &(pms[m]->SpMat()));
-      }
-      pM = pmMat->CreateMonolithic();
+      pmMat->SetBlock(m, m, &(pms[m]->SpMat()));
    }
+   pM = pmMat->CreateMonolithic();
 }
 
 void StokesSolver::AssembleInterfaceMatrixes()