From ec9f02712adc14894c9a14984c721fcb9d5e9257 Mon Sep 17 00:00:00 2001 From: "Kevin\" Seung Whan Chung" Date: Wed, 19 Jun 2024 12:17:01 -0700 Subject: [PATCH] Unsteady sampling (#46) * MultiBlockSolver::Solve takes sample generator as input. * MultiBlockSolver::SaveSnapshots takes sample generator. * ROM linear elements * templates for ROMTensorElement and ROMEQPElement --- CMakeLists.txt | 5 +- include/advdiff_solver.hpp | 4 +- include/linelast_solver.hpp | 8 +- include/multiblock_solver.hpp | 44 ++--- include/poisson_solver.hpp | 8 +- include/rom_element_collection.hpp | 118 +++++++++++ include/sample_generator.hpp | 4 +- include/steady_ns_solver.hpp | 3 +- include/stokes_solver.hpp | 8 +- include/unsteady_ns_solver.hpp | 2 +- src/advdiff_solver.cpp | 20 +- src/linelast_solver.cpp | 46 ++--- src/main_workflow.cpp | 16 +- src/multiblock_solver.cpp | 306 ++++------------------------- src/poisson_solver.cpp | 46 ++--- src/rom_element_collection.cpp | 244 +++++++++++++++++++++++ src/sample_generator.cpp | 37 +++- src/steady_ns_solver.cpp | 24 +-- src/stokes_solver.cpp | 62 +++--- src/unsteady_ns_solver.cpp | 9 +- 20 files changed, 569 insertions(+), 445 deletions(-) create mode 100644 include/rom_element_collection.hpp create mode 100644 src/rom_element_collection.cpp diff --git a/CMakeLists.txt b/CMakeLists.txt index e908c1a7..f9d5c1c5 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -158,8 +158,9 @@ set(scaleupROMObj_SOURCES include/block_smoother.hpp src/block_smoother.cpp - # include/navier_solver.hpp - # src/navier_solver.cpp + include/rom_element_collection.hpp + src/rom_element_collection.cpp + include/unsteady_ns_solver.hpp src/unsteady_ns_solver.cpp diff --git a/include/advdiff_solver.hpp b/include/advdiff_solver.hpp index b0b4e9fc..62a39533 100644 --- a/include/advdiff_solver.hpp +++ b/include/advdiff_solver.hpp @@ -46,9 +46,9 @@ friend class ParameterizedProblem; void BuildDomainOperators() override; // Component-wise assembly - void BuildCompROMLinElems(Array &fes_comp) override; + void BuildCompROMLinElems() override; - bool Solve() override; + bool Solve(SampleGenerator *sample_generator = NULL) override; void SetFlowAtSubdomain(std::function F, const int m=-1); diff --git a/include/linelast_solver.hpp b/include/linelast_solver.hpp index 8b5aa3e5..2ea41596 100644 --- a/include/linelast_solver.hpp +++ b/include/linelast_solver.hpp @@ -81,7 +81,7 @@ class LinElastSolver : public MultiBlockSolver // system-specific. virtual void AssembleInterfaceMatrices(); - virtual bool Solve(); + bool Solve(SampleGenerator *sample_generator = NULL) override; virtual void SetupBCVariables() override; virtual void SetupIC(std::function F); @@ -93,9 +93,9 @@ class LinElastSolver : public MultiBlockSolver virtual void SetupDomainBCOperators(); // Component-wise assembly - virtual void BuildCompROMLinElems(Array &fes_comp); - virtual void BuildBdrROMLinElems(Array &fes_comp); - virtual void BuildItfaceROMLinElems(Array &fes_comp); + void BuildCompROMLinElems() override; + void BuildBdrROMLinElems() override; + void BuildItfaceROMLinElems() override; virtual void ProjectOperatorOnReducedBasis(); diff --git a/include/multiblock_solver.hpp b/include/multiblock_solver.hpp index 2c824582..8e0dfe0b 100644 --- a/include/multiblock_solver.hpp +++ b/include/multiblock_solver.hpp @@ -13,6 +13,7 @@ #include "rom_handler.hpp" #include "hdf5_utils.hpp" #include "sample_generator.hpp" +#include "rom_element_collection.hpp" // By convention we only use mfem namespace as default, not CAROM. using namespace mfem; @@ -112,13 +113,8 @@ friend class ParameterizedProblem; bool separate_variable_basis = false; // Used for bottom-up building, only with ComponentTopologyHandler. - // For now, assumes ROM basis represents the entire vector solution. - Array comp_mats; // Size(num_components); - // boundary condition is enforced via forcing term. - Array *> bdr_mats; - Array port_mats; // reference ports. - // DenseTensor objects from nonlinear operators - // will be defined per each derived MultiBlockSolver. + Array comp_fes; + ROMLinearElement *rom_elems = NULL; public: MultiBlockSolver(); @@ -202,31 +198,23 @@ friend class ParameterizedProblem; // Component-wise assembly void GetComponentFESpaces(Array &comp_fes); - virtual void AllocateROMLinElems(); + // virtual void AllocateROMLinElems(); void BuildROMLinElems(); - virtual void BuildCompROMLinElems(Array &fes_comp) = 0; - virtual void BuildBdrROMLinElems(Array &fes_comp) = 0; + virtual void BuildCompROMLinElems() = 0; + virtual void BuildBdrROMLinElems() = 0; // TODO(kevin): part of this can be transferred to InterfaceForm. - virtual void BuildItfaceROMLinElems(Array &fes_comp) = 0; - - void SaveROMLinElems(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. - virtual void SaveCompBdrROMLinElems(hid_t &file_id); - void SaveBdrROMLinElems(hid_t &comp_grp_id, const int &comp_idx); - void SaveItfaceROMLinElems(hid_t &file_id); - - void LoadROMLinElems(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. - virtual void LoadCompBdrROMLinElems(hid_t &file_id); - void LoadBdrROMLinElems(hid_t &comp_grp_id, const int &comp_idx); - void LoadItfaceROMLinElems(hid_t &file_id); + virtual void BuildItfaceROMLinElems() = 0; + + void SaveROMLinElems(const std::string &filename) + { assert(rom_elems); rom_elems->Save(filename); } + + void LoadROMLinElems(const std::string &filename) + { assert(rom_elems); rom_elems->Load(filename); } void AssembleROMMat(); - virtual bool Solve() = 0; + virtual bool Solve(SampleGenerator *sample_generator = NULL) = 0; virtual void InitVisualization(const std::string& output_dir = ""); virtual void InitUnifiedParaview(const std::string &file_prefix); @@ -271,7 +259,9 @@ friend class ParameterizedProblem; void InitROMHandler(); void GetBasisTags(std::vector &basis_tags); - virtual void PrepareSnapshots(BlockVector* &U_snapshots, std::vector &basis_tags); + virtual BlockVector* PrepareSnapshots(std::vector &basis_tags); + void SaveSnapshots(SampleGenerator *sample_generator); + virtual void LoadReducedBasis() { rom_handler->LoadReducedBasis(); } virtual void ProjectOperatorOnReducedBasis() = 0; virtual void ProjectRHSOnReducedBasis(); diff --git a/include/poisson_solver.hpp b/include/poisson_solver.hpp index 8eeaaf81..9aa72f86 100644 --- a/include/poisson_solver.hpp +++ b/include/poisson_solver.hpp @@ -88,11 +88,11 @@ friend class ParameterizedProblem; virtual void AssembleInterfaceMatrices(); // Component-wise assembly - virtual void BuildCompROMLinElems(Array &fes_comp); - virtual void BuildBdrROMLinElems(Array &fes_comp); - virtual void BuildItfaceROMLinElems(Array &fes_comp); + virtual void BuildCompROMLinElems() override; + void BuildBdrROMLinElems() override; + void BuildItfaceROMLinElems() override; - virtual bool Solve(); + virtual bool Solve(SampleGenerator *sample_generator = NULL); virtual void ProjectOperatorOnReducedBasis(); diff --git a/include/rom_element_collection.hpp b/include/rom_element_collection.hpp new file mode 100644 index 00000000..fd7b04a3 --- /dev/null +++ b/include/rom_element_collection.hpp @@ -0,0 +1,118 @@ +// Copyright 2023 Lawrence Livermore National Security, LLC. See the top-level LICENSE file for details. +// +// SPDX-License-Identifier: MIT + +#ifndef SCALEUPROM_ROM_ELEMENT_COLLECTION_HPP +#define SCALEUPROM_ROM_ELEMENT_COLLECTION_HPP + +#include "topology_handler.hpp" +#include "rom_nonlinearform.hpp" +#include "rom_interfaceform.hpp" +#include "mfem.hpp" +#include "hdf5_utils.hpp" + +// By convention we only use mfem namespace as default, not CAROM. +using namespace mfem; + +class ROMElementCollection +{ +protected: + const int num_var; + const int num_comp; + const int num_ref_ports; + const bool separate_variable; + + TopologyHandler *topol_handler = NULL; // not owned + Array fes; // not owned + +public: + ROMElementCollection(TopologyHandler *topol_handler_, const Array &fes_, + const bool separate_variable_) + : topol_handler(topol_handler_), fes(fes_), + num_comp(topol_handler_->GetNumComponents()), + num_var(fes_.Size() / topol_handler_->GetNumComponents()), + num_ref_ports(topol_handler_->GetNumRefPorts()), + separate_variable(separate_variable_) + { + assert(num_comp * num_var == fes.Size()); + assert(topol_handler->GetType() == TopologyHandlerMode::COMPONENT); + } + + virtual ~ROMElementCollection() {} + + virtual void Save(const std::string &filename) = 0; + virtual void Load(const std::string &filename) = 0; +}; + +class ROMLinearElement : public ROMElementCollection +{ +public: + Array comp; // Size(num_components); + // boundary condition is enforced via forcing term. + Array *> bdr; + Array port; // reference ports. + +public: + ROMLinearElement(TopologyHandler *topol_handler_, + const Array &fes_, + const bool separate_variable_); + + virtual ~ROMLinearElement(); + + void Save(const std::string &filename) override; + void Load(const std::string &filename) override; + +private: + void SaveCompBdrElems(hid_t &file_id); + void SaveBdrElems(hid_t &comp_grp_id, const int &comp_idx); + void SaveItfaceElems(hid_t &file_id); + + void LoadCompBdrElems(hid_t &file_id); + void LoadBdrElems(hid_t &comp_grp_id, const int &comp_idx); + void LoadItfaceElems(hid_t &file_id); +}; + +class ROMTensorElement : public ROMElementCollection +{ +public: + Array comp; // Size(num_components); + + /* boundary/interface is not implemented yet.. should consider */ + // Array *> bdr; + // Array port; // reference ports. + +public: + ROMTensorElement(TopologyHandler *topol_handler_, + const Array &fes_, + const bool separate_variable_) + : ROMElementCollection(topol_handler_, fes_, separate_variable_) {} + + virtual ~ROMTensorElement(); + + void Save(const std::string &filename) override {} + void Load(const std::string &filename) override {} + +}; + +class ROMEQPElement : public ROMElementCollection +{ +public: + Array comp; // Size(num_components); + // boundary condition is enforced via forcing term. + Array *> bdr; + Array port; // reference ports. + +public: + ROMEQPElement(TopologyHandler *topol_handler_, + const Array &fes_, + const bool separate_variable_) + : ROMElementCollection(topol_handler_, fes_, separate_variable_) {} + + virtual ~ROMEQPElement() {} + + void Save(const std::string &filename) override {} + void Load(const std::string &filename) override {} + +}; + +#endif diff --git a/include/sample_generator.hpp b/include/sample_generator.hpp index 639f5c89..3fabdceb 100644 --- a/include/sample_generator.hpp +++ b/include/sample_generator.hpp @@ -64,10 +64,10 @@ class SampleGenerator std::vector basis_tags; std::map basis_tag2idx; - /* snapshot pairs per interface port, for nonlinear EQP */ + /* snapshot pairs per interface port, for nonlinear interface EQP */ std::vector port_tags; std::map port_tag2idx; - Array *> port_colidxs; + Array *> port_colidxs; public: SampleGenerator(MPI_Comm comm); diff --git a/include/steady_ns_solver.hpp b/include/steady_ns_solver.hpp index e1ae4dcb..dc3374a2 100644 --- a/include/steady_ns_solver.hpp +++ b/include/steady_ns_solver.hpp @@ -137,7 +137,6 @@ friend class SteadyNSOperator; // component ROM element for nonlinear convection. Array comp_tensors, subdomain_tensors; Array comp_eqps, subdomain_eqps; - Array comp_fes; // pointers to existing fespace, no need to delete Solver *J_solver = NULL; GMRESSolver *J_gmres = NULL; @@ -161,7 +160,7 @@ friend class SteadyNSOperator; void SaveROMOperator(const std::string input_prefix="") override; void LoadROMOperatorFromFile(const std::string input_prefix="") override; - bool Solve() override; + bool Solve(SampleGenerator *sample_generator = NULL) override; void ProjectOperatorOnReducedBasis() override; diff --git a/include/stokes_solver.hpp b/include/stokes_solver.hpp index 73244397..e1171956 100644 --- a/include/stokes_solver.hpp +++ b/include/stokes_solver.hpp @@ -184,11 +184,11 @@ friend class ParameterizedProblem; virtual void SetupPressureMassMatrix(); // Component-wise assembly - virtual void BuildCompROMLinElems(Array &fes_comp); - virtual void BuildBdrROMLinElems(Array &fes_comp); - virtual void BuildItfaceROMLinElems(Array &fes_comp); + void BuildCompROMLinElems() override; + void BuildBdrROMLinElems() override; + void BuildItfaceROMLinElems() override; - virtual bool Solve(); + virtual bool Solve(SampleGenerator *sample_generator = NULL); virtual void Solve_obsolete(); virtual void LoadReducedBasis() override; diff --git a/include/unsteady_ns_solver.hpp b/include/unsteady_ns_solver.hpp index 75571b7f..ddfde72d 100644 --- a/include/unsteady_ns_solver.hpp +++ b/include/unsteady_ns_solver.hpp @@ -75,7 +75,7 @@ friend class SteadyNSOperator; void LoadROMOperatorFromFile(const std::string input_prefix="") override { mfem_error("UnsteadyNSSolver::LoadROMOperatorFromFile is not implemented yet!\n"); } - bool Solve() override; + bool Solve(SampleGenerator *sample_generator = NULL) override; using MultiBlockSolver::SaveVisualization; void SaveVisualization(const int step, const double time) override; diff --git a/src/advdiff_solver.cpp b/src/advdiff_solver.cpp index 7823c080..94b486f9 100644 --- a/src/advdiff_solver.cpp +++ b/src/advdiff_solver.cpp @@ -47,20 +47,18 @@ void AdvDiffSolver::BuildDomainOperators() } } -void AdvDiffSolver::BuildCompROMLinElems(Array &fes_comp) +void AdvDiffSolver::BuildCompROMLinElems() { mfem_error("AdvDiffSolver::BuildCompROMLinElems is not implemented yet!\n"); assert(train_mode == UNIVERSAL); assert(rom_handler->BasisLoaded()); + assert(rom_elems); - const int num_comp = fes_comp.Size(); - assert(comp_mats.Size() == num_comp); - - for (int c = 0; c < num_comp; c++) + for (int c = 0; c < topol_handler->GetNumComponents(); c++) { Mesh *comp = topol_handler->GetComponentMesh(c); - BilinearForm a_comp(fes_comp[c]); + BilinearForm a_comp(comp_fes[c]); a_comp.AddDomainIntegrator(new DiffusionIntegrator); if (full_dg) @@ -72,12 +70,12 @@ void AdvDiffSolver::BuildCompROMLinElems(Array &fes_comp) a_comp.Finalize(); // Poisson equation has only one solution variable. - comp_mats[c]->SetSize(1, 1); - (*comp_mats[c])(0, 0) = rom_handler->ProjectToRefBasis(c, c, &(a_comp.SpMat())); + rom_elems->comp[c]->SetSize(1, 1); + (*rom_elems->comp[c])(0, 0) = rom_handler->ProjectToRefBasis(c, c, &(a_comp.SpMat())); } } -bool AdvDiffSolver::Solve() +bool AdvDiffSolver::Solve(SampleGenerator *sample_generator) { // If using direct solver, returns always true. bool converged = true; @@ -151,6 +149,10 @@ bool AdvDiffSolver::Solve() delete solver; } + /* save solution if sample generator is provided */ + if (converged && sample_generator) + SaveSnapshots(sample_generator); + return converged; } diff --git a/src/linelast_solver.cpp b/src/linelast_solver.cpp index 63b462bc..a36b6193 100644 --- a/src/linelast_solver.cpp +++ b/src/linelast_solver.cpp @@ -305,7 +305,7 @@ void LinElastSolver::AssembleInterfaceMatrices() a_itf->AssembleInterfaceMatrices(mats); } -bool LinElastSolver::Solve() +bool LinElastSolver::Solve(SampleGenerator *sample_generator) { // If using direct solver, returns always true. bool converged = true; @@ -381,6 +381,10 @@ bool LinElastSolver::Solve() delete solver; } + /* save solution if sample generator is provided */ + if (converged && sample_generator) + SaveSnapshots(sample_generator); + return converged; } @@ -494,18 +498,16 @@ void LinElastSolver::ProjectOperatorOnReducedBasis() } // Component-wise assembly -void LinElastSolver::BuildCompROMLinElems(Array &fes_comp) +void LinElastSolver::BuildCompROMLinElems() { assert(train_mode == UNIVERSAL); assert(rom_handler->BasisLoaded()); + assert(rom_elems); - const int num_comp = fes_comp.Size(); - assert(comp_mats.Size() == num_comp); - - for (int c = 0; c < num_comp; c++) + for (int c = 0; c < topol_handler->GetNumComponents(); c++) { Mesh *comp = topol_handler->GetComponentMesh(c); - BilinearForm a_comp(fes_comp[c]); + BilinearForm a_comp(comp_fes[c]); a_comp.AddDomainIntegrator(new ElasticityIntegrator(*(lambda_c[c]), *(mu_c[c]))); if (full_dg) @@ -515,23 +517,21 @@ void LinElastSolver::BuildCompROMLinElems(Array &fes_comp) a_comp.Finalize(); // Elasticity equation has only one solution variable. - comp_mats[c]->SetSize(1, 1); - (*comp_mats[c])(0, 0) = rom_handler->ProjectToRefBasis(c, c, &(a_comp.SpMat())); + rom_elems->comp[c]->SetSize(1, 1); + (*rom_elems->comp[c])(0, 0) = rom_handler->ProjectToRefBasis(c, c, &(a_comp.SpMat())); } } -void LinElastSolver::BuildBdrROMLinElems(Array &fes_comp) +void LinElastSolver::BuildBdrROMLinElems() { assert(train_mode == UNIVERSAL); assert(rom_handler->BasisLoaded()); + assert(rom_elems); - const int num_comp = fes_comp.Size(); - assert(bdr_mats.Size() == num_comp); - - for (int c = 0; c < num_comp; c++) + for (int c = 0; c < topol_handler->GetNumComponents(); c++) { Mesh *comp = topol_handler->GetComponentMesh(c); - assert(bdr_mats[c]->Size() == comp->bdr_attributes.Size()); + assert(rom_elems->bdr[c]->Size() == comp->bdr_attributes.Size()); MatrixBlocks *bdr_mat; for (int b = 0; b < comp->bdr_attributes.Size(); b++) @@ -539,31 +539,31 @@ void LinElastSolver::BuildBdrROMLinElems(Array &fes_comp) Array bdr_marker(comp->bdr_attributes.Max()); bdr_marker = 0; bdr_marker[comp->bdr_attributes[b] - 1] = 1; - BilinearForm a_comp(fes_comp[c]); + BilinearForm a_comp(comp_fes[c]); a_comp.AddBdrFaceIntegrator(new DGElasticityIntegrator(*(lambda_c[c]), *(mu_c[c]), alpha, kappa), bdr_marker); a_comp.Assemble(); a_comp.Finalize(); - bdr_mat = (*bdr_mats[c])[b]; + bdr_mat = (*rom_elems->bdr[c])[b]; bdr_mat->SetSize(1, 1); (*bdr_mat)(0, 0) = rom_handler->ProjectToRefBasis(c, c, &(a_comp.SpMat())); } } } -void LinElastSolver::BuildItfaceROMLinElems(Array &fes_comp) +void LinElastSolver::BuildItfaceROMLinElems() { assert(topol_mode == TopologyHandlerMode::COMPONENT); assert(train_mode == UNIVERSAL); assert(rom_handler->BasisLoaded()); + assert(rom_elems); const int num_ref_ports = topol_handler->GetNumRefPorts(); - assert(port_mats.Size() == num_ref_ports); for (int p = 0; p < num_ref_ports; p++) { - assert(port_mats[p]->nrows == 2); - assert(port_mats[p]->ncols == 2); + assert(rom_elems->port[p]->nrows == 2); + assert(rom_elems->port[p]->ncols == 2); int c1, c2; topol_handler->GetComponentPair(p, c1, c2); @@ -577,11 +577,11 @@ void LinElastSolver::BuildItfaceROMLinElems(Array &fes_com // NOTE: If comp1 == comp2, using comp1 and comp2 directly leads to an incorrect penalty matrix. // Need to use two copied instances. - a_itf->AssembleInterfaceMatrixAtPort(p, fes_comp, spmats); + a_itf->AssembleInterfaceMatrixAtPort(p, comp_fes, spmats); for (int i = 0; i < 2; i++) for (int j = 0; j < 2; j++) - (*port_mats[p])(i, j) = rom_handler->ProjectToRefBasis(c_idx[i], c_idx[j], spmats(i,j)); + (*rom_elems->port[p])(i, j) = rom_handler->ProjectToRefBasis(c_idx[i], c_idx[j], spmats(i,j)); for (int i = 0; i < 2; i++) for (int j = 0; j < 2; j++) delete spmats(i, j); diff --git a/src/main_workflow.cpp b/src/main_workflow.cpp index cb475b2c..e2c5996f 100644 --- a/src/main_workflow.cpp +++ b/src/main_workflow.cpp @@ -185,7 +185,7 @@ void GenerateSamples(MPI_Comm comm) test->SetupBCOperators(); test->Assemble(); - bool converged = test->Solve(); + bool converged = test->Solve(sample_generator); if (!converged) { // If deterministic, terminate the sampling here. @@ -203,19 +203,8 @@ void GenerateSamples(MPI_Comm comm) test->SaveSolution(sol_file); test->SaveVisualization(); - // test->SaveSnapshot(file_idx); - // View-Vector of the entire solution of test, splitted according to ROM basis setup. - BlockVector *U_snapshots = NULL; - // Basis tags for each block of U_snapshots. - std::vector basis_tags; - Array col_idxs; - test->PrepareSnapshots(U_snapshots, basis_tags); - sample_generator->SaveSnapshot(U_snapshots, basis_tags, col_idxs); - sample_generator->SaveSnapshotPorts(test->GetTopologyHandler(), test->GetTrainMode(), col_idxs); - sample_generator->ReportStatus(s); - delete U_snapshots; delete test; s++; @@ -434,7 +423,6 @@ void BuildROM(MPI_Comm comm) if (topol_mode == TopologyHandlerMode::SUBMESH) mfem_error("Submesh does not support component rom building level!\n"); - test->AllocateROMLinElems(); test->BuildROMLinElems(); test->SaveROMLinElems(filename); @@ -525,8 +513,6 @@ double SingleRun(MPI_Comm comm, const std::string output_file) { if (topol_mode == TopologyHandlerMode::SUBMESH) mfem_error("Submesh does not support component rom building level!\n"); - - test->AllocateROMLinElems(); printf("Loading ROM projected elements.. "); test->LoadROMLinElems(filename); diff --git a/src/multiblock_solver.cpp b/src/multiblock_solver.cpp index b3c0668d..8eff3288 100644 --- a/src/multiblock_solver.cpp +++ b/src/multiblock_solver.cpp @@ -72,16 +72,8 @@ MultiBlockSolver::~MultiBlockSolver() delete rom_handler; delete topol_handler; - DeletePointers(comp_mats); - - for (int c = 0; c < bdr_mats.Size(); c++) - { - DeletePointers((*bdr_mats[c])); - - delete bdr_mats[c]; - } - - DeletePointers(port_mats); + DeletePointers(comp_fes); + delete rom_elems; } void MultiBlockSolver::ParseInputs() @@ -266,277 +258,26 @@ void MultiBlockSolver::GetComponentFESpaces(Array &comp_fe } } -void MultiBlockSolver::AllocateROMLinElems() -{ - assert(topol_mode == TopologyHandlerMode::COMPONENT); - assert(train_mode == UNIVERSAL); - - const int num_comp = topol_handler->GetNumComponents(); - const int num_ref_ports = topol_handler->GetNumRefPorts(); - - comp_mats.SetSize(num_comp); - const int block_size = (separate_variable_basis) ? num_var : 1; - for (int c = 0; c < num_comp; c++) - comp_mats[c] = new MatrixBlocks(block_size, block_size); - - bdr_mats.SetSize(num_comp); - for (int c = 0; c < num_comp; c++) - { - Mesh *comp = topol_handler->GetComponentMesh(c); - bdr_mats[c] = new Array(comp->bdr_attributes.Size()); - for (int b = 0; b < bdr_mats[c]->Size(); b++) - (*bdr_mats[c])[b] = new MatrixBlocks(block_size, block_size); - } - port_mats.SetSize(num_ref_ports); - for (int p = 0; p < num_ref_ports; p++) - port_mats[p] = new MatrixBlocks(2 * block_size, 2 * block_size); -} - void MultiBlockSolver::BuildROMLinElems() { assert(topol_mode == TopologyHandlerMode::COMPONENT); assert(train_mode == UNIVERSAL); assert(rom_handler->BasisLoaded()); - // Component domain system - const int num_comp = topol_handler->GetNumComponents(); - Array fes_comp; - GetComponentFESpaces(fes_comp); - - BuildCompROMLinElems(fes_comp); + BuildCompROMLinElems(); // Boundary penalty matrices - BuildBdrROMLinElems(fes_comp); + BuildBdrROMLinElems(); // Port penalty matrices - BuildItfaceROMLinElems(fes_comp); - - for (int k = 0 ; k < fes_comp.Size(); k++) delete fes_comp[k]; -} - -void MultiBlockSolver::SaveROMLinElems(const std::string &filename) -{ - assert(topol_mode == TopologyHandlerMode::COMPONENT); - assert(train_mode == UNIVERSAL); - - hid_t file_id; - herr_t errf = 0; - file_id = H5Fcreate(filename.c_str(), H5F_ACC_TRUNC, H5P_DEFAULT, H5P_DEFAULT); - assert(file_id >= 0); - - // components + boundary - SaveCompBdrROMLinElems(file_id); - - // (reference) ports - SaveItfaceROMLinElems(file_id); - - errf = H5Fclose(file_id); - assert(errf >= 0); - return; -} - -void MultiBlockSolver::SaveCompBdrROMLinElems(hid_t &file_id) -{ - assert(file_id >= 0); - hid_t grp_id; - herr_t errf; - - grp_id = H5Gcreate(file_id, "components", H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT); - assert(grp_id >= 0); - - const int num_comp = topol_handler->GetNumComponents(); - assert(comp_mats.Size() == num_comp); - assert(bdr_mats.Size() == num_comp); - - hdf5_utils::WriteAttribute(grp_id, "number_of_components", num_comp); - - std::string dset_name; - for (int c = 0; c < num_comp; c++) - { - dset_name = topol_handler->GetComponentName(c); - - hid_t comp_grp_id; - comp_grp_id = H5Gcreate(grp_id, dset_name.c_str(), H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT); - assert(comp_grp_id >= 0); - - hdf5_utils::WriteDataset(comp_grp_id, "domain", *comp_mats[c]); - - // boundaries are saved for each component group. - SaveBdrROMLinElems(comp_grp_id, c); - - errf = H5Gclose(comp_grp_id); - assert(errf >= 0); - } // for (int c = 0; c < num_comp; c++) - - errf = H5Gclose(grp_id); - assert(errf >= 0); - return; -} - -void MultiBlockSolver::SaveBdrROMLinElems(hid_t &comp_grp_id, const int &comp_idx) -{ - assert(comp_grp_id >= 0); - herr_t errf; - hid_t bdr_grp_id; - bdr_grp_id = H5Gcreate(comp_grp_id, "boundary", H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT); - assert(bdr_grp_id >= 0); - - const int num_bdr = bdr_mats[comp_idx]->Size(); - Mesh *comp = topol_handler->GetComponentMesh(comp_idx); - assert(num_bdr == comp->bdr_attributes.Size()); - - hdf5_utils::WriteAttribute(bdr_grp_id, "number_of_boundaries", num_bdr); - - Array *bdr_mat_c = bdr_mats[comp_idx]; - for (int b = 0; b < num_bdr; b++) - hdf5_utils::WriteDataset(bdr_grp_id, std::to_string(b), *((*bdr_mat_c)[b])); - - errf = H5Gclose(bdr_grp_id); - assert(errf >= 0); - return; -} - -void MultiBlockSolver::SaveItfaceROMLinElems(hid_t &file_id) -{ - assert(file_id >= 0); - herr_t errf; - hid_t grp_id; - grp_id = H5Gcreate(file_id, "ports", H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT); - assert(grp_id >= 0); - - const int num_ref_ports = topol_handler->GetNumRefPorts(); - assert(port_mats.Size() == num_ref_ports); - - hdf5_utils::WriteAttribute(grp_id, "number_of_ports", num_ref_ports); - - std::string dset_name; - int c1, c2, a1, a2; - for (int p = 0; p < num_ref_ports; p++) - { - topol_handler->GetRefPortInfo(p, c1, c2, a1, a2); - dset_name = topol_handler->GetComponentName(c1) + ":" + topol_handler->GetComponentName(c2); - dset_name += "-" + std::to_string(a1) + ":" + std::to_string(a2); - hdf5_utils::WriteDataset(grp_id, dset_name, *port_mats[p]); - } - - errf = H5Gclose(grp_id); - assert(errf >= 0); - return; -} - -void MultiBlockSolver::LoadROMLinElems(const std::string &filename) -{ - assert(topol_mode == TopologyHandlerMode::COMPONENT); - assert(train_mode == UNIVERSAL); - - hid_t file_id; - herr_t errf = 0; - file_id = H5Fopen(filename.c_str(), H5F_ACC_RDONLY, H5P_DEFAULT); - assert(file_id >= 0); - - // components - LoadCompBdrROMLinElems(file_id); - - // (reference) ports - LoadItfaceROMLinElems(file_id); - - errf = H5Fclose(file_id); - assert(errf >= 0); - - return; -} - -void MultiBlockSolver::LoadCompBdrROMLinElems(hid_t &file_id) -{ - assert(file_id >= 0); - herr_t errf; - hid_t grp_id; - grp_id = H5Gopen2(file_id, "components", H5P_DEFAULT); - assert(grp_id >= 0); - - int num_comp; - hdf5_utils::ReadAttribute(grp_id, "number_of_components", num_comp); - assert(num_comp >= topol_handler->GetNumComponents()); - assert(comp_mats.Size() == topol_handler->GetNumComponents()); - assert(bdr_mats.Size() == topol_handler->GetNumComponents()); - - std::string dset_name; - for (int c = 0; c < topol_handler->GetNumComponents(); c++) - { - dset_name = topol_handler->GetComponentName(c); - - hid_t comp_grp_id; - comp_grp_id = H5Gopen2(grp_id, dset_name.c_str(), H5P_DEFAULT); - assert(comp_grp_id >= 0); - - hdf5_utils::ReadDataset(comp_grp_id, "domain", *comp_mats[c]); - - // boundary - LoadBdrROMLinElems(comp_grp_id, c); - - errf = H5Gclose(comp_grp_id); - assert(errf >= 0); - } // for (int c = 0; c < num_comp; c++) - - errf = H5Gclose(grp_id); - assert(errf >= 0); -} - -void MultiBlockSolver::LoadBdrROMLinElems(hid_t &comp_grp_id, const int &comp_idx) -{ - assert(comp_grp_id >= 0); - herr_t errf; - hid_t bdr_grp_id; - bdr_grp_id = H5Gopen2(comp_grp_id, "boundary", H5P_DEFAULT); - assert(bdr_grp_id >= 0); - - int num_bdr; - hdf5_utils::ReadAttribute(bdr_grp_id, "number_of_boundaries", num_bdr); - - Mesh *comp = topol_handler->GetComponentMesh(comp_idx); - assert(num_bdr == comp->bdr_attributes.Size()); - assert(num_bdr = bdr_mats[comp_idx]->Size()); - - Array *bdr_mat_c = bdr_mats[comp_idx]; - for (int b = 0; b < num_bdr; b++) - hdf5_utils::ReadDataset(bdr_grp_id, std::to_string(b), *(*bdr_mat_c)[b]); - - errf = H5Gclose(bdr_grp_id); - assert(errf >= 0); - return; -} - -void MultiBlockSolver::LoadItfaceROMLinElems(hid_t &file_id) -{ - assert(file_id >= 0); - herr_t errf; - hid_t grp_id; - grp_id = H5Gopen2(file_id, "ports", H5P_DEFAULT); - assert(grp_id >= 0); - - int num_ref_ports; - hdf5_utils::ReadAttribute(grp_id, "number_of_ports", num_ref_ports); - assert(num_ref_ports >= topol_handler->GetNumRefPorts()); - assert(port_mats.Size() == topol_handler->GetNumRefPorts()); - - std::string dset_name; - int c1, c2, a1, a2; - for (int p = 0; p < topol_handler->GetNumRefPorts(); p++) - { - topol_handler->GetRefPortInfo(p, c1, c2, a1, a2); - dset_name = topol_handler->GetComponentName(c1) + ":" + topol_handler->GetComponentName(c2); - dset_name += "-" + std::to_string(a1) + ":" + std::to_string(a2); - hdf5_utils::ReadDataset(grp_id, dset_name, *port_mats[p]); - } - - errf = H5Gclose(grp_id); - assert(errf >= 0); + BuildItfaceROMLinElems(); } void MultiBlockSolver::AssembleROMMat() { assert(topol_mode == TopologyHandlerMode::COMPONENT); assert(train_mode == UNIVERSAL); + assert(rom_elems); const Array *rom_block_offsets = rom_handler->GetBlockOffsets(); BlockMatrix *romMat = new BlockMatrix(*rom_block_offsets); @@ -558,7 +299,7 @@ void MultiBlockSolver::AssembleROMMat() for (int k = 0; k < num_block; k++) midx[k] = midx0 + k; - MatrixBlocks *comp_mat = comp_mats[c_type]; + MatrixBlocks *comp_mat = rom_elems->comp[c_type]; AddToBlockMatrix(midx, midx, *comp_mat, *romMat); // boundary matrices of each component. @@ -574,7 +315,7 @@ void MultiBlockSolver::AssembleROMMat() if (bdr_type[global_idx] == BoundaryType::NEUMANN) continue; - MatrixBlocks *bdr_mat = (*bdr_mats[c_type])[b]; + MatrixBlocks *bdr_mat = (*(rom_elems->bdr[c_type]))[b]; AddToBlockMatrix(midx, midx, *bdr_mat, *romMat); } // for (int b = 0; b < bdr_c2g->Size(); b++) } // for (int m = 0; m < numSub; m++) @@ -584,7 +325,7 @@ void MultiBlockSolver::AssembleROMMat() { const PortInfo *pInfo = topol_handler->GetPortInfo(p); const int p_type = topol_handler->GetPortType(p); - MatrixBlocks *port_mat = port_mats[p_type]; + MatrixBlocks *port_mat = rom_elems->port[p_type]; const int m1 = pInfo->Mesh1; const int m2 = pInfo->Mesh2; @@ -918,6 +659,12 @@ void MultiBlockSolver::AssembleROMNlinOper() void MultiBlockSolver::InitROMHandler() { rom_handler = new MFEMROMHandler(train_mode, topol_handler, var_offsets, var_names, separate_variable_basis); + + if (!((topol_mode == TopologyHandlerMode::COMPONENT) && (train_mode == UNIVERSAL))) + return; + + GetComponentFESpaces(comp_fes); + rom_elems = new ROMLinearElement(topol_handler, comp_fes, separate_variable_basis); } void MultiBlockSolver::GetBasisTags(std::vector &basis_tags) @@ -937,8 +684,10 @@ void MultiBlockSolver::GetBasisTags(std::vector &basis_tags) } } -void MultiBlockSolver::PrepareSnapshots(BlockVector* &U_snapshots, std::vector &basis_tags) +BlockVector* MultiBlockSolver::PrepareSnapshots(std::vector &basis_tags) { + BlockVector *U_snapshots = NULL; + // View vector for U. if (separate_variable_basis) U_snapshots = new BlockVector(U->GetData(), var_offsets); @@ -947,6 +696,25 @@ void MultiBlockSolver::PrepareSnapshots(BlockVector* &U_snapshots, std::vectorNumBlocks() == basis_tags.size()); + + return U_snapshots; +} + +void MultiBlockSolver::SaveSnapshots(SampleGenerator *sample_generator) +{ + assert(sample_generator); + + /* split the solution into each component with the corresponding tag */ + std::vector basis_tags; + BlockVector *U_snapshots = PrepareSnapshots(basis_tags); + + Array col_idxs; + sample_generator->SaveSnapshot(U_snapshots, basis_tags, col_idxs); + sample_generator->SaveSnapshotPorts(topol_handler, train_mode, col_idxs); + + /* delete only the view vector, not the data itself. */ + delete U_snapshots; + return; } void MultiBlockSolver::ProjectRHSOnReducedBasis() diff --git a/src/poisson_solver.cpp b/src/poisson_solver.cpp index 2d0e9090..b4f23a5c 100644 --- a/src/poisson_solver.cpp +++ b/src/poisson_solver.cpp @@ -342,18 +342,16 @@ void PoissonSolver::AssembleInterfaceMatrices() a_itf->AssembleInterfaceMatrices(mats); } -void PoissonSolver::BuildCompROMLinElems(Array &fes_comp) +void PoissonSolver::BuildCompROMLinElems() { assert(train_mode == UNIVERSAL); assert(rom_handler->BasisLoaded()); + assert(rom_elems); - const int num_comp = fes_comp.Size(); - assert(comp_mats.Size() == num_comp); - - for (int c = 0; c < num_comp; c++) + for (int c = 0; c < topol_handler->GetNumComponents(); c++) { Mesh *comp = topol_handler->GetComponentMesh(c); - BilinearForm a_comp(fes_comp[c]); + BilinearForm a_comp(comp_fes[c]); a_comp.AddDomainIntegrator(new DiffusionIntegrator); if (full_dg) @@ -363,23 +361,21 @@ void PoissonSolver::BuildCompROMLinElems(Array &fes_comp) a_comp.Finalize(); // Poisson equation has only one solution variable. - comp_mats[c]->SetSize(1, 1); - (*comp_mats[c])(0, 0) = rom_handler->ProjectToRefBasis(c, c, &(a_comp.SpMat())); + rom_elems->comp[c]->SetSize(1, 1); + (*rom_elems->comp[c])(0, 0) = rom_handler->ProjectToRefBasis(c, c, &(a_comp.SpMat())); } } -void PoissonSolver::BuildBdrROMLinElems(Array &fes_comp) +void PoissonSolver::BuildBdrROMLinElems() { assert(train_mode == UNIVERSAL); assert(rom_handler->BasisLoaded()); + assert(rom_elems); - const int num_comp = fes_comp.Size(); - assert(bdr_mats.Size() == num_comp); - - for (int c = 0; c < num_comp; c++) + for (int c = 0; c < topol_handler->GetNumComponents(); c++) { Mesh *comp = topol_handler->GetComponentMesh(c); - assert(bdr_mats[c]->Size() == comp->bdr_attributes.Size()); + assert(rom_elems->bdr[c]->Size() == comp->bdr_attributes.Size()); MatrixBlocks *bdr_mat; for (int b = 0; b < comp->bdr_attributes.Size(); b++) @@ -387,31 +383,31 @@ void PoissonSolver::BuildBdrROMLinElems(Array &fes_comp) Array bdr_marker(comp->bdr_attributes.Max()); bdr_marker = 0; bdr_marker[comp->bdr_attributes[b] - 1] = 1; - BilinearForm a_comp(fes_comp[c]); + BilinearForm a_comp(comp_fes[c]); a_comp.AddBdrFaceIntegrator(new DGDiffusionIntegrator(sigma, kappa), bdr_marker); a_comp.Assemble(); a_comp.Finalize(); - bdr_mat = (*bdr_mats[c])[b]; + bdr_mat = (*rom_elems->bdr[c])[b]; bdr_mat->SetSize(1, 1); (*bdr_mat)(0, 0) = rom_handler->ProjectToRefBasis(c, c, &(a_comp.SpMat())); } } } -void PoissonSolver::BuildItfaceROMLinElems(Array &fes_comp) +void PoissonSolver::BuildItfaceROMLinElems() { assert(topol_mode == TopologyHandlerMode::COMPONENT); assert(train_mode == UNIVERSAL); assert(rom_handler->BasisLoaded()); + assert(rom_elems); const int num_ref_ports = topol_handler->GetNumRefPorts(); - assert(port_mats.Size() == num_ref_ports); for (int p = 0; p < num_ref_ports; p++) { - assert(port_mats[p]->nrows == 2); - assert(port_mats[p]->ncols == 2); + assert(rom_elems->port[p]->nrows == 2); + assert(rom_elems->port[p]->ncols == 2); int c1, c2; topol_handler->GetComponentPair(p, c1, c2); @@ -425,18 +421,18 @@ void PoissonSolver::BuildItfaceROMLinElems(Array &fes_comp // NOTE: If comp1 == comp2, using comp1 and comp2 directly leads to an incorrect penalty matrix. // Need to use two copied instances. - a_itf->AssembleInterfaceMatrixAtPort(p, fes_comp, spmats); + a_itf->AssembleInterfaceMatrixAtPort(p, comp_fes, spmats); for (int i = 0; i < 2; i++) for (int j = 0; j < 2; j++) - (*port_mats[p])(i, j) = rom_handler->ProjectToRefBasis(c_idx[i], c_idx[j], spmats(i,j)); + (*rom_elems->port[p])(i, j) = rom_handler->ProjectToRefBasis(c_idx[i], c_idx[j], spmats(i,j)); for (int i = 0; i < 2; i++) for (int j = 0; j < 2; j++) delete spmats(i, j); } // for (int p = 0; p < num_ref_ports; p++) } -bool PoissonSolver::Solve() +bool PoissonSolver::Solve(SampleGenerator *sample_generator) { // If using direct solver, returns always true. bool converged = true; @@ -510,6 +506,10 @@ bool PoissonSolver::Solve() delete solver; } + /* save solution if sample generator is provided */ + if (converged && sample_generator) + SaveSnapshots(sample_generator); + return converged; } diff --git a/src/rom_element_collection.cpp b/src/rom_element_collection.cpp new file mode 100644 index 00000000..5c2e4fb6 --- /dev/null +++ b/src/rom_element_collection.cpp @@ -0,0 +1,244 @@ +// Copyright 2023 Lawrence Livermore National Security, LLC. See the top-level LICENSE file for details. +// +// SPDX-License-Identifier: MIT + +#include "rom_element_collection.hpp" +#include "etc.hpp" + +using namespace mfem; + +ROMLinearElement::ROMLinearElement( + TopologyHandler *topol_handler_, const Array &fes_, const bool separate_variable_) + : ROMElementCollection(topol_handler_, fes_, separate_variable_) +{ + comp.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); + + bdr.SetSize(num_comp); + for (int c = 0; c < num_comp; c++) + { + Mesh *comp = topol_handler->GetComponentMesh(c); + bdr[c] = new Array(comp->bdr_attributes.Size()); + for (int b = 0; b < bdr[c]->Size(); b++) + (*bdr[c])[b] = new MatrixBlocks(block_size, block_size); + } + port.SetSize(num_ref_ports); + for (int p = 0; p < num_ref_ports; p++) + port[p] = new MatrixBlocks(2 * block_size, 2 * block_size); +} + +ROMLinearElement::~ROMLinearElement() +{ + DeletePointers(comp); + for (int c = 0; c < bdr.Size(); c++) + { + DeletePointers((*bdr[c])); + delete bdr[c]; + } + DeletePointers(port); +} + +void ROMLinearElement::Save(const std::string &filename) +{ + hid_t file_id; + herr_t errf = 0; + file_id = H5Fcreate(filename.c_str(), H5F_ACC_TRUNC, H5P_DEFAULT, H5P_DEFAULT); + assert(file_id >= 0); + + // components + boundary + SaveCompBdrElems(file_id); + + // (reference) ports + SaveItfaceElems(file_id); + + errf = H5Fclose(file_id); + assert(errf >= 0); + return; +} + +void ROMLinearElement::SaveCompBdrElems(hid_t &file_id) +{ + assert(file_id >= 0); + hid_t grp_id; + herr_t errf; + + grp_id = H5Gcreate(file_id, "components", H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT); + assert(grp_id >= 0); + + hdf5_utils::WriteAttribute(grp_id, "number_of_components", num_comp); + + std::string dset_name; + for (int c = 0; c < num_comp; c++) + { + dset_name = topol_handler->GetComponentName(c); + + hid_t comp_grp_id; + comp_grp_id = H5Gcreate(grp_id, dset_name.c_str(), H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT); + assert(comp_grp_id >= 0); + + hdf5_utils::WriteDataset(comp_grp_id, "domain", *comp[c]); + + // boundaries are saved for each component group. + SaveBdrElems(comp_grp_id, c); + + errf = H5Gclose(comp_grp_id); + assert(errf >= 0); + } // for (int c = 0; c < num_comp; c++) + + errf = H5Gclose(grp_id); + assert(errf >= 0); + return; +} + +void ROMLinearElement::SaveBdrElems(hid_t &comp_grp_id, const int &comp_idx) +{ + assert(comp_grp_id >= 0); + herr_t errf; + hid_t bdr_grp_id; + bdr_grp_id = H5Gcreate(comp_grp_id, "boundary", H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT); + assert(bdr_grp_id >= 0); + + const int num_bdr = bdr[comp_idx]->Size(); + Mesh *comp = topol_handler->GetComponentMesh(comp_idx); + assert(num_bdr == comp->bdr_attributes.Size()); + + hdf5_utils::WriteAttribute(bdr_grp_id, "number_of_boundaries", num_bdr); + + Array *bdr_c = bdr[comp_idx]; + for (int b = 0; b < num_bdr; b++) + hdf5_utils::WriteDataset(bdr_grp_id, std::to_string(b), *((*bdr_c)[b])); + + errf = H5Gclose(bdr_grp_id); + assert(errf >= 0); + return; +} + +void ROMLinearElement::SaveItfaceElems(hid_t &file_id) +{ + assert(file_id >= 0); + herr_t errf; + hid_t grp_id; + grp_id = H5Gcreate(file_id, "ports", H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT); + assert(grp_id >= 0); + + hdf5_utils::WriteAttribute(grp_id, "number_of_ports", num_ref_ports); + + std::string dset_name; + int c1, c2, a1, a2; + for (int p = 0; p < num_ref_ports; p++) + { + topol_handler->GetRefPortInfo(p, c1, c2, a1, a2); + dset_name = topol_handler->GetComponentName(c1) + ":" + topol_handler->GetComponentName(c2); + dset_name += "-" + std::to_string(a1) + ":" + std::to_string(a2); + hdf5_utils::WriteDataset(grp_id, dset_name, *port[p]); + } + + errf = H5Gclose(grp_id); + assert(errf >= 0); + return; +} + +void ROMLinearElement::Load(const std::string &filename) +{ + hid_t file_id; + herr_t errf = 0; + file_id = H5Fopen(filename.c_str(), H5F_ACC_RDONLY, H5P_DEFAULT); + assert(file_id >= 0); + + // components + LoadCompBdrElems(file_id); + + // (reference) ports + LoadItfaceElems(file_id); + + errf = H5Fclose(file_id); + assert(errf >= 0); + + return; +} + +void ROMLinearElement::LoadCompBdrElems(hid_t &file_id) +{ + assert(file_id >= 0); + herr_t errf; + hid_t grp_id; + grp_id = H5Gopen2(file_id, "components", H5P_DEFAULT); + assert(grp_id >= 0); + + int num_comp_; + hdf5_utils::ReadAttribute(grp_id, "number_of_components", num_comp_); + assert(num_comp_ >= num_comp); + + std::string dset_name; + for (int c = 0; c < num_comp; c++) + { + dset_name = topol_handler->GetComponentName(c); + + hid_t comp_grp_id; + comp_grp_id = H5Gopen2(grp_id, dset_name.c_str(), H5P_DEFAULT); + assert(comp_grp_id >= 0); + + hdf5_utils::ReadDataset(comp_grp_id, "domain", *comp[c]); + + // boundary + LoadBdrElems(comp_grp_id, c); + + errf = H5Gclose(comp_grp_id); + assert(errf >= 0); + } // for (int c = 0; c < num_comp; c++) + + errf = H5Gclose(grp_id); + assert(errf >= 0); +} + +void ROMLinearElement::LoadBdrElems(hid_t &comp_grp_id, const int &comp_idx) +{ + assert(comp_grp_id >= 0); + herr_t errf; + hid_t bdr_grp_id; + bdr_grp_id = H5Gopen2(comp_grp_id, "boundary", H5P_DEFAULT); + assert(bdr_grp_id >= 0); + + int num_bdr; + hdf5_utils::ReadAttribute(bdr_grp_id, "number_of_boundaries", num_bdr); + + Mesh *comp = topol_handler->GetComponentMesh(comp_idx); + assert(num_bdr == comp->bdr_attributes.Size()); + assert(num_bdr = bdr[comp_idx]->Size()); + + Array *bdr_c = bdr[comp_idx]; + for (int b = 0; b < num_bdr; b++) + hdf5_utils::ReadDataset(bdr_grp_id, std::to_string(b), *(*bdr_c)[b]); + + errf = H5Gclose(bdr_grp_id); + assert(errf >= 0); + return; +} + +void ROMLinearElement::LoadItfaceElems(hid_t &file_id) +{ + assert(file_id >= 0); + herr_t errf; + hid_t grp_id; + grp_id = H5Gopen2(file_id, "ports", H5P_DEFAULT); + assert(grp_id >= 0); + + int num_ref_ports_; + hdf5_utils::ReadAttribute(grp_id, "number_of_ports", num_ref_ports_); + assert(num_ref_ports_ >= num_ref_ports); + + std::string dset_name; + int c1, c2, a1, a2; + for (int p = 0; p < num_ref_ports; p++) + { + topol_handler->GetRefPortInfo(p, c1, c2, a1, a2); + dset_name = topol_handler->GetComponentName(c1) + ":" + topol_handler->GetComponentName(c2); + dset_name += "-" + std::to_string(a1) + ":" + std::to_string(a2); + hdf5_utils::ReadDataset(grp_id, dset_name, *port[p]); + } + + errf = H5Gclose(grp_id); + assert(errf >= 0); +} \ No newline at end of file diff --git a/src/sample_generator.cpp b/src/sample_generator.cpp index 967a2406..3a6d5658 100644 --- a/src/sample_generator.cpp +++ b/src/sample_generator.cpp @@ -182,18 +182,22 @@ void SampleGenerator::SaveSnapshot(BlockVector *U_snapshots, std::vectorNumBlocks() == snapshot_basis_tags.size()); col_idxs.SetSize(U_snapshots->NumBlocks()); + /* add snapshots according to their tags */ for (int s = 0; s < snapshot_basis_tags.size(); s++) { + /* if the tag was never seen before, create a new snapshot generator */ if (!basis_tag2idx.count(snapshot_basis_tags[s])) { const int fom_vdofs = U_snapshots->BlockSize(s); AddSnapshotGenerator(fom_vdofs, GetSamplePrefix(), snapshot_basis_tags[s]); } + /* add the snapshot into the corresponding snapshot generator */ int index = basis_tag2idx[snapshot_basis_tags[s]]; bool addSample = snapshot_generators[index]->takeSample(U_snapshots->GetBlock(s).GetData()); assert(addSample); + /* save the column index in each snapshot matrix, for port data. */ col_idxs[s] = snapshot_generators[index]->getNumSamples(); } } @@ -207,29 +211,42 @@ void SampleGenerator::SaveSnapshotPorts(TopologyHandler *topol_handler, const Tr assert(col_idxs.Size() % numSub == 0); const int num_var = col_idxs.Size() / numSub; - /* - We save the snapshot ports by dom%d (INDIVIDUAL) or component names (UNIVERSAL). - We does not consider basis separation at this point. - Though we use basis tags here, it is only for the sake of - getting consistent domain/component mesh names. - */ + /* iterate every port */ for (int p = 0; p < numPorts; p++) { const PortInfo *port = topol_handler->GetPortInfo(p); + + /* + We save the snapshot ports by dom%d (INDIVIDUAL) or component names (UNIVERSAL). + Variable separation does not matter at this point. + Though we use basis tags here, it is only for the sake of + getting consistent domain/component mesh names. + */ std::string mesh1, mesh2; mesh1 = GetBasisTag(port->Mesh1, train_mode, topol_handler); mesh2 = GetBasisTag(port->Mesh2, train_mode, topol_handler); + /* + note the mesh names are not necessarily the same as the subdomain names + depending on the train mode. + */ PortTag tag = {.Mesh1 = mesh1, .Mesh2 = mesh2, .Attr1 = port->Attr1, .Attr2 = port->Attr2}; + + /* if the port was never seen, create a new collector */ if (!port_tag2idx.count(tag)) { port_tag2idx[tag] = port_tags.size(); port_tags.push_back(tag); - port_colidxs.Append(new Array(0)); + port_colidxs.Append(new Array2D); } const int index = port_tag2idx[tag]; int col1, col2; + /* + column indices in snapshot matrices must be identical for all variables, + whether variables are separated for training or not. + We thus pick the column index of the first variable. + */ col1 = col_idxs[num_var * (port->Mesh1)]; col2 = col_idxs[num_var * (port->Mesh2)]; for (int v = 0; v < num_var; v++) @@ -238,8 +255,10 @@ void SampleGenerator::SaveSnapshotPorts(TopologyHandler *topol_handler, const Tr assert(col2 == col_idxs[v + num_var * (port->Mesh2)]); } - port_colidxs[index]->Append(col1); - port_colidxs[index]->Append(col2); + int row_idx = port_colidxs[index]->NumRows(); + port_colidxs[index]->SetSize(row_idx + 1, 2); + (*port_colidxs[index])(row_idx, 0) = col1; + (*port_colidxs[index])(row_idx, 1) = col2; } } diff --git a/src/steady_ns_solver.cpp b/src/steady_ns_solver.cpp index a37efe7d..2fc0efeb 100644 --- a/src/steady_ns_solver.cpp +++ b/src/steady_ns_solver.cpp @@ -479,7 +479,7 @@ void SteadyNSSolver::LoadROMOperatorFromFile(const std::string input_prefix) } } -bool SteadyNSSolver::Solve() +bool SteadyNSSolver::Solve(SampleGenerator *sample_generator) { int maxIter = config.GetOption("solver/max_iter", 100); double rtol = config.GetOption("solver/relative_tolerance", 1.e-10); @@ -566,6 +566,10 @@ bool SteadyNSSolver::Solve() SortBySubdomains(sol_byvar, *U); + /* save solution if sample generator is provided */ + if (converged && sample_generator) + SaveSnapshots(sample_generator); + return converged; } @@ -769,23 +773,7 @@ void SteadyNSSolver::AllocateROMEQPElems() const int num_comp = topol_handler->GetNumComponents(); comp_eqps.SetSize(num_comp); - comp_fes.SetSize(num_comp); comp_eqps = NULL; - comp_fes = NULL; - - for (int c = 0; c < num_comp; c++) - { - int midx = -1; - for (int m = 0; m < numSub; m++) - if (rom_handler->GetRefIndexForSubdomain(m) == c) - { - midx = m; - break; - } - assert((midx >= 0) && (midx < numSub)); - - comp_fes[c] = ufes[midx]; - } DenseMatrix *basis; for (int c = 0; c < num_comp; c++) @@ -796,7 +784,7 @@ void SteadyNSSolver::AllocateROMEQPElems() auto nl_integ_tmp = new VectorConvectionTrilinearFormIntegrator(*zeta_coeff); nl_integ_tmp->SetIntRule(ir_nl); - comp_eqps[c] = new ROMNonlinearForm(basis->NumCols(), comp_fes[c]); + comp_eqps[c] = new ROMNonlinearForm(basis->NumCols(), comp_fes[c * num_var]); comp_eqps[c]->AddDomainIntegrator(nl_integ_tmp); comp_eqps[c]->SetBasis(*basis); comp_eqps[c]->SetPrecomputeMode(precompute); diff --git a/src/stokes_solver.cpp b/src/stokes_solver.cpp index a18b0afa..7e6579df 100644 --- a/src/stokes_solver.cpp +++ b/src/stokes_solver.cpp @@ -521,16 +521,14 @@ void StokesSolver::AssembleInterfaceMatrices() b_itf->AssembleInterfaceMatrices(b_mats); } -void StokesSolver::BuildCompROMLinElems(Array &fes_comp) +void StokesSolver::BuildCompROMLinElems() { assert(topol_mode == TopologyHandlerMode::COMPONENT); assert(train_mode == UNIVERSAL); assert(rom_handler->BasisLoaded()); + assert(rom_elems); const int num_comp = topol_handler->GetNumComponents(); - assert(comp_mats.Size() == num_comp); - assert(fes_comp.Size() == num_comp * num_var); - assert(nu_coeff); int num_blocks = (separate_variable_basis) ? num_var: 1; @@ -539,8 +537,8 @@ void StokesSolver::BuildCompROMLinElems(Array &fes_comp) const int fidx = c * num_var; Mesh *comp = topol_handler->GetComponentMesh(c); - BilinearForm m_comp(fes_comp[fidx]); - MixedBilinearFormDGExtension b_comp(fes_comp[fidx], fes_comp[fidx+1]); + BilinearForm m_comp(comp_fes[fidx]); + MixedBilinearFormDGExtension b_comp(comp_fes[fidx], comp_fes[fidx+1]); m_comp.AddDomainIntegrator(new VectorDiffusionIntegrator(*nu_coeff)); if (full_dg) @@ -559,18 +557,18 @@ void StokesSolver::BuildCompROMLinElems(Array &fes_comp) SparseMatrix *b_mat = &(b_comp.SpMat()); SparseMatrix *bt_mat = Transpose(*b_mat); - comp_mats[c]->SetSize(num_blocks, num_blocks); + rom_elems->comp[c]->SetSize(num_blocks, num_blocks); if (separate_variable_basis) { - (*comp_mats[c])(0, 0) = rom_handler->ProjectToRefBasis(fidx, fidx, m_mat); - (*comp_mats[c])(1, 0) = rom_handler->ProjectToRefBasis(fidx+1, fidx, b_mat); - (*comp_mats[c])(0, 1) = rom_handler->ProjectToRefBasis(fidx, fidx+1, bt_mat); + (*rom_elems->comp[c])(0, 0) = rom_handler->ProjectToRefBasis(fidx, fidx, m_mat); + (*rom_elems->comp[c])(1, 0) = rom_handler->ProjectToRefBasis(fidx+1, fidx, b_mat); + (*rom_elems->comp[c])(0, 1) = rom_handler->ProjectToRefBasis(fidx, fidx+1, bt_mat); } else { Array dummy1, dummy2; BlockMatrix *sys_comp = FormBlockMatrix(m_mat, b_mat, bt_mat, dummy1, dummy2); - (*comp_mats[c])(0, 0) = rom_handler->ProjectToRefBasis(c, c, sys_comp); + (*rom_elems->comp[c])(0, 0) = rom_handler->ProjectToRefBasis(c, c, sys_comp); delete sys_comp; } @@ -578,16 +576,14 @@ void StokesSolver::BuildCompROMLinElems(Array &fes_comp) } } -void StokesSolver::BuildBdrROMLinElems(Array &fes_comp) +void StokesSolver::BuildBdrROMLinElems() { assert(topol_mode == TopologyHandlerMode::COMPONENT); assert(train_mode == UNIVERSAL); assert(rom_handler->BasisLoaded()); + assert(rom_elems); const int num_comp = topol_handler->GetNumComponents(); - assert(bdr_mats.Size() == num_comp); - assert(fes_comp.Size() == num_comp * num_var); - assert(nu_coeff); int num_blocks = (separate_variable_basis) ? num_var: 1; @@ -595,7 +591,7 @@ void StokesSolver::BuildBdrROMLinElems(Array &fes_comp) { const int fidx = c * num_var; Mesh *comp = topol_handler->GetComponentMesh(c); - assert(bdr_mats[c]->Size() == comp->bdr_attributes.Size()); + assert(rom_elems->bdr[c]->Size() == comp->bdr_attributes.Size()); for (int b = 0; b < comp->bdr_attributes.Size(); b++) { @@ -603,8 +599,8 @@ void StokesSolver::BuildBdrROMLinElems(Array &fes_comp) bdr_marker = 0; bdr_marker[comp->bdr_attributes[b] - 1] = 1; - BilinearForm m_comp(fes_comp[fidx]); - MixedBilinearFormDGExtension b_comp(fes_comp[fidx], fes_comp[fidx+1]); + BilinearForm m_comp(comp_fes[fidx]); + MixedBilinearFormDGExtension b_comp(comp_fes[fidx], comp_fes[fidx+1]); m_comp.AddBdrFaceIntegrator(new DGVectorDiffusionIntegrator(*nu_coeff, sigma, kappa), bdr_marker); b_comp.AddBdrFaceIntegrator(new DGNormalFluxIntegrator, bdr_marker); @@ -618,7 +614,7 @@ void StokesSolver::BuildBdrROMLinElems(Array &fes_comp) SparseMatrix *b_mat = &(b_comp.SpMat()); SparseMatrix *bt_mat = Transpose(*b_mat); - MatrixBlocks *bdr_mat = (*bdr_mats[c])[b]; + MatrixBlocks *bdr_mat = (*rom_elems->bdr[c])[b]; bdr_mat->SetSize(num_blocks, num_blocks); if (separate_variable_basis) { @@ -639,25 +635,24 @@ void StokesSolver::BuildBdrROMLinElems(Array &fes_comp) } } -void StokesSolver::BuildItfaceROMLinElems(Array &fes_comp) +void StokesSolver::BuildItfaceROMLinElems() { assert(topol_mode == TopologyHandlerMode::COMPONENT); assert(train_mode == UNIVERSAL); assert(rom_handler->BasisLoaded()); + assert(rom_elems); const int num_comp = topol_handler->GetNumComponents(); - assert(fes_comp.Size() == num_comp * num_var); Array ufes_comp(num_comp), pfes_comp(num_comp); for (int c = 0; c < num_comp; c++) { - ufes_comp[c] = fes_comp[c * num_var]; - pfes_comp[c] = fes_comp[c * num_var + 1]; + ufes_comp[c] = comp_fes[c * num_var]; + pfes_comp[c] = comp_fes[c * num_var + 1]; } int num_blocks = (separate_variable_basis) ? 2 * num_var: 2; const int num_ref_ports = topol_handler->GetNumRefPorts(); - assert(port_mats.Size() == num_ref_ports); for (int p = 0; p < num_ref_ports; p++) { int c1, c2; @@ -686,22 +681,25 @@ void StokesSolver::BuildItfaceROMLinElems(Array &fes_comp) // NOTE: the index also should be transposed. bt_mats_p(j, i) = Transpose(*b_mats_p(i, j)); - port_mats[p]->SetSize(num_blocks, num_blocks); + rom_elems->port[p]->SetSize(num_blocks, num_blocks); for (int i = 0; i < 2; i++) for (int j = 0; j < 2; j++) if (separate_variable_basis) { - (*port_mats[p])(i * num_var, j * num_var) = rom_handler->ProjectToRefBasis(c_idx[i], c_idx[j], m_mats_p(i, j)); - (*port_mats[p])(i * num_var + 1, j * num_var) = rom_handler->ProjectToRefBasis(c_idx[i] + 1, c_idx[j], b_mats_p(i, j)); - (*port_mats[p])(i * num_var, j * num_var + 1) = rom_handler->ProjectToRefBasis(c_idx[i], c_idx[j] + 1, bt_mats_p(i, j)); + (*rom_elems->port[p])(i * num_var, j * num_var) = + rom_handler->ProjectToRefBasis(c_idx[i], c_idx[j], m_mats_p(i, j)); + (*rom_elems->port[p])(i * num_var + 1, j * num_var) = + rom_handler->ProjectToRefBasis(c_idx[i] + 1, c_idx[j], b_mats_p(i, j)); + (*rom_elems->port[p])(i * num_var, j * num_var + 1) = + rom_handler->ProjectToRefBasis(c_idx[i], c_idx[j] + 1, bt_mats_p(i, j)); } else { Array dummy1, dummy2; BlockMatrix *tmp_mat = FormBlockMatrix(m_mats_p(i,j), b_mats_p(i,j), bt_mats_p(i,j), dummy1, dummy2); - (*port_mats[p])(i, j) = rom_handler->ProjectToRefBasis(c_idx[i], c_idx[j], tmp_mat); + (*rom_elems->port[p])(i, j) = rom_handler->ProjectToRefBasis(c_idx[i], c_idx[j], tmp_mat); delete tmp_mat; } @@ -711,7 +709,7 @@ void StokesSolver::BuildItfaceROMLinElems(Array &fes_comp) } // for (int p = 0; p < num_ref_ports; p++) } -bool StokesSolver::Solve() +bool StokesSolver::Solve(SampleGenerator *sample_generator) { // If using direct solver, returns always true. bool converged = true; @@ -814,6 +812,10 @@ bool StokesSolver::Solve() SortBySubdomains(sol_byvar, *U); + /* save solution if sample generator is provided */ + if (converged && sample_generator) + SaveSnapshots(sample_generator); + return converged; } diff --git a/src/unsteady_ns_solver.cpp b/src/unsteady_ns_solver.cpp index d7b41234..8d67c7c2 100644 --- a/src/unsteady_ns_solver.cpp +++ b/src/unsteady_ns_solver.cpp @@ -77,7 +77,7 @@ void UnsteadyNSSolver::AssembleOperator() massMat->SetBlock(i, i, &(mass[i]->SpMat())); } -bool UnsteadyNSSolver::Solve() +bool UnsteadyNSSolver::Solve(SampleGenerator *sample_generator) { bool converged = true; int initial_step = 0; @@ -92,6 +92,9 @@ bool UnsteadyNSSolver::Solve() LoadSolutionWithTime(restart_file, 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); + InitializeTimeIntegration(); SortByVariables(*U, *U_step); @@ -116,6 +119,10 @@ bool UnsteadyNSSolver::Solve() 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)) + SaveSnapshots(sample_generator); } SortBySubdomains(*U_step, *U);