diff --git a/include/hdf5_utils.hpp b/include/hdf5_utils.hpp index d89c71c5..8394ce4e 100644 --- a/include/hdf5_utils.hpp +++ b/include/hdf5_utils.hpp @@ -8,6 +8,22 @@ #include "mfem.hpp" #include "hdf5.h" #include "linalg_utils.hpp" +#include "rom_handler.hpp" +#include "hyperreduction_integ.hpp" + +namespace mfem +{ + +enum IntegratorType +{ + DOMAIN, + INTERIORFACE, + BDRFACE, + INTERFACE, + NUM_INTEG_TYPE +}; + +} using namespace mfem; @@ -25,6 +41,9 @@ int GetDatasetSize(hid_t &source, std::string dataset, hsize_t* &dims); void ReadAttribute(hid_t &source, std::string attribute, std::string &value); void WriteAttribute(hid_t &source, std::string attribute, const std::string &value); +void ReadAttribute(hid_t &source, std::string attribute, BasisTag &value); +void WriteAttribute(hid_t &source, std::string attribute, const BasisTag &value); + template void ReadAttribute(hid_t &source, std::string attribute, T &value) { herr_t status; @@ -68,8 +87,11 @@ void ReadDataset(hid_t &source, std::string dataset, Array &value) assert(errf >= 0); value.SetSize(dims[0]); - hid_t dataType = hdf5_utils::GetType(value[0]); - H5Dread(dset_id, dataType, H5S_ALL, H5S_ALL, H5P_DEFAULT, value.Write()); + if (dims[0] > 0) + { + hid_t dataType = hdf5_utils::GetType(value[0]); + H5Dread(dset_id, dataType, H5S_ALL, H5S_ALL, H5P_DEFAULT, value.Write()); + } errf = H5Dclose(dset_id); assert(errf >= 0); @@ -115,7 +137,8 @@ void WriteDataset(hid_t &source, std::string dataset, const Array &value) { herr_t errf = 0; - hid_t dataType = hdf5_utils::GetType(value[0]); + T tmp; + hid_t dataType = hdf5_utils::GetType(tmp); hsize_t dims[1]; dims[0] = value.Size(); @@ -125,8 +148,11 @@ void WriteDataset(hid_t &source, std::string dataset, const Array &value) hid_t dset_id = H5Dcreate2(source, dataset.c_str(), dataType, dspace_id, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT); assert(dset_id >= 0); - errf = H5Dwrite(dset_id, dataType, H5S_ALL, H5S_ALL, H5P_DEFAULT, value.Read()); - assert(errf >= 0); + if (dims[0] > 0) + { + errf = H5Dwrite(dset_id, dataType, H5S_ALL, H5S_ALL, H5P_DEFAULT, value.Read()); + assert(errf >= 0); + } errf = H5Dclose(dset_id); assert(errf >= 0); @@ -174,6 +200,9 @@ void WriteDataset(hid_t &source, std::string dataset, const DenseTensor &value); void ReadDataset(hid_t &source, std::string dataset, MatrixBlocks &value); void WriteDataset(hid_t &source, std::string dataset, const MatrixBlocks &value); +void ReadDataset(hid_t &source, std::string dataset, const IntegratorType type, Array &value); +void WriteDataset(hid_t &source, std::string dataset, const IntegratorType type, const Array &value); + inline bool pathExists(hid_t id, const std::string& path) { return H5Lexists(id, path.c_str(), H5P_DEFAULT) > 0; diff --git a/include/hyperreduction_integ.hpp b/include/hyperreduction_integ.hpp index b29fe841..75698591 100644 --- a/include/hyperreduction_integ.hpp +++ b/include/hyperreduction_integ.hpp @@ -11,12 +11,13 @@ namespace mfem { struct SampleInfo { - // TODO(kevin): all integrators could use the same integer variable, reducing the memory size of this struct. - int el; // element index (used for DomainIntegrator) - int face; // face index (used for InteriorFaceIntegrator) - int be; // boundary element index (used for BdrFaceIntegrator) - int itf; // interface info index (used for InterfaceIntegrator) - // can add dofs for other hyper reductions. + /* + - For DomainIntegrator: element index + - For BdrFaceIntegrator: boundary element index + - For InteriorFaceIntegrator: face index + - For InterfaceIntegrator: interface info index + */ + int el; int qp; // quadrature point double qw; // quadrature weight diff --git a/include/main_workflow.hpp b/include/main_workflow.hpp index 8fcf9043..946701a2 100644 --- a/include/main_workflow.hpp +++ b/include/main_workflow.hpp @@ -15,10 +15,12 @@ void RunExample(); MultiBlockSolver* InitSolver(); SampleGenerator* InitSampleGenerator(MPI_Comm comm); -std::vector GetGlobalBasisTagList(const TopologyHandlerMode &topol_mode, bool separate_variable_basis); +std::vector GetGlobalBasisTagList(const TopologyHandlerMode &topol_mode, bool separate_variable_basis); void GenerateSamples(MPI_Comm comm); void CollectSamples(SampleGenerator *sample_generator); +void CollectSamplesByPort(SampleGenerator *sample_generator, const std::string &basis_prefix); +void CollectSamplesByBasis(SampleGenerator *sample_generator, const std::string &basis_prefix); void BuildROM(MPI_Comm comm); void TrainROM(MPI_Comm comm); // supremizer-enrichment etc.. @@ -26,7 +28,7 @@ void AuxiliaryTrainROM(MPI_Comm comm, SampleGenerator *sample_generator); // EQP training, could include hypre-reduction optimization. void TrainEQP(MPI_Comm comm); // Input parsing routine to list out all snapshot files for training a basis. -void FindSnapshotFilesForBasis(const std::string &basis_tag, const std::string &default_filename, std::vector &file_list); +void FindSnapshotFilesForBasis(const BasisTag &basis_tag, const std::string &default_filename, std::vector &file_list); // return relative error if comparing solution. double SingleRun(MPI_Comm comm, const std::string output_file = ""); diff --git a/include/multiblock_solver.hpp b/include/multiblock_solver.hpp index 9311ccc6..6901cb26 100644 --- a/include/multiblock_solver.hpp +++ b/include/multiblock_solver.hpp @@ -228,36 +228,23 @@ friend class ParameterizedProblem; void LoadSolutionWithTime(const std::string &filename, int &step, double &time); void CopySolution(BlockVector *input_sol); - virtual void AllocateROMTensorElems() - { mfem_error("Abstract method MultiBlockSolver::AllocateROMTensorElems!\n"); } + virtual void AllocateROMNlinElems() + { mfem_error("Abstract method MultiBlockSolver::AllocateROMNlinElems!\n"); } virtual void BuildROMTensorElems() { mfem_error("Abstract method MultiBlockSolver::BuildROMTensorElems!\n"); } - virtual void SaveROMTensorElems(const std::string &filename) - { mfem_error("Abstract method MultiBlockSolver::SaveROMTensorElems!\n"); } - virtual void LoadROMTensorElems(const std::string &filename) - { mfem_error("Abstract method MultiBlockSolver::LoadROMTensorElems!\n"); } - virtual void AssembleROMTensorOper() - { mfem_error("Abstract method MultiBlockSolver::AssembleROMTensorOper!\n"); } - - virtual void AllocateROMEQPElems() - { mfem_error("Abstract method MultiBlockSolver::AllocateROMEQPElems!\n"); } - virtual void TrainEQPElems(SampleGenerator *sample_generator) - { mfem_error("Abstract method MultiBlockSolver::TrainEQPElems!\n"); } - virtual void SaveEQPElems(const std::string &filename) - { mfem_error("Abstract method MultiBlockSolver::SaveEQP!\n"); } - virtual void LoadEQPElems(const std::string &filename) - { mfem_error("Abstract method MultiBlockSolver::LoadEQP!\n"); } - virtual void AssembleROMEQPOper() - { mfem_error("Abstract method MultiBlockSolver::AssembleROMEQPOper!\n"); } - - void AllocateROMNlinElems(); - void LoadROMNlinElems(const std::string &input_prefix); - void AssembleROMNlinOper(); - - void InitROMHandler(); - void GetBasisTags(std::vector &basis_tags); - - virtual BlockVector* PrepareSnapshots(std::vector &basis_tags); + virtual void TrainROMEQPElems(SampleGenerator *sample_generator) + { mfem_error("Abstract method MultiBlockSolver::TrainROMEQPElems!\n"); } + virtual void SaveROMNlinElems(const std::string &input_prefix) + { mfem_error("Abstract method MultiBlockSolver::SaveROMNlinElems!\n"); } + virtual void LoadROMNlinElems(const std::string &input_prefix) + { mfem_error("Abstract method MultiBlockSolver::LoadROMNlinElems!\n"); } + virtual void AssembleROMNlinOper() + { mfem_error("Abstract method MultiBlockSolver::AssembleROMNlinOper!\n"); } + + virtual void InitROMHandler(); + void GetBasisTags(std::vector &basis_tags); + + virtual BlockVector* PrepareSnapshots(std::vector &basis_tags); void SaveSnapshots(SampleGenerator *sample_generator); virtual void LoadReducedBasis() { rom_handler->LoadReducedBasis(); } diff --git a/include/parameterized_problem.hpp b/include/parameterized_problem.hpp index bbb55e70..70d287d9 100644 --- a/include/parameterized_problem.hpp +++ b/include/parameterized_problem.hpp @@ -82,6 +82,13 @@ namespace backward_facing_step void pic(const Vector &x, double t, Vector &y); } +namespace lid_driven_cavity +{ + extern double u0, L; + + void ubdr(const Vector &x, double t, Vector &y); +} + namespace periodic_flow_past_array { extern Vector f; @@ -339,12 +346,26 @@ class BackwardFacingStep : public FlowProblem virtual ~BackwardFacingStep() {}; }; +class LidDrivenCavity : public FlowProblem +{ +public: + LidDrivenCavity(); + virtual ~LidDrivenCavity() {}; +}; + class PeriodicFlowPastArray : public FlowProblem { public: PeriodicFlowPastArray(); }; +class ForceDrivenCorner : public PeriodicFlowPastArray +{ +public: + ForceDrivenCorner(); + virtual ~ForceDrivenCorner() {}; +}; + class LinElastProblem : public ParameterizedProblem { friend class LinElastSolver; diff --git a/include/rom_element_collection.hpp b/include/rom_element_collection.hpp index fd7b04a3..99471bce 100644 --- a/include/rom_element_collection.hpp +++ b/include/rom_element_collection.hpp @@ -100,7 +100,7 @@ class ROMEQPElement : public ROMElementCollection Array comp; // Size(num_components); // boundary condition is enforced via forcing term. Array *> bdr; - Array port; // reference ports. + ROMInterfaceForm *port = NULL; // reference ports. public: ROMEQPElement(TopologyHandler *topol_handler_, diff --git a/include/rom_handler.hpp b/include/rom_handler.hpp index 98df6134..1dd66d96 100644 --- a/include/rom_handler.hpp +++ b/include/rom_handler.hpp @@ -30,8 +30,50 @@ enum NonlinearHandling NUM_NLNHNDL }; -const std::string GetBasisTagForComponent(const int &comp_idx, const TopologyHandler *topol_handler, const std::string var_name=""); -const std::string GetBasisTag(const int &subdomain_index, const TopologyHandler *topol_handler, const std::string var_name=""); +struct BasisTag +{ + /* component mesh name */ + std::string comp = ""; + /* variable name, if separate basis is used */ + std::string var = ""; + + BasisTag() {} + + BasisTag(const std::string &comp_, const std::string &var_="") + : comp(comp_), var(var_) {} + + const std::string print() const + { + std::string tag = comp; + if (var != "") + tag += "_" + var; + return tag; + } + + bool operator==(const BasisTag &tag) const + { + return ((comp == tag.comp) && (var == tag.var)); + } + + bool operator<(const BasisTag &tag) const + { + if (comp == tag.comp) + return (var < tag.var); + + return (comp < tag.comp); + } + + BasisTag& operator=(const BasisTag &tag) + { + comp = tag.comp; + var = tag.var; + return *this; + } +}; + + +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=""); class ROMHandlerBase { @@ -81,7 +123,7 @@ class ROMHandlerBase Array dim_ref_basis; Array carom_ref_basis; Array rom_comp_block_offsets; - std::vector basis_tags; + std::vector basis_tags; bool basis_loaded; bool operator_loaded; @@ -126,7 +168,7 @@ class ROMHandlerBase const bool SeparateVariable() { return separate_variable; } const std::string GetOperatorPrefix() { return operator_prefix; } const std::string GetBasisPrefix() { return basis_prefix; } - const std::string GetRefBasisTag(const int ref_idx) { return basis_tags[ref_idx]; } + 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; diff --git a/include/rom_interfaceform.hpp b/include/rom_interfaceform.hpp index 2d641da8..0f20ef07 100644 --- a/include/rom_interfaceform.hpp +++ b/include/rom_interfaceform.hpp @@ -7,6 +7,7 @@ #include "interface_form.hpp" #include "rom_handler.hpp" +#include "hdf5_utils.hpp" namespace mfem { @@ -15,13 +16,24 @@ class ROMInterfaceForm : public InterfaceForm { protected: const int numPorts; + const int numRefPorts; + const int num_comp; + + /* component finite element space */ + Array comp_fes; // not owned + + /* size of num_comp */ + Array comp_basis; // not owned + Array comp_basis_dof_offsets; /* size of numSub */ + /* view Arrays of component arrays */ Array basis; // not owned Array basis_dof_offsets; /* EQP sample point info */ /* + View Array of fnfi_ref_sample. Array of size (fnfi.Size() * topol_handler->GetNumPorts()), where each element is another Array of EQP samples at the given port p and the given integrator i. @@ -30,15 +42,25 @@ class ROMInterfaceForm : public InterfaceForm */ Array *> fnfi_sample; + /* + Array of size (fnfi.Size() * topol_handler->GetNumRefPorts()), + where each element is another Array of EQP samples + at the given reference port p and the given integrator i. + For the reference port p and the integrator i, + fnfi_sample[p + i * numRefPorts] + */ + Array *> fnfi_ref_sample; + /// @brief Flag for precomputing necessary coefficients for fast computation. bool precompute = false; public: - ROMInterfaceForm(Array &meshes_, Array &fes_, TopologyHandler *topol_); + ROMInterfaceForm(Array &meshes_, Array &fes_, + Array &comp_fes_, TopologyHandler *topol_); virtual ~ROMInterfaceForm(); - void SetBasisAtSubdomain(const int m, DenseMatrix &basis_, const int offset=0); + void SetBasisAtComponent(const int c, DenseMatrix &basis_, const int offset=0); void UpdateBlockOffsets(); /// Adds new Interior Face Integrator. @@ -46,32 +68,43 @@ class ROMInterfaceForm : public InterfaceForm { InterfaceForm::AddInterfaceIntegrator(nlfi); + for (int p = 0; p < numRefPorts; p++) + fnfi_ref_sample.Append(NULL); for (int p = 0; p < numPorts; p++) fnfi_sample.Append(NULL); } - void UpdateInterFaceIntegratorSampling(const int i, const int p, const Array &itf, - const Array &qp, const Array &qw) + void UpdateInterFaceIntegratorSampling(const int i, const int rp, + const Array &samples) { assert((i >= 0) && (i < fnfi.Size())); - assert((p >= 0) && (p < numPorts)); - assert(fnfi_sample.Size() == fnfi.Size() * numPorts); + assert((rp >= 0) && (rp < numRefPorts)); + assert(fnfi_ref_sample.Size() == fnfi.Size() * numRefPorts); + + const int idx = rp + i * numRefPorts; + if (fnfi_ref_sample[idx]) + delete fnfi_ref_sample[idx]; - const int idx = p + i * numPorts; - if (fnfi_sample[idx]) - delete fnfi_sample[idx]; + fnfi_ref_sample[idx] = new Array(samples); - fnfi_sample[idx] = new Array(0); - for (int s = 0; s < itf.Size(); s++) - fnfi_sample[idx]->Append({.el=-1, .face=-1, .be=-1, .itf=itf[s], .qp=qp[s], .qw=qw[s]}); + for (int p = 0; p < numPorts; p++) + { + if (topol_handler->GetPortType(p) != rp) + continue; + + fnfi_sample[p + i * numPorts] = fnfi_ref_sample[idx]; + } } void InterfaceAddMult(const Vector &x, Vector &y) const override; void InterfaceGetGradient(const Vector &x, Array2D &mats) const override; + void TrainEQPForRefPort(const int p, const CAROM::Matrix &snapshot1, const CAROM::Matrix &snapshot2, + const Array2D &snap_pair_idx, const double eqp_tol); + void SetupEQPSystem(const CAROM::Matrix &snapshot1, const CAROM::Matrix &snapshot2, - const Array &snap1_idx, const Array &snap2_idx, + const Array2D &snap_pair_idx, DenseMatrix &basis1, DenseMatrix &basis2, const int &basis1_offset, const int &basis2_offset, FiniteElementSpace *fes1, FiniteElementSpace *fes2, @@ -81,7 +114,10 @@ class ROMInterfaceForm : public InterfaceForm void TrainEQPForIntegrator(const int nqe, const CAROM::Matrix &Gt, const CAROM::Vector &rhs_Gw, const double eqp_tol, - Array &sample_el, Array &sample_qp, Array &sample_qw); + Array &samples); + + void SaveEQPForIntegrator(const int k, hid_t file_id, const std::string &dsetname); + void LoadEQPForIntegrator(const int k, hid_t file_id, const std::string &dsetname); private: /* These methods are not available in this class */ diff --git a/include/rom_nonlinearform.hpp b/include/rom_nonlinearform.hpp index bdc64fbd..bec0c609 100644 --- a/include/rom_nonlinearform.hpp +++ b/include/rom_nonlinearform.hpp @@ -13,14 +13,6 @@ namespace mfem { -enum IntegratorType -{ - DOMAIN, - INTERIORFACE, - BDRFACE, - NUM_INTEG_TYPE -}; - class ROMNonlinearForm : public NonlinearForm { private: @@ -71,7 +63,7 @@ class ROMNonlinearForm : public NonlinearForm void TrainEQP(const CAROM::Matrix &snapshots, const double eqp_tol = 1.0e-2); void TrainEQPForIntegrator(HyperReductionIntegrator *nlfi, const CAROM::Matrix &Gt, const CAROM::Vector &rhs_Gw, const double eqp_tol, - Array &sample_el, Array &sample_qp, Array &sample_qw); + Array &samples); void SetupEQPSystemForDomainIntegrator(const CAROM::Matrix &snapshots, HyperReductionIntegrator *nlfi, CAROM::Matrix &Gt, CAROM::Vector &rhs_Gw); void SetupEQPSystemForInteriorFaceIntegrator(const CAROM::Matrix &snapshots, HyperReductionIntegrator *nlfi, @@ -79,7 +71,7 @@ class ROMNonlinearForm : public NonlinearForm void SetupEQPSystemForBdrFaceIntegrator(const CAROM::Matrix &snapshots, HyperReductionIntegrator *nlfi, const Array &bdr_attr_marker, CAROM::Matrix &Gt, CAROM::Vector &rhs_Gw, Array &bidxs); - void GetEQPForIntegrator(const IntegratorType type, const int k, Array &sample_el, Array &sample_qp, Array &sample_qw); + Array* GetEQPForIntegrator(const IntegratorType type, const int k); void SaveEQPForIntegrator(const IntegratorType type, const int k, hid_t file_id, const std::string &dsetname); void LoadEQPForIntegrator(const IntegratorType type, const int k, hid_t file_id, const std::string &dsetname); @@ -90,14 +82,12 @@ class ROMNonlinearForm : public NonlinearForm dnfi_sample.Append(NULL); } - void UpdateDomainIntegratorSampling(const int i, const Array &el, const Array &qp, const Array &qw) + void UpdateDomainIntegratorSampling(const int i, const Array &samples) { assert((i >= 0) && (i < dnfi.Size())); assert(dnfi.Size() == dnfi_sample.Size()); - dnfi_sample[i] = new Array(0); - for (int s = 0; s < el.Size(); s++) - dnfi_sample[i]->Append({.el=el[s], .face=-1, .be=-1, .qp=qp[s], .qw=qw[s]}); + dnfi_sample[i] = new Array(samples); } /// Access all integrators added with AddDomainIntegrator(). @@ -111,14 +101,12 @@ class ROMNonlinearForm : public NonlinearForm fnfi_sample.Append(NULL); } - void UpdateInteriorFaceIntegratorSampling(const int i, const Array &face, const Array &qp, const Array &qw) + void UpdateInteriorFaceIntegratorSampling(const int i, const Array &samples) { assert((i >= 0) && (i < fnfi.Size())); assert(fnfi.Size() == fnfi_sample.Size()); - fnfi_sample[i] = new Array(0); - for (int s = 0; s < face.Size(); s++) - fnfi_sample[i]->Append({.el=-1, .face=face[s], .be=-1, .qp=qp[s], .qw=qw[s]}); + fnfi_sample[i] = new Array(samples); } /** @brief Access all interior face integrators added with @@ -140,18 +128,17 @@ class ROMNonlinearForm : public NonlinearForm Array &bdr_marker) { bfnfi.Append(nfi); - bfnfi_marker.Append(&bdr_marker); + // Decided to own it, as opposed to the parent class. + bfnfi_marker.Append(new Array(bdr_marker)); bfnfi_sample.Append(NULL); } - void UpdateBdrFaceIntegratorSampling(const int i, const Array &be, const Array &qp, const Array &qw) + void UpdateBdrFaceIntegratorSampling(const int i, const Array &samples) { assert((i >= 0) && (i < bfnfi.Size())); assert(bfnfi.Size() == bfnfi_sample.Size()); - bfnfi_sample[i] = new Array(0); - for (int s = 0; s < be.Size(); s++) - bfnfi_sample[i]->Append({.el=-1, .face=-1, .be=be[s], .qp=qp[s], .qw=qw[s]}); + bfnfi_sample[i] = new Array(samples); } /** @brief Access all boundary face integrators added with diff --git a/include/sample_generator.hpp b/include/sample_generator.hpp index 5bca52d2..91a20bcd 100644 --- a/include/sample_generator.hpp +++ b/include/sample_generator.hpp @@ -61,8 +61,8 @@ class SampleGenerator Array snapshot_options; Array snapshot_generators; // each snapshot will be sorted out by its basis tag. - std::vector basis_tags; - std::map basis_tag2idx; + std::vector basis_tags; + std::map basis_tag2idx; /* snapshot pairs per interface port, for nonlinear interface EQP */ std::vector port_tags; @@ -84,8 +84,8 @@ class SampleGenerator Parameter* GetParam(const int &k) { return params[k]; } const std::string GetSamplePrefix() { return sample_dir + "/" + sample_prefix + "_sample"; } - const std::string GetBaseFilename(const std::string &prefix, const std::string &basis_tag) - { return prefix + "_" + basis_tag; } + const std::string GetBaseFilename(const std::string &prefix, const BasisTag &basis_tag) + { return prefix + "_" + basis_tag.print(); } // Generate parameter space as listed in sample_generation/problem_name. virtual void SetParamSpaceSizes(); @@ -110,21 +110,27 @@ class SampleGenerator Number of blocks in U_snapshots must be equal to the size of snapshot_basis_tags. The appended column indices of each basis tag are stored in col_idxs. */ - void SaveSnapshot(BlockVector *U_snapshots, std::vector &snapshot_basis_tags, Array &col_idxs); + void SaveSnapshot(BlockVector *U_snapshots, std::vector &snapshot_basis_tags, Array &col_idxs); void SaveSnapshotPorts(TopologyHandler *topol_handler, const Array &col_idxs); - void AddSnapshotGenerator(const int &fom_vdofs, const std::string &prefix, const std::string &basis_tag); + void AddSnapshotGenerator(const int &fom_vdofs, const std::string &prefix, const BasisTag &basis_tag); void WriteSnapshots(); void WriteSnapshotPorts(); - const CAROM::Matrix* LookUpSnapshot(const std::string &basis_tag); + const CAROM::Matrix* LookUpSnapshot(const BasisTag &basis_tag); + Array2D* LookUpSnapshotPortColOffsets(const PortTag &port_tag); void ReportStatus(const int &sample_idx); /* Collect snapshot matrices from the file list to the specified basis tag. */ - void CollectSnapshots(const std::string &basis_prefix, - const std::string &basis_tag, - const std::vector &file_list); + void CollectSnapshotsByBasis(const std::string &basis_prefix, + const BasisTag &basis_tag, + const std::vector &file_list); + /* + Collect snapshot matrices from the file list to the specified port tag file. + */ + void CollectSnapshotsByPort(const std::string &basis_prefix, + const std::string &port_tag_file); /* Perform SVD over snapshot for basis_tag. Calculate the energy fraction for num_basis. @@ -137,6 +143,8 @@ class SampleGenerator // Save all singular value spectrum. Calculate the coverage for ref_num_basis (optional). void SaveSV(CAROM::BasisGenerator *basis_generator, const std::string& prefix, const int &ref_num_basis = -1); + const int GetSnapshotOffset(const std::string &comp); + }; #endif diff --git a/include/steady_ns_solver.hpp b/include/steady_ns_solver.hpp index dc3374a2..a75cf5f3 100644 --- a/include/steady_ns_solver.hpp +++ b/include/steady_ns_solver.hpp @@ -97,15 +97,25 @@ class SteadyNSEQPROM : public SteadyNSROM { protected: Array hs; // not owned by SteadyNSEQPROM. + 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_, const Array &block_offsets_, const bool direct_solve_=true) - : SteadyNSROM(linearOp_, hs_.Size(), block_offsets_, direct_solve_), hs(hs_) {} + SteadyNSEQPROM(SparseMatrix *linearOp_, Array &hs_, ROMInterfaceForm *itf_, + const Array &block_offsets_, const bool direct_solve_=true); virtual ~SteadyNSEQPROM() {} virtual void Mult(const Vector &x, Vector &y) const; virtual Operator &GetGradient(const Vector &x) const; + +private: + void GetVel(const Vector &x, Vector &x_u) const; + void AddVel(const Vector &y_u, Vector &y) const; }; class SteadyNSSolver : public StokesSolver @@ -123,12 +133,14 @@ friend class SteadyNSOperator; double zeta = 1.0; ConstantCoefficient *zeta_coeff = NULL, *minus_zeta = NULL, *minus_half_zeta = NULL; + VectorConstantCoefficient *zero = NULL; // operator for nonlinear convection. Array hs; // integration rule for nonliear operator. const IntegrationRule *ir_nl = NULL; + const IntegrationRule *ir_face = NULL; // interface integrator InterfaceForm *nl_itf = NULL; @@ -137,6 +149,7 @@ friend class SteadyNSOperator; // component ROM element for nonlinear convection. Array comp_tensors, subdomain_tensors; Array comp_eqps, subdomain_eqps; + ROMInterfaceForm *itf_eqp = NULL; Solver *J_solver = NULL; GMRESSolver *J_gmres = NULL; @@ -149,8 +162,6 @@ friend class SteadyNSOperator; using StokesSolver::GetVariableNames; - void InitVariables() override; - void BuildDomainOperators() override; void SetupDomainBCOperators() override; @@ -166,20 +177,27 @@ friend class SteadyNSOperator; void SolveROM() override; - void AllocateROMTensorElems() override; - void BuildROMTensorElems() override; - void SaveROMTensorElems(const std::string &filename) override; - void LoadROMTensorElems(const std::string &filename) override; - void AssembleROMTensorOper() override; + void InitROMHandler() override; - void AllocateROMEQPElems() override; - void TrainEQPElems(SampleGenerator *sample_generator) override; - void SaveEQPElems(const std::string &filename) override; - void LoadEQPElems(const std::string &filename) override; - void AssembleROMEQPOper() override; + virtual void AllocateROMNlinElems() override; + virtual void BuildROMTensorElems() override; + virtual void TrainROMEQPElems(SampleGenerator *sample_generator) override; + virtual void SaveROMNlinElems(const std::string &input_prefix) override; + virtual void LoadROMNlinElems(const std::string &input_prefix) override; + virtual void AssembleROMNlinOper() override; private: DenseTensor* GetReducedTensor(DenseMatrix *basis, FiniteElementSpace *fespace); + + void AllocateROMTensorElems(); + void SaveROMTensorElems(const std::string &filename); + void LoadROMTensorElems(const std::string &filename); + void AssembleROMTensorOper(); + + void AllocateROMEQPElems(); + void SaveEQPElems(const std::string &filename); + void LoadEQPElems(const std::string &filename); + void AssembleROMEQPOper(); }; #endif diff --git a/include/topology_handler.hpp b/include/topology_handler.hpp index 285f1f2d..71ee1238 100644 --- a/include/topology_handler.hpp +++ b/include/topology_handler.hpp @@ -168,8 +168,27 @@ class SubMeshTopologyHandler : public TopologyHandler virtual Mesh* GetMesh(const int k) { return &(*meshes[k]); } virtual Mesh* GetGlobalMesh() { return pmesh; } - // SubMeshTopologyHandler assumes only one component. - virtual Mesh* GetComponentMesh(const int &c) { return GetMesh(0); } + // In SubMeshTopologyHandler, each subdomain is unique component. + virtual Mesh* GetComponentMesh(const int &c) { return GetMesh(c); } + + virtual const int GetPortType(const int &port_idx) { return port_idx; } + virtual const int GetNumRefPorts() { return num_ports; } + // return component indexes for a reference port + virtual void GetComponentPair(const int &ref_port_idx, int &comp1, int &comp2) + { + comp1 = port_infos[ref_port_idx].Mesh1; + comp2 = port_infos[ref_port_idx].Mesh2; + return; + } + virtual void GetRefPortInfo(const int &ref_port_idx, int &comp1, int &comp2, int &attr1, int &attr2) + { + comp1 = port_infos[ref_port_idx].Mesh1; + comp2 = port_infos[ref_port_idx].Mesh2; + attr1 = port_infos[ref_port_idx].Attr1; + attr2 = port_infos[ref_port_idx].Attr2; + } + virtual Array* const GetRefInterfaceInfos(const int &k) + { return GetInterfaceInfos(k); } // Export mesh pointers and interface info. virtual void ExportInfo(Array &mesh_ptrs, TopologyData &topol_data); diff --git a/include/unsteady_ns_solver.hpp b/include/unsteady_ns_solver.hpp index ddfde72d..939f2cc2 100644 --- a/include/unsteady_ns_solver.hpp +++ b/include/unsteady_ns_solver.hpp @@ -65,8 +65,7 @@ friend class SteadyNSOperator; virtual ~UnsteadyNSSolver(); using SteadyNSSolver::GetVariableNames; - - void InitVariables() override; + void BuildDomainOperators() override; void AssembleOperator() override; @@ -88,28 +87,16 @@ friend class SteadyNSOperator; void SolveROM() override { mfem_error("UnsteadyNSSolver::SolveROM is not implemented yet!\n"); } - void AllocateROMEQPElems() override - { mfem_error("UnsteadyNSSolver::AllocateROMEQPElems is not implemented yet!\n"); } - void TrainEQPElems(SampleGenerator *sample_generator) override - { mfem_error("UnsteadyNSSolver::TrainEQPElems is not implemented yet!\n"); } - void SaveEQPElems(const std::string &filename) override - { mfem_error("UnsteadyNSSolver::SaveEQPElems is not implemented yet!\n"); } - void LoadEQPElems(const std::string &filename) override - { mfem_error("UnsteadyNSSolver::LoadEQPElems is not implemented yet!\n"); } - void AssembleROMEQPOper() override - { mfem_error("UnsteadyNSSolver::AssembleROMEQPOper is not implemented yet!\n"); } - - /* Tensorial ROM is not supported for unsteady NS */ - void AllocateROMTensorElems() override - { mfem_error("UnsteadyNSSolver::AllocateROMTensorElems is not implemented yet!\n"); } void BuildROMTensorElems() override { mfem_error("UnsteadyNSSolver::BuildROMTensorElems is not implemented yet!\n"); } - void SaveROMTensorElems(const std::string &filename) override - { mfem_error("UnsteadyNSSolver::SaveROMTensorElems is not implemented yet!\n"); } - void LoadROMTensorElems(const std::string &filename) override - { mfem_error("UnsteadyNSSolver::LoadROMTensorElems is not implemented yet!\n"); } - void AssembleROMTensorOper() override - { mfem_error("UnsteadyNSSolver::AssembleROMTensorOper 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(); diff --git a/sketches/ns_dg_mms.cpp b/sketches/ns_dg_mms.cpp index dd901acf..cd30cd1b 100644 --- a/sketches/ns_dg_mms.cpp +++ b/sketches/ns_dg_mms.cpp @@ -413,6 +413,8 @@ int main(int argc, char *argv[]) if (use_dg) nVarf->AddInteriorFaceIntegrator(lf_integ2); nVarf->AddBdrFaceIntegrator(lf_bdr_integ1, u_ess_attr); + /* having linear-only term causes sub-optimal convergence */ + // fform->AddBdrFaceIntegrator(new DGBdrLaxFriedrichsLFIntegrator(ucoeff, &zeta_coeff), u_ess_attr); } break; case Scheme::TEMAM: diff --git a/sketches/ns_rom.cpp b/sketches/ns_rom.cpp index f51a70bf..5841a2c0 100644 --- a/sketches/ns_rom.cpp +++ b/sketches/ns_rom.cpp @@ -1435,18 +1435,16 @@ int main(int argc, char *argv[]) rom_nlinf->SetBasis(u_basis); rom_nlinf->TrainEQP(*snapshots, eqp_tol); - Array sample_el, sample_qp; - Array sample_qw; - rom_nlinf->GetEQPForIntegrator(IntegratorType::DOMAIN, 0, sample_el, sample_qp, sample_qw); + Array *samples = rom_nlinf->GetEQPForIntegrator(IntegratorType::DOMAIN, 0); { // save empirical quadrature point locations. ElementTransformation *T; Vector loc(dim); - DenseMatrix qp_loc(sample_el.Size(), dim); - for (int p = 0; p < sample_el.Size(); p++) + DenseMatrix qp_loc(samples->Size(), dim); + for (int p = 0; p < samples->Size(); p++) { - T = ufes->GetElementTransformation(sample_el[p]); - const IntegrationPoint &ip = ir->IntPoint(sample_qp[p]); + T = ufes->GetElementTransformation((*samples)[p].el); + const IntegrationPoint &ip = ir->IntPoint((*samples)[p].qp); T->Transform(ip, loc); for (int d = 0; d < dim; d++) diff --git a/src/advdiff_solver.cpp b/src/advdiff_solver.cpp index f4c21dd3..9b16800c 100644 --- a/src/advdiff_solver.cpp +++ b/src/advdiff_solver.cpp @@ -242,6 +242,7 @@ void AdvDiffSolver::GetFlowField(ParameterizedProblem *flow_problem) stokes_solver = new StokesSolver; stokes_solver->InitVariables(); + if (use_rom) stokes_solver->InitROMHandler(); stokes_solver->SetSolutionSaveMode(save_flow); bool flow_loaded = false; diff --git a/src/component_topology_handler.cpp b/src/component_topology_handler.cpp index 660af315..25a2a863 100644 --- a/src/component_topology_handler.cpp +++ b/src/component_topology_handler.cpp @@ -12,7 +12,7 @@ using namespace std; using namespace mfem; ComponentTopologyHandler::ComponentTopologyHandler() - : TopologyHandler(COMPONENT) + : TopologyHandler(TopologyHandlerMode::COMPONENT) { verbose = config.GetOption("mesh/component-wise/verbose", false); write_ports = config.GetOption("mesh/component-wise/write_ports", false); @@ -732,6 +732,10 @@ void ComponentTopologyHandler::SetupPorts() assert(port_infos[p].Attr1 == ref_port->Attr1); assert(port_infos[p].Attr2 == ref_port->Attr2); + if ((ref_port->Component1 == ref_port->Component2) && + (port_infos[p].Mesh1 == port_infos[p].Mesh2)) + mfem_error("Currently interface to a mesh itself is not supported!\n"); + interface_infos[p] = ref_interfaces[port_types[p]]; } diff --git a/src/dg_linear.cpp b/src/dg_linear.cpp index 7f987598..afafd481 100644 --- a/src/dg_linear.cpp +++ b/src/dg_linear.cpp @@ -474,7 +474,7 @@ void DGBdrTemamLFIntegrator::AssembleRHSElementVect( } w = ip.weight; - if (Z) w *= Z->Eval(*(Tr.Face), eip); + if (Z) w *= Z->Eval(Tr, ip); Q.Eval(Qvec, *Tr.Elem1, eip); diff --git a/src/hdf5_utils.cpp b/src/hdf5_utils.cpp index 29cfdc52..11092412 100644 --- a/src/hdf5_utils.cpp +++ b/src/hdf5_utils.cpp @@ -83,6 +83,18 @@ void WriteAttribute(hid_t &dest, std::string attribute, const std::string &value H5Tclose(attrType); } +void ReadAttribute(hid_t &source, std::string attribute, BasisTag &value) +{ + ReadAttribute(source, attribute + "_comp", value.comp); + ReadAttribute(source, attribute + "_var", value.var); +} + +void WriteAttribute(hid_t &source, std::string attribute, const BasisTag &value) +{ + WriteAttribute(source, attribute + "_comp", value.comp); + WriteAttribute(source, attribute + "_var", value.var); +} + void ReadDataset(hid_t &source, std::string dataset, DenseMatrix &value) { herr_t errf = 0; @@ -479,4 +491,89 @@ void WriteDataset(hid_t &source, std::string dataset, const MatrixBlocks &value) assert(errf >= 0); } +void ReadDataset(hid_t &source, std::string dataset, const IntegratorType type, Array &value) +{ + std::string eldset; + switch (type) + { + case IntegratorType::DOMAIN: eldset = "elem"; break; + case IntegratorType::INTERIORFACE: eldset = "face"; break; + case IntegratorType::BDRFACE: eldset = "be"; break; + case IntegratorType::INTERFACE: eldset = "itf"; break; + default: + mfem_error("hdf5_utils::ReadDataset- Unknown IntegratorType!\n"); + } + + Array el, qp; + Array qw; + + assert(source >= 0); + hid_t grp_id; + herr_t errf; + + grp_id = H5Gopen2(source, dataset.c_str(), H5P_DEFAULT); + assert(grp_id >= 0); + + hdf5_utils::ReadDataset(grp_id, eldset, el); + hdf5_utils::ReadDataset(grp_id, "quad-pt", qp); + hdf5_utils::ReadDataset(grp_id, "quad-wt", qw); + + errf = H5Gclose(grp_id); + assert(errf >= 0); + + value.SetSize(el.Size()); + for (int k = 0; k < el.Size(); k++) + { + value[k].qp = qp[k]; + value[k].qw = qw[k]; + value[k].el = el[k]; + } + + return; +} + +void WriteDataset(hid_t &source, std::string dataset, const IntegratorType type, const Array &value) +{ + std::string eldset; + switch (type) + { + case IntegratorType::DOMAIN: eldset = "elem"; break; + case IntegratorType::INTERIORFACE: eldset = "face"; break; + case IntegratorType::BDRFACE: eldset = "be"; break; + case IntegratorType::INTERFACE: eldset = "itf"; break; + default: + mfem_error("hdf5_utils::WriteDataset- Unknown IntegratorType!\n"); + } + + Array el, qp; + Array qw; + + el.SetSize(0); + qp.SetSize(0); + qw.SetSize(0); + + for (int s = 0; s < value.Size(); s++) + { + el.Append(value[s].el); + qp.Append(value[s].qp); + qw.Append(value[s].qw); + } + + assert(source >= 0); + hid_t grp_id; + herr_t errf; + + grp_id = H5Gcreate(source, dataset.c_str(), H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT); + assert(grp_id >= 0); + + hdf5_utils::WriteDataset(grp_id, eldset, el); + hdf5_utils::WriteDataset(grp_id, "quad-pt", qp); + hdf5_utils::WriteDataset(grp_id, "quad-wt", qw); + + errf = H5Gclose(grp_id); + assert(errf >= 0); + + return; +} + } diff --git a/src/interfaceinteg.cpp b/src/interfaceinteg.cpp index e793760f..89e28048 100644 --- a/src/interfaceinteg.cpp +++ b/src/interfaceinteg.cpp @@ -1712,7 +1712,7 @@ void DGLaxFriedrichsFluxIntegrator::AssembleQuadratureGrad( void DGLaxFriedrichsFluxIntegrator::AppendPrecomputeInteriorFaceCoeffs( const FiniteElementSpace *fes, DenseMatrix &basis, const SampleInfo &sample) { - const int face = sample.face; + const int face = sample.el; FaceElementTransformations *T = fes->GetMesh()->GetInteriorFaceTransformations(face); assert(T != NULL); assert(T->Elem2No >= 0); @@ -1723,7 +1723,7 @@ void DGLaxFriedrichsFluxIntegrator::AppendPrecomputeInteriorFaceCoeffs( void DGLaxFriedrichsFluxIntegrator::AppendPrecomputeBdrFaceCoeffs( const FiniteElementSpace *fes, DenseMatrix &basis, const SampleInfo &sample) { - const int be = sample.be; + const int be = sample.el; FaceElementTransformations *T = fes->GetMesh()->GetBdrFaceTransformations(be); assert(T != NULL); assert(T->Elem2No < 0); diff --git a/src/linalg_utils.cpp b/src/linalg_utils.cpp index 123ee196..ea40d1fd 100644 --- a/src/linalg_utils.cpp +++ b/src/linalg_utils.cpp @@ -137,8 +137,8 @@ namespace mfem void AddToBlockMatrix(const Array &ridx, const Array &cidx, const MatrixBlocks &mats, BlockMatrix &bmat) { assert(bmat.owns_blocks); - assert(bmat.NumRowBlocks() >= mats.nrows); - assert(bmat.NumColBlocks() >= mats.ncols); + // assert(bmat.NumRowBlocks() >= mats.nrows); + // assert(bmat.NumColBlocks() >= mats.ncols); assert(ridx.Size() == mats.nrows); assert(cidx.Size() == mats.ncols); assert((ridx.Min() >= 0) && (ridx.Max() < bmat.NumRowBlocks())); diff --git a/src/linelast_solver.cpp b/src/linelast_solver.cpp index 70e89e74..b78b0f4c 100644 --- a/src/linelast_solver.cpp +++ b/src/linelast_solver.cpp @@ -133,8 +133,6 @@ void LinElastSolver::InitVariables() // Does this make any difference? us[m]->SetTrueVector(); } - if (use_rom) - MultiBlockSolver::InitROMHandler(); } void LinElastSolver::BuildOperators() diff --git a/src/main_workflow.cpp b/src/main_workflow.cpp index 0d71af25..5abb72e0 100644 --- a/src/main_workflow.cpp +++ b/src/main_workflow.cpp @@ -90,9 +90,9 @@ SampleGenerator* InitSampleGenerator(MPI_Comm comm) return generator; } -std::vector GetGlobalBasisTagList(const TopologyHandlerMode &topol_mode, bool separate_variable_basis) +std::vector GetGlobalBasisTagList(const TopologyHandlerMode &topol_mode, bool separate_variable_basis) { - std::vector basis_tags(0); + std::vector basis_tags(0); std::vector component_list(0); if (topol_mode == TopologyHandlerMode::SUBMESH) @@ -132,9 +132,9 @@ std::vector GetGlobalBasisTagList(const TopologyHandlerMode &topol_ { if (separate_variable_basis) for (int v = 0; v < var_list.size(); v++) - basis_tags.push_back(component_list[c] + "_" + var_list[v]); + basis_tags.push_back(BasisTag(component_list[c], var_list[v])); else - basis_tags.push_back(component_list[c]); + basis_tags.push_back(BasisTag(component_list[c])); } return basis_tags; @@ -160,6 +160,8 @@ void GenerateSamples(MPI_Comm comm) test = InitSolver(); test->InitVariables(); + if (test->UseRom()) + test->InitROMHandler(); problem->SetSingleRun(); test->SetParameterizedProblem(problem); @@ -207,6 +209,40 @@ void GenerateSamples(MPI_Comm comm) } void CollectSamples(SampleGenerator *sample_generator) +{ + std::string mode = config.GetOption("sample_collection/mode", "basis"); + std::string basis_prefix = config.GetOption("basis/prefix", "basis"); + + if (mode == "basis") + CollectSamplesByBasis(sample_generator, basis_prefix); + else if (mode == "port") + CollectSamplesByPort(sample_generator, basis_prefix); + else + mfem_error("CollectSamples: unknown sample collection mode!\n"); +} + +void CollectSamplesByPort(SampleGenerator *sample_generator, const std::string &basis_prefix) +{ + // parse the sample snapshot file list. + std::vector file_list = config.GetOption>( + "sample_collection/port_files", std::vector(0)); + YAML::Node port_format = config.FindNode("sample_collection/port_fileformat"); + // if file list is specified with a format, parse through the format. + if (port_format) + { + FilenameParam port_param("", port_format); + port_param.ParseFilenames(file_list); + } + + // if additional inputs are not specified for port files, set default port file name. + if (file_list.size() == 0) + file_list.push_back(sample_generator->GetSamplePrefix() + ".port.h5"); + + for (int f = 0; f < file_list.size(); f++) + sample_generator->CollectSnapshotsByPort(basis_prefix, file_list[f]); +} + +void CollectSamplesByBasis(SampleGenerator *sample_generator, const std::string &basis_prefix) { assert(sample_generator); @@ -214,13 +250,11 @@ void CollectSamples(SampleGenerator *sample_generator) bool separate_variable_basis = config.GetOption("model_reduction/separate_variable_basis", false); // Find the all required basis tags. - std::vector basis_tags = GetGlobalBasisTagList(topol_mode, separate_variable_basis); + std::vector basis_tags = GetGlobalBasisTagList(topol_mode, separate_variable_basis); // tag-specific optional inputs. YAML::Node basis_list = config.FindNode("basis/tags"); - std::string basis_prefix = config.GetOption("basis/prefix", "basis"); - // loop over the required basis tag list. for (int p = 0; p < basis_tags.size(); p++) { @@ -230,7 +264,7 @@ void CollectSamples(SampleGenerator *sample_generator) FindSnapshotFilesForBasis(basis_tags[p], default_filename, file_list); assert(file_list.size() > 0); - sample_generator->CollectSnapshots(basis_prefix, basis_tags[p], file_list); + sample_generator->CollectSnapshotsByBasis(basis_prefix, basis_tags[p], file_list); } // for (int p = 0; p < basis_tags.size(); p++) } @@ -264,6 +298,7 @@ void AuxiliaryTrainROM(MPI_Comm comm, SampleGenerator *sample_generator) if (!solver->UseRom()) mfem_error("ROM must be enabled for supremizer enrichment!\n"); solver->InitVariables(); + solver->InitROMHandler(); // This time needs to be ROMHandler, in order not to run StokesSolver::LoadSupremizer. solver->GetROMHandler()->LoadReducedBasis(); @@ -289,6 +324,7 @@ void TrainEQP(MPI_Comm comm) MultiBlockSolver *test = NULL; test = InitSolver(); test->InitVariables(); + test->InitROMHandler(); if (!test->IsNonlinear()) { @@ -300,6 +336,7 @@ void TrainEQP(MPI_Comm comm) if (!test->UseRom()) mfem_error("ROM must be enabled for EQP training!\n"); test->LoadReducedBasis(); + test->AllocateROMNlinElems(); ROMHandlerBase *rom = test->GetROMHandler(); ROMBuildingLevel save_operator = rom->GetBuildingLevel(); @@ -312,7 +349,7 @@ void TrainEQP(MPI_Comm comm) else mfem_error("Unknown TopologyHandler Mode!\n"); - std::string filename = rom->GetOperatorPrefix() + ".eqp.h5"; + std::string oper_prefix = rom->GetOperatorPrefix(); switch (save_operator) { case ROMBuildingLevel::COMPONENT: @@ -320,9 +357,8 @@ void TrainEQP(MPI_Comm comm) if (topol_mode == TopologyHandlerMode::SUBMESH) mfem_error("Submesh does not support component rom building level!\n"); - test->AllocateROMEQPElems(); - test->TrainEQPElems(sample_generator); - test->SaveEQPElems(filename); + test->TrainROMEQPElems(sample_generator); + test->SaveROMNlinElems(oper_prefix); break; } case ROMBuildingLevel::GLOBAL: @@ -341,7 +377,7 @@ void TrainEQP(MPI_Comm comm) delete test; } -void FindSnapshotFilesForBasis(const std::string &basis_tag, const std::string &default_filename, std::vector &file_list) +void FindSnapshotFilesForBasis(const BasisTag &basis_tag, const std::string &default_filename, std::vector &file_list) { file_list.clear(); @@ -352,7 +388,7 @@ void FindSnapshotFilesForBasis(const std::string &basis_tag, const std::string & if (basis_list) { // Find if additional inputs are specified for basis_tag. - YAML::Node basis_tag_input = config.LookUpFromDict("name", basis_tag, basis_list); + YAML::Node basis_tag_input = config.LookUpFromDict("name", basis_tag.print(), basis_list); // If basis_tag has additional inputs, parse them. if (basis_tag_input) @@ -383,6 +419,7 @@ void BuildROM(MPI_Comm comm) test = InitSolver(); if (!test->UseRom()) mfem_error("ROM must be enabled for BuildROM!\n"); test->InitVariables(); + test->InitROMHandler(); // test->InitVisualization(); // The ROM operator will be built based on the parameter specified for single-run. @@ -393,6 +430,9 @@ void BuildROM(MPI_Comm comm) test->BuildOperators(); test->SetupBCOperators(); test->LoadReducedBasis(); + + if (test->IsNonlinear()) + test->AllocateROMNlinElems(); TopologyHandlerMode topol_mode = test->GetTopologyMode(); ROMHandlerBase *rom = test->GetROMHandler(); @@ -402,7 +442,7 @@ void BuildROM(MPI_Comm comm) if (save_operator == ROMBuildingLevel::GLOBAL) test->Assemble(); - std::string filename = rom->GetOperatorPrefix() + ".h5"; + std::string oper_prefix = rom->GetOperatorPrefix(); switch (save_operator) { case ROMBuildingLevel::COMPONENT: @@ -411,20 +451,19 @@ void BuildROM(MPI_Comm comm) mfem_error("Submesh does not support component rom building level!\n"); test->BuildROMLinElems(); - test->SaveROMLinElems(filename); + test->SaveROMLinElems(oper_prefix + ".h5"); if ((test->IsNonlinear()) && (rom->GetNonlinearHandling() == NonlinearHandling::TENSOR)) { - test->AllocateROMTensorElems(); test->BuildROMTensorElems(); - test->SaveROMTensorElems(filename); + test->SaveROMNlinElems(oper_prefix); } break; } case ROMBuildingLevel::GLOBAL: { test->ProjectOperatorOnReducedBasis(); - test->SaveROMOperator(filename); + test->SaveROMOperator(oper_prefix + ".h5"); break; } case ROMBuildingLevel::NONE: @@ -455,6 +494,7 @@ double SingleRun(MPI_Comm comm, const std::string output_file) ParameterizedProblem *problem = InitParameterizedProblem(); MultiBlockSolver *test = InitSolver(); test->InitVariables(); + if (test->UseRom()) test->InitROMHandler(); test->InitVisualization(); StopWatch solveTimer; @@ -479,6 +519,9 @@ double SingleRun(MPI_Comm comm, const std::string output_file) { rom = test->GetROMHandler(); test->LoadReducedBasis(); + + if (test->IsNonlinear()) + test->AllocateROMNlinElems(); } solveTimer.Start(); @@ -511,7 +554,6 @@ double SingleRun(MPI_Comm comm, const std::string output_file) if (test->IsNonlinear()) { - test->AllocateROMNlinElems(); test->LoadROMNlinElems(rom->GetOperatorPrefix()); test->AssembleROMNlinOper(); } diff --git a/src/multiblock_solver.cpp b/src/multiblock_solver.cpp index fdd24bc3..acb387d8 100644 --- a/src/multiblock_solver.cpp +++ b/src/multiblock_solver.cpp @@ -383,7 +383,7 @@ void MultiBlockSolver::InitIndividualParaview(const std::string& file_prefix) { assert(var_names.size() == num_var); assert((visual.domain_offset >= 0) && (visual.domain_offset < numSub)); - assert((visual.domain_interval > 0) && (visual.domain_interval < numSub)); + assert((visual.domain_interval > 0) && (visual.domain_interval <= numSub)); paraviewColls.SetSize(numSub); paraviewColls = NULL; @@ -485,7 +485,7 @@ void MultiBlockSolver::SaveVisualization() else { assert((visual.domain_offset >= 0) && (visual.domain_offset < numSub)); - assert((visual.domain_interval > 0) && (visual.domain_interval < numSub)); + assert((visual.domain_interval > 0) && (visual.domain_interval <= numSub)); } for (int m = 0; m < paraviewColls.Size(); m++) @@ -622,36 +622,6 @@ void MultiBlockSolver::CopySolution(BlockVector *input_sol) *U = *input_sol; } -void MultiBlockSolver::AllocateROMNlinElems() -{ - assert(nonlinear_mode); - assert(rom_handler); - NonlinearHandling rom_nlin_mode = rom_handler->GetNonlinearHandling(); - - if (rom_nlin_mode == NonlinearHandling::TENSOR) AllocateROMTensorElems(); - else if (rom_nlin_mode == NonlinearHandling::EQP) AllocateROMEQPElems(); -} - -void MultiBlockSolver::LoadROMNlinElems(const std::string &input_prefix) -{ - assert(nonlinear_mode); - assert(rom_handler); - NonlinearHandling rom_nlin_mode = rom_handler->GetNonlinearHandling(); - - if (rom_nlin_mode == NonlinearHandling::TENSOR) LoadROMTensorElems(input_prefix + ".h5"); - else if (rom_nlin_mode == NonlinearHandling::EQP) LoadEQPElems(input_prefix + ".eqp.h5"); -} - -void MultiBlockSolver::AssembleROMNlinOper() -{ - assert(nonlinear_mode); - assert(rom_handler); - NonlinearHandling rom_nlin_mode = rom_handler->GetNonlinearHandling(); - - if (rom_nlin_mode == NonlinearHandling::TENSOR) AssembleROMTensorOper(); - else if (rom_nlin_mode == NonlinearHandling::EQP) AssembleROMEQPOper(); -} - void MultiBlockSolver::InitROMHandler() { rom_handler = new MFEMROMHandler(topol_handler, var_offsets, var_names, separate_variable_basis); @@ -663,7 +633,7 @@ void MultiBlockSolver::InitROMHandler() rom_elems = new ROMLinearElement(topol_handler, comp_fes, separate_variable_basis); } -void MultiBlockSolver::GetBasisTags(std::vector &basis_tags) +void MultiBlockSolver::GetBasisTags(std::vector &basis_tags) { if (separate_variable_basis) { @@ -680,7 +650,7 @@ void MultiBlockSolver::GetBasisTags(std::vector &basis_tags) } } -BlockVector* MultiBlockSolver::PrepareSnapshots(std::vector &basis_tags) +BlockVector* MultiBlockSolver::PrepareSnapshots(std::vector &basis_tags) { BlockVector *U_snapshots = NULL; @@ -701,7 +671,7 @@ void MultiBlockSolver::SaveSnapshots(SampleGenerator *sample_generator) assert(sample_generator); /* split the solution into each component with the corresponding tag */ - std::vector basis_tags; + std::vector basis_tags; BlockVector *U_snapshots = PrepareSnapshots(basis_tags); Array col_idxs; diff --git a/src/parameterized_problem.cpp b/src/parameterized_problem.cpp index dbaf86a6..10699fc0 100644 --- a/src/parameterized_problem.cpp +++ b/src/parameterized_problem.cpp @@ -197,6 +197,22 @@ void pic(const Vector &x, double t, Vector &y) } // namespace backward_facing_step +namespace lid_driven_cavity +{ + +double u0, L; + +void ubdr(const Vector &x, double t, Vector &y) +{ + const int dim = x.Size(); + y.SetSize(dim); + y = 0.0; + + y(0) = 4.0 / L / L * u0 * x(0) * (L - x(0)); +} + +} // namespace lid_driven_cavity + namespace periodic_flow_past_array { @@ -497,6 +513,14 @@ ParameterizedProblem* InitParameterizedProblem() { problem = new BackwardFacingStep(); } + else if (problem_name == "lid_driven_cavity") + { + problem = new LidDrivenCavity(); + } + else if (problem_name == "force_driven_corner") + { + problem = new ForceDrivenCorner(); + } else if (problem_name == "periodic_flow_past_array") { problem = new PeriodicFlowPastArray(); @@ -677,7 +701,7 @@ ChannelFlow::ChannelFlow() battr.SetSize(5); for (int b = 0; b < 5; b++) battr[b] = b+1; - bdr_type.SetSize(4); + bdr_type.SetSize(5); bdr_type = BoundaryType::ZERO; bdr_type[1] = BoundaryType::NEUMANN; bdr_type[3] = BoundaryType::DIRICHLET; @@ -884,6 +908,49 @@ BackwardFacingStep::BackwardFacingStep() } } +/* + LidDrivenCavity +*/ + +LidDrivenCavity::LidDrivenCavity() + : FlowProblem() +{ + /* Assume there are only five boundary attributes: 1, 2, 3, 4, 5 */ + battr.SetSize(5); + for (int b = 0; b < 5; b++) + battr[b] = b+1; + + /* + All boundaries are zero except boundary attribute 3 + */ + bdr_type.SetSize(5); + bdr_type = BoundaryType::ZERO; + bdr_type[2] = BoundaryType::DIRICHLET; + bdr_type[0] = BoundaryType::DIRICHLET; + + // pointer to static function. + vector_bdr_ptr.SetSize(5); + vector_rhs_ptr = NULL; + /* technically only vector_bdr_ptr[2] will be used. */ + vector_bdr_ptr = &(function_factory::flow_problem::lid_driven_cavity::ubdr); + + param_num = 3; + + // Default values. + function_factory::flow_problem::nu = 1.0; + function_factory::flow_problem::lid_driven_cavity::u0 = 1.0; + function_factory::flow_problem::lid_driven_cavity::L = 1.0; + + param_map["nu"] = 0; + param_map["u0"] = 1; + param_map["L"] = 2; + + param_ptr.SetSize(param_num); + param_ptr[0] = &(function_factory::flow_problem::nu); + param_ptr[1] = &(function_factory::flow_problem::lid_driven_cavity::u0); + param_ptr[2] = &(function_factory::flow_problem::lid_driven_cavity::L); +} + /* PeriodicFlowPastArray */ @@ -920,6 +987,26 @@ PeriodicFlowPastArray::PeriodicFlowPastArray() } +/* + ForceDrivenCorner +*/ + +ForceDrivenCorner::ForceDrivenCorner() + : PeriodicFlowPastArray() +{ + battr.SetSize(5); + for (int b = 0; b < 5; b++) + battr[b] = b+1; + bdr_type.SetSize(5); + bdr_type = BoundaryType::ZERO; + bdr_type[2] = BoundaryType::NEUMANN; + bdr_type[3] = BoundaryType::NEUMANN; + + // pointer to static function. + vector_bdr_ptr.SetSize(5); + vector_bdr_ptr = NULL; +} + /* LinElastDisp */ diff --git a/src/poisson_solver.cpp b/src/poisson_solver.cpp index 7c4284fc..948ed837 100644 --- a/src/poisson_solver.cpp +++ b/src/poisson_solver.cpp @@ -152,8 +152,6 @@ void PoissonSolver::InitVariables() } rhs_coeffs.SetSize(0); - - if (use_rom) MultiBlockSolver::InitROMHandler(); } void PoissonSolver::BuildOperators() diff --git a/src/rom_handler.cpp b/src/rom_handler.cpp index e6311c08..3c400994 100644 --- a/src/rom_handler.cpp +++ b/src/rom_handler.cpp @@ -16,16 +16,17 @@ using namespace std; namespace mfem { -const std::string GetBasisTagForComponent( +const BasisTag GetBasisTagForComponent( const int &comp_idx, const TopologyHandler *topol_handler, const std::string var_name) { - std::string tag = topol_handler->GetComponentName(comp_idx); + BasisTag tag; + tag.comp = topol_handler->GetComponentName(comp_idx); if (var_name != "") - tag += "_" + var_name; + tag.var = var_name; return tag; } -const std::string GetBasisTag( +const BasisTag GetBasisTag( const int &subdomain_index, const TopologyHandler *topol_handler, const std::string var_name) { int c_type = topol_handler->GetMeshType(subdomain_index); @@ -105,18 +106,18 @@ void ROMHandlerBase::ParseInputs() for (int b = 0; b < num_rom_ref_blocks; b++) { /* determine basis tag */ - std::string basis_tag; + BasisTag basis_tag; if (separate_variable) basis_tag = GetBasisTagForComponent(b / num_var, topol_handler, fom_var_names[b % num_var]); else basis_tag = GetBasisTagForComponent(b, topol_handler); /* determine number of basis */ - num_ref_basis[b] = config.LookUpFromDict("name", basis_tag, "number_of_basis", num_ref_basis_default, basis_list); + num_ref_basis[b] = config.LookUpFromDict("name", basis_tag.print(), "number_of_basis", num_ref_basis_default, basis_list); if (num_ref_basis[b] < 0) { - printf("Cannot find the number of basis for %s!\n", basis_tag.c_str()); + printf("Cannot find the number of basis for %s!\n", basis_tag.print().c_str()); mfem_error("Or specify the default number of basis!\n"); } @@ -167,12 +168,12 @@ void ROMHandlerBase::ParseSupremizerInput(Array &num_ref_supreme, ArraygetSpatialBasis(num_ref_basis[k])); carom_ref_basis[k]->gather(); @@ -572,7 +573,10 @@ void MFEMROMHandler::NonlinearSolve(Operator &oper, BlockVector* U, Solver *prec if (use_restart) ProjectGlobalToDomainBasis(U, reduced_sol); else - (*reduced_sol) = 0.0; + { + for (int k = 0; k < reduced_sol->Size(); k++) + (*reduced_sol)(k) = 1.0e-1 * UniformRandom(); + } int maxIter = config.GetOption("solver/max_iter", 100); double rtol = config.GetOption("solver/relative_tolerance", 1.e-10); diff --git a/src/rom_interfaceform.cpp b/src/rom_interfaceform.cpp index 876c3fb5..30f024e0 100644 --- a/src/rom_interfaceform.cpp +++ b/src/rom_interfaceform.cpp @@ -13,9 +13,14 @@ namespace mfem { ROMInterfaceForm::ROMInterfaceForm( - Array &meshes_, Array &fes_, TopologyHandler *topol_) - : InterfaceForm(meshes_, fes_, topol_), numPorts(topol_->GetNumPorts()) + Array &meshes_, Array &fes_, Array &comp_fes_, TopologyHandler *topol_) + : InterfaceForm(meshes_, fes_, topol_), numPorts(topol_->GetNumPorts()), numRefPorts(topol_->GetNumRefPorts()), + num_comp(topol_->GetNumComponents()), comp_fes(comp_fes_) { + comp_basis.SetSize(num_comp); + comp_basis = NULL; + comp_basis_dof_offsets.SetSize(num_comp); + basis.SetSize(numSub); basis = NULL; basis_dof_offsets.SetSize(numSub); @@ -26,17 +31,26 @@ ROMInterfaceForm::ROMInterfaceForm( ROMInterfaceForm::~ROMInterfaceForm() { - DeletePointers(fnfi_sample); + DeletePointers(fnfi_ref_sample); } -void ROMInterfaceForm::SetBasisAtSubdomain(const int m, DenseMatrix &basis_, const int offset) +void ROMInterfaceForm::SetBasisAtComponent(const int c, DenseMatrix &basis_, const int offset) { // assert(basis_.NumCols() == height); - assert(basis_.NumRows() >= fes[m]->GetTrueVSize() + offset); - assert(basis.Size() == numSub); + assert(basis_.NumRows() >= comp_fes[c]->GetTrueVSize() + offset); + assert(comp_basis.Size() == num_comp); + + comp_basis[c] = &basis_; + comp_basis_dof_offsets[c] = offset; - basis[m] = &basis_; - basis_dof_offsets[m] = offset; + for (int m = 0; m < numSub; m++) + { + if (topol_handler->GetMeshType(m) != c) + continue; + + basis[m] = &basis_; + basis_dof_offsets[m] = offset; + } } void ROMInterfaceForm::UpdateBlockOffsets() @@ -53,7 +67,7 @@ void ROMInterfaceForm::UpdateBlockOffsets() } void ROMInterfaceForm::InterfaceAddMult(const Vector &x, Vector &y) const -{ +{ assert(block_offsets.Min() >= 0); assert(x.Size() == block_offsets.Last()); assert(y.Size() == block_offsets.Last()); @@ -118,7 +132,7 @@ void ROMInterfaceForm::InterfaceAddMult(const Vector &x, Vector &y) const SampleInfo *sample = sample_info->GetData(); for (int i = 0; i < sample_info->Size(); i++, sample++) { - int itf = sample->itf; + int itf = sample->el; InterfaceInfo *if_info = &((*interface_infos)[itf]); topol_handler->GetInterfaceTransformations(mesh1, mesh2, if_info, tr1, tr2); const IntegrationPoint &ip = ir->IntPoint(sample->qp); @@ -241,7 +255,7 @@ void ROMInterfaceForm::InterfaceGetGradient(const Vector &x, Array2DGetData(); for (int i = 0; i < sample_info->Size(); i++, sample++) { - int itf = sample->itf; + int itf = sample->el; InterfaceInfo *if_info = &((*interface_infos)[itf]); topol_handler->GetInterfaceTransformations(mesh1, mesh2, if_info, tr1, tr2); const IntegrationPoint &ip = ir->IntPoint(sample->qp); @@ -291,9 +305,63 @@ void ROMInterfaceForm::InterfaceGetGradient(const Vector &x, Array2D &snap_pair_idx, const double eqp_tol) +{ + int c1, c2, a1, a2; + // TODO(kevin): at least component topology handler maintain attrs the same for both reference and subdomain. + // Need to check submesh topology handler. + topol_handler->GetRefPortInfo(p, c1, c2, a1, a2); + Array *itf_info = topol_handler->GetRefInterfaceInfos(p); + + FiniteElementSpace *fes1 = comp_fes[c1]; + FiniteElementSpace *fes2 = comp_fes[c2]; + assert(fes1 && fes2); + + /* + TODO(kevin): this is a boilerplate for parallel POD/EQP training. + Full parallelization will have to consider local matrix construction from local snapshot/basis matrix. + */ + assert(snapshot1.distributed()); + assert(snapshot2.distributed()); + CAROM::Matrix snapshot1_work(snapshot1), snapshot2_work(snapshot2); + snapshot1_work.gather(); + snapshot2_work.gather(); + assert(snapshot1_work.numRows() >= fes1->GetTrueVSize()); + assert(snapshot2_work.numRows() >= fes2->GetTrueVSize()); + + DenseMatrix *basis1 = comp_basis[c1]; + DenseMatrix *basis2 = comp_basis[c2]; + const int basis1_offset = comp_basis_dof_offsets[c1]; + const int basis2_offset = comp_basis_dof_offsets[c2]; + assert(basis1 && basis2); + + // NOTE(kevin): these will be resized within the routines as needed. + // just initializing with distribute option. + CAROM::Matrix Gt(1,1, true); + CAROM::Vector rhs_Gw(1, false); + + Array samples; + Array fidxs; + + for (int it = 0; it < fnfi.Size(); it++) + { + const IntegrationRule *ir = fnfi[it]->GetIntegrationRule(); + assert(ir); + const int nqe = ir->GetNPoints(); + + SetupEQPSystem(snapshot1_work, snapshot2_work, snap_pair_idx, + *basis1, *basis2, basis1_offset, basis2_offset, + fes1, fes2, itf_info, fnfi[it], Gt, rhs_Gw); + TrainEQPForIntegrator(nqe, Gt, rhs_Gw, eqp_tol, samples); + UpdateInterFaceIntegratorSampling(it, p, samples); + } +} + void ROMInterfaceForm::SetupEQPSystem( const CAROM::Matrix &snapshot1, const CAROM::Matrix &snapshot2, - const Array &snap1_idx, const Array &snap2_idx, + const Array2D &snap_pair_idx, DenseMatrix &basis1, DenseMatrix &basis2, const int &basis1_offset, const int &basis2_offset, FiniteElementSpace *fes1, FiniteElementSpace *fes2, @@ -335,12 +403,12 @@ void ROMInterfaceForm::SetupEQPSystem( We take the snapshot pairs given by snap1/2_idx. These two arrays should have the same size. */ - const int nsnap = snap1_idx.Size(); - assert(nsnap == snap2_idx.Size()); + const int nsnap = snap_pair_idx.NumRows(); + assert(snap_pair_idx.NumCols() == 2); const int nsnap1 = snapshot1.numColumns(); const int nsnap2 = snapshot2.numColumns(); - assert((snap1_idx.Min() >= 0) && (snap1_idx.Max() < nsnap1)); - assert((snap2_idx.Min() >= 0) && (snap2_idx.Max() < nsnap2)); + // assert((snap1_idx.Min() >= 0) && (snap1_idx.Max() < nsnap1)); + // assert((snap2_idx.Min() >= 0) && (snap2_idx.Max() < nsnap2)); const int ne_global = itf_infos->Size(); const int ne = CAROM::split_dimension(ne_global, MPI_COMM_WORLD); @@ -376,9 +444,9 @@ void ROMInterfaceForm::SetupEQPSystem( { // NOTE(kevin): have to copy the vector since libROM matrix is row-major. for (int k = 0; k < fes1->GetTrueVSize(); ++k) - v1_i[k] = snapshot1(k, snap1_idx[i]); + v1_i[k] = snapshot1(k, snap_pair_idx(i, 0)); for (int k = 0; k < fes2->GetTrueVSize(); ++k) - v2_i[k] = snapshot2(k, snap2_idx[i]); + v2_i[k] = snapshot2(k, snap_pair_idx(i, 1)); for (int e = elem_offsets[rank], eidx = 0; e < elem_offsets[rank+1]; e++, eidx++) { @@ -457,7 +525,7 @@ void ROMInterfaceForm::SetupEQPSystem( void ROMInterfaceForm::TrainEQPForIntegrator( const int nqe, const CAROM::Matrix &Gt, const CAROM::Vector &rhs_Gw, - const double eqp_tol, Array &sample_el, Array &sample_qp, Array &sample_qw) + const double eqp_tol, Array &samples) { // void SolveNNLS(const int rank, const double nnls_tol, const int maxNNLSnnz, // CAROM::Vector const& w, CAROM::Matrix & Gt, @@ -536,22 +604,85 @@ void ROMInterfaceForm::TrainEQPForIntegrator( MPI_Allgatherv(eqpSol.getData(), eqpSol.dim(), MPI_DOUBLE, eqpSol_global.getData(), eqp_sol_cnts.data(), eqp_sol_offsets.data(), MPI_DOUBLE, MPI_COMM_WORLD); - sample_el.SetSize(0); - sample_qp.SetSize(0); - sample_qw.SetSize(0); + samples.SetSize(0); for (int i = 0; i < eqpSol_global.dim(); ++i) { if (eqpSol_global(i) > 1.0e-12) { const int e = i / nqe; // Element index - sample_el.Append(i / nqe); - sample_qp.Append(i % nqe); - sample_qw.Append(eqpSol_global(i)); + samples.Append({.el = i / nqe, .qp = i % nqe, .qw = eqpSol_global(i)}); } } - printf("Size of sampled qp: %d\n", sample_el.Size()); - if (nnz != sample_el.Size()) + printf("Size of sampled qp: %d\n", samples.Size()); + if (nnz != samples.Size()) printf("Sample quadrature points with weight < 1.0e-12 are neglected.\n"); } +void ROMInterfaceForm::SaveEQPForIntegrator(const int k, hid_t file_id, const std::string &dsetname) +{ + assert(file_id >= 0); + hid_t grp_id; + herr_t errf; + + grp_id = H5Gcreate(file_id, dsetname.c_str(), H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT); + assert(grp_id >= 0); + + Array *ref_sample = NULL; + int c1, c2, a1, a2; + std::string port_dset; + for (int p = 0; p < numRefPorts; p++) + { + topol_handler->GetRefPortInfo(p, c1, c2, a1, a2); + port_dset = topol_handler->GetComponentName(c1) + ":" + topol_handler->GetComponentName(c2); + port_dset += "-" + std::to_string(a1) + ":" + std::to_string(a2); + ref_sample = fnfi_ref_sample[p + k * numRefPorts]; + + hdf5_utils::WriteDataset(grp_id, port_dset, IntegratorType::INTERFACE, *ref_sample); + } + + errf = H5Gclose(grp_id); + assert(errf >= 0); + return; +} + +void ROMInterfaceForm::LoadEQPForIntegrator(const int k, hid_t file_id, const std::string &dsetname) +{ + assert(file_id >= 0); + hid_t grp_id; + herr_t errf; + + grp_id = H5Gopen2(file_id, dsetname.c_str(), H5P_DEFAULT); + assert(grp_id >= 0); + + Array *ref_sample = NULL; + int c1, c2, a1, a2; + std::string port_dset; + for (int p = 0; p < numRefPorts; p++) + { + topol_handler->GetRefPortInfo(p, c1, c2, a1, a2); + port_dset = topol_handler->GetComponentName(c1) + ":" + topol_handler->GetComponentName(c2); + port_dset += "-" + std::to_string(a1) + ":" + std::to_string(a2); + + const int port_idx = p + k * numRefPorts; + if (fnfi_ref_sample[port_idx]) + delete fnfi_ref_sample[port_idx]; + fnfi_ref_sample[port_idx] = new Array(0); + + hdf5_utils::ReadDataset(grp_id, port_dset, IntegratorType::INTERFACE, + *fnfi_ref_sample[port_idx]); + } + + /* update subdomain ports */ + for (int p = 0; p < numPorts; p++) + { + int rp = topol_handler->GetPortType(p); + + fnfi_sample[p + k * numPorts] = fnfi_ref_sample[rp + k * numRefPorts]; + } + + errf = H5Gclose(grp_id); + assert(errf >= 0); + return; +} + } diff --git a/src/rom_nonlinearform.cpp b/src/rom_nonlinearform.cpp index 614ab3c3..f4c73cf2 100644 --- a/src/rom_nonlinearform.cpp +++ b/src/rom_nonlinearform.cpp @@ -40,6 +40,8 @@ ROMNonlinearForm::~ROMNonlinearForm() { delete bfnfi[i]; delete bfnfi_sample[i]; + // Decided to own it, as opposed to the parent class. + delete bfnfi_marker[i]; } printf("%10s\t%10s\t%10s\t%10s\t%10s\t%10s\t%10s\n", @@ -160,7 +162,7 @@ void ROMNonlinearForm::Mult(const Vector &x, Vector &y) const SampleInfo *sample = sample_info->GetData(); for (int i = 0; i < sample_info->Size(); i++, sample++) { - int face = sample->face; + int face = sample->el; tr = mesh->GetInteriorFaceTransformations(face); const IntegrationPoint &ip = ir->IntPoint(sample->qp); @@ -236,7 +238,7 @@ void ROMNonlinearForm::Mult(const Vector &x, Vector &y) const SampleInfo *sample = sample_info->GetData(); for (int i = 0; i < sample_info->Size(); i++, sample++) { - int be = sample->be; + int be = sample->el; const int bdr_attr = mesh->GetBdrAttribute(be); if (bfnfi_marker[k] && (*bfnfi_marker[k])[bdr_attr-1] == 0) @@ -417,7 +419,7 @@ Operator& ROMNonlinearForm::GetGradient(const Vector &x) const SampleInfo *sample = sample_info->GetData(); for (int i = 0; i < sample_info->Size(); i++, sample++) { - int face = sample->face; + int face = sample->el; tr = mesh->GetInteriorFaceTransformations(face); const IntegrationPoint &ip = ir->IntPoint(sample->qp); @@ -493,7 +495,7 @@ Operator& ROMNonlinearForm::GetGradient(const Vector &x) const SampleInfo *sample = sample_info->GetData(); for (int i = 0; i < sample_info->Size(); i++, sample++) { - int be = sample->be; + int be = sample->el; const int bdr_attr = mesh->GetBdrAttribute(be); if (bfnfi_marker[k] && (*bfnfi_marker[k])[bdr_attr-1] == 0) @@ -637,7 +639,7 @@ void ROMNonlinearForm::PrecomputeCoefficients() SampleInfo *sample = sample_info->GetData(); for (int i = 0; i < sample_info->Size(); i++, sample++) { - int be = sample->be; + int be = sample->el; const int bdr_attr = mesh->GetBdrAttribute(be); if (bfnfi_marker[k] && @@ -687,25 +689,24 @@ void ROMNonlinearForm::TrainEQP(const CAROM::Matrix &snapshots, const double eqp CAROM::Matrix Gt(1,1, true); CAROM::Vector rhs_Gw(1, false); - Array el, qp; - Array qw; + Array samples; Array fidxs; Mesh *mesh = fes->GetMesh(); for (int k = 0; k < dnfi.Size(); k++) { SetupEQPSystemForDomainIntegrator(snapshots_work, dnfi[k], Gt, rhs_Gw); - TrainEQPForIntegrator(dnfi[k], Gt, rhs_Gw, eqp_tol, el, qp, qw); - UpdateDomainIntegratorSampling(k, el, qp, qw); + TrainEQPForIntegrator(dnfi[k], Gt, rhs_Gw, eqp_tol, samples); + UpdateDomainIntegratorSampling(k, samples); } for (int k = 0; k < fnfi.Size(); k++) { SetupEQPSystemForInteriorFaceIntegrator(snapshots_work, fnfi[k], Gt, rhs_Gw, fidxs); - TrainEQPForIntegrator(fnfi[k], Gt, rhs_Gw, eqp_tol, el, qp, qw); - for (int s = 0; s < el.Size(); s++) - el[s] = fidxs[el[s]]; - UpdateInteriorFaceIntegratorSampling(k, el, qp, qw); + TrainEQPForIntegrator(fnfi[k], Gt, rhs_Gw, eqp_tol, samples); + for (int s = 0; s < samples.Size(); s++) + samples[s].el = fidxs[samples[s].el]; + UpdateInteriorFaceIntegratorSampling(k, samples); } // Which boundary attributes need to be processed? @@ -731,16 +732,16 @@ void ROMNonlinearForm::TrainEQP(const CAROM::Matrix &snapshots, const double eqp } SetupEQPSystemForBdrFaceIntegrator(snapshots_work, bfnfi[k], bdr_attr_marker, Gt, rhs_Gw, fidxs); - TrainEQPForIntegrator(bfnfi[k], Gt, rhs_Gw, eqp_tol, el, qp, qw); - for (int s = 0; s < el.Size(); s++) - el[s] = fidxs[el[s]]; - UpdateBdrFaceIntegratorSampling(k, el, qp, qw); + TrainEQPForIntegrator(bfnfi[k], Gt, rhs_Gw, eqp_tol, samples); + for (int s = 0; s < samples.Size(); s++) + samples[s].el = fidxs[samples[s].el]; + UpdateBdrFaceIntegratorSampling(k, samples); } } void ROMNonlinearForm::TrainEQPForIntegrator( HyperReductionIntegrator *nlfi, const CAROM::Matrix &Gt, const CAROM::Vector &rhs_Gw, - const double eqp_tol, Array &sample_el, Array &sample_qp, Array &sample_qw) + const double eqp_tol, Array &samples) { const IntegrationRule *ir = nlfi->GetIntegrationRule(); @@ -751,6 +752,7 @@ void ROMNonlinearForm::TrainEQPForIntegrator( // void SolveNNLS(const int rank, const double nnls_tol, const int maxNNLSnnz, // CAROM::Vector const& w, CAROM::Matrix & Gt, // CAROM::Vector & sol) + const double normRHS = rhs_Gw.norm(); double nnls_tol = 1.0e-11; int maxNNLSnnz = 0; CAROM::Vector eqpSol(Gt.numRows(), true); @@ -803,7 +805,6 @@ void ROMNonlinearForm::TrainEQPForIntegrator( Gt.transposeMult(eqpSol, res); const double normGsol = res.norm(); - const double normRHS = rhs_Gw.norm(); res -= rhs_Gw; const double relNorm = res.norm() / std::max(normGsol, normRHS); @@ -825,21 +826,17 @@ void ROMNonlinearForm::TrainEQPForIntegrator( MPI_Allgatherv(eqpSol.getData(), eqpSol.dim(), MPI_DOUBLE, eqpSol_global.getData(), eqp_sol_cnts.data(), eqp_sol_offsets.data(), MPI_DOUBLE, MPI_COMM_WORLD); - sample_el.SetSize(0); - sample_qp.SetSize(0); - sample_qw.SetSize(0); + samples.SetSize(0); for (int i = 0; i < eqpSol_global.dim(); ++i) { if (eqpSol_global(i) > 1.0e-12) { const int e = i / nqe; // Element index - sample_el.Append(i / nqe); - sample_qp.Append(i % nqe); - sample_qw.Append(eqpSol_global(i)); + samples.Append({.el = i / nqe, .qp = i % nqe, .qw = eqpSol_global(i)}); } } - printf("Size of sampled qp: %d\n", sample_el.Size()); - if (nnz != sample_el.Size()) + printf("Size of sampled qp: %d\n", samples.Size()); + if (nnz != samples.Size()) printf("Sample quadrature points with weight < 1.0e-12 are neglected.\n"); } @@ -1221,8 +1218,7 @@ void ROMNonlinearForm::SetupEQPSystemForBdrFaceIntegrator( return; } -void ROMNonlinearForm::GetEQPForIntegrator( - const IntegratorType type, const int k, Array &el, Array &qp, Array &qw) +Array* ROMNonlinearForm::GetEQPForIntegrator(const IntegratorType type, const int k) { Array *sample = NULL; switch (type) @@ -1252,94 +1248,90 @@ void ROMNonlinearForm::GetEQPForIntegrator( mfem_error("Unknown Integrator type!\n"); } - el.SetSize(0); - qp.SetSize(0); - qw.SetSize(0); - - for (int s = 0; s < sample->Size(); s++) - { - switch (type) - { - case IntegratorType::DOMAIN: el.Append((*sample)[s].el); break; - case IntegratorType::INTERIORFACE: el.Append((*sample)[s].face); break; - case IntegratorType::BDRFACE: el.Append((*sample)[s].be); break; - } - qp.Append((*sample)[s].qp); - qw.Append((*sample)[s].qw); - } + return sample; } void ROMNonlinearForm::SaveEQPForIntegrator( const IntegratorType type, const int k, hid_t file_id, const std::string &dsetname) { - std::string eldset; + Array *sample = NULL; switch (type) { - case IntegratorType::DOMAIN: eldset = "elem"; break; - case IntegratorType::INTERIORFACE: eldset = "face"; break; - case IntegratorType::BDRFACE: eldset = "be"; break; + case IntegratorType::DOMAIN: + { + assert((k >= 0) && (k < dnfi.Size())); + assert(dnfi.Size() == dnfi_sample.Size()); + sample = dnfi_sample[k]; + } + break; + case IntegratorType::INTERIORFACE: + { + assert((k >= 0) && (k < fnfi.Size())); + assert(fnfi.Size() == fnfi_sample.Size()); + sample = fnfi_sample[k]; + } + break; + case IntegratorType::BDRFACE: + { + assert((k >= 0) && (k < bfnfi.Size())); + assert(bfnfi.Size() == bfnfi_sample.Size()); + sample = bfnfi_sample[k]; + } + break; default: - mfem_error("ROMNonlinearForm::SaveEQPForIntegrator- Unknown IntegratorType!\n"); + mfem_error("Unknown Integrator type!\n"); } + assert(sample); - Array el, qp; - Array qw; - GetEQPForIntegrator(type, k, el, qp, qw); - - assert(file_id >= 0); - hid_t grp_id; - herr_t errf; - - grp_id = H5Gcreate(file_id, dsetname.c_str(), H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT); - assert(grp_id >= 0); - - hdf5_utils::WriteDataset(grp_id, eldset, el); - hdf5_utils::WriteDataset(grp_id, "quad-pt", qp); - hdf5_utils::WriteDataset(grp_id, "quad-wt", qw); - - errf = H5Gclose(grp_id); - assert(errf >= 0); + hdf5_utils::WriteDataset(file_id, dsetname, type, *sample); return; } void ROMNonlinearForm::LoadEQPForIntegrator( const IntegratorType type, const int k, hid_t file_id, const std::string &dsetname) { - std::string eldset; + Array *sample = NULL; switch (type) { - case IntegratorType::DOMAIN: eldset = "elem"; break; - case IntegratorType::INTERIORFACE: eldset = "face"; break; - case IntegratorType::BDRFACE: eldset = "be"; break; - default: - mfem_error("ROMNonlinearForm::LoadEQPForIntegrator- Unknown IntegratorType!\n"); - } - - Array el, qp; - Array qw; - - assert(file_id >= 0); - hid_t grp_id; - herr_t errf; - - grp_id = H5Gopen2(file_id, dsetname.c_str(), H5P_DEFAULT); - assert(grp_id >= 0); - - hdf5_utils::ReadDataset(grp_id, eldset, el); - hdf5_utils::ReadDataset(grp_id, "quad-pt", qp); - hdf5_utils::ReadDataset(grp_id, "quad-wt", qw); + case IntegratorType::DOMAIN: + { + assert((k >= 0) && (k < dnfi.Size())); + assert(dnfi.Size() == dnfi_sample.Size()); + if (dnfi_sample[k]) + delete dnfi_sample[k]; - errf = H5Gclose(grp_id); - assert(errf >= 0); + dnfi_sample[k] = new Array(0); + sample = dnfi_sample[k]; + } + break; + case IntegratorType::INTERIORFACE: + { + assert((k >= 0) && (k < fnfi.Size())); + assert(fnfi.Size() == fnfi_sample.Size()); + if (fnfi_sample[k]) + delete fnfi_sample[k]; - switch (type) - { - case IntegratorType::DOMAIN: UpdateDomainIntegratorSampling(k, el, qp, qw); break; - case IntegratorType::INTERIORFACE: UpdateInteriorFaceIntegratorSampling(k, el, qp, qw); break; - case IntegratorType::BDRFACE: UpdateBdrFaceIntegratorSampling(k, el, qp, qw); break; + fnfi_sample[k] = new Array(0); + sample = fnfi_sample[k]; + } + break; + case IntegratorType::BDRFACE: + { + assert((k >= 0) && (k < bfnfi.Size())); + assert(bfnfi.Size() == bfnfi_sample.Size()); + if (bfnfi_sample[k]) + delete bfnfi_sample[k]; + + bfnfi_sample[k] = new Array(0); + sample = bfnfi_sample[k]; + } + break; default: - mfem_error("ROMNonlinearForm::LoadEQPForIntegrator- Unknown IntegratorType!\n"); + mfem_error("Unknown Integrator type!\n"); } + assert(sample); + + hdf5_utils::ReadDataset(file_id, dsetname, type, *sample); return; } diff --git a/src/sample_generator.cpp b/src/sample_generator.cpp index b2c22b57..4e51a26e 100644 --- a/src/sample_generator.cpp +++ b/src/sample_generator.cpp @@ -177,7 +177,7 @@ const std::string SampleGenerator::GetSamplePath(const int &idx, const std::stri return full_path; } -void SampleGenerator::SaveSnapshot(BlockVector *U_snapshots, std::vector &snapshot_basis_tags, Array &col_idxs) +void SampleGenerator::SaveSnapshot(BlockVector *U_snapshots, std::vector &snapshot_basis_tags, Array &col_idxs) { assert(U_snapshots->NumBlocks() == snapshot_basis_tags.size()); col_idxs.SetSize(U_snapshots->NumBlocks()); @@ -218,13 +218,13 @@ void SampleGenerator::SaveSnapshotPorts(TopologyHandler *topol_handler, const Ar /* We save the snapshot ports by component names. - 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. + Variable separation does not matter at this point.. */ + int c1 = topol_handler->GetMeshType(port->Mesh1); + int c2 = topol_handler->GetMeshType(port->Mesh2); std::string mesh1, mesh2; - mesh1 = GetBasisTag(port->Mesh1, topol_handler); - mesh2 = GetBasisTag(port->Mesh2, topol_handler); + mesh1 = topol_handler->GetComponentName(c1); + mesh2 = topol_handler->GetComponentName(c2); /* note the mesh names are not necessarily the same as the subdomain names @@ -262,7 +262,7 @@ void SampleGenerator::SaveSnapshotPorts(TopologyHandler *topol_handler, const Ar } } -void SampleGenerator::AddSnapshotGenerator(const int &fom_vdofs, const std::string &prefix, const std::string &basis_tag) +void SampleGenerator::AddSnapshotGenerator(const int &fom_vdofs, const std::string &prefix, const BasisTag &basis_tag) { const std::string filename = GetBaseFilename(prefix, basis_tag); @@ -298,9 +298,13 @@ void SampleGenerator::WriteSnapshotPorts() file_id = H5Fcreate(filename.c_str(), H5F_ACC_TRUNC, H5P_DEFAULT, H5P_DEFAULT); assert(file_id >= 0); - hdf5_utils::WriteAttribute(file_id, "sample_prefix", GetSamplePrefix()); + /* this is the path to all associated snapshot matrices */ + char sample_path[PATH_MAX]; + realpath(GetSamplePrefix().c_str(), sample_path); + hdf5_utils::WriteAttribute(file_id, "sample_prefix", std::string(sample_path)); + + /* write port information */ hdf5_utils::WriteAttribute(file_id, "number_of_ports", (int) port_tags.size()); - for (int p = 0; p < port_tags.size(); p++) { hid_t grp_id; @@ -317,12 +321,17 @@ void SampleGenerator::WriteSnapshotPorts() assert(errf >= 0); } + /* write basis tag list */ + hdf5_utils::WriteAttribute(file_id, "number_of_basistags", (int) basis_tags.size()); + for (int b = 0; b < basis_tags.size(); b++) + hdf5_utils::WriteAttribute(file_id, std::string("basistag" + std::to_string(b)).c_str(), basis_tags[b]); + errf = H5Fclose(file_id); assert(errf >= 0); return; } -const CAROM::Matrix* SampleGenerator::LookUpSnapshot(const std::string &basis_tag) +const CAROM::Matrix* SampleGenerator::LookUpSnapshot(const BasisTag &basis_tag) { assert(snapshot_generators.Size() > 0); assert(snapshot_generators.Size() == basis_tags.size()); @@ -337,13 +346,28 @@ const CAROM::Matrix* SampleGenerator::LookUpSnapshot(const std::string &basis_ta if (idx < 0) { - printf("basis tag: %s\n", basis_tag.c_str()); + printf("basis tag: %s\n", basis_tag.print().c_str()); mfem_error("SampleGenerator::LookUpSnapshot- basis tag does not exist in snapshot list!\n"); } return snapshot_generators[idx]->getSnapshotMatrix(); } +Array2D* SampleGenerator::LookUpSnapshotPortColOffsets(const PortTag &port_tag) +{ + assert(port_colidxs.Size() > 0); + assert(port_colidxs.Size() == port_tags.size()); + + int idx = -1; + if (port_tag2idx.count(port_tag)) + idx = port_tag2idx[port_tag]; + + if (idx < 0) + mfem_error("SampleGenerator::LookUpSnapshotPortColOffsets- port tag does not exist in port list!\n"); + + return port_colidxs[idx]; +} + void SampleGenerator::ReportStatus(const int &sample_idx) { if (sample_idx % report_freq != 0) return; @@ -354,12 +378,12 @@ void SampleGenerator::ReportStatus(const int &sample_idx) printf("Basis tags: %ld\n", basis_tags.size()); printf("%20.20s\t# of snapshots\n", "Basis tags"); for (int k = 0; k < basis_tags.size(); k++) - printf("%20.20s\t%d\n", basis_tags[k].c_str(), snapshot_generators[k]->getNumSamples()); + printf("%20.20s\t%d\n", basis_tags[k].print().c_str(), snapshot_generators[k]->getNumSamples()); printf("==============================================\n"); } -void SampleGenerator::CollectSnapshots(const std::string &basis_prefix, - const std::string &basis_tag, +void SampleGenerator::CollectSnapshotsByBasis(const std::string &basis_prefix, + const BasisTag &basis_tag, const std::vector &file_list) { // Get dimension from the first snapshot file. @@ -373,14 +397,91 @@ void SampleGenerator::CollectSnapshots(const std::string &basis_prefix, */ int local_num_vdofs = CAROM::split_dimension(fom_num_vdof, MPI_COMM_WORLD); - // Append snapshot generator. - AddSnapshotGenerator(local_num_vdofs, basis_prefix, basis_tag); - CAROM::BasisGenerator *basis_generator = snapshot_generators.Last(); + /* if the tag was never seen before, append a new snapshot generator */ + if (!basis_tag2idx.count(basis_tag)) + AddSnapshotGenerator(local_num_vdofs, basis_prefix, basis_tag); + int index = basis_tag2idx[basis_tag]; + CAROM::BasisGenerator *basis_generator = snapshot_generators[index]; for (int s = 0; s < file_list.size(); s++) basis_generator->loadSamples(file_list[s], "snapshot", 1e9, CAROM::Database::formats::HDF5_MPIO); } +void SampleGenerator::CollectSnapshotsByPort( + const std::string &basis_prefix, const std::string &port_tag_file) +{ + // TODO(kevin): we need to see how this impacts the parallel I/O.. + + hid_t file_id; + herr_t errf = 0; + file_id = H5Fopen(port_tag_file.c_str(), H5F_ACC_RDONLY, H5P_DEFAULT); + assert(file_id >= 0); + + std::string sample_path; + hdf5_utils::ReadAttribute(file_id, "sample_prefix", sample_path); + printf("SampleGenerator: snapshot port prefix=%s\n", sample_path.c_str()); + + int num_ports = -1; + hdf5_utils::ReadAttribute(file_id, "number_of_ports", num_ports); + + for (int p = 0; p < num_ports; p++) + { + hid_t grp_id; + grp_id = H5Gopen2(file_id, std::to_string(p).c_str(), H5P_DEFAULT); + assert(grp_id >= 0); + + PortTag port_tag; + hdf5_utils::ReadAttribute(grp_id, "Mesh1", port_tag.Mesh1); + hdf5_utils::ReadAttribute(grp_id, "Mesh2", port_tag.Mesh2); + hdf5_utils::ReadAttribute(grp_id, "Attr1", port_tag.Attr1); + hdf5_utils::ReadAttribute(grp_id, "Attr2", port_tag.Attr2); + + /* if the tag was never seen before, create a new port tag */ + if (!port_tag2idx.count(port_tag)) + { + port_tag2idx[port_tag] = port_tags.size(); + port_tags.push_back(port_tag); + port_colidxs.Append(new Array2D); + } + int idx = port_tag2idx[port_tag]; + + /* if other snapshots are already collected, column indices must be offseted */ + int col_offset1 = GetSnapshotOffset(port_tag.Mesh1); + int col_offset2 = GetSnapshotOffset(port_tag.Mesh2); + + Array2D tmp_colidx; + hdf5_utils::ReadDataset(grp_id, "col_idxs", tmp_colidx); + + int row_offset = port_colidxs[idx]->NumRows(); + port_colidxs[idx]->SetSize(row_offset + tmp_colidx.NumRows(), 2); + for (int r = 0, r0 = row_offset; r < tmp_colidx.NumRows(); r++, r0++) + { + (*port_colidxs[idx])(r0, 0) = tmp_colidx(r, 0) + col_offset1; + (*port_colidxs[idx])(r0, 1) = tmp_colidx(r, 1) + col_offset2; + } + + errf = H5Gclose(grp_id); + assert(errf >= 0); + } // for (int p = 0; p < num_ports; p++) + + /* read snapshots from basis tag list */ + int num_basistag = -1; + hdf5_utils::ReadAttribute(file_id, "number_of_basistags", num_basistag); + assert(num_basistag > 0); + BasisTag basis_tag; + std::vector snapshot_file(1); + for (int b = 0; b < num_basistag; b++) + { + hdf5_utils::ReadAttribute(file_id, std::string("basistag" + std::to_string(b)).c_str(), basis_tag); + snapshot_file[0] = GetBaseFilename(sample_path, basis_tag) + "_snapshot"; + CollectSnapshotsByBasis(basis_prefix, basis_tag, snapshot_file); + } + + errf = H5Fclose(file_id); + assert(errf >= 0); + return; +} + void SampleGenerator::FormReducedBasis(const std::string &basis_prefix) { assert(snapshot_generators.Size() > 0); @@ -402,7 +503,7 @@ void SampleGenerator::FormReducedBasis(const std::string &basis_prefix) if (basis_list) { // Find if additional inputs are specified for basis_tags[p]. - YAML::Node basis_tag_input = config.LookUpFromDict("name", basis_tags[k], basis_list); + YAML::Node basis_tag_input = config.LookUpFromDict("name", basis_tags[k].print(), basis_list); // If basis_tags[p] has additional inputs, parse them. // parse tag-specific number of basis. @@ -480,4 +581,25 @@ void SampleGenerator::SaveSV(CAROM::BasisGenerator *basis_generator, const std:: // TODO: hdf5 format + parallel case. std::string filename = prefix + "_sv.txt"; CAROM::PrintVector(*rom_sv, filename); +} + +const int SampleGenerator::GetSnapshotOffset(const std::string &comp) +{ + int offset = 0, tmp = -1; + + assert(basis_tags.size() == snapshot_generators.Size()); + for (int b = 0; b < basis_tags.size(); b++) + { + if (basis_tags[b].comp != comp) continue; + + offset = snapshot_generators[b]->getNumSamples(); + + if (tmp < 0) + tmp = offset; + + if (offset != tmp) + mfem_error("SampleGenerator::GetSnapshotOffset- Number of samples over variables does not match!\n"); + } + + return offset; } \ No newline at end of file diff --git a/src/steady_ns_solver.cpp b/src/steady_ns_solver.cpp index d9196095..f4ef8119 100644 --- a/src/steady_ns_solver.cpp +++ b/src/steady_ns_solver.cpp @@ -214,6 +214,37 @@ Operator& SteadyNSTensorROM::GetGradient(const Vector &x) const SteadyNSTensorROM */ +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_) +{ + 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; + } + + u_offsets.SetSize(numSub+1); + u_offsets = 0; + for (int k = 0; k < numSub; k++) + { + u_offsets[k+1] = block_offsets[1 + k * num_var] - block_offsets[k * num_var]; + } + u_offsets.PartialSum(); + + 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 { y = 0.0; @@ -227,6 +258,13 @@ void SteadyNSEQPROM::Mult(const Vector &x, Vector &y) const hs[m]->AddMult(x_comp, y_comp); } + + if (!itf) return; + + GetVel(x, x_u); + y_u = 0.0; + itf->InterfaceAddMult(x_u, y_u); + AddVel(y_u, y); } Operator& SteadyNSEQPROM::GetGradient(const Vector &x) const @@ -234,8 +272,36 @@ Operator& SteadyNSEQPROM::GetGradient(const Vector &x) const delete jac_mono; delete jac_hypre; jac_mono = new SparseMatrix(*linearOp); - DenseMatrix *jac_comp; + if (itf) + { + BlockMatrix jac(block_offsets); + jac.owns_blocks = true; + + Array2D mats(numSub, numSub); + 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; + mats(i, j) = new SparseMatrix(block_offsets[iidx+1] - block_offsets[iidx], + block_offsets[jidx+1] - block_offsets[jidx]); + } + + GetVel(x, x_u); + itf->InterfaceGetGradient(x_u, mats); + + MatrixBlocks u_mats(mats); + + AddToBlockMatrix(u_idxs, u_idxs, u_mats, jac); + jac.Finalize(); + SparseMatrix *tmp = jac.CreateMonolithic(); + (*jac_mono) += *tmp; + + delete tmp; + } + + DenseMatrix *jac_comp; for (int m = 0; m < numSub; m++) { int midx = (separate_variable) ? m * num_var : m; @@ -256,6 +322,43 @@ Operator& SteadyNSEQPROM::GetGradient(const Vector &x) const return *jac_mono; } +void SteadyNSEQPROM::GetVel(const Vector &x, Vector &x_u) const +{ + if (!separate_variable) + { + x_u = x; + return; + } + x_u.SetSize(u_offsets.Last()); + + BlockVector x_block(const_cast(x).GetData(), block_offsets); + BlockVector u_block(x_u.GetData(), u_offsets); + for (int m = 0; m < numSub; m++) + { + Vector tmp; + u_block.GetBlockView(m, tmp); + tmp = x_block.GetBlock(m * num_var); + } +} + +void SteadyNSEQPROM::AddVel(const Vector &y_u, Vector &y) const +{ + if (!separate_variable) + { + y += y_u; + return; + } + + BlockVector u_block(const_cast(y_u).GetData(), u_offsets); + BlockVector y_block(y.GetData(), block_offsets); + for (int m = 0; m < numSub; m++) + { + Vector tmp; + y_block.GetBlockView(m * num_var, tmp); + tmp += u_block.GetBlock(m); + } +} + /* SteadyNSSolver */ @@ -265,6 +368,10 @@ SteadyNSSolver::SteadyNSSolver() { nonlinear_mode = true; + Vector zero_vec(dim); + zero_vec = 0.0; + zero = new VectorConstantCoefficient(zero_vec); + // StokesSolver reads viscosity from stokes/nu. nu = config.GetOption("stokes/nu", 1.0); delete nu_coeff; @@ -281,7 +388,10 @@ SteadyNSSolver::SteadyNSSolver() else mfem_error("SteadyNSSolver: unknown operator type!\n"); - ir_nl = &(IntRules.Get(ufes[0]->GetFE(0)->GetGeomType(), (int)(ceil(1.5 * (2 * ufes[0]->GetMaxElementOrder() - 1))))); + ir_nl = &(IntRules.Get(ufes[0]->GetFE(0)->GetGeomType(), + (int)(ceil(1.5 * (2 * ufes[0]->GetMaxElementOrder() - 1))))); + ir_face = &(IntRules.Get(meshes[0]->GetBdrElementGeometry(0), + (int)(ceil(1.5 * (2 * ufes[0]->GetMaxElementOrder() - 1))))); /* SteadyNSSolver requires all the meshes to have the same element type. */ int num_comp = topol_handler->GetNumComponents(); @@ -293,6 +403,7 @@ SteadyNSSolver::SteadyNSSolver() SteadyNSSolver::~SteadyNSSolver() { + delete zero; delete zeta_coeff; delete minus_zeta; delete minus_half_zeta; @@ -312,20 +423,10 @@ SteadyNSSolver::~SteadyNSSolver() DeletePointers(subdomain_tensors); } else if (rom_handler->GetNonlinearHandling() == NonlinearHandling::EQP) + { DeletePointers(comp_eqps); - } -} - -void SteadyNSSolver::InitVariables() -{ - StokesSolver::InitVariables(); - if (use_rom) - { - rom_handler->SetNonlinearMode(nonlinear_mode); - subdomain_tensors.SetSize(numSub); - subdomain_tensors = NULL; - subdomain_eqps.SetSize(numSub); - subdomain_eqps = NULL; + delete itf_eqp; + } } } @@ -351,10 +452,12 @@ void SteadyNSSolver::BuildDomainOperators() { auto *lf_integ1 = new IncompressibleInviscidFluxNLFIntegrator(*minus_zeta); lf_integ1->SetIntRule(ir_nl); + auto *lf_integ2 = new DGLaxFriedrichsFluxIntegrator(*minus_zeta); + lf_integ2->SetIntRule(ir_face); hs[m]->AddDomainIntegrator(lf_integ1); if (full_dg) - hs[m]->AddInteriorFaceIntegrator(new DGLaxFriedrichsFluxIntegrator(*minus_zeta)); + hs[m]->AddInteriorFaceIntegrator(lf_integ2); } break; default: @@ -366,8 +469,9 @@ void SteadyNSSolver::BuildDomainOperators() if (oper_type == OperType::LF) { nl_itf = new InterfaceForm(meshes, ufes, topol_handler); - nl_itf->AddInterfaceIntegrator(new DGLaxFriedrichsFluxIntegrator(*minus_zeta)); - // nl_interface->SetIntRule(ir_nl); + auto *lf_integ2 = new DGLaxFriedrichsFluxIntegrator(*minus_zeta); + lf_integ2->SetIntRule(ir_face); + nl_itf->AddInterfaceIntegrator(lf_integ2); } } @@ -377,6 +481,8 @@ void SteadyNSSolver::SetupDomainBCOperators() if (oper_type != OperType::LF) return; + HyperReductionIntegrator *lf_integ2 = NULL; + assert(hs.Size() == numSub); for (int m = 0; m < numSub; m++) { @@ -388,12 +494,15 @@ void SteadyNSSolver::SetupDomainBCOperators() // TODO: Non-homogeneous Neumann stress bc if (bdr_type[b] == BoundaryType::NEUMANN) - hs[m]->AddBdrFaceIntegrator(new DGLaxFriedrichsFluxIntegrator(*minus_zeta), *bdr_markers[b]); + lf_integ2 = new DGLaxFriedrichsFluxIntegrator(*minus_zeta); else { assert(BCExistsOnBdr(b)); - hs[m]->AddBdrFaceIntegrator(new DGLaxFriedrichsFluxIntegrator(*minus_zeta, ud_coeffs[b]), *bdr_markers[b]); + lf_integ2 = new DGLaxFriedrichsFluxIntegrator(*minus_zeta, ud_coeffs[b]); } + + lf_integ2->SetIntRule(ir_face); + hs[m]->AddBdrFaceIntegrator(lf_integ2, *bdr_markers[b]); } } @@ -417,12 +526,11 @@ void SteadyNSSolver::SaveROMOperator(const std::string input_prefix) if (rom_handler->GetNonlinearHandling() != NonlinearHandling::TENSOR) return; - std::string filename = rom_handler->GetOperatorPrefix() + ".h5"; - assert(FileExists(filename)); + std::string filename = rom_handler->GetOperatorPrefix() + ".tensor.h5"; hid_t file_id; herr_t errf = 0; - file_id = H5Fopen(filename.c_str(), H5F_ACC_RDWR, H5P_DEFAULT); + file_id = H5Fcreate(filename.c_str(), H5F_ACC_TRUNC, H5P_DEFAULT, H5P_DEFAULT); assert(file_id >= 0); hid_t grp_id; @@ -451,7 +559,7 @@ void SteadyNSSolver::LoadROMOperatorFromFile(const std::string input_prefix) subdomain_tensors.SetSize(numSub); subdomain_tensors = NULL; { - std::string filename = rom_handler->GetOperatorPrefix() + ".h5"; + std::string filename = rom_handler->GetOperatorPrefix() + ".tensor.h5"; assert(FileExists(filename)); hid_t file_id; @@ -621,7 +729,7 @@ void SteadyNSSolver::SolveROM() { 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, *(rom_handler->GetBlockOffsets())); + rom_oper = new SteadyNSEQPROM(rom_handler->GetOperator(), subdomain_eqps, itf_eqp, *(rom_handler->GetBlockOffsets())); } rom_handler->NonlinearSolve(*rom_oper, U_domain); @@ -629,10 +737,93 @@ void SteadyNSSolver::SolveROM() delete rom_oper; } -void SteadyNSSolver::AllocateROMTensorElems() +void SteadyNSSolver::InitROMHandler() { - assert(topol_mode == TopologyHandlerMode::COMPONENT); + StokesSolver::InitROMHandler(); + + rom_handler->SetNonlinearMode(nonlinear_mode); + subdomain_tensors.SetSize(numSub); + subdomain_tensors = NULL; + subdomain_eqps.SetSize(numSub); + subdomain_eqps = NULL; + + if (oper_type == OperType::LF) + { + Array comp_ufes(topol_handler->GetNumComponents()); + for (int c = 0; c < comp_ufes.Size(); c++) + comp_ufes[c] = comp_fes[c * num_var]; + itf_eqp = new ROMInterfaceForm(meshes, ufes, comp_ufes, topol_handler); + } +} + +void SteadyNSSolver::AllocateROMNlinElems() +{ + assert(rom_handler); + + switch (rom_handler->GetNonlinearHandling()) + { + case NonlinearHandling::TENSOR: + AllocateROMTensorElems(); + break; + case NonlinearHandling::EQP: + AllocateROMEQPElems(); + break; + default: + mfem_error("SteadyNSSolver::InitROMHandler- cannot initiate ROM elements!"); + break; + } +} + +void SteadyNSSolver::SaveROMNlinElems(const std::string &input_prefix) +{ + switch (rom_handler->GetNonlinearHandling()) + { + case NonlinearHandling::TENSOR: + SaveROMTensorElems(input_prefix + ".tensor.h5"); + break; + case NonlinearHandling::EQP: + SaveEQPElems(input_prefix + ".eqp.h5"); + break; + default: + mfem_error("SteadyNSSolver::SaveROMNlinElems- cannot initiate ROM elements!"); + break; + } +} + +void SteadyNSSolver::LoadROMNlinElems(const std::string &input_prefix) +{ + switch (rom_handler->GetNonlinearHandling()) + { + case NonlinearHandling::TENSOR: + LoadROMTensorElems(input_prefix + ".tensor.h5"); + break; + case NonlinearHandling::EQP: + LoadEQPElems(input_prefix + ".eqp.h5"); + break; + default: + mfem_error("SteadyNSSolver::SaveROMNlinElems- cannot initiate ROM elements!"); + break; + } +} + +void SteadyNSSolver::AssembleROMNlinOper() +{ + switch (rom_handler->GetNonlinearHandling()) + { + case NonlinearHandling::TENSOR: + AssembleROMTensorOper(); + break; + case NonlinearHandling::EQP: + AssembleROMEQPOper(); + break; + default: + mfem_error("SteadyNSSolver::SaveROMNlinElems- cannot initiate ROM elements!"); + break; + } +} +void SteadyNSSolver::AllocateROMTensorElems() +{ const int num_comp = topol_handler->GetNumComponents(); comp_tensors.SetSize(num_comp); comp_tensors = NULL; @@ -645,8 +836,6 @@ void SteadyNSSolver::BuildROMTensorElems() // Component domain system const int num_comp = topol_handler->GetNumComponents(); - Array fes_comp; - GetComponentFESpaces(fes_comp); DenseMatrix *basis = NULL; assert(comp_tensors.Size() == num_comp); @@ -656,27 +845,24 @@ void SteadyNSSolver::BuildROMTensorElems() const int fidx = c * num_var; const int cidx = (separate_variable_basis) ? fidx : c; rom_handler->GetReferenceBasis(cidx, basis); - comp_tensors[c] = GetReducedTensor(basis, fes_comp[fidx]); + comp_tensors[c] = GetReducedTensor(basis, comp_fes[fidx]); } // for (int c = 0; c < num_comp; c++) - - for (int k = 0 ; k < fes_comp.Size(); k++) delete fes_comp[k]; } void SteadyNSSolver::SaveROMTensorElems(const std::string &filename) { assert(topol_mode == TopologyHandlerMode::COMPONENT); - assert(FileExists(filename)); hid_t file_id; herr_t errf = 0; - file_id = H5Fopen(filename.c_str(), H5F_ACC_RDWR, H5P_DEFAULT); + file_id = H5Fcreate(filename.c_str(), H5F_ACC_TRUNC, H5P_DEFAULT, H5P_DEFAULT); assert(file_id >= 0); const int num_comp = topol_handler->GetNumComponents(); assert(comp_tensors.Size() == num_comp); hid_t grp_id; - grp_id = H5Gopen2(file_id, "components", H5P_DEFAULT); + grp_id = H5Gcreate(file_id, "components", H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT); assert(grp_id >= 0); std::string dset_name; @@ -686,7 +872,7 @@ void SteadyNSSolver::SaveROMTensorElems(const std::string &filename) dset_name = topol_handler->GetComponentName(c); hid_t comp_grp_id; - comp_grp_id = H5Gopen2(grp_id, dset_name.c_str(), H5P_DEFAULT); + 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, "tensor", *comp_tensors[c]); @@ -758,8 +944,6 @@ void SteadyNSSolver::AssembleROMTensorOper() void SteadyNSSolver::AllocateROMEQPElems() { - assert(topol_mode == TopologyHandlerMode::COMPONENT); - assert(rom_handler); assert(rom_handler->BasisLoaded()); @@ -770,22 +954,86 @@ void SteadyNSSolver::AllocateROMEQPElems() comp_eqps = NULL; DenseMatrix *basis; + HyperReductionIntegrator *nl_integ_tmp = NULL; + InterfaceNonlinearFormIntegrator *lf_integ2 = NULL; for (int c = 0; c < num_comp; c++) { int idx = (separate_variable_basis) ? c * num_var : c; rom_handler->GetReferenceBasis(idx, basis); - auto nl_integ_tmp = new VectorConvectionTrilinearFormIntegrator(*zeta_coeff); - nl_integ_tmp->SetIntRule(ir_nl); - comp_eqps[c] = new ROMNonlinearForm(basis->NumCols(), comp_fes[c * num_var]); - comp_eqps[c]->AddDomainIntegrator(nl_integ_tmp); + + switch (oper_type) + { + case (OperType::BASE): + nl_integ_tmp = new VectorConvectionTrilinearFormIntegrator(*zeta_coeff); + nl_integ_tmp->SetIntRule(ir_nl); + comp_eqps[c]->AddDomainIntegrator(nl_integ_tmp); + break; + + case (OperType::LF): + nl_integ_tmp = new IncompressibleInviscidFluxNLFIntegrator(*minus_zeta); + nl_integ_tmp->SetIntRule(ir_nl); + lf_integ2 = new DGLaxFriedrichsFluxIntegrator(*minus_zeta); + lf_integ2->SetIntRule(ir_face); + + comp_eqps[c]->AddDomainIntegrator(nl_integ_tmp); + if (full_dg) + comp_eqps[c]->AddInteriorFaceIntegrator(lf_integ2); + break; + + default: + break; + } + comp_eqps[c]->SetBasis(*basis); comp_eqps[c]->SetPrecomputeMode(precompute); - } + } // for (int c = 0; c < num_comp; c++) + + if (oper_type != OperType::LF) + return; + + /* allocate rom nonlinear form with boundary integrators */ + InterfaceNonlinearFormIntegrator *lf_bdr = NULL; + for (int c = 0; c < num_comp; c++) + { + /* For LF integrator, boundary eqp integrators are needed */ + Mesh *comp = topol_handler->GetComponentMesh(c); + for (int b = 0; b < comp->bdr_attributes.Size(); b++) + { + Array bdr_marker(comp->bdr_attributes.Max()); + bdr_marker = 0; + bdr_marker[comp->bdr_attributes[b] - 1] = 1; + + /* non-homogeneous Dirichlet + Neumann */ + lf_bdr = new DGLaxFriedrichsFluxIntegrator(*minus_zeta); + lf_bdr->SetIntRule(ir_face); + comp_eqps[c]->AddBdrFaceIntegrator(lf_bdr, bdr_marker); + + /* homogeneous Dirichlet */ + lf_bdr = new DGLaxFriedrichsFluxIntegrator(*minus_zeta, zero); + lf_bdr->SetIntRule(ir_face); + comp_eqps[c]->AddBdrFaceIntegrator(lf_bdr, bdr_marker); + } // for (int b = 0; b < comp->bdr_attributes.Size(); b++) + } // for (int c = 0; c < num_comp; c++) + + /* allocate rom interface form */ + assert(itf_eqp); + lf_integ2 = new DGLaxFriedrichsFluxIntegrator(*minus_zeta); + lf_integ2->SetIntRule(ir_face); + + itf_eqp->AddInterfaceIntegrator(lf_integ2); + + for (int c = 0; c < num_comp; c++) + { + int idx = (separate_variable_basis) ? c * num_var : c; + rom_handler->GetReferenceBasis(idx, basis); + itf_eqp->SetBasisAtComponent(c, *basis); + } + itf_eqp->UpdateBlockOffsets(); } -void SteadyNSSolver::TrainEQPElems(SampleGenerator *sample_generator) +void SteadyNSSolver::TrainROMEQPElems(SampleGenerator *sample_generator) { assert(topol_mode == TopologyHandlerMode::COMPONENT); @@ -799,7 +1047,7 @@ void SteadyNSSolver::TrainEQPElems(SampleGenerator *sample_generator) double eqp_tol = config.GetOption("model_reduction/eqp/relative_tolerance", 1.0e-2); /* EQP NNLS for each reference ROM component */ - std::string basis_tag; + BasisTag basis_tag; for (int c = 0; c < num_comp; c++) { int idx = (separate_variable_basis) ? c * num_var : c; @@ -808,6 +1056,50 @@ void SteadyNSSolver::TrainEQPElems(SampleGenerator *sample_generator) const CAROM::Matrix *snapshots = sample_generator->LookUpSnapshot(basis_tag); comp_eqps[c]->TrainEQP(*snapshots, eqp_tol); } + + if (oper_type != OperType::LF) return; + + /* EQP NNLS for interface ROM, for each reference port */ + for (int p = 0; p < topol_handler->GetNumRefPorts(); p++) + { + int c1, c2, a1, a2; + // TODO(kevin): at least component topology handler maintain attrs the same for both reference and subdomain. + // Need to check submesh topology handler. + topol_handler->GetRefPortInfo(p, c1, c2, a1, a2); + + PortTag tag = {.Mesh1 = topol_handler->GetComponentName(c1), + .Mesh2 = topol_handler->GetComponentName(c2), + .Attr1 = a1, .Attr2 = a2}; + + /* Load snapshot matrices for the reference port */ + int idx1 = (separate_variable_basis) ? c1 * num_var : c1; + int idx2 = (separate_variable_basis) ? c2 * num_var : c2; + BasisTag basis_tag1 = rom_handler->GetRefBasisTag(idx1); + BasisTag basis_tag2 = rom_handler->GetRefBasisTag(idx2); + + /* + BasisGenerator::getSnapshotMatrix deletes the existing snapshot matrix, + and creates a new snapshot matrix. + If LookUpSnapshot happens to find the same basis twice, + it will nullify the first pointer. + */ + const CAROM::Matrix *snapshots1 = sample_generator->LookUpSnapshot(basis_tag1); + const CAROM::Matrix *snapshots2 = NULL; + if (basis_tag1 == basis_tag2) + snapshots2 = snapshots1; + else + snapshots2 = sample_generator->LookUpSnapshot(basis_tag2); + + /* Load bases for the reference port */ + DenseMatrix *basis1, *basis2; + rom_handler->GetReferenceBasis(idx1, basis1); + rom_handler->GetReferenceBasis(idx2, basis2); + + /* Load column indices for the reference port */ + Array2D *port_colidx = sample_generator->LookUpSnapshotPortColOffsets(tag); + + itf_eqp->TrainEQPForRefPort(p, *snapshots1, *snapshots2, *port_colidx, eqp_tol); + } // for (int p = 0; p < topol_handler->GetNumRefPorts(); p++) } void SteadyNSSolver::SaveEQPElems(const std::string &filename) @@ -841,13 +1133,30 @@ void SteadyNSSolver::SaveEQPElems(const std::string &filename) assert(comp_eqps[c]); dset_name = topol_handler->GetComponentName(c); - // only one integrator exists in each nonlinear form. - comp_eqps[c]->SaveEQPForIntegrator(IntegratorType::DOMAIN, 0, grp_id, dset_name); + comp_eqps[c]->SaveEQPForIntegrator(IntegratorType::DOMAIN, 0, grp_id, dset_name + "_integ0"); + if (oper_type == OperType::LF) + { + if (full_dg) + comp_eqps[c]->SaveEQPForIntegrator(IntegratorType::INTERIORFACE, 0, grp_id, dset_name + "_integ1"); + + Mesh *comp = topol_handler->GetComponentMesh(c); + int num_bdr = comp->bdr_attributes.Size(); + for (int b = 0; b < num_bdr; b++) + { + comp_eqps[c]->SaveEQPForIntegrator(IntegratorType::BDRFACE, 2 * b, grp_id, + dset_name + "_bdr" + std::to_string(b)); + comp_eqps[c]->SaveEQPForIntegrator(IntegratorType::BDRFACE, 2 * b + 1, grp_id, + dset_name + "_zero" + std::to_string(b)); + } // for (int b = 0; b < num_bdr; b++) + } // if (oper_type == OperType::LF) } // for (int c = 0; c < num_comp; c++) errf = H5Gclose(grp_id); assert(errf >= 0); + if (oper_type == OperType::LF) + itf_eqp->SaveEQPForIntegrator(0, file_id, "interface_integ0"); + errf = H5Fclose(file_id); assert(errf >= 0); } @@ -881,7 +1190,22 @@ void SteadyNSSolver::LoadEQPElems(const std::string &filename) dset_name = topol_handler->GetComponentName(c); // only one integrator exists in each nonlinear form. - comp_eqps[c]->LoadEQPForIntegrator(IntegratorType::DOMAIN, 0, grp_id, dset_name); + comp_eqps[c]->LoadEQPForIntegrator(IntegratorType::DOMAIN, 0, grp_id, dset_name + "_integ0"); + if (oper_type == OperType::LF) + { + if (full_dg) + comp_eqps[c]->LoadEQPForIntegrator(IntegratorType::INTERIORFACE, 0, grp_id, dset_name + "_integ1"); + + Mesh *comp = topol_handler->GetComponentMesh(c); + int num_bdr = comp->bdr_attributes.Size(); + for (int b = 0; b < num_bdr; b++) + { + comp_eqps[c]->LoadEQPForIntegrator(IntegratorType::BDRFACE, 2 * b, grp_id, + dset_name + "_bdr" + std::to_string(b)); + comp_eqps[c]->LoadEQPForIntegrator(IntegratorType::BDRFACE, 2 * b + 1, grp_id, + dset_name + "_zero" + std::to_string(b)); + } // for (int b = 0; b < num_bdr; b++) + } // if (oper_type == OperType::LF) if (comp_eqps[c]->PrecomputeMode()) comp_eqps[c]->PrecomputeCoefficients(); @@ -890,6 +1214,9 @@ void SteadyNSSolver::LoadEQPElems(const std::string &filename) errf = H5Gclose(grp_id); assert(errf >= 0); + if (oper_type == OperType::LF) + itf_eqp->LoadEQPForIntegrator(0, file_id, "interface_integ0"); + errf = H5Fclose(file_id); assert(errf >= 0); } @@ -903,8 +1230,85 @@ void SteadyNSSolver::AssembleROMEQPOper() subdomain_eqps.SetSize(numSub); subdomain_eqps = NULL; + + DenseMatrix *basis; + HyperReductionIntegrator *nl_integ_tmp = NULL; + InterfaceNonlinearFormIntegrator *lf_integ2 = NULL, *lf_bdr = NULL; for (int m = 0; m < numSub; m++) - subdomain_eqps[m] = comp_eqps[rom_handler->GetRefIndexForSubdomain(m)]; + { + int c_type = topol_handler->GetMeshType(m); + int idx = (separate_variable_basis) ? m * num_var : m; + rom_handler->GetDomainBasis(idx, basis); + + subdomain_eqps[m] = new ROMNonlinearForm(basis->NumCols(), ufes[m]); + + switch (oper_type) + { + case (OperType::BASE): + nl_integ_tmp = new VectorConvectionTrilinearFormIntegrator(*zeta_coeff); + nl_integ_tmp->SetIntRule(ir_nl); + subdomain_eqps[m]->AddDomainIntegrator(nl_integ_tmp); + + subdomain_eqps[m]->UpdateDomainIntegratorSampling(0, + *(comp_eqps[c_type]->GetEQPForIntegrator(IntegratorType::DOMAIN, 0))); + break; + + case (OperType::LF): + { + nl_integ_tmp = new IncompressibleInviscidFluxNLFIntegrator(*minus_zeta); + nl_integ_tmp->SetIntRule(ir_nl); + lf_integ2 = new DGLaxFriedrichsFluxIntegrator(*minus_zeta); + lf_integ2->SetIntRule(ir_face); + + subdomain_eqps[m]->AddDomainIntegrator(nl_integ_tmp); + if (full_dg) + subdomain_eqps[m]->AddInteriorFaceIntegrator(lf_integ2); + + subdomain_eqps[m]->UpdateDomainIntegratorSampling(0, + *(comp_eqps[c_type]->GetEQPForIntegrator(IntegratorType::DOMAIN, 0))); + if (full_dg) + subdomain_eqps[m]->UpdateInteriorFaceIntegratorSampling(0, + *(comp_eqps[c_type]->GetEQPForIntegrator(IntegratorType::INTERIORFACE, 0))); + + Array *bdr_c2g = topol_handler->GetBdrAttrComponentToGlobalMap(m); + int idx = 0; + for (int b = 0; b < bdr_c2g->Size(); b++) + { + int global_idx = global_bdr_attributes.Find((*bdr_c2g)[b]); + if (global_idx < 0) continue; + + // TODO: Non-homogeneous Neumann stress bc + if (bdr_type[global_idx] == BoundaryType::NEUMANN) + lf_integ2 = new DGLaxFriedrichsFluxIntegrator(*minus_zeta); + else + { + assert(BCExistsOnBdr(global_idx)); + lf_integ2 = new DGLaxFriedrichsFluxIntegrator(*minus_zeta, ud_coeffs[global_idx]); + } + + lf_integ2->SetIntRule(ir_face); + subdomain_eqps[m]->AddBdrFaceIntegrator(lf_integ2, *bdr_markers[global_idx]); + + int kind = (bdr_type[global_idx] == BoundaryType::ZERO) ? 1 : 0; + subdomain_eqps[m]->UpdateBdrFaceIntegratorSampling(idx, + *(comp_eqps[c_type]->GetEQPForIntegrator(IntegratorType::BDRFACE, 2 * b + kind))); + + idx++; + } // for (int b = 0; b < bdr_c2g->Size(); b++) + } break; + + default: + break; + } // switch (oper_type) + + subdomain_eqps[m]->SetBasis(*basis); + // TODO(kevin): load these coefficients, not re-computing. + subdomain_eqps[m]->SetPrecomputeMode(comp_eqps[c_type]->PrecomputeMode()); + if (subdomain_eqps[m]->PrecomputeMode()) + subdomain_eqps[m]->PrecomputeCoefficients(); + } // for (int m = 0; m < numSub; m++) + + /* ROMInterfaceForm is already loaded */ } DenseTensor* SteadyNSSolver::GetReducedTensor(DenseMatrix *basis, FiniteElementSpace *fespace) diff --git a/src/stokes_solver.cpp b/src/stokes_solver.cpp index 679f37b3..c6d0fc7d 100644 --- a/src/stokes_solver.cpp +++ b/src/stokes_solver.cpp @@ -217,8 +217,6 @@ void StokesSolver::InitVariables() } f_coeffs.SetSize(0); - - if (use_rom) MultiBlockSolver::InitROMHandler(); } void StokesSolver::DeterminePressureDirichlet() @@ -967,7 +965,7 @@ void StokesSolver::LoadSupremizer() for (int m = 0; m < size; m++) { // Load the supremizer basis. - std::string basis_tag = GetBasisTagForComponent(m, topol_handler, "sup"); + BasisTag basis_tag = GetBasisTagForComponent(m, topol_handler, "sup"); /* TODO(kevin): this is a boilerplate for parallel POD/EQP training. In full parallelization, all processes will participate in this supremizer reading procedure. @@ -976,7 +974,7 @@ void StokesSolver::LoadSupremizer() { hid_t file_id; herr_t errf = 0; - std::string filename(basis_prefix + "_" + basis_tag + ".h5"); + std::string filename(basis_prefix + "_" + basis_tag.print() + ".h5"); printf("\nOpening file %s.. ", filename.c_str()); file_id = H5Fopen(filename.c_str(), H5F_ACC_RDONLY, H5P_DEFAULT); assert(file_id >= 0); @@ -1333,7 +1331,7 @@ void StokesSolver::EnrichSupremizer() Orthonormalize(*ubasis, *supreme); // Save the supremizer basis. - std::string basis_tag = GetBasisTagForComponent(m, topol_handler, "sup"); + BasisTag basis_tag = GetBasisTagForComponent(m, topol_handler, "sup"); /* TODO(kevin): this is a boilerplate for parallel POD/EQP training. In full parallelization, all processes will participate in this supremizer writing procedure. @@ -1342,7 +1340,7 @@ void StokesSolver::EnrichSupremizer() { hid_t file_id; herr_t errf = 0; - std::string filename(basis_prefix + "_" + basis_tag + ".h5"); + std::string filename(basis_prefix + "_" + basis_tag.print() + ".h5"); file_id = H5Fcreate(filename.c_str(), H5F_ACC_TRUNC, H5P_DEFAULT, H5P_DEFAULT); assert(file_id >= 0); diff --git a/src/unsteady_ns_solver.cpp b/src/unsteady_ns_solver.cpp index 8d67c7c2..a61e2d1f 100644 --- a/src/unsteady_ns_solver.cpp +++ b/src/unsteady_ns_solver.cpp @@ -45,11 +45,6 @@ UnsteadyNSSolver::~UnsteadyNSSolver() delete Hop; } -void UnsteadyNSSolver::InitVariables() -{ - StokesSolver::InitVariables(); -} - void UnsteadyNSSolver::BuildDomainOperators() { SteadyNSSolver::BuildDomainOperators(); diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index 5722d39d..4c9ba78b 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -62,6 +62,9 @@ file(COPY meshes/test.1x1x1.tet.mesh DESTINATION ${CMAKE_BINARY_DIR}/test/meshes file(COPY meshes/test.2x2x2.tet.mesh DESTINATION ${CMAKE_BINARY_DIR}/test/meshes) file(COPY meshes/test_topol.3d.h5 DESTINATION ${CMAKE_BINARY_DIR}/test/meshes) +file(COPY inputs/test.interface.yml DESTINATION ${CMAKE_BINARY_DIR}/test/inputs) +file(COPY meshes/box-channel.1x1.x_periodic.h5 DESTINATION ${CMAKE_BINARY_DIR}/test/meshes) + add_executable(test_block_smoother test_block_smoother.cpp $) add_executable(nonlinear_integ_grad nonlinear_integ_grad.cpp $) diff --git a/test/gmsh/CMakeLists.txt b/test/gmsh/CMakeLists.txt index c1222e06..c43e48c1 100644 --- a/test/gmsh/CMakeLists.txt +++ b/test/gmsh/CMakeLists.txt @@ -6,8 +6,12 @@ file(COPY square-circle.geo DESTINATION ${CMAKE_BINARY_DIR}/test/gmsh/) file(COPY square.mesh DESTINATION ${CMAKE_BINARY_DIR}/test/gmsh/) file(COPY square.tri.mesh DESTINATION ${CMAKE_BINARY_DIR}/test/gmsh/) file(COPY test.multi_comp.h5 DESTINATION ${CMAKE_BINARY_DIR}/test/gmsh/) +file(COPY box-channel.2x2.periodic.h5 DESTINATION ${CMAKE_BINARY_DIR}/test/gmsh/) +file(COPY box-channel.1x2.h5 DESTINATION ${CMAKE_BINARY_DIR}/test/gmsh/) 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/) ADD_CUSTOM_COMMAND( OUTPUT ${CMAKE_BINARY_DIR}/test/gmsh/square-circle.msh diff --git a/test/gmsh/box-channel.1x1.periodic.h5 b/test/gmsh/box-channel.1x1.periodic.h5 new file mode 100644 index 00000000..e9564942 Binary files /dev/null and b/test/gmsh/box-channel.1x1.periodic.h5 differ diff --git a/test/gmsh/box-channel.1x2.h5 b/test/gmsh/box-channel.1x2.h5 new file mode 100644 index 00000000..a342cf03 Binary files /dev/null and b/test/gmsh/box-channel.1x2.h5 differ diff --git a/test/gmsh/box-channel.2x2.periodic.h5 b/test/gmsh/box-channel.2x2.periodic.h5 new file mode 100644 index 00000000..abaa0555 Binary files /dev/null and b/test/gmsh/box-channel.2x2.periodic.h5 differ diff --git a/test/gmsh/steadyns.interface_eqp.yml b/test/gmsh/steadyns.interface_eqp.yml new file mode 100644 index 00000000..ab2e6e04 --- /dev/null +++ b/test/gmsh/steadyns.interface_eqp.yml @@ -0,0 +1,85 @@ +main: +#mode: run_example/sample_generation/build_rom/single_run + mode: single_run + use_rom: true + solver: steady-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 + +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: "stokes" + parameters: + - key: single_run/periodic_flow_past_array/nu + type: double + sample_size: 2 + minimum: 1.0 + maximum: 1.1 + +sample_collection: + mode: port + port_fileformat: + format: stokes%d_sample.port.h5 + minimum: 0 + maximum: 1 + +basis: + prefix: "stokes" + number_of_basis: 6 + svd: + save_spectrum: true + update_right_sv: false + visualization: + enabled: false + +model_reduction: + separate_variable_basis: true + rom_handler_type: mfem + # individual/universal + subdomain_training: universal + 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: us diff --git a/test/gmsh/steadyns.lf.yml b/test/gmsh/steadyns.lf.yml new file mode 100644 index 00000000..c1d09646 --- /dev/null +++ b/test/gmsh/steadyns.lf.yml @@ -0,0 +1,80 @@ +main: +#mode: run_example/sample_generation/build_rom/single_run + mode: single_run + use_rom: true + solver: steady-ns + +navier-stokes: + operator-type: lf + +mesh: + type: component-wise + component-wise: + global_config: "box-channel.1x2.h5" + components: + - name: "square-circle" + file: "square-circle.msh.mfem" + +domain-decomposition: + type: interior_penalty + +discretization: + order: 1 + full-discrete-galerkin: true + +solver: + direct_solve: true + +visualization: + enable: false + output_dir: dd_mms_output + +parameterized_problem: + name: force_driven_corner + +single_run: + force_driven_corner: + nu: 1.1 + fx: 0. + fy: -0.5 + +sample_generation: + maximum_number_of_snapshots: 400 + component_sampling: false + file_path: + prefix: "stokes" + parameters: + - key: single_run/force_driven_corner/nu + type: double + sample_size: 3 + minimum: 1.0 + maximum: 1.2 + +sample_collection: + mode: port + +basis: + prefix: "stokes" + number_of_basis: 6 + svd: + save_spectrum: true + update_right_sv: false + visualization: + enabled: false + +model_reduction: + separate_variable_basis: true + rom_handler_type: mfem + # individual/universal + subdomain_training: universal + 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: us diff --git a/test/gmsh/test_multi_comp_workflow.cpp b/test/gmsh/test_multi_comp_workflow.cpp index 7a0e605a..0e7355c0 100644 --- a/test/gmsh/test_multi_comp_workflow.cpp +++ b/test/gmsh/test_multi_comp_workflow.cpp @@ -11,7 +11,7 @@ using namespace mfem; static const double threshold = 1.0e-14; static const double stokes_threshold = 2.0e-12; -static const double ns_threshold = 1.0e-9; +static const double ns_threshold = 1.0e-7; /** * Simple smoke test to make sure Google Test is properly linked @@ -219,6 +219,73 @@ TEST(ComponentWiseTest, SteadyNSTest_SeparateVariable_EQP) return; } +TEST(SteadyNS_Workflow, LF) +{ + config = InputParser("steadyns.lf.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"; + 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; +} + +TEST(SteadyNS_Workflow, InterfaceEQP) +{ + config = InputParser("steadyns.interface_eqp.yml"); + + printf("\nSample Generation \n\n"); + + config.dict_["main"]["mode"] = "sample_generation"; + config.dict_["sample_generation"]["file_path"]["prefix"] = "stokes0"; + GenerateSamples(MPI_COMM_WORLD); + + config.dict_["sample_generation"]["file_path"]["prefix"] = "stokes1"; + config.dict_["sample_generation"]["parameters"][0]["sample_size"] = 1; + config.dict_["sample_generation"]["parameters"][0]["minimum"] = 1.2; + config.dict_["sample_generation"]["parameters"][0]["maximum"] = 1.2; + 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 < ns_threshold); + + return; +} + TEST(MultiComponentGlobalROM, StokesTest) { config = InputParser("stokes.component.yml"); diff --git a/test/inputs/test.interface.yml b/test/inputs/test.interface.yml new file mode 100644 index 00000000..ff89ff3f --- /dev/null +++ b/test/inputs/test.interface.yml @@ -0,0 +1,14 @@ +mesh: + type: component-wise + component-wise: + global_config: "meshes/box-channel.1x1.x_periodic.h5" + components: + - name: "empty" + file: "meshes/square.tri.mesh" + +domain-decomposition: + type: interior_penalty + +discretization: + order: 1 + full-discrete-galerkin: true diff --git a/test/meshes/box-channel.1x1.x_periodic.h5 b/test/meshes/box-channel.1x1.x_periodic.h5 new file mode 100644 index 00000000..d80ef454 Binary files /dev/null and b/test/meshes/box-channel.1x1.x_periodic.h5 differ diff --git a/test/test_rom_interfaceform.cpp b/test/test_rom_interfaceform.cpp index dd907955..d48ad7cf 100644 --- a/test/test_rom_interfaceform.cpp +++ b/test/test_rom_interfaceform.cpp @@ -6,6 +6,7 @@ #include "interfaceinteg.hpp" #include "rom_interfaceform.hpp" #include "etc.hpp" +#include "component_topology_handler.hpp" using namespace std; using namespace mfem; @@ -98,10 +99,11 @@ TEST(ROMInterfaceForm, InterfaceAddMult) for (int m = 0; m < numSub; m++) basis[m]->MultTranspose(y.GetBlock(m), Pty.GetBlock(m)); - ROMInterfaceForm *rform = new ROMInterfaceForm(meshes, fes, submesh); + ROMInterfaceForm *rform = new ROMInterfaceForm(meshes, fes, fes, submesh); rform->AddInterfaceIntegrator(integ2); + /* For submesh, mesh index is the same as component index */ for (int m = 0; m < numSub; m++) - rform->SetBasisAtSubdomain(m, *basis[m]); + rform->SetBasisAtComponent(m, *basis[m]); rform->UpdateBlockOffsets(); // we set the full elements/quadrature points, @@ -112,17 +114,14 @@ TEST(ROMInterfaceForm, InterfaceAddMult) for (int p = 0; p < nport; p++) { Array *interface_infos = submesh->GetInterfaceInfos(p); - Array sample_itf(0), sample_qp(0); - Array sample_qw(0); + Array samples(0); for (int itf = 0; itf < interface_infos->Size(); itf++) for (int q = 0; q < nqe; q++) { - sample_itf.Append(itf); - sample_qp.Append(q); - sample_qw.Append(w_el[q]); + samples.Append({.el=itf, .qp=q, .qw=w_el[q]}); } - rform->UpdateInterFaceIntegratorSampling(0, p, sample_itf, sample_qp, sample_qw); + rform->UpdateInterFaceIntegratorSampling(0, p, samples); } rom_y = 0.0; @@ -204,10 +203,11 @@ TEST(ROMInterfaceForm, InterfaceGetGradient) BlockVector rom_y(rom_block_offsets); - ROMInterfaceForm *rform = new ROMInterfaceForm(meshes, fes, submesh); + ROMInterfaceForm *rform = new ROMInterfaceForm(meshes, fes, fes, submesh); rform->AddInterfaceIntegrator(integ2); + /* For submesh, mesh index is the same as component index */ for (int m = 0; m < numSub; m++) - rform->SetBasisAtSubdomain(m, *basis[m]); + rform->SetBasisAtComponent(m, *basis[m]); rform->UpdateBlockOffsets(); // we set the full elements/quadrature points, @@ -218,17 +218,14 @@ TEST(ROMInterfaceForm, InterfaceGetGradient) for (int p = 0; p < nport; p++) { Array *interface_infos = submesh->GetInterfaceInfos(p); - Array sample_itf(0), sample_qp(0); - Array sample_qw(0); + Array samples(0); for (int itf = 0; itf < interface_infos->Size(); itf++) for (int q = 0; q < nqe; q++) { - sample_itf.Append(itf); - sample_qp.Append(q); - sample_qw.Append(w_el[q]); + samples.Append({.el=itf, .qp=q, .qw=w_el[q]}); } - rform->UpdateInterFaceIntegratorSampling(0, p, sample_itf, sample_qp, sample_qw); + rform->UpdateInterFaceIntegratorSampling(0, p, samples); } rom_y = 0.0; @@ -349,11 +346,11 @@ TEST(ROMInterfaceForm, SetupEQPSystem_for_a_port) snapshots2(i, j) = 2.0 * UniformRandom() - 1.0; /* fictitious pair indices */ - Array snap1_idx(num_snap), snap2_idx(num_snap); + Array2D snap_pair_idx(num_snap, 2); for (int i = 0; i < num_snap; i++) { - snap1_idx[i] = UniformRandom(0, nsnap1-1); - snap2_idx[i] = UniformRandom(0, nsnap2-1); + snap_pair_idx(i, 0) = UniformRandom(0, nsnap1-1); + snap_pair_idx(i, 1) = UniformRandom(0, nsnap2-1); } /* bases generated from fictitious snapshots */ @@ -421,18 +418,19 @@ TEST(ROMInterfaceForm, SetupEQPSystem_for_a_port) nform->AddInterfaceIntegrator(integ1); /* EQP interface operator */ - ROMInterfaceForm *rform(new ROMInterfaceForm(meshes, fes, submesh)); + ROMInterfaceForm *rform(new ROMInterfaceForm(meshes, fes, fes, submesh)); rform->AddInterfaceIntegrator(integ2); Array bases(numSub); + /* For submesh, mesh index is the same as component index */ for (int m = 0; m < bases.Size(); m++) { bases[m] = new DenseMatrix(fes[m]->GetTrueVSize(), nb_def); *bases[m] = 0.0; - rform->SetBasisAtSubdomain(m, *bases[m]); + rform->SetBasisAtComponent(m, *bases[m]); } - rform->SetBasisAtSubdomain(midx[0], basis1); - rform->SetBasisAtSubdomain(midx[1], basis2); + rform->SetBasisAtComponent(midx[0], basis1); + rform->SetBasisAtComponent(midx[1], basis2); rform->UpdateBlockOffsets(); // rform->SetPrecomputeMode(false); @@ -446,8 +444,8 @@ TEST(ROMInterfaceForm, SetupEQPSystem_for_a_port) Vector basis_col; for (int s = 0; s < num_snap; s++) { - snapshots1.GetColumnReference(snap1_idx[s], snap1); - snapshots2.GetColumnReference(snap2_idx[s], snap2); + snapshots1.GetColumnReference(snap_pair_idx(s, 0), snap1); + snapshots2.GetColumnReference(snap_pair_idx(s, 1), snap2); rhs_vec1 = 0.0; rhs_vec2 = 0.0; nform->AssembleInterfaceVector(mesh1, mesh2, fes1, fes2, itf_infos, snap1, snap2, rhs_vec1, rhs_vec2); @@ -475,27 +473,26 @@ TEST(ROMInterfaceForm, SetupEQPSystem_for_a_port) /* equivalent operation must happen within this routine */ rform->SetupEQPSystem(carom_snapshots1_work, carom_snapshots2_work, - snap1_idx, snap2_idx, basis1, basis2, 0, 0, + snap_pair_idx, basis1, basis2, 0, 0, fes1, fes2, itf_infos, integ2, Gt, rhs2); for (int k = 0; k < rhs1.dim(); k++) EXPECT_NEAR(rhs1(k), rhs2(k), threshold); double eqp_tol = 1.0e-10; - Array sample_el(0), sample_qp(0); - Array sample_qw(0); + Array samples(0); const int nqe = ir->GetNPoints(); - rform->TrainEQPForIntegrator(nqe, Gt, rhs2, eqp_tol, sample_el, sample_qp, sample_qw); + rform->TrainEQPForIntegrator(nqe, Gt, rhs2, eqp_tol, samples); // if (rform->PrecomputeMode()) rform->PrecomputeCoefficients(); - rform->UpdateInterFaceIntegratorSampling(0, pidx, sample_el, sample_qp, sample_qw); + rform->UpdateInterFaceIntegratorSampling(0, pidx, samples); DenseMatrix rom_rhs1(rhs1.getData(), NB, num_snap), rom_rhs2(NB, num_snap); Array rom_blocks = rform->GetBlockOffsets(); BlockVector rom_sol(rom_blocks), rom_rhs2_vec(rom_blocks); for (int s = 0; s < num_snap; s++) { - snapshots1.GetColumnReference(snap1_idx[s], snap1); - snapshots2.GetColumnReference(snap2_idx[s], snap2); + snapshots1.GetColumnReference(snap_pair_idx(s, 0), snap1); + snapshots2.GetColumnReference(snap_pair_idx(s, 1), snap2); rom_sol = 0.0; basis1.MultTranspose(snap1, rom_sol.GetBlock(midx[0])); @@ -522,6 +519,219 @@ TEST(ROMInterfaceForm, SetupEQPSystem_for_a_port) delete basis2_generator; } +// // TODO(kevin): figure out how to set up interface for self. +// TEST(ROMInterfaceForm, Self_Interface) +// { +// config = InputParser("inputs/test.interface.yml"); +// const int order = UniformRandom(1, 1); + +// ComponentTopologyHandler *topol = new ComponentTopologyHandler(); + +// Array meshes; +// TopologyData topol_data; +// topol->ExportInfo(meshes, topol_data); +// const int dim = topol_data.dim; +// const int numSub = topol_data.numSub; + +// const int pidx = UniformRandom(0, topol->GetNumPorts()-1); +// const PortInfo *pInfo = topol->GetPortInfo(pidx); +// Array* const itf_infos = topol->GetInterfaceInfos(pidx); + +// Array midx(2); +// Mesh *mesh1, *mesh2; + +// midx[0] = pInfo->Mesh1; +// midx[1] = pInfo->Mesh2; +// assert(midx[0] == midx[1]); + +// mesh1 = meshes[midx[0]]; +// mesh2 = meshes[midx[1]]; + +// FiniteElementCollection *dg_coll(new DG_FECollection(order, dim)); +// Array fes(numSub); +// for (int m = 0; m < numSub; m++) +// fes[m] = new FiniteElementSpace(meshes[m], dg_coll, dim); + +// FiniteElementSpace *fes1, *fes2; +// fes1 = fes[midx[0]]; +// fes2 = fes[midx[1]]; + +// const int ndofs1 = fes1->GetTrueVSize(); +// const int ndofs2 = fes2->GetTrueVSize(); +// assert(ndofs1 == ndofs2); + +// /* number of snapshots on the domain */ +// /* we only have only one domain */ +// const int nsnap = UniformRandom(3, 5); +// /* same number of basis to fully represent snapshots */ +// const int num_basis = nsnap; +// const int NB = num_basis + num_basis; +// const int num_snap = nsnap; + +// /* fictitious snapshots */ +// DenseMatrix snapshots(ndofs1, nsnap); +// for (int i = 0; i < ndofs1; i++) +// for (int j = 0; j < nsnap; j++) +// snapshots(i, j) = 2.0 * UniformRandom() - 1.0; + +// /* fictitious pair indices */ +// Array2D snap_pair_idx(num_snap, 2); +// for (int i = 0; i < num_snap; i++) +// { +// snap_pair_idx(i, 0) = i; +// snap_pair_idx(i, 1) = i; +// } + +// /* bases generated from fictitious snapshots */ +// DenseMatrix basis(ndofs1, num_basis); +// const CAROM::Matrix *carom_snapshots; +// CAROM::BasisGenerator *basis_generator; +// { +// CAROM::Options options(ndofs1, nsnap, 1, true); +// options.static_svd_preserve_snapshot = true; +// basis_generator = new CAROM::BasisGenerator(options, false, "test_basis"); +// Vector snapshot(ndofs1); +// for (int s = 0; s < nsnap; s++) +// { +// snapshots.GetColumnReference(s, snapshot); +// basis_generator->takeSample(snapshot.GetData()); +// } +// basis_generator->endSamples(); +// carom_snapshots = basis_generator->getSnapshotMatrix(); + +// CAROM::BasisReader basis_reader("test_basis"); +// const CAROM::Matrix *carom_basis = basis_reader.getSpatialBasis(num_basis); +// CAROM::CopyMatrix(*carom_basis, basis); +// } + +// /* get integration rule for eqp */ +// const IntegrationRule *ir = NULL; +// for (int be = 0; be < fes1->GetNBE(); be++) +// { +// FaceElementTransformations *tr = mesh1->GetBdrFaceTransformations(be); +// if (tr != NULL) +// { +// ir = &IntRules.Get(tr->GetGeometryType(), +// (int)(ceil(1.5 * (2 * fes1->GetMaxElementOrder() - 1)))); +// break; +// } +// } +// assert(ir); + +// ConstantCoefficient pi(3.141592); +// auto *integ1 = new DGLaxFriedrichsFluxIntegrator(pi); +// integ1->SetIntRule(ir); +// auto *integ2 = new DGLaxFriedrichsFluxIntegrator(pi); +// integ2->SetIntRule(ir); + +// /* FOM interface operator */ +// InterfaceForm *nform(new InterfaceForm(meshes, fes, topol)); +// nform->AddInterfaceIntegrator(integ1); + +// /* EQP interface operator */ +// ROMInterfaceForm *rform(new ROMInterfaceForm(meshes, fes, fes, topol)); +// rform->AddInterfaceIntegrator(integ2); + +// rform->SetBasisAtComponent(midx[0], basis); +// rform->UpdateBlockOffsets(); +// // rform->SetPrecomputeMode(false); + +// CAROM::Vector rhs1(num_snap * NB, false); +// CAROM::Vector rhs2(num_snap * NB, false); +// CAROM::Matrix Gt(1, 1, true); + +// /* exact right-hand side by inner product of basis and fom vectors */ +// Vector rhs_vec1(ndofs1), rhs_vec2(ndofs2); +// Vector snap1(ndofs1), snap2(ndofs2); +// Vector basis_col; +// for (int s = 0; s < num_snap; s++) +// { +// snapshots.GetColumnReference(snap_pair_idx(s, 0), snap1); +// snapshots.GetColumnReference(snap_pair_idx(s, 1), snap2); + +// rhs_vec1 = 0.0; rhs_vec2 = 0.0; +// nform->AssembleInterfaceVector(mesh1, mesh2, fes1, fes2, itf_infos, snap1, snap2, rhs_vec1, rhs_vec2); + +// for (int b = 0; b < num_basis; b++) +// { +// basis.GetColumnReference(b, basis_col); +// rhs1(b + s * NB) = basis_col * rhs_vec1; +// } +// for (int b = 0; b < num_basis; b++) +// { +// basis.GetColumnReference(b, basis_col); +// rhs1(b + num_basis + s * NB) = basis_col * rhs_vec2; +// } +// } + +// /* +// TODO(kevin): this is a boilerplate for parallel POD/EQP training. +// Will have to consider parallel-compatible test. +// */ +// CAROM::Matrix carom_snapshots1_work(*carom_snapshots); +// carom_snapshots1_work.gather(); +// CAROM::Matrix carom_snapshots2_work(*carom_snapshots); +// carom_snapshots2_work.gather(); + +// /* equivalent operation must happen within this routine */ +// rform->SetupEQPSystem(carom_snapshots1_work, carom_snapshots2_work, +// snap_pair_idx, basis, basis, 0, 0, +// fes1, fes2, itf_infos, integ2, Gt, rhs2); + +// for (int k = 0; k < rhs1.dim(); k++) +// EXPECT_NEAR(rhs1(k), rhs2(k), threshold); + +// double eqp_tol = 1.0e-10; +// Array samples(0); +// const int nqe = ir->GetNPoints(); +// rform->TrainEQPForIntegrator(nqe, Gt, rhs2, eqp_tol, samples); +// // if (rform->PrecomputeMode()) rform->PrecomputeCoefficients(); +// rform->UpdateInterFaceIntegratorSampling(0, pidx, samples); + +// DenseMatrix rom_rhs1(rhs1.getData(), NB, num_snap), rom_rhs2(NB, num_snap); +// Array rom_blocks = rform->GetBlockOffsets(); +// BlockVector rom_sol(rom_blocks), rom_rhs2_vec(rom_blocks); +// for (int s = 0; s < num_snap; s++) +// { +// snapshots.GetColumnReference(snap_pair_idx(s, 0), snap1); +// snapshots.GetColumnReference(snap_pair_idx(s, 1), snap2); + +// rom_sol = 0.0; +// basis.MultTranspose(snap1, rom_sol.GetBlock(midx[0])); +// basis.MultTranspose(snap2, rom_sol.GetBlock(midx[1])); +// { +// printf("rom_sol:\n"); +// for (int k = 0; k < rom_sol.Size(); k++) +// printf("%.3e\t", rom_sol(k)); +// printf("\n"); +// } +// rom_rhs2_vec = 0.0; +// rform->InterfaceAddMult(rom_sol, rom_rhs2_vec); +// { +// printf("rom_rhs2_vec:\n"); +// for (int k = 0; k < rom_rhs2_vec.Size(); k++) +// printf("%.3e\t", rom_rhs2_vec(k)); +// printf("\n"); +// } + +// for (int b = 0; b < num_basis; b++) +// rom_rhs2(b, s) = rom_rhs2_vec.GetBlock(midx[0])(b); +// for (int b = 0; b < num_basis; b++) +// rom_rhs2(b + num_basis, s) = rom_rhs2_vec.GetBlock(midx[1])(b); +// } + +// for (int i = 0; i < NB; i++) +// for (int j = 0; j < num_snap; j++) +// EXPECT_NEAR(rom_rhs1(i, j), rom_rhs2(i, j), threshold); + +// delete topol; +// delete dg_coll; +// DeletePointers(fes); +// delete nform; +// delete rform; +// delete basis_generator; +// } + int main(int argc, char* argv[]) { MPI_Init(&argc, &argv); diff --git a/test/test_rom_nonlinearform.cpp b/test/test_rom_nonlinearform.cpp index f5f51343..b32f8f57 100644 --- a/test/test_rom_nonlinearform.cpp +++ b/test/test_rom_nonlinearform.cpp @@ -58,16 +58,15 @@ TEST(ROMNonlinearForm, VectorConvectionTrilinearFormIntegrator) const int nqe = ir.GetNPoints(); const int ne = fes->GetNE(); Array const& w_el = ir.GetWeights(); - Array sample_el(ne * nqe), sample_qp(ne * nqe); - Array sample_qw(ne * nqe); + Array samples(ne * nqe); for (int e = 0, idx = 0; e < ne; e++) for (int q = 0; q < nqe; q++, idx++) { - sample_el[idx] = e; - sample_qp[idx] = q; - sample_qw[idx] = w_el[q]; + samples[idx].el = e; + samples[idx].qp = q; + samples[idx].qw = w_el[q]; } - rform->UpdateDomainIntegratorSampling(0, sample_el, sample_qp, sample_qw); + rform->UpdateDomainIntegratorSampling(0, samples); Vector rom_u(num_basis), u(fes->GetTrueVSize()); for (int k = 0; k < rom_u.Size(); k++) @@ -128,16 +127,15 @@ TEST(ROMNonlinearForm, IncompressibleInviscidFluxNLFIntegrator) const int nqe = ir.GetNPoints(); const int ne = fes->GetNE(); Array const& w_el = ir.GetWeights(); - Array sample_el(ne * nqe), sample_qp(ne * nqe); - Array sample_qw(ne * nqe); + Array samples(ne * nqe); for (int e = 0, idx = 0; e < ne; e++) for (int q = 0; q < nqe; q++, idx++) { - sample_el[idx] = e; - sample_qp[idx] = q; - sample_qw[idx] = w_el[q]; + samples[idx].el = e; + samples[idx].qp = q; + samples[idx].qw = w_el[q]; } - rform->UpdateDomainIntegratorSampling(0, sample_el, sample_qp, sample_qw); + rform->UpdateDomainIntegratorSampling(0, samples); Vector rom_u(num_basis), u(fes->GetTrueVSize()); for (int k = 0; k < rom_u.Size(); k++) @@ -199,8 +197,7 @@ TEST(ROMNonlinearForm, DGLaxFriedrichsFluxIntegrator) const int nqe = ir.GetNPoints(); const int nf = mesh->GetNumFaces(); Array const& w_el = ir.GetWeights(); - Array sample_el(0), sample_qp(0); - Array sample_qw(0); + Array samples(0); for (int f = 0; f < nf; f++) { tr = mesh->GetInteriorFaceTransformations(f); @@ -208,12 +205,10 @@ TEST(ROMNonlinearForm, DGLaxFriedrichsFluxIntegrator) for (int q = 0; q < nqe; q++) { - sample_el.Append(f); - sample_qp.Append(q); - sample_qw.Append(w_el[q]); + samples.Append({.el = f, .qp = q, .qw = w_el[q]}); } } - rform->UpdateInteriorFaceIntegratorSampling(0, sample_el, sample_qp, sample_qw); + rform->UpdateInteriorFaceIntegratorSampling(0, samples); Vector rom_u(num_basis), u(fes->GetTrueVSize()); for (int k = 0; k < rom_u.Size(); k++) @@ -270,15 +265,14 @@ TEST(ROMNonlinearForm_gradient, VectorConvectionTrilinearFormIntegrator) const int nqe = ir.GetNPoints(); const int ne = fes->GetNE(); Array const& w_el = ir.GetWeights(); - Array sample_el(nsample), sample_qp(nsample); - Array sample_qw(nsample); + Array samples(nsample); for (int s = 0; s < nsample; s++) { - sample_el[s] = UniformRandom(0, ne-1); - sample_qp[s] = UniformRandom(0, nqe-1); - sample_qw[s] = UniformRandom(); + samples[s].el = UniformRandom(0, ne-1); + samples[s].qp = UniformRandom(0, nqe-1); + samples[s].qw = UniformRandom(); } - rform->UpdateDomainIntegratorSampling(0, sample_el, sample_qp, sample_qw); + rform->UpdateDomainIntegratorSampling(0, samples); Vector rom_u(num_basis); for (int k = 0; k < rom_u.Size(); k++) @@ -359,15 +353,14 @@ TEST(ROMNonlinearForm_gradient, IncompressibleInviscidFluxNLFIntegrator) const int nqe = ir.GetNPoints(); const int ne = fes->GetNE(); Array const& w_el = ir.GetWeights(); - Array sample_el(nsample), sample_qp(nsample); - Array sample_qw(nsample); + Array samples(nsample); for (int s = 0; s < nsample; s++) { - sample_el[s] = UniformRandom(0, ne-1); - sample_qp[s] = UniformRandom(0, nqe-1); - sample_qw[s] = UniformRandom(); + samples[s].el = UniformRandom(0, ne-1); + samples[s].qp = UniformRandom(0, nqe-1); + samples[s].qw = UniformRandom(); } - rform->UpdateDomainIntegratorSampling(0, sample_el, sample_qp, sample_qw); + rform->UpdateDomainIntegratorSampling(0, samples); Vector rom_u(num_basis); for (int k = 0; k < rom_u.Size(); k++) @@ -448,8 +441,7 @@ TEST(ROMNonlinearForm_gradient, DGLaxFriedrichsFluxIntegrator) const int nqe = ir.GetNPoints(); const int nf = mesh->GetNumFaces(); Array const& w_el = ir.GetWeights(); - Array sample_el(nsample), sample_qp(nsample); - Array sample_qw(nsample); + Array samples(nsample); int s = 0; FaceElementTransformations *tr; @@ -459,13 +451,13 @@ TEST(ROMNonlinearForm_gradient, DGLaxFriedrichsFluxIntegrator) tr = mesh->GetInteriorFaceTransformations(f); if (tr == NULL) continue; - sample_el[s] = f; - sample_qp[s] = UniformRandom(0, nqe-1); - sample_qw[s] = UniformRandom(); + samples[s].el = f; + samples[s].qp = UniformRandom(0, nqe-1); + samples[s].qw = UniformRandom(); s++; } - rform->UpdateInteriorFaceIntegratorSampling(0, sample_el, sample_qp, sample_qw); + rform->UpdateInteriorFaceIntegratorSampling(0, samples); Vector rom_u(num_basis); for (int k = 0; k < rom_u.Size(); k++) @@ -547,15 +539,14 @@ TEST(ROMNonlinearForm_fast, VectorConvectionTrilinearFormIntegrator) const int nqe = ir.GetNPoints(); const int ne = fes->GetNE(); Array const& w_el = ir.GetWeights(); - Array sample_el(nsample), sample_qp(nsample); - Array sample_qw(nsample); + Array samples(nsample); for (int s = 0; s < nsample; s++) { - sample_el[s] = UniformRandom(0, ne-1); - sample_qp[s] = UniformRandom(0, nqe-1); - sample_qw[s] = UniformRandom(); + samples[s].el = UniformRandom(0, ne-1); + samples[s].qp = UniformRandom(0, nqe-1); + samples[s].qw = UniformRandom(); } - rform->UpdateDomainIntegratorSampling(0, sample_el, sample_qp, sample_qw); + rform->UpdateDomainIntegratorSampling(0, samples); rform->PrecomputeCoefficients(); Vector rom_u(num_basis); @@ -644,15 +635,14 @@ TEST(ROMNonlinearForm_fast, IncompressibleInviscidFluxNLFIntegrator) const int nqe = ir.GetNPoints(); const int ne = fes->GetNE(); Array const& w_el = ir.GetWeights(); - Array sample_el(nsample), sample_qp(nsample); - Array sample_qw(nsample); + Array samples(nsample); for (int s = 0; s < nsample; s++) { - sample_el[s] = UniformRandom(0, ne-1); - sample_qp[s] = UniformRandom(0, nqe-1); - sample_qw[s] = UniformRandom(); + samples[s].el = UniformRandom(0, ne-1); + samples[s].qp = UniformRandom(0, nqe-1); + samples[s].qw = UniformRandom(); } - rform->UpdateDomainIntegratorSampling(0, sample_el, sample_qp, sample_qw); + rform->UpdateDomainIntegratorSampling(0, samples); rform->PrecomputeCoefficients(); Vector rom_u(num_basis); @@ -741,8 +731,11 @@ TEST(ROMNonlinearForm_fast, DGLaxFriedrichsFluxIntegrator) assert(ir); ConstantCoefficient pi(3.141592); + Vector ud(dim); + for (int d = 0; d < dim; d++) ud(d) = 2.0 * UniformRandom() - 1.0; + VectorConstantCoefficient ud_coeff(ud); auto *integ = new DGLaxFriedrichsFluxIntegrator(pi); - auto *integ_bdr = new DGLaxFriedrichsFluxIntegrator(pi); + auto *integ_bdr = new DGLaxFriedrichsFluxIntegrator(pi, &ud_coeff); integ->SetIntRule(ir); integ_bdr->SetIntRule(ir); @@ -758,8 +751,7 @@ TEST(ROMNonlinearForm_fast, DGLaxFriedrichsFluxIntegrator) const int nf = mesh->GetNumFaces(); const int nbe = fes->GetNBE(); Array const& w_el = ir->GetWeights(); - Array sample_el(nsample), sample_qp(nsample); - Array sample_qw(nsample); + Array samples(nsample); FaceElementTransformations *tr; @@ -771,13 +763,13 @@ TEST(ROMNonlinearForm_fast, DGLaxFriedrichsFluxIntegrator) tr = mesh->GetInteriorFaceTransformations(f); if (tr == NULL) continue; - sample_el[s] = f; - sample_qp[s] = UniformRandom(0, nqe-1); - sample_qw[s] = UniformRandom(); + samples[s].el = f; + samples[s].qp = UniformRandom(0, nqe-1); + samples[s].qw = UniformRandom(); s++; } - rform->UpdateInteriorFaceIntegratorSampling(0, sample_el, sample_qp, sample_qw); + rform->UpdateInteriorFaceIntegratorSampling(0, samples); /* fill up boundary face samples */ s = 0; @@ -787,13 +779,13 @@ TEST(ROMNonlinearForm_fast, DGLaxFriedrichsFluxIntegrator) tr = mesh->GetBdrFaceTransformations(be); if (tr == NULL) continue; - sample_el[s] = be; - sample_qp[s] = UniformRandom(0, nqe-1); - sample_qw[s] = UniformRandom(); + samples[s].el = be; + samples[s].qp = UniformRandom(0, nqe-1); + samples[s].qw = UniformRandom(); s++; } - rform->UpdateBdrFaceIntegratorSampling(0, sample_el, sample_qp, sample_qw); + rform->UpdateBdrFaceIntegratorSampling(0, samples); rform->PrecomputeCoefficients(); @@ -1107,9 +1099,13 @@ TEST(ROMNonlinearForm, SetupEQPSystemForBdrFaceIntegrator) IntegrationRule ir = IntRules.Get(fes->GetFE(0)->GetGeomType(), (int)(ceil(1.5 * (2 * fes->GetMaxElementOrder() - 1)))); ConstantCoefficient pi(3.141592); - auto *integ1 = new DGLaxFriedrichsFluxIntegrator(pi); + Vector ud(dim); + for (int d = 0; d < dim; d++) ud(d) = 2.0 * UniformRandom() - 1.0; + VectorConstantCoefficient ud_coeff(ud); + + auto *integ1 = new DGLaxFriedrichsFluxIntegrator(pi, &ud_coeff); integ1->SetIntRule(&ir); - auto *integ2 = new DGLaxFriedrichsFluxIntegrator(pi); + auto *integ2 = new DGLaxFriedrichsFluxIntegrator(pi, &ud_coeff); integ2->SetIntRule(&ir); NonlinearForm *nform(new NonlinearForm(fes)); diff --git a/utils/python/box_channel_config.py b/utils/python/box_channel_config.py index b036e226..ce3db3b9 100644 --- a/utils/python/box_channel_config.py +++ b/utils/python/box_channel_config.py @@ -33,9 +33,9 @@ class BoxChannelConfig(Configuration): (0, 1): 3, (-1, 0): 4} comp_used = [] - periodic = False + periodic = [False, False] - def __init__(self, nx_, ny_, periodic=False): + def __init__(self, nx_, ny_, periodic=[False, False]): Configuration.__init__(self) self.nx, self.ny = nx_, ny_ self.nmesh = self.nx * self.ny @@ -98,19 +98,20 @@ def interfaceForPair(self, midx1, midx2): # respective position between two meshes. # if periodic, more possible respective positions are available. disp12 = [loc2 - loc1] - if (self.periodic): - from copy import deepcopy - disp0 = deepcopy(disp12[0]) + from copy import deepcopy + disp0 = deepcopy(disp12[0]) + if (self.periodic[0]): if (disp0[0] == (self.nx-1)): disp12 += [deepcopy(disp0)] disp12[-1][0] -= self.nx - if (disp0[0] == -(self.nx-1)): + if ((disp0[0] == -(self.nx-1)) and not (self.nx == 1)): disp12 += [deepcopy(disp0)] disp12[-1][0] += self.nx + if (self.periodic[1]): if (disp0[1] == (self.ny-1)): disp12 += [deepcopy(disp0)] disp12[-1][1] -= self.ny - if (disp0[1] == -(self.ny-1)): + if ((disp0[1] == -(self.ny-1)) and not (self.ny == 1)): disp12 += [deepcopy(disp0)] disp12[-1][1] += self.ny