From 97b53555f8735636bb6eb3e0da216305384f7427 Mon Sep 17 00:00:00 2001 From: "Kevin\" Seung Whan Chung" Date: Tue, 9 Jul 2024 09:16:42 -0700 Subject: [PATCH] EQP ROM workflow for unsteady NS (#50) * ROMHandler: support block ordering by variable. * MultiBlockSolver::AssembleROMMat kernel * UnsteadyNSSolver: assemble ROM Mat. * UnsteadyNSSolver::SolveROM verified. * minor fix * UnsteadyNSSolver::SetupInitialCondition * minor fix * SampleGenerator::SaveSnapshot- indexing bug hot fix * fix modulo by zero bug.. --- include/multiblock_solver.hpp | 3 + include/rom_element_collection.hpp | 2 + include/rom_handler.hpp | 48 +++-- include/steady_ns_solver.hpp | 14 +- include/unsteady_ns_solver.hpp | 28 ++- src/main_workflow.cpp | 4 +- src/multiblock_solver.cpp | 67 +++---- src/rom_element_collection.cpp | 8 + src/rom_handler.cpp | 121 +++++++++--- src/sample_generator.cpp | 3 +- src/steady_ns_solver.cpp | 51 +++-- src/stokes_solver.cpp | 12 +- src/unsteady_ns_solver.cpp | 253 ++++++++++++++++++++++--- test/gmsh/CMakeLists.txt | 1 + test/gmsh/test_multi_comp_workflow.cpp | 32 ++++ test/gmsh/usns.periodic.yml | 95 ++++++++++ test/test_workflow.cpp | 73 +++++++ 17 files changed, 657 insertions(+), 158 deletions(-) create mode 100644 test/gmsh/usns.periodic.yml diff --git a/include/multiblock_solver.hpp b/include/multiblock_solver.hpp index 6901cb26..c215ef49 100644 --- a/include/multiblock_solver.hpp +++ b/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 fom_sols, Array rom_sols, Vector &error); void CompareSolution(BlockVector &test_U, Vector &error); + +protected: + virtual void AssembleROMMat(BlockMatrix &romMat); }; #endif diff --git a/include/rom_element_collection.hpp b/include/rom_element_collection.hpp index 814a217e..14b4ce88 100644 --- a/include/rom_element_collection.hpp +++ b/include/rom_element_collection.hpp @@ -50,6 +50,8 @@ class ROMLinearElement : public ROMElementCollection Array *> bdr; Array port; // reference ports. + Array mass; // Size(num_components); + public: ROMLinearElement(TopologyHandler *topol_handler_, const Array &fes_, diff --git a/include/rom_handler.hpp b/include/rom_handler.hpp index 1eb373e7..490ce95e 100644 --- a/include/rom_handler.hpp +++ b/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 num_basis; Array 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 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 &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* GetBlockOffsets() { return &rom_block_offsets; } - const Array* 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 &num_ref_supreme, Array &num_supreme); @@ -157,9 +168,9 @@ class ROMHandlerBase virtual void ProjectToDomainBasis(const Array &idx_i, const Array &idx_j, const Array2D &mats, Array2D &rom_mats) = 0; virtual void ProjectComponentToRefBasis(const int &c, const Array2D &mats, Array2D &rom_mats) = 0; - virtual void ProjectComponentToDomainBasis(const int &c, const Array2D &mats, Array2D &rom_mats) = 0; + virtual void ProjectComponentToDomainBasis(const int &m, const Array2D &mats, Array2D &rom_mats) = 0; virtual void ProjectInterfaceToRefBasis(const int &c1, const int &c2, const Array2D &mats, Array2D &rom_mats) = 0; - virtual void ProjectInterfaceToDomainBasis(const int &c1, const int &c2, const Array2D &mats, Array2D &rom_mats) = 0; + virtual void ProjectInterfaceToDomainBasis(const int &m1, const int &m2, const Array2D &mats, Array2D &rom_mats) = 0; virtual void ProjectVariableToDomainBasis(const int &vi, const int &vj, const Array2D &mats, Array2D &rom_mats) = 0; virtual void ProjectGlobalToDomainBasis(const Array2D &mats, Array2D &rom_mats) = 0; virtual void ProjectOperatorOnReducedBasis(const Array2D &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 &fes, const std::vector &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 &input_var_offsets, @@ -241,9 +250,9 @@ class MFEMROMHandler : public ROMHandlerBase virtual void ProjectToDomainBasis(const Array &idx_i, const Array &idx_j, const Array2D &mats, Array2D &rom_mats); virtual void ProjectComponentToRefBasis(const int &c, const Array2D &mats, Array2D &rom_mats); - virtual void ProjectComponentToDomainBasis(const int &c, const Array2D &mats, Array2D &rom_mats); + virtual void ProjectComponentToDomainBasis(const int &m, const Array2D &mats, Array2D &rom_mats); virtual void ProjectInterfaceToRefBasis(const int &c1, const int &c2, const Array2D &mats, Array2D &rom_mats); - virtual void ProjectInterfaceToDomainBasis(const int &c1, const int &c2, const Array2D &mats, Array2D &rom_mats); + virtual void ProjectInterfaceToDomainBasis(const int &m1, const int &m2, const Array2D &mats, Array2D &rom_mats); virtual void ProjectVariableToDomainBasis(const int &vi, const int &vj, const Array2D &mats, Array2D &rom_mats); virtual void ProjectGlobalToDomainBasis(const Array2D &mats, Array2D &rom_mats); virtual void ProjectOperatorOnReducedBasis(const Array2D &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 &fes, const std::vector &var_names); diff --git a/include/steady_ns_solver.hpp b/include/steady_ns_solver.hpp index 88cf719c..13ee520d 100644 --- a/include/steady_ns_solver.hpp +++ b/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 block_offsets; + /* block index of velocity on FOM variable */ + Array midxs; Array *> 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 &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 hs; // not owned by SteadyNSTensorROM. public: - SteadyNSTensorROM(SparseMatrix *linearOp_, Array &hs_, const Array &block_offsets_, const bool direct_solve_=true) - : SteadyNSROM(linearOp_, hs_.Size(), block_offsets_, direct_solve_), hs(hs_) {} + SteadyNSTensorROM(ROMHandlerBase *rom_handler, Array &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 u_offsets; - Array u_idxs; mutable Vector x_u, y_u; public: - SteadyNSEQPROM(SparseMatrix *linearOp_, Array &hs_, ROMInterfaceForm *itf_, - const Array &block_offsets_, const bool direct_solve_=true); + SteadyNSEQPROM(ROMHandlerBase *rom_handler, Array &hs_, + ROMInterfaceForm *itf_, const bool direct_solve_=true); virtual ~SteadyNSEQPROM() {} diff --git a/include/unsteady_ns_solver.hpp b/include/unsteady_ns_solver.hpp index 939f2cc2..fd3e367d 100644 --- a/include/unsteady_ns_solver.hpp +++ b/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 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 &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 diff --git a/src/main_workflow.cpp b/src/main_workflow.cpp index 5abb72e0..6a1202c5 100644 --- a/src/main_workflow.cpp +++ b/src/main_workflow.cpp @@ -288,12 +288,14 @@ void AuxiliaryTrainROM(MPI_Comm comm, SampleGenerator *sample_generator) bool separate_variable_basis = config.GetOption("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"); diff --git a/src/multiblock_solver.cpp b/src/multiblock_solver.cpp index acb387d8..d946c710 100644 --- a/src/multiblock_solver.cpp +++ b/src/multiblock_solver.cpp @@ -272,31 +272,33 @@ void MultiBlockSolver::BuildROMLinElems() void MultiBlockSolver::AssembleROMMat() { - assert(topol_mode == TopologyHandlerMode::COMPONENT); - assert(rom_elems); - const Array *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 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 *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 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 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) diff --git a/src/rom_element_collection.cpp b/src/rom_element_collection.cpp index 5c2e4fb6..47ffa6d5 100644 --- a/src/rom_element_collection.cpp +++ b/src/rom_element_collection.cpp @@ -12,9 +12,13 @@ ROMLinearElement::ROMLinearElement( : ROMElementCollection(topol_handler_, fes_, separate_variable_) { comp.SetSize(num_comp); + mass.SetSize(num_comp); const int block_size = (separate_variable) ? num_var : 1; for (int c = 0; c < num_comp; c++) + { comp[c] = new MatrixBlocks(block_size, block_size); + mass[c] = new MatrixBlocks(block_size, block_size); + } bdr.SetSize(num_comp); for (int c = 0; c < num_comp; c++) @@ -80,6 +84,8 @@ void ROMLinearElement::SaveCompBdrElems(hid_t &file_id) hdf5_utils::WriteDataset(comp_grp_id, "domain", *comp[c]); + hdf5_utils::WriteDataset(comp_grp_id, "mass", *mass[c]); + // boundaries are saved for each component group. SaveBdrElems(comp_grp_id, c); @@ -182,6 +188,8 @@ void ROMLinearElement::LoadCompBdrElems(hid_t &file_id) hdf5_utils::ReadDataset(comp_grp_id, "domain", *comp[c]); + hdf5_utils::ReadDataset(comp_grp_id, "mass", *mass[c]); + // boundary LoadBdrElems(comp_grp_id, c); diff --git a/src/rom_handler.cpp b/src/rom_handler.cpp index 3c400994..bfeed644 100644 --- a/src/rom_handler.cpp +++ b/src/rom_handler.cpp @@ -63,6 +63,8 @@ ROMHandlerBase::ROMHandlerBase( ROMHandlerBase::~ROMHandlerBase() { DeletePointers(carom_ref_basis); + delete reduced_rhs; + delete reduced_sol; } void ROMHandlerBase::ParseInputs() @@ -156,6 +158,12 @@ void ROMHandlerBase::ParseInputs() if (nlin_handle_str == "tensor") nlin_handle = NonlinearHandling::TENSOR; else if (nlin_handle_str == "eqp") nlin_handle = NonlinearHandling::EQP; assert((!nonlinear_mode) || (nlin_handle != NonlinearHandling::NUM_NLNHNDL)); + + std::string ordering_str = config.GetOption("model_reduction/ordering", "domain"); + if (ordering_str == "domain") ordering = ROMOrderBy::DOMAIN; + else if (ordering_str == "variable") ordering = ROMOrderBy::VARIABLE; + else + mfem_error("ROMHandlerBase: unknown ordering!\n"); } void ROMHandlerBase::ParseSupremizerInput(Array &num_ref_supreme, Array &num_supreme) @@ -183,6 +191,46 @@ void ROMHandlerBase::ParseSupremizerInput(Array &num_ref_supreme, ArrayGetMeshType(k); - for (int v = 0; v < num_var; v++, idx++) + for (int v = 0; v < num_var; v++) { + idx = GetBlockIndex(k, v); + rom_block_offsets[idx+1] = num_ref_basis[c * num_var + v]; - rom_varblock_offsets[1 + k + v * numSub] = num_ref_basis[c * num_var + v]; - } - } + } // for (int v = 0; v < num_var; v++) + } // for (int k = 0; k < numSub; k++) } else { @@ -264,14 +312,12 @@ void ROMHandlerBase::SetBlockSizes() int c = topol_handler->GetMeshType(k); rom_block_offsets[k+1] = num_ref_basis[c]; } - rom_varblock_offsets = rom_block_offsets; } rom_block_offsets.GetSubArray(1, num_rom_blocks, num_basis); rom_block_offsets.PartialSum(); rom_comp_block_offsets.PartialSum(); - rom_varblock_offsets.PartialSum(); } /* @@ -316,8 +362,6 @@ MFEMROMHandler::~MFEMROMHandler() DeletePointers(ref_basis); delete romMat; delete romMat_mono; - delete reduced_rhs; - delete reduced_sol; delete romMat_hypre; delete mumps; } @@ -336,12 +380,13 @@ void MFEMROMHandler::LoadReducedBasis() } dom_basis.SetSize(num_rom_blocks); + int m, c, v, idx; for (int k = 0; k < num_rom_blocks; k++) { - int m = (separate_variable) ? k / num_var : k; - int c = topol_handler->GetMeshType(m); - int v = (separate_variable) ? k % num_var : 0; - int idx = (separate_variable) ? c * num_var + v : c; + GetDomainAndVariableIndex(k, m, v); + + c = topol_handler->GetMeshType(m); + idx = (separate_variable) ? c * num_var + v : c; dom_basis[k] = ref_basis[idx]; } @@ -432,8 +477,14 @@ void MFEMROMHandler::ProjectGlobalToDomainBasis(const BlockVector* vec, BlockVec rom_vec = new BlockVector(rom_block_offsets); + int m, v, fom_idx; for (int i = 0; i < num_rom_blocks; i++) - ProjectToDomainBasis(i, vec->GetBlock(i), rom_vec->GetBlock(i)); + { + GetDomainAndVariableIndex(i, m, v); + fom_idx = (separate_variable)? v + m * num_var : m; + + ProjectToDomainBasis(i, vec->GetBlock(fom_idx), rom_vec->GetBlock(i)); + } } void MFEMROMHandler::LiftUpFromRefBasis(const int &i, const Vector &rom_vec, Vector &vec) @@ -467,19 +518,20 @@ void MFEMROMHandler::LiftUpGlobal(const BlockVector &rom_vec, BlockVector &vec) assert(rom_vec.NumBlocks() == num_rom_blocks); assert(vec.NumBlocks() == num_rom_blocks); + int m, v, fom_idx; for (int i = 0; i < num_rom_blocks; i++) - LiftUpFromDomainBasis(i, rom_vec.GetBlock(i), vec.GetBlock(i)); + { + GetDomainAndVariableIndex(i, m, v); + fom_idx = (separate_variable)? v + m * num_var : m; + + LiftUpFromDomainBasis(i, rom_vec.GetBlock(i), vec.GetBlock(fom_idx)); + } } -void MFEMROMHandler::Solve(BlockVector* U) +void MFEMROMHandler::Solve(BlockVector &rhs, BlockVector &sol) { - assert(U->NumBlocks() == num_rom_blocks); assert(operator_loaded); - printf("Solve ROM.\n"); - reduced_sol = new BlockVector(rom_block_offsets); - (*reduced_sol) = 0.0; - int maxIter = config.GetOption("solver/max_iter", 10000); double rtol = config.GetOption("solver/relative_tolerance", 1.e-15); double atol = config.GetOption("solver/absolute_tolerance", 1.e-15); @@ -490,7 +542,7 @@ void MFEMROMHandler::Solve(BlockVector* U) { assert(mumps); mumps->SetPrintLevel(print_level); - mumps->Mult(*reduced_rhs, *reduced_sol); + mumps->Mult(rhs, sol); } else { @@ -548,7 +600,7 @@ void MFEMROMHandler::Solve(BlockVector* U) // StopWatch solveTimer; // solveTimer.Start(); - solver->Mult(*reduced_rhs, *reduced_sol); + solver->Mult(rhs, sol); // solveTimer.Stop(); // printf("ROM-solve-only time: %f seconds.\n", solveTimer.RealTime()); @@ -558,6 +610,18 @@ void MFEMROMHandler::Solve(BlockVector* U) delete M; delete solver; } +} + +void MFEMROMHandler::Solve(BlockVector* U) +{ + assert(U->NumBlocks() == num_rom_blocks); + assert(reduced_rhs); + + printf("Solve ROM.\n"); + reduced_sol = new BlockVector(rom_block_offsets); + (*reduced_sol) = 0.0; + + Solve(*reduced_rhs, *reduced_sol); // 23. reconstruct FOM state LiftUpGlobal(*reduced_sol, *U); @@ -847,7 +911,7 @@ void MFEMROMHandler::LoadOperatorFromFile(const std::string filename) assert(errf >= 0); } -void MFEMROMHandler::SetRomMat(BlockMatrix *input_mat) +void MFEMROMHandler::SetRomMat(BlockMatrix *input_mat, const bool init_direct_solver) { if (romMat != input_mat) { @@ -858,7 +922,8 @@ void MFEMROMHandler::SetRomMat(BlockMatrix *input_mat) delete romMat_mono; romMat_mono = romMat->CreateMonolithic(); - if (linsol_type == SolverType::DIRECT) SetupDirectSolver(); + if ((linsol_type == SolverType::DIRECT) && (init_direct_solver)) + SetupDirectSolver(); operator_loaded = true; } @@ -976,7 +1041,7 @@ IterativeSolver* MFEMROMHandler::SetIterativeSolver(const MFEMROMHandler::Solver void MFEMROMHandler::SetupDirectSolver() { // If nonlinear mode, Jacobian will keep changing within Solve, thus no need of initial LU factorization. - if ((linsol_type != MFEMROMHandler::SolverType::DIRECT) || nonlinear_mode) + if ((linsol_type != MFEMROMHandler::SolverType::DIRECT)) return; assert(romMat_mono); diff --git a/src/sample_generator.cpp b/src/sample_generator.cpp index 4e51a26e..8703772b 100644 --- a/src/sample_generator.cpp +++ b/src/sample_generator.cpp @@ -198,7 +198,8 @@ void SampleGenerator::SaveSnapshot(BlockVector *U_snapshots, std::vectorgetNumSamples(); + /* 0-based index */ + col_idxs[s] = snapshot_generators[index]->getNumSamples() - 1; } } diff --git a/src/steady_ns_solver.cpp b/src/steady_ns_solver.cpp index c99300fb..024d6443 100644 --- a/src/steady_ns_solver.cpp +++ b/src/steady_ns_solver.cpp @@ -134,9 +134,11 @@ Operator& SteadyNSOperator::GetGradient(const Vector &x) const */ SteadyNSROM::SteadyNSROM( - SparseMatrix *linearOp_, const int numSub_, const Array &block_offsets_, const bool direct_solve_) - : Operator(linearOp_->Height(), linearOp_->Width()), linearOp(linearOp_), - numSub(numSub_), block_offsets(block_offsets_), direct_solve(direct_solve_) + const int numSub_, ROMHandlerBase *rom_handler, const bool direct_solve_) + : Operator(rom_handler->GetOperator()->Height(), rom_handler->GetOperator()->Width()), + linearOp(rom_handler->GetOperator()), + numSub(numSub_), block_offsets(*rom_handler->GetBlockOffsets()), + direct_solve(direct_solve_) { separate_variable = (block_offsets.Size() == num_var * numSub + 1); assert(separate_variable || (block_offsets.Size() == numSub + 1)); @@ -147,9 +149,12 @@ SteadyNSROM::SteadyNSROM( sys_row_starts[1] = height; block_idxs.SetSize(numSub); + midxs.SetSize(numSub); for (int m = 0; m < numSub; m++) { - int midx = (separate_variable) ? m * num_var : m; + int midx = rom_handler->GetBlockIndex(m, 0); + midxs[m] = midx; + block_idxs[m] = new Array(block_offsets[midx+1] - block_offsets[midx]); for (int k = 0, idx = block_offsets[midx]; k < block_idxs[m]->Size(); k++, idx++) (*block_idxs[m])[k] = idx; @@ -172,7 +177,7 @@ void SteadyNSTensorROM::Mult(const Vector &x, Vector &y) const for (int m = 0; m < numSub; m++) { - int midx = (separate_variable) ? m * num_var : m; + int midx = midxs[m]; x_comp.MakeRef(const_cast(x), block_offsets[midx], block_offsets[midx+1] - block_offsets[midx]); y_comp.MakeRef(y, block_offsets[midx], block_offsets[midx+1] - block_offsets[midx]); @@ -189,7 +194,7 @@ Operator& SteadyNSTensorROM::GetGradient(const Vector &x) const for (int m = 0; m < numSub; m++) { - int midx = (separate_variable) ? m * num_var : m; + int midx = midxs[m]; x_comp.MakeRef(const_cast(x), block_offsets[midx], block_offsets[midx+1] - block_offsets[midx]); jac_comp.SetSize(block_offsets[midx+1] - block_offsets[midx]); @@ -215,17 +220,15 @@ Operator& SteadyNSTensorROM::GetGradient(const Vector &x) const */ SteadyNSEQPROM::SteadyNSEQPROM( - SparseMatrix *linearOp_, Array &hs_, ROMInterfaceForm *itf_, - const Array &block_offsets_, const bool direct_solve_) - : SteadyNSROM(linearOp_, hs_.Size(), block_offsets_, direct_solve_), hs(hs_), itf(itf_) + ROMHandlerBase *rom_handler, Array &hs_, + ROMInterfaceForm *itf_, const bool direct_solve_) + : SteadyNSROM(hs_.Size(), rom_handler, direct_solve_), hs(hs_), itf(itf_) { if (!separate_variable) { u_offsets = block_offsets; x_u.SetSize(block_offsets.Last()); y_u.SetSize(block_offsets.Last()); - u_idxs.SetSize(numSub); - for (int m = 0; m < numSub; m++) u_idxs[m] = m; return; } @@ -239,10 +242,6 @@ SteadyNSEQPROM::SteadyNSEQPROM( x_u.SetSize(u_offsets.Last()); y_u.SetSize(u_offsets.Last()); - - u_idxs.SetSize(numSub); - for (int m = 0; m < numSub; m++) - u_idxs[m] = m * num_var; } void SteadyNSEQPROM::Mult(const Vector &x, Vector &y) const @@ -252,7 +251,7 @@ void SteadyNSEQPROM::Mult(const Vector &x, Vector &y) const for (int m = 0; m < numSub; m++) { - int midx = (separate_variable) ? m * num_var : m; + int midx = midxs[m]; x_comp.MakeRef(const_cast(x), block_offsets[midx], block_offsets[midx+1] - block_offsets[midx]); y_comp.MakeRef(y, block_offsets[midx], block_offsets[midx+1] - block_offsets[midx]); @@ -282,8 +281,8 @@ Operator& SteadyNSEQPROM::GetGradient(const Vector &x) const for (int i = 0; i < numSub; i++) for (int j = 0; j < numSub; j++) { - int iidx = (separate_variable) ? i * num_var : i; - int jidx = (separate_variable) ? j * num_var : j; + int iidx = midxs[i]; + int jidx = midxs[j]; mats(i, j) = new SparseMatrix(block_offsets[iidx+1] - block_offsets[iidx], block_offsets[jidx+1] - block_offsets[jidx]); } @@ -293,7 +292,7 @@ Operator& SteadyNSEQPROM::GetGradient(const Vector &x) const MatrixBlocks u_mats(mats); - AddToBlockMatrix(u_idxs, u_idxs, u_mats, jac); + AddToBlockMatrix(midxs, midxs, u_mats, jac); jac.Finalize(); SparseMatrix *tmp = jac.CreateMonolithic(); (*jac_mono) += *tmp; @@ -304,7 +303,7 @@ Operator& SteadyNSEQPROM::GetGradient(const Vector &x) const DenseMatrix *jac_comp; for (int m = 0; m < numSub; m++) { - int midx = (separate_variable) ? m * num_var : m; + int midx = midxs[m]; x_comp.MakeRef(const_cast(x), block_offsets[midx], block_offsets[midx+1] - block_offsets[midx]); // NOTE(kevin): jac_comp is owned by hs[m]. No need of deleting it. @@ -337,7 +336,7 @@ void SteadyNSEQPROM::GetVel(const Vector &x, Vector &x_u) const { Vector tmp; u_block.GetBlockView(m, tmp); - tmp = x_block.GetBlock(m * num_var); + tmp = x_block.GetBlock(midxs[m]); } } @@ -354,7 +353,7 @@ void SteadyNSEQPROM::AddVel(const Vector &y_u, Vector &y) const for (int m = 0; m < numSub; m++) { Vector tmp; - y_block.GetBlockView(m * num_var, tmp); + y_block.GetBlockView(midxs[m], tmp); tmp += u_block.GetBlock(m); } } @@ -694,7 +693,7 @@ void SteadyNSSolver::ProjectOperatorOnReducedBasis() DenseMatrix *basis = NULL; for (int m = 0; m < numSub; m++) { - int idx = (separate_variable_basis) ? m * num_var : m; + int idx = rom_handler->GetBlockIndex(m, 0); rom_handler->GetDomainBasis(idx, basis); subdomain_tensors[m] = GetReducedTensor(basis, ufes[m]); } @@ -723,13 +722,13 @@ void SteadyNSSolver::SolveROM() { assert(subdomain_tensors.Size() == numSub); for (int m = 0; m < numSub; m++) assert(subdomain_tensors[m]); - rom_oper = new SteadyNSTensorROM(rom_handler->GetOperator(), subdomain_tensors, *(rom_handler->GetBlockOffsets())); + rom_oper = new SteadyNSTensorROM(rom_handler, subdomain_tensors); } else if (rom_handler->GetNonlinearHandling() == NonlinearHandling::EQP) { assert(subdomain_eqps.Size() == numSub); for (int m = 0; m < numSub; m++) assert(subdomain_eqps[m]); - rom_oper = new SteadyNSEQPROM(rom_handler->GetOperator(), subdomain_eqps, itf_eqp, *(rom_handler->GetBlockOffsets())); + rom_oper = new SteadyNSEQPROM(rom_handler, subdomain_eqps, itf_eqp); } rom_handler->NonlinearSolve(*rom_oper, U_domain); @@ -1244,7 +1243,7 @@ void SteadyNSSolver::AssembleROMEQPOper() for (int m = 0; m < numSub; m++) { int c_type = topol_handler->GetMeshType(m); - int idx = (separate_variable_basis) ? m * num_var : m; + int idx = rom_handler->GetBlockIndex(m, 0); rom_handler->GetDomainBasis(idx, basis); subdomain_eqps[m] = new ROMNonlinearForm(basis->NumCols(), ufes[m], false); diff --git a/src/stokes_solver.cpp b/src/stokes_solver.cpp index c6d0fc7d..2ff664a1 100644 --- a/src/stokes_solver.cpp +++ b/src/stokes_solver.cpp @@ -1018,6 +1018,7 @@ void StokesSolver::ProjectOperatorOnReducedBasis() Array2D *> ioffsets(numSub, numSub), joffsets(numSub, numSub); tmp = NULL; ioffsets = NULL; joffsets = NULL; + int ui, pi, uj, pj; for (int i = 0; i < numSub; i++) for (int j = 0; j < numSub; j++) { @@ -1026,9 +1027,14 @@ void StokesSolver::ProjectOperatorOnReducedBasis() if (separate_variable_basis) { - tmp(i * num_var, j * num_var) = m_mats(i, j); - tmp(i * num_var + 1, j * num_var) = b_mats(i, j); - tmp(i * num_var, j * num_var + 1) = bt_mats(i, j); + ui = rom_handler->GetBlockIndex(i, 0); + pi = rom_handler->GetBlockIndex(i, 1); + uj = rom_handler->GetBlockIndex(j, 0); + pj = rom_handler->GetBlockIndex(j, 1); + + tmp(ui, uj) = m_mats(i, j); + tmp(pi, uj) = b_mats(i, j); + tmp(ui, pj) = bt_mats(i, j); } else { diff --git a/src/unsteady_ns_solver.cpp b/src/unsteady_ns_solver.cpp index a61e2d1f..6cdd7ec6 100644 --- a/src/unsteady_ns_solver.cpp +++ b/src/unsteady_ns_solver.cpp @@ -31,6 +31,9 @@ UnsteadyNSSolver::UnsteadyNSSolver() if (time_order != 1) mfem_error("UnsteadyNSSolver supports only first-order time integration for now!\n"); + + if (use_rom && !separate_variable_basis) + mfem_error("UnsteadyNSSolver does not allow unified basis for all variables!\n"); } UnsteadyNSSolver::~UnsteadyNSSolver() @@ -43,6 +46,10 @@ UnsteadyNSSolver::~UnsteadyNSSolver() delete U_stepview; delete RHS_stepview; delete Hop; + delete rom_mass; + + delete u_ic; + delete p_ic; } void UnsteadyNSSolver::BuildDomainOperators() @@ -78,14 +85,10 @@ bool UnsteadyNSSolver::Solve(SampleGenerator *sample_generator) int initial_step = 0; double time = 0.0; - bool use_restart = config.GetOption("solver/use_restart", false); std::string restart_file, file_fmt; file_fmt = "%s/%s_%08d.h5"; - if (use_restart) - { - restart_file = config.GetRequiredOption("solver/restart_file"); - LoadSolutionWithTime(restart_file, initial_step, time); - } + + SetupInitialCondition(initial_step, time); int sample_interval = config.GetOption("sample_generation/time-integration/sample_interval", 0); int bootstrap = config.GetOption("sample_generation/time-integration/bootstrap", 0); @@ -103,20 +106,24 @@ bool UnsteadyNSSolver::Solve(SampleGenerator *sample_generator) cfl = ComputeCFL(dt); SanityCheck(step); - if (((step+1) % report_interval) == 0) + if (report_interval && + ((step+1) % report_interval) == 0) printf("Time step: %05d, CFL: %.3e\n", step+1, cfl); - if (((step+1) % visual.time_interval) == 0) + if (visual.time_interval && + ((step+1) % visual.time_interval) == 0) SaveVisualization(step+1, time); - if ((save_sol) && (((step+1) % restart_interval) == 0)) + if (save_sol && restart_interval && + (((step+1) % restart_interval) == 0)) { restart_file = string_format(file_fmt, sol_dir.c_str(), sol_prefix.c_str(), step+1); SaveSolutionWithTime(restart_file, step+1, time); } /* save solution if sample generator is provided */ - if (sample_generator && ((step+1) > bootstrap) && (((step+1) % sample_interval) == 0)) + if (sample_generator && sample_interval && + ((step+1) > bootstrap) && (((step+1) % sample_interval) == 0)) SaveSnapshots(sample_generator); } @@ -149,6 +156,8 @@ void UnsteadyNSSolver::InitializeTimeIntegration() for (int m = 0; m < numSub; m++) Hop->SetDiagonalBlock(m, hs[m]); + /* the routines above can be merged to AssembleOperator */ + offsets_byvar.SetSize(num_var * numSub + 1); offsets_byvar = 0; for (int k = 0; k < numSub; k++) @@ -166,6 +175,37 @@ void UnsteadyNSSolver::InitializeTimeIntegration() Cu1.SetSize(U_stepview->BlockSize(0)); } +void UnsteadyNSSolver::SetupInitialCondition(int &initial_step, double &time) +{ + bool use_restart = config.GetOption("solver/use_restart", false); + std::string restart_file, file_fmt; + file_fmt = "%s/%s_%08d.h5"; + + if (use_restart) + { + restart_file = config.GetRequiredOption("solver/restart_file"); + LoadSolutionWithTime(restart_file, initial_step, time); + } + else + { + if (u_ic && p_ic) + for (int m = 0; m < numSub; m++) + { + vels[m]->ProjectCoefficient(*u_ic); + ps[m]->ProjectCoefficient(*p_ic); + } + + initial_step = 0; + time = 0.0; + } + + if ((!use_restart) && save_sol) + { + restart_file = string_format(file_fmt, sol_dir.c_str(), sol_prefix.c_str(), initial_step); + SaveSolutionWithTime(restart_file, initial_step, time); + } +} + void UnsteadyNSSolver::Step(double &time, int step) { /* set time for forcing/boundary. At this point, time remains at the previous timestep. */ @@ -221,15 +261,20 @@ void UnsteadyNSSolver::SetParameterizedProblem(ParameterizedProblem *problem) *ps[m] = 0.0; } + Vector zero_vel(dim), zero_pres(1); + zero_vel = 0.0; zero_pres = 0.0; + + if (u_ic) delete u_ic; + if (p_ic) delete p_ic; + /* set up initial condition */ bool ic_setup = (problem->ic_ptr.Size() >= 2); if (!ic_setup) + { + u_ic = new VectorConstantCoefficient(zero_vel); + p_ic = new VectorConstantCoefficient(zero_pres); return; - - Vector zero_vel(dim), zero_pres(1); - zero_vel = 0.0; zero_pres = 0.0; - - VectorCoefficient *u_ic, *p_ic; + } if (problem->ic_ptr[0]) u_ic = new VectorFunctionCoefficient(vdim[0], problem->ic_ptr[0]); @@ -240,15 +285,14 @@ void UnsteadyNSSolver::SetParameterizedProblem(ParameterizedProblem *problem) p_ic = new VectorFunctionCoefficient(1, problem->ic_ptr[1]); else p_ic = new VectorConstantCoefficient(zero_pres); +} - for (int m = 0; m < numSub; m++) - { - vels[m]->ProjectCoefficient(*u_ic); - ps[m]->ProjectCoefficient(*p_ic); - } +BlockVector* UnsteadyNSSolver::PrepareSnapshots(std::vector &basis_tags) +{ + /* copy to original solution variable */ + SortBySubdomains(*U_step, *U); - delete u_ic; - delete p_ic; + return MultiBlockSolver::PrepareSnapshots(basis_tags); } double UnsteadyNSSolver::ComputeCFL(const double dt_) @@ -314,4 +358,169 @@ void UnsteadyNSSolver::SetTime(const double time) if (ud_coeffs[k]) ud_coeffs[k]->SetTime(time); for (int k = 0; k < sn_coeffs.Size(); k++) if (sn_coeffs[k]) sn_coeffs[k]->SetTime(time); +} + +void UnsteadyNSSolver::InitROMHandler() +{ + SteadyNSSolver::InitROMHandler(); + + if (rom_handler->GetOrdering() != ROMOrderBy::VARIABLE) + mfem_error("UnsteadyNSSolver::InitROMHandler- unsteady NS solver only allows ROM ordering by variable!\n"); +} + +void UnsteadyNSSolver::BuildCompROMLinElems() +{ + SteadyNSSolver::BuildCompROMLinElems(); + + if (!separate_variable_basis) + mfem_error("UnsteadyNSSolver::BuildCompROMLinElems- unified basis is not supported!\n"); + + /* We only need velocity mass matrix */ + const int num_comp = topol_handler->GetNumComponents(); + for (int c = 0; c < num_comp; c++) + { + const int fidx = c * num_var; + Mesh *comp = topol_handler->GetComponentMesh(c); + + BilinearForm m_comp(comp_fes[fidx]); + + m_comp.AddDomainIntegrator(new VectorMassIntegrator); + m_comp.Assemble(); + m_comp.Finalize(); + + SparseMatrix *m_mat = &(m_comp.SpMat()); + + assert(rom_elems->mass[c]->nrows > 0); + assert(rom_elems->mass[c]->ncols > 0); + + if (separate_variable_basis) + (*rom_elems->mass[c])(0, 0) = rom_handler->ProjectToRefBasis(fidx, fidx, m_mat); + } +} + +void UnsteadyNSSolver::AssembleROMMat(BlockMatrix &romMat) +{ + SteadyNSSolver::AssembleROMMat(romMat); + + assert(rom_handler->GetOrdering() == ROMOrderBy::VARIABLE); + + rom_handler->GetBlockOffsets()->GetSubArray(0, numSub+1, rom_u_offsets); + rom_mass = new BlockMatrix(rom_u_offsets); + rom_mass->owns_blocks = true; + + for (int m = 0; m < numSub; m++) + { + int c_type = topol_handler->GetMeshType(m); + + Array midx(1); + midx[0] = rom_handler->GetBlockIndex(m, 0); + + MatrixBlocks mass_mat(1, 1); + mass_mat(0, 0) = new SparseMatrix(*(*rom_elems->mass[c_type])(0,0)); + + /* mass matrix itself for RHS */ + AddToBlockMatrix(midx, midx, mass_mat, *rom_mass); + + /* add mass matrices on LHS */ + *mass_mat(0, 0) *= bd0 / dt; + AddToBlockMatrix(midx, midx, mass_mat, romMat); + } // for (int m = 0; m < numSub; m++) +} + +void UnsteadyNSSolver::SolveROM() +{ + assert(rom_handler->GetOrdering() == ROMOrderBy::VARIABLE); + + int initial_step = 0; + double time = 0.0; + + std::string restart_file, file_fmt; + file_fmt = "%s/%s_%08d.h5"; + + SetupInitialCondition(initial_step, time); + + const Array *rom_block_offsets = rom_handler->GetBlockOffsets(); + + Array rom_u_offsets(numSub + 1); + Array rom_p_offsets(numSub + 1); + rom_block_offsets->GetSubArray(0, numSub + 1, rom_u_offsets); + rom_block_offsets->GetSubArray(numSub, numSub + 1, rom_p_offsets); + for (int k = numSub; k >= 0; k--) + rom_p_offsets[k] -= rom_p_offsets[0]; + + Array rom_var_offsets(num_var + 1); + rom_var_offsets[0] = 0; + rom_var_offsets[1] = (*rom_block_offsets)[numSub]; + rom_var_offsets[2] = (*rom_block_offsets)[2 * numSub]; + + BlockVector *reduced_sol = NULL; + rom_handler->ProjectGlobalToDomainBasis(U, reduced_sol); + BlockVector reduced_rhs(*rom_handler->GetReducedRHS()); + + BlockVector *rsol_view = new BlockVector(reduced_sol->GetData(), rom_var_offsets); + BlockVector *rrhs_view = new BlockVector(reduced_rhs.GetData(), rom_var_offsets); + + u1.SetSize(rsol_view->BlockSize(0)); + Cu1.SetSize(rsol_view->BlockSize(0)); + + Hop = new BlockOperator(rom_u_offsets); + for (int m = 0; m < numSub; m++) + Hop->SetDiagonalBlock(m, subdomain_eqps[m]); + + int pN = -1; + BlockVector rom_ones(rom_p_offsets); + rom_ones = 0.0; + if (!pres_dbc) + { + pN = p_offsets.Last(); + BlockVector fom_ones(p_offsets); + fom_ones = 1.0; + for (int m = 0; m < numSub; m++) + { + int idx = rom_handler->GetBlockIndex(m, 1); + rom_handler->ProjectToDomainBasis(idx, fom_ones.GetBlock(m), rom_ones.GetBlock(m)); + } + } + + for (int step = initial_step; step < nt; step++) + { + /* set time for forcing/boundary. At this point, time remains at the previous timestep. */ + SetTime(time); + + /* copy velocity */ + u1 = rsol_view->GetBlock(0); + + /* evaluate nonlinear advection at previous time step */ + Hop->Mult(u1, Cu1); + itf_eqp->InterfaceAddMult(u1, Cu1); + + /* Base right-hand side for boundary conditions and forcing */ + reduced_rhs = *(rom_handler->GetReducedRHS()); + + /* Add nonlinear convection */ + rrhs_view->GetBlock(0).Add(-ab1, Cu1); + + /* Add time derivative term */ + // TODO: extend for high order bdf schemes + rom_mass->AddMult(u1, rrhs_view->GetBlock(0), -bd1 / dt); + + /* Solve for the next step */ + rom_handler->Solve(reduced_rhs, *reduced_sol); + + /* remove pressure scalar if all dirichlet bc */ + if (!pres_dbc) + { + double p_const = (rom_ones * (rsol_view->GetBlock(1))) / pN; + + rsol_view->GetBlock(1).Add(-p_const, rom_ones); + } + + time += dt; + } + + rom_handler->LiftUpGlobal(*reduced_sol, *U); + + delete rsol_view; + delete reduced_sol; + return; } \ No newline at end of file diff --git a/test/gmsh/CMakeLists.txt b/test/gmsh/CMakeLists.txt index c43e48c1..4e297876 100644 --- a/test/gmsh/CMakeLists.txt +++ b/test/gmsh/CMakeLists.txt @@ -12,6 +12,7 @@ file(COPY test.component.yml DESTINATION ${CMAKE_BINARY_DIR}/test/gmsh/) file(COPY stokes.component.yml DESTINATION ${CMAKE_BINARY_DIR}/test/gmsh/) file(COPY steadyns.lf.yml DESTINATION ${CMAKE_BINARY_DIR}/test/gmsh/) file(COPY steadyns.interface_eqp.yml DESTINATION ${CMAKE_BINARY_DIR}/test/gmsh/) +file(COPY usns.periodic.yml DESTINATION ${CMAKE_BINARY_DIR}/test/gmsh/) ADD_CUSTOM_COMMAND( OUTPUT ${CMAKE_BINARY_DIR}/test/gmsh/square-circle.msh diff --git a/test/gmsh/test_multi_comp_workflow.cpp b/test/gmsh/test_multi_comp_workflow.cpp index 0e7355c0..9e76d863 100644 --- a/test/gmsh/test_multi_comp_workflow.cpp +++ b/test/gmsh/test_multi_comp_workflow.cpp @@ -420,6 +420,38 @@ TEST(MultiComponentGlobalROM, SteadyNSTest_SeparateVariable) return; } +TEST(UnsteadyNS_Workflow, Periodic) +{ + config = InputParser("usns.periodic.yml"); + + printf("\nSample Generation \n\n"); + + config.dict_["main"]["mode"] = "sample_generation"; + GenerateSamples(MPI_COMM_WORLD); + + config.dict_["main"]["mode"] = "train_rom"; + TrainROM(MPI_COMM_WORLD); + + config.dict_["main"]["mode"] = "train_eqp"; + TrainEQP(MPI_COMM_WORLD); + + printf("\nBuild ROM \n\n"); + + config.dict_["main"]["mode"] = "build_rom"; + BuildROM(MPI_COMM_WORLD); + + config.dict_["main"]["mode"] = "single_run"; + // config.dict_["solver"]["use_restart"] = true; + // config.dict_["solver"]["restart_file"] = "usns_restart_00000000.h5"; + double error = SingleRun(MPI_COMM_WORLD, "test_output.h5"); + + // This reproductive case must have a very small error at the level of finite-precision. + printf("Error: %.15E\n", error); + EXPECT_TRUE(error < ns_threshold); + + return; +} + int main(int argc, char* argv[]) { ::testing::InitGoogleTest(&argc, argv); diff --git a/test/gmsh/usns.periodic.yml b/test/gmsh/usns.periodic.yml new file mode 100644 index 00000000..92fa56df --- /dev/null +++ b/test/gmsh/usns.periodic.yml @@ -0,0 +1,95 @@ +main: +#mode: run_example/sample_generation/build_rom/single_run + mode: single_run + use_rom: true + solver: unsteady-ns + +navier-stokes: + operator-type: lf + +mesh: + type: component-wise + component-wise: + global_config: "box-channel.2x2.periodic.h5" + components: + - name: "empty" + file: "square.tri.mesh" + - name: "square-circle" + file: "square-circle.msh.mfem" + +domain-decomposition: + type: interior_penalty + +discretization: + order: 1 + full-discrete-galerkin: true + +solver: + direct_solve: true + +time-integration: + timestep_size: 0.01 + number_of_timesteps: 3 + +save_solution: + enabled: false + file_path: + prefix: usns_restart + +visualization: + enable: false + output_dir: dd_mms_output + +parameterized_problem: + name: periodic_flow_past_array + +single_run: + periodic_flow_past_array: + nu: 1.1 + fx: 0.5 + fy: -0.5 + +sample_generation: + maximum_number_of_snapshots: 400 + file_path: + prefix: "usns0" + parameters: + - key: single_run/periodic_flow_past_array/nu + type: double + sample_size: 1 + minimum: 1.1 + maximum: 1.1 + time-integration: + sample_interval: 1 + bootstrap: 0 + +sample_collection: + mode: port + port_fileformat: + format: usns%d_sample.port.h5 + minimum: 0 + maximum: 0 + +basis: + prefix: "usns" + number_of_basis: 6 + svd: + save_spectrum: true + update_right_sv: false + visualization: + enabled: false + +model_reduction: + separate_variable_basis: true + ordering: variable + nonlinear_handling: eqp + eqp: + relative_tolerance: 1.0e-11 + precompute: true + save_operator: + level: component + prefix: "test.rom_elem" + compare_solution: + enabled: true + linear_solver_type: direct + linear_system_type: sid diff --git a/test/test_workflow.cpp b/test/test_workflow.cpp index 8b44c0b8..f0c59c09 100644 --- a/test/test_workflow.cpp +++ b/test/test_workflow.cpp @@ -280,6 +280,40 @@ TEST(Stokes_Workflow, ComponentSeparateVariable) return; } +TEST(Stokes_Workflow, ROM_OrderingByVariable) +{ + config = InputParser("inputs/stokes.component.yml"); + + config.dict_["model_reduction"]["ordering"] = "variable"; + + config.dict_["model_reduction"]["separate_variable_basis"] = true; + config.dict_["solver"]["direct_solve"] = true; + config.dict_["model_reduction"]["linear_solver_type"] = "direct"; + config.dict_["model_reduction"]["linear_system_type"] = "us"; + + printf("\nSample Generation \n\n"); + + config.dict_["main"]["mode"] = "sample_generation"; + GenerateSamples(MPI_COMM_WORLD); + + config.dict_["main"]["mode"] = "train_rom"; + TrainROM(MPI_COMM_WORLD); + + printf("\nBuild ROM \n\n"); + + config.dict_["main"]["mode"] = "build_rom"; + BuildROM(MPI_COMM_WORLD); + + config.dict_["main"]["mode"] = "single_run"; + double error = SingleRun(MPI_COMM_WORLD); + + // This reproductive case must have a very small error at the level of finite-precision. + printf("Error: %.15E\n", error); + EXPECT_TRUE(error < stokes_threshold); + + return; +} + TEST(SteadyNS_Workflow, SubmeshTest) { config = InputParser("inputs/steady_ns.base.yml"); @@ -477,6 +511,45 @@ TEST(SteadyNS_Workflow, ComponentSeparateVariable_EQP) return; } +TEST(SteadyNS_Workflow, ROM_OrderingByVariable) +{ + config = InputParser("inputs/steady_ns.component.yml"); + + config.dict_["model_reduction"]["ordering"] = "variable"; + + config.dict_["model_reduction"]["separate_variable_basis"] = true; + config.dict_["model_reduction"]["linear_solver_type"] = "direct"; + config.dict_["model_reduction"]["linear_system_type"] = "us"; + config.dict_["model_reduction"]["nonlinear_handling"] = "eqp"; + config.dict_["model_reduction"]["eqp"]["relative_tolerance"] = 1.0e-11; + config.dict_["model_reduction"]["eqp"]["precompute"] = true; + + printf("\nSample Generation \n\n"); + + config.dict_["main"]["mode"] = "sample_generation"; + GenerateSamples(MPI_COMM_WORLD); + + config.dict_["main"]["mode"] = "train_rom"; + TrainROM(MPI_COMM_WORLD); + + config.dict_["main"]["mode"] = "train_eqp"; + TrainEQP(MPI_COMM_WORLD); + + printf("\nBuild ROM \n\n"); + + config.dict_["main"]["mode"] = "build_rom"; + BuildROM(MPI_COMM_WORLD); + + config.dict_["main"]["mode"] = "single_run"; + double error = SingleRun(MPI_COMM_WORLD, "test_output.h5"); + + // This reproductive case must have a very small error at the level of finite-precision. + printf("Error: %.15E\n", error); + EXPECT_TRUE(error < stokes_threshold); + + return; +} + TEST(LinElast_Workflow, SubmeshTest) { config = InputParser("inputs/linelast.base.yml");