From 704e18125c56801bae4ac5370c824e31ecdcff4e Mon Sep 17 00:00:00 2001 From: Arne Morten Kvarving Date: Wed, 11 Nov 2020 11:47:52 +0100 Subject: [PATCH] changed: cache L2 projection matrix --- src/ASM/ASMbase.C | 9 ++++ src/ASM/ASMbase.h | 2 + src/ASM/GlbL2projector.C | 110 +++++++++++++++++++++------------------ src/ASM/GlbL2projector.h | 13 +++-- src/LinAlg/PETScMatrix.C | 43 ++++++++------- src/LinAlg/PETScMatrix.h | 1 + 6 files changed, 101 insertions(+), 77 deletions(-) diff --git a/src/ASM/ASMbase.C b/src/ASM/ASMbase.C index b4f917aa5..442e45ac1 100644 --- a/src/ASM/ASMbase.C +++ b/src/ASM/ASMbase.C @@ -16,6 +16,7 @@ #include "ASM3D.h" #include "IFEM.h" #include "MPC.h" +#include "SparseMatrix.h" #include "Tensor.h" #include "Vec3.h" #include "Vec3Oper.h" @@ -63,6 +64,7 @@ ASMbase::ASMbase (unsigned char n_p, unsigned char n_s, unsigned char n_f) nel = nnod = 0; idx = 0; firstIp = 0; + glbL2_A = nullptr; } @@ -80,6 +82,7 @@ ASMbase::ASMbase (const ASMbase& patch, unsigned char n_f) nnod = patch.nnod; idx = patch.idx; firstIp = patch.firstIp; + glbL2_A = nullptr; // Note: Properties are _not_ copied } @@ -96,6 +99,7 @@ ASMbase::ASMbase (const ASMbase& patch) nnod = patch.nnod; idx = patch.idx; firstIp = patch.firstIp; + glbL2_A = nullptr; // Only copy the regular part of the FE data, leave out any extraordinaries @@ -127,6 +131,8 @@ ASMbase::~ASMbase () { for (MPC* mpc : mpcs) delete mpc; + + delete glbL2_A; } @@ -174,6 +180,9 @@ void ASMbase::clear (bool retainGeometry) BCode.clear(); dCode.clear(); mpcs.clear(); + + delete glbL2_A; + glbL2_A = nullptr; } diff --git a/src/ASM/ASMbase.h b/src/ASM/ASMbase.h index c6bba4e4c..dfef35463 100644 --- a/src/ASM/ASMbase.h +++ b/src/ASM/ASMbase.h @@ -918,6 +918,8 @@ class ASMbase static int gEl; //!< Global element counter static int gNod; //!< Global node counter + mutable SparseMatrix* glbL2_A = nullptr; //!< Global L2-projection (mass) matrix + private: std::vector myLMTypes; //!< Type of Lagrange multiplier ('L' or 'G') std::set myLMs; //!< Nodal indices of the Lagrange multipliers diff --git a/src/ASM/GlbL2projector.C b/src/ASM/GlbL2projector.C index b25fca2ae..6ce53dc8a 100644 --- a/src/ASM/GlbL2projector.C +++ b/src/ASM/GlbL2projector.C @@ -166,7 +166,7 @@ public: L2Mats(GlbL2& p, size_t nen, size_t nf, LocalIntegral* q = nullptr) : gl2Int(p), elmData(q) { - this->resize(1,nf); + this->resize(p.pA ? 1 : 0,nf); this->redim(nen); } @@ -193,7 +193,8 @@ class L2GlobalInt : public GlobalIntegral { public: //! \brief The constructor initializes the system matrix references. - L2GlobalInt(SparseMatrix& A_, StdVector& B_) : A(A_), B(B_) {} + L2GlobalInt(int dim_, SparseMatrix* A_, StdVector& B_) + : dim(dim_), A(A_), B(B_) {} //! \brief Empty destructor. virtual ~L2GlobalInt() {} @@ -205,41 +206,47 @@ public: for (size_t i = 0; i < elm->mnpc.size(); i++) { int inod = elm->mnpc[i]+1; - for (size_t j = 0; j < elm->mnpc.size(); j++) - { - int jnod = elm->mnpc[j]+1; - A(inod,jnod) += elm->A.front()(i+1,j+1); + if (A) { + for (size_t j = 0; j < elm->mnpc.size(); j++) + { + int jnod = elm->mnpc[j]+1; + (*A)(inod,jnod) += elm->A.front()(i+1,j+1); + } } for (const Vector& b : elm->b) { B(inod) += b[i]; - inod += A.dim(); + inod += dim; } } return true; } private: - SparseMatrix& A; //!< Reference to left-hand-side matrix + int dim; //!< Dimension of matrix + SparseMatrix* A; //!< Reference to left-hand-side matrix StdVector& B; //!< Reference to right-hand-side vector }; -GlbL2::GlbL2 (IntegrandBase* p, size_t n) : problem(p) +GlbL2::GlbL2 (IntegrandBase* p, size_t n, SparseMatrix* A) : + pA(A), problem(p) { nrhs = p->getNoFields(2); this->allocate(n); } -GlbL2::GlbL2 (FunctionBase* f, size_t n) : problem(nullptr), functions({f}) +GlbL2::GlbL2 (FunctionBase* f, size_t n, SparseMatrix* A) : + pA(A), problem(nullptr), functions({f}) { nrhs = f->dim(); this->allocate(n); } -GlbL2::GlbL2 (const FunctionVec& f, size_t n) : problem(nullptr), functions(f) +GlbL2::GlbL2 (const FunctionVec& f, size_t n, SparseMatrix* A) : + pA(A), problem(nullptr), functions(f) { nrhs = 0; for (FunctionBase* func : f) @@ -250,32 +257,32 @@ GlbL2::GlbL2 (const FunctionVec& f, size_t n) : problem(nullptr), functions(f) GlbL2::~GlbL2() { - delete pA; delete pB; -#ifdef HAS_PETSC - delete adm; -#endif } void GlbL2::allocate (size_t n) { #ifdef HAS_PETSC - adm = nullptr; if (GlbL2::MatrixType == LinAlg::PETSC && GlbL2::SolverParams) { - adm = new ProcessAdm(); - pA = new PETScMatrix(*adm,*GlbL2::SolverParams); - pB = new PETScVector(*adm,n*nrhs); + static ProcessAdm adm; + if (!pA) { + pA = new PETScMatrix(adm,*GlbL2::SolverParams); + pA->redim(n,n); + } + + pB = new PETScVector(adm,n*nrhs); } else #endif { - pA = new SparseMatrix(SparseMatrix::SUPERLU); + if (!pA) { + pA = new SparseMatrix(SparseMatrix::SUPERLU); + pA->redim(n,n); + } pB = new StdVector(n*nrhs); } - - pA->redim(n,n); } @@ -353,7 +360,8 @@ bool GlbL2::evalInt (LocalIntegral& elmInt, solPt.insert(solPt.end(),funcPt.begin(),funcPt.end()); } - gl2.A.front().outer_product(fe.N,fe.N,true,fe.detJxW); + if (!gl2.A.empty()) + gl2.A.front().outer_product(fe.N,fe.N,true,fe.detJxW); for (size_t j = 0; j < solPt.size(); j++) gl2.b[j].add(fe.N,solPt[j]*fe.detJxW); @@ -376,7 +384,8 @@ bool GlbL2::evalIntMx (LocalIntegral& elmInt, if (!problem->diverged(fe.iGP+1)) return false; - gl2.A.front().outer_product(fe.N,fe.N,true,fe.detJxW); + if (!gl2.A.empty()) + gl2.A.front().outer_product(fe.N,fe.N,true,fe.detJxW); for (size_t j = 0; j < solPt.size(); j++) gl2.b[j].add(fe.N,solPt[j]*fe.detJxW); @@ -418,6 +427,7 @@ bool GlbL2::solve (Matrix& sField) #if SP_DEBUG > 1 std::cout <<"\nSolution:"<< sField; #endif + return true; } @@ -473,11 +483,11 @@ bool ASMbase::L2projection (Matrix& sField, const TimeDomain& time) { PROFILE2("ASMbase::L2projection"); - - GlbL2 gl2(integrand,this->getNoNodes(1)); - L2GlobalInt dummy(*gl2.pA,*gl2.pB); - - gl2.preAssemble(MNPC,this->getNoElms(true)); + GlbL2 gl2(integrand,this->getNoNodes(1),glbL2_A); + L2GlobalInt dummy(this->getNoNodes(1),glbL2_A ? nullptr : gl2.pA,*gl2.pB); + if (!glbL2_A) + gl2.preAssemble(MNPC,this->getNoElms(true)); + glbL2_A = gl2.pA; return this->integrate(gl2,dummy,time) && gl2.solve(sField); } @@ -485,12 +495,12 @@ bool ASMbase::L2projection (Matrix& sField, bool ASMbase::L2projection (Matrix& sField, FunctionBase* function, double t) { PROFILE2("ASMbase::L2projection"); - - GlbL2 gl2(function,this->getNoNodes(1)); - L2GlobalInt dummy(*gl2.pA,*gl2.pB); + GlbL2 gl2(function,this->getNoNodes(1),glbL2_A); + L2GlobalInt dummy(this->getNoNodes(1),glbL2_A ? nullptr : gl2.pA,*gl2.pB); + if (!glbL2_A) + gl2.preAssemble(MNPC,this->getNoElms(true)); + glbL2_A = gl2.pA; TimeDomain time; time.t = t; - - gl2.preAssemble(MNPC,this->getNoElms(true)); return this->integrate(gl2,dummy,time) && gl2.solve(sField); } @@ -499,12 +509,12 @@ bool ASMbase::L2projection (const std::vector& sField, const FunctionVec& function, double t) { PROFILE2("ASMbase::L2projection"); - - GlbL2 gl2(function,this->getNoNodes(1)); - L2GlobalInt dummy(*gl2.pA,*gl2.pB); + GlbL2 gl2(function,this->getNoNodes(1),glbL2_A); + L2GlobalInt dummy(this->getNoNodes(1),glbL2_A ? nullptr : gl2.pA,*gl2.pB); + if (!glbL2_A) + gl2.preAssemble(MNPC,this->getNoElms(true)); + glbL2_A = gl2.pA; TimeDomain time; time.t = t; - - gl2.preAssemble(MNPC,this->getNoElms(true)); return this->integrate(gl2,dummy,time) && gl2.solve(sField); } @@ -520,44 +530,45 @@ bool ASMbase::globalL2projection (Matrix& sField, // Assemble the projection matrices size_t i, nnod = this->getNoProjectionNodes(); size_t j, ncomp = integrand.dim(); - SparseMatrix* A; StdVector* B; switch (GlbL2::MatrixType) { case LinAlg::UMFPACK: - A = new SparseMatrix(SparseMatrix::UMFPACK); + if (!glbL2_A) + glbL2_A = new SparseMatrix(SparseMatrix::UMFPACK); B = new StdVector(nnod*ncomp); break; #ifdef HAS_PETSC case LinAlg::PETSC: if (GlbL2::SolverParams) { - A = new PETScMatrix(ProcessAdm(), *GlbL2::SolverParams); - B = new PETScVector(ProcessAdm(), nnod*ncomp); + static ProcessAdm adm; + if (!glbL2_A) + glbL2_A = new PETScMatrix(adm, *GlbL2::SolverParams); + B = new PETScVector(adm, nnod*ncomp); break; } #endif default: - A = new SparseMatrix(SparseMatrix::SUPERLU); + glbL2_A = new SparseMatrix(SparseMatrix::SUPERLU); B = new StdVector(nnod*ncomp); } - A->redim(nnod,nnod); + glbL2_A->redim(nnod,nnod); - if (!this->assembleL2matrices(*A,*B,integrand,continuous)) + if (!this->assembleL2matrices(*glbL2_A,*B,integrand,continuous)) { - delete A; delete B; return false; } #if SP_DEBUG > 1 - std::cout <<"---- Matrix A -----\n"<< *A + std::cout <<"---- Matrix A -----\n"<< *glbL2_A <<"-------------------"<< std::endl; std::cout <<"---- Vector B -----"<< *B <<"-------------------"<< std::endl; #endif // Solve the patch-global equation system - if (!A->solve(*B)) return false; + if (!glbL2_A->solve(*B)) return false; // Store the control-point values of the projected field sField.resize(ncomp,nnod); @@ -569,7 +580,6 @@ bool ASMbase::globalL2projection (Matrix& sField, std::cout <<"- Solution Vector -"<< sField <<"-------------------"<< std::endl; #endif - delete A; delete B; return true; } diff --git a/src/ASM/GlbL2projector.h b/src/ASM/GlbL2projector.h index e851fab72..fa4604988 100644 --- a/src/ASM/GlbL2projector.h +++ b/src/ASM/GlbL2projector.h @@ -24,7 +24,6 @@ class FunctionBase; class LinSolParams; class SparseMatrix; class StdVector; -class ProcessAdm; typedef std::vector IntVec; //!< Vector of integers typedef std::vector uIntVec; //!< Vector of unsigned integers @@ -113,15 +112,18 @@ class GlbL2 : public Integrand //! \brief The constructor initializes the projection matrices. //! \param[in] p The main problem integrand //! \param[in] n Dimension of the L2-projection matrices (number of nodes) - GlbL2(IntegrandBase* p, size_t n); + //! \param[in] A Sparse matrix to use + GlbL2(IntegrandBase* p, size_t n, SparseMatrix* A); //! \brief Alternative constructor for projection of an explicit function. //! \param[in] f The function to do L2-projection on //! \param[in] n Dimension of the L2-projection matrices (number of nodes) - GlbL2(FunctionBase* f, size_t n); + //! \param[in] A Sparse matrix to use + GlbL2(FunctionBase* f, size_t n, SparseMatrix* A); //! \brief Alternative constructor for projection of explicit functions. //! \param[in] f The functions to do L2-projection on //! \param[in] n Dimension of the L2-projection matrices (number of nodes) - GlbL2(const FunctionVec& f, size_t n); + //! \param[in] A Sparse matrix to use + GlbL2(const FunctionVec& f, size_t n, SparseMatrix* A); //! \brief The destructor frees the system matrix and system vector. virtual ~GlbL2(); @@ -203,9 +205,6 @@ class GlbL2 : public Integrand IntegrandBase* problem; //!< The main problem integrand FunctionVec functions; //!< Explicit functions to L2-project size_t nrhs; //!< Number of right-hand-size vectors -#ifdef HAS_PETSC - ProcessAdm* adm; //!< Process administrator for PETSc -#endif }; #endif diff --git a/src/LinAlg/PETScMatrix.C b/src/LinAlg/PETScMatrix.C index 81f798215..6a2d5bd8f 100644 --- a/src/LinAlg/PETScMatrix.C +++ b/src/LinAlg/PETScMatrix.C @@ -504,7 +504,7 @@ bool PETScMatrix::solve (SystemVector& B, bool newLHS, Real*) if (!Bptr) return false; - if (A.empty() || !assembled) + if (A.empty() || !assembled || B.dim() != this->rows()) return this->solveDirect(*Bptr); Vec x; @@ -582,27 +582,30 @@ bool PETScMatrix::solveDirect(PETScVector& B) { // the sparsity pattern has been grown in-place, we need to init PETsc state. // this is currently only used for patch-global L2 systems. - if (A.empty() && !this->optimiseSLU()) - return false; - - // Set correct number of rows and columns for matrix. size_t nrow = IA.size()-1; - MatSetSizes(pA, nrow, nrow, PETSC_DECIDE, PETSC_DECIDE); - MatSetFromOptions(pA); - PetscInt max = 0; - for (size_t i = 0; i < nrow; ++i) // symmetric so row/column sizes should be the same - if (IA[i+1]-IA[i] > max) - max = IA[i+1]-IA[i]; - MatSeqAIJSetPreallocation(pA, max, nullptr); - MatSetOption(pA, MAT_NEW_NONZERO_LOCATION_ERR, PETSC_FALSE); - MatSetUp(pA); - - for (size_t j = 0; j < nrow; ++j) - for (int i = IA[j]; i < IA[j+1]; ++i) - MatSetValue(pA, JA[i], j, A[i], INSERT_VALUES); + if (!assembled) { + if (A.empty() && !this->optimiseSLU()) + return false; + // Set correct number of rows and columns for matrix. + nrow = IA.size()-1; + MatSetSizes(pA, nrow, nrow, PETSC_DECIDE, PETSC_DECIDE); + MatSetFromOptions(pA); + PetscInt max = 0; + for (size_t i = 0; i < nrow; ++i) // symmetric so row/column sizes should be the same + if (IA[i+1]-IA[i] > max) + max = IA[i+1]-IA[i]; + MatSeqAIJSetPreallocation(pA, max, nullptr); + MatSetOption(pA, MAT_NEW_NONZERO_LOCATION_ERR, PETSC_FALSE); + MatSetUp(pA); - MatAssemblyBegin(pA,MAT_FINAL_ASSEMBLY); - MatAssemblyEnd(pA,MAT_FINAL_ASSEMBLY); + for (size_t j = 0; j < nrow; ++j) + for (int i = IA[j]; i < IA[j+1]; ++i) + MatSetValue(pA, JA[i], j, A[i], INSERT_VALUES); + + MatAssemblyBegin(pA,MAT_FINAL_ASSEMBLY); + MatAssemblyEnd(pA,MAT_FINAL_ASSEMBLY); + assembled = true; + } Vec B1, x; VecCreate(PETSC_COMM_SELF, &B1); diff --git a/src/LinAlg/PETScMatrix.h b/src/LinAlg/PETScMatrix.h index ca5a1fe57..f48170f59 100644 --- a/src/LinAlg/PETScMatrix.h +++ b/src/LinAlg/PETScMatrix.h @@ -82,6 +82,7 @@ class PETScVector : public StdVector //! \brief Return associated process administrator const ProcessAdm& getAdm() const { return adm; } + protected: Vec x; //!< The actual PETSc vector const ProcessAdm& adm; //!< Process administrator