From f628f3dd53191e9dfdc9d1e0bdcb54073cd42844 Mon Sep 17 00:00:00 2001 From: Arne Morten Kvarving Date: Fri, 14 Oct 2016 15:30:37 +0200 Subject: [PATCH] added: support for K-refinement cycles --- Apps/Common/SIMKCyclePC.h | 302 ++++++++++++++++++++++++++++++++++++ Apps/Common/SIMMxV.h | 69 ++++++++ Apps/Common/SIMSolverKRef.h | 164 ++++++++++++++++++++ src/ASM/ASMs2D.h | 12 +- src/ASM/ASMs3D.h | 12 +- src/ASM/SAMpatchPETSc.C | 2 +- src/LinAlg/LinSolParams.C | 2 + src/LinAlg/PETScMatrix.C | 108 ++++++++++--- src/LinAlg/PETScMatrix.h | 25 ++- src/LinAlg/PETScSolParams.C | 24 +++ src/LinAlg/PETScSolParams.h | 35 +++++ src/SIM/SIM2D.C | 1 + src/SIM/SIM3D.C | 1 + src/SIM/SIMinput.h | 5 + 14 files changed, 728 insertions(+), 34 deletions(-) create mode 100644 Apps/Common/SIMKCyclePC.h create mode 100644 Apps/Common/SIMMxV.h create mode 100644 Apps/Common/SIMSolverKRef.h diff --git a/Apps/Common/SIMKCyclePC.h b/Apps/Common/SIMKCyclePC.h new file mode 100644 index 000000000..cdc2c1a8e --- /dev/null +++ b/Apps/Common/SIMKCyclePC.h @@ -0,0 +1,302 @@ +//============================================================================== +//! +//! \file SIMKCyclePC.h +//! +//! \date Nov 28 2012 +//! +//! \author Arne Morten Kvarving / SINTEF +//! +//! \brief Use a simulator as a preconditioner for K-refinement cycles. +//! +//============================================================================== +#ifndef SIM_K_CYCLE_PC_H_ +#define SIM_K_CYCLE_PC_H_ + +#ifdef HAS_PETSC + +#include "AlgEqSystem.h" +#include "ASMs2D.h" +#include "ASMs3D.h" +#include "SAM.h" +#include +#include + + +/*! + \brief Template for K-cycle preconditioning. + \details Currently only single patch models. +*/ + +template +class SIMKCyclePC : public PETScPC +{ +public: + //! \brief Default constructor. + //! \param fromSim The model for the solution + //! \param pcSim The model for the preconditioner + SIMKCyclePC(Sim& fromSim, Sim& pcSim) : + S1(pcSim), fSim(fromSim), + B(SparseMatrix::SUPERLU), BT(SparseMatrix::SUPERLU) + { + setup(); + } + + //! \brief Setup the restriction/prolongation operators. + //! \details Operators are constructed through Greville point projection. + void setup() + { + if (fSim.getPatch(1)->getNoSpaceDim() == 2) { + const ASMs2D* fpch = dynamic_cast(fSim.getPatch(1)); + const ASMs2D* pch = dynamic_cast(S1.getPatch(1)); + std::array gPrm; + fpch->getGrevilleParameters(gPrm[0], 0); + fpch->getGrevilleParameters(gPrm[1], 1); + N.resize(gPrm[0].size()*gPrm[1].size(), S1.getNoDOFs()); + NT.resize(S1.getNoDOFs(), gPrm[0].size()*gPrm[1].size()); + std::vector spline; + pch->getBasis(1)->computeBasisGrid(gPrm[0], gPrm[1], spline); + int p1 = pch->getBasis(1)->order_u(); + int p2 = pch->getBasis(1)->order_v(); + int n1 = pch->getBasis(1)->numCoefs_u(); + int n2 = pch->getBasis(1)->numCoefs_v(); + for (size_t ip = 0; ip < gPrm[0].size()*gPrm[1].size(); ++ip) { + IntVec idx; + ASMs2D::scatterInd(n1,n2,p1,p2,spline[ip].left_idx,idx); + for (size_t k = 0; k < spline[ip].basisValues.size(); ++k) { + N(ip+1, idx[k]+1) += spline[ip].basisValues[k]; + NT(idx[k]+1, ip+1) += spline[ip].basisValues[k]; + } + } + fpch->getBasis(1)->computeBasisGrid(gPrm[0], gPrm[1], spline); + p1 = fpch->getBasis(1)->order_u(); + p2 = fpch->getBasis(1)->order_v(); + n1 = fpch->getBasis(1)->numCoefs_u(); + n2 = fpch->getBasis(1)->numCoefs_v(); + B.resize(fSim.getNoDOFs(), fSim.getNoDOFs()); + BT.resize(fSim.getNoDOFs(), fSim.getNoDOFs()); + + for (size_t ip = 0; ip < gPrm[0].size()*gPrm[1].size(); ++ip) { + IntVec idx; + ASMs2D::scatterInd(n1,n2,p1,p2,spline[ip].left_idx,idx); + for (size_t k = 0; k < spline[ip].basisValues.size(); ++k) { + B(ip+1, idx[k]+1) += spline[ip].basisValues[k]; + BT(idx[k]+1, ip+1) += spline[ip].basisValues[k]; + } + } + } else { + const ASMs3D* fpch = dynamic_cast(fSim.getPatch(1)); + const ASMs3D* pch = dynamic_cast(S1.getPatch(1)); + std::array gPrm; + fpch->getGrevilleParameters(gPrm[0], 0); + fpch->getGrevilleParameters(gPrm[1], 1); + fpch->getGrevilleParameters(gPrm[2], 2); + N.resize(gPrm[0].size()*gPrm[1].size()*gPrm[2].size(), S1.getNoDOFs()); + NT.resize(S1.getNoDOFs(), gPrm[0].size()*gPrm[1].size()*gPrm[2].size()); + std::vector spline; + pch->getBasis(1)->computeBasisGrid(gPrm[0], gPrm[1], gPrm[2], spline); + int p1 = pch->getBasis(1)->order(0); + int p2 = pch->getBasis(1)->order(1); + int p3 = pch->getBasis(1)->order(2); + int n1 = pch->getBasis(1)->numCoefs(0); + int n2 = pch->getBasis(1)->numCoefs(1); + int n3 = pch->getBasis(1)->numCoefs(2); + for (size_t ip = 0; ip < gPrm[0].size()*gPrm[1].size()*gPrm[2].size(); ++ip) { + IntVec idx; + ASMs3D::scatterInd(n1,n2,n3,p1,p2,p3,spline[ip].left_idx,idx); + for (size_t k = 0; k < spline[ip].basisValues.size(); ++k) { + N(ip+1, idx[k]+1) += spline[ip].basisValues[k]; + NT(idx[k]+1, ip+1) += spline[ip].basisValues[k]; + } + } + fpch->getBasis(1)->computeBasisGrid(gPrm[0], gPrm[1], gPrm[2], spline); + p1 = fpch->getBasis(1)->order(0); + p2 = fpch->getBasis(1)->order(1); + p3 = fpch->getBasis(1)->order(2); + n1 = fpch->getBasis(1)->numCoefs(0); + n2 = fpch->getBasis(1)->numCoefs(1); + n3 = fpch->getBasis(1)->numCoefs(2); + B.resize(fSim.getNoDOFs(), fSim.getNoDOFs()); + BT.resize(fSim.getNoDOFs(), fSim.getNoDOFs()); + + for (size_t ip = 0; ip < gPrm[0].size()*gPrm[1].size()*gPrm[2].size(); ++ip) { + IntVec idx; + ASMs3D::scatterInd(n1,n2,n3,p1,p2,p3,spline[ip].left_idx,idx); + for (size_t k = 0; k < spline[ip].basisValues.size(); ++k) { + B(ip+1, idx[k]+1) += spline[ip].basisValues[k]; + BT(idx[k]+1, ip+1) += spline[ip].basisValues[k]; + } + } + } + } + + //! \brief Evaluate preconditioner. + bool eval(Vec& x, Vec& y) override + { + // copy vector into a standard vector for expansion + int len; + VecGetSize(x, &len); + StdVector X(len); + std::vector idx(len); + std::iota(idx.begin(), idx.end(), 0); + VecGetValues(x, len, idx.data(), X.getPtr()); + + // Expand solution vector x to DOF vector + StdVector dofSol(fSim.getNoDOFs()); + fSim.getSAM()->expandSolution(X, dofSol); + + // Evaluate restriction - N^T * (B^T)^-1 * dofSol + BT.solve(dofSol); + StdVector restrict(NT.rows()); + NT.multiply(dofSol, restrict); + + // copy to solution vector + StdVector* eqsVec = S1.getVector(); + eqsVec->init(); + vectorToSolVec(S1, *eqsVec, restrict); + eqsVec->beginAssembly(); + eqsVec->endAssembly(); + + // Perform inner solve, store in restrict + double cond; + int oldLev = this->S1.msgLevel; + this->S1.msgLevel = 0; + if (!this->S1.solveSystem(restrict, 0, &cond, "displacement", false)) + return false; + this->S1.msgLevel = oldLev; + + // Evaluate prolonged solution - dofSol = B^-1 * N * restrict + N.multiply(restrict, dofSol); + B.solve(dofSol); + + // Inject into solution vector y + const int* meqn = fSim.getSAM()->getMEQN(); + for (size_t i = 1; i <= fSim.getPatch(1)->getNoNodes(); ++i) { + int eq = meqn[i-1]; + if (eq > 0) + VecSetValue(y, eq-1, dofSol(i), INSERT_VALUES); + } + + return true; + } + + //! \brief Evaluate the right-hand-side system - the prolongiation of the coarse system solution. + //! \param rhs The resulting vector + void evalRHS(PETScVector& rhs) + { + // Solve coarse system, store in solCoarse + StdVector solCoarse(S1.getNoDOFs()); + S1.solveSystem(solCoarse); + + // Evaluate prolongiation - rhs = B^-1 * N * solCoarse + StdVector prolong(N.rows()); + N.multiply(solCoarse, prolong); + B.solve(prolong); + vectorToSolVec(fSim, rhs, prolong); + + rhs.beginAssembly(); + rhs.endAssembly(); + } + + //! \brief Returns name of this preconditioner. + const char* getName() const override { return "K-cycle preconditioner"; } + +protected: + //! \brief Helper template for copying a equation vector to a DOF vector. + //! \param sim The simulator to use + //! \param V1 The resulting DOF vector + //! \param V2 The initial equation vecto + template + void vectorToSolVec(Sim& sim, V1& dof, const V2& eqn) + { + const int* meqn = sim.getSAM()->getMEQN(); + for (size_t i = 1; i <= sim.getPatch(1)->getNoNodes(); ++i) { + int eq = meqn[i-1]; + if (eq > 0) + dof(eq) = eqn(i); + } + } + + Sim& S1; //!< Preconditioner simulator + Sim& fSim; //!< Solution model + SparseMatrix B, BT; //!< Greville mass matrix for fSim + SparseMatrix N; //!< Basis functions of S1 in greville of fSim + SparseMatrix NT; //!< Stupid no-transpose multiply +}; + + +/*! + \brief Driver class for K-refined NURBS-based FEM analysis of linear problems. +*/ + +template +class SIMKRefWrap : public T1 +{ +public: + //! \brief Default constructor. + explicit SIMKRefWrap(const typename T1::SetupProps& props) : T1(props) {} + + //! \brief Clears out patches, FE model and linear system. + void clearModel() + { + this->clearProperties(); + delete this->mySam; + this->mySam=nullptr; + for (ASMbase* patch : this->myModel) + delete patch; + this->myModel.clear(); + delete this->myEqSys; + this->myEqSys = nullptr; + } + + //! \brief Set a linear solver parameter. + //! \param key Parameter to set + //! \param value Value for parameter + //! \param old If non-null, previous value is stored here + void setLinSolParam(const std::string& key, const std::string& value, std::string* old = nullptr) + { + if (!this->mySolParams) + this->mySolParams = new LinSolParams; + + if (old) + *old = this->mySolParams->getStringValue(key); + + this->mySolParams->addValue(key, value); + } + + //! \brief Remove shell matrix-vector product. + //! \param oldSolver Solver used before it was overridden with PETSc for shell solver + //! \param oldRtol Old relative tolerance before it was overridden with shell solver tolerance + void clearMxV(int oldSolver, const std::string& oldRtol) + { + this->mySolParams->addValue("matrixfree", "0"); + this->mySolParams->addValue("rtol", oldRtol); + this->getMatrix()->clearMxV(); + + // We have to flag that we want to solve using sparse direct solver for preconditioner. + // We have assembled into a PETSc matrix previously. + if (oldSolver != SystemMatrix::PETSC) + static_cast(this->myEqSys->getMatrix(0))->setSolveSparse(true); + } + + //! \brief Obtain a pointer to current PETSc system matrix. + PETScMatrix* getMatrix() + { + return dynamic_cast(this->myEqSys->getMatrix(0)); + } + + //! \brief Obtain a pointer to current system vector. + StdVector* getVector() + { + return dynamic_cast(this->myEqSys->getVector(0)); + } + + //! \brief Assemble system. + bool assembleStep() + { + return this->assembleSystem(); + } +}; + +#endif + +#endif diff --git a/Apps/Common/SIMMxV.h b/Apps/Common/SIMMxV.h new file mode 100644 index 000000000..d61fee55f --- /dev/null +++ b/Apps/Common/SIMMxV.h @@ -0,0 +1,69 @@ +//============================================================================== +//! +//! \file SIMMxV.h +//! +//! \date Nov 28 2012 +//! +//! \author Arne Morten Kvarving / SINTEF +//! +//! \brief Use a simulator as a matrix-vector product for iterative solvers. +//! +//============================================================================== +#ifndef SIM_MXV_H_ +#define SIM_MXV_H_ + +#ifdef HAS_PETSC + +#include "PETScMatrix.h" + + +/*! + \brief Template for applying the matrix assembled in a simulator to a vector + for use in iterative solvers. +*/ + +template +class SIMMxV : public PETScMxV { +public: + //! \brief Default constructor. + //! \param sim The simulator to wrap + SIMMxV(Sim& sim) : S1(sim) + { + S1.getMatrix()->setMxV(this,true); + } + + //! \brief Evaluate the matrix-vector product y = A*y + bool evalMxV(Vec& x, Vec& y) override + { + MatMult(S1.getMatrix()->getBlockMatrices()[0], x, y); + return true; + } + + //! \brief Set a matrix-free preconditioner. + void setPC(PETScPC* pc) + { + S1.getMatrix()->setPC(pc); + } + + //! \brief Solve at time level. + bool solveStep() + { + TimeStep tp; + return S1.solveStep(tp); + } + + //! \brief The matrix used for building the preconditioner. + //! \details Defaults to the system matrix + bool evalPC(Mat& P) override + { + MatDuplicate(S1.getMatrix()->getMatrix(), MAT_COPY_VALUES, &P); + return true; + } + +protected: + Sim& S1; //!< Reference to wrapped simulator +}; + +#endif + +#endif diff --git a/Apps/Common/SIMSolverKRef.h b/Apps/Common/SIMSolverKRef.h new file mode 100644 index 000000000..a1c37c901 --- /dev/null +++ b/Apps/Common/SIMSolverKRef.h @@ -0,0 +1,164 @@ +// $Id$ +//============================================================================== +//! +//! \file SIMSolverKRef.h +//! +//! \date Feb 2 2016 +//! +//! \author Arne Morten Kvarving / SINTEF +//! +//! \brief Stationary solver class template for K-refinement. +//! +//============================================================================== + +#ifndef _SIM_SOLVER_KREF_H_ +#define _SIM_SOLVER_KREF_H_ + +#include "SIMSolver.h" +#include "SIMenums.h" +#include "SIMMxV.h" +#include "SIMKCyclePC.h" +#include "Utilities.h" + + +/*! + \brief Template class for stationary K-refinement simulator drivers. + \details This template can be instanciated over any type implementing the + ISolver interface. It provides a solver loop with data output. +*/ + +template class SIMSolverKRef : public SIMSolverStat> +{ +public: + //! \brief The constructor forwards to the parent class constructor. + SIMSolverKRef(SIMKRefWrap& s1, SIMKRefWrap& s2) : + SIMSolverStat>(s1), S2(s2) + { + } + + //! \brief Empty destructor. + virtual ~SIMSolverKRef() = default; + + //! \brief Read k-refinement settings. + bool read(const char* infile) override { return this->SIMadmin::read(infile); } + + //! \brief Solves the problem (all k-refinement cycles). + int solveProblem(char* infile, + const char* heading = nullptr, bool = false) override + { + int geoBlk=0, nBlock=0; + + this->printHeading(heading); + + this->S1.getProcessAdm().cout <<"\nK-refinement cycle 1" << std::endl; + + TimeStep tp; + if (!this->initSystem(infile, this->S2, false)) + return 1; + + SIMKRefWrap* s1 = &this->S1; + SIMKRefWrap* s2 = &this->S2; + + // First assemble first preconditioner system. + s2->initDirichlet(); + s2->assembleStep(); + + for (tp.iter = 1; tp.iter <= cycles; ++tp.iter) { + if (tp.iter != 1) + this->S1.getProcessAdm().cout <<"\nK-refinement cycle "<< tp.iter << std::endl; + + // Reload main simulator + s1->increaseAdditionalRefinement(); + if (!this->initSystem(infile,*s1)) + return 3; + + SIMMxV> mxv(*s1); + SIMKCyclePC> pc(*s1,*s2); + mxv.setPC(&pc); + + // First RHS is solution of coarse system prolonged onto fine mesh. + PETScVector rhs(s1->getProcessAdm(), + s1->getSAM()->getNoEquations()); + s2->initDirichlet(); + pc.evalRHS(rhs); + s1->getMatrix()->setInitialGuess(&rhs); + s2->updateDirichlet(0.0); + + s1->initDirichlet(); + if (!mxv.solveStep()) + return 1; + + if (!s1->saveModel(tp.iter == 1 ? infile : nullptr, geoBlk, nBlock)) + return 2; + else if (!s1->saveStep(tp, nBlock)) + return 2; + + // Reset linear solver parameters + s1->clearMxV(oldSolver, oldRtol); + + // Increase additional refinement for s2 (s1 in next step) and clear out model. + s2->increaseAdditionalRefinement(); + s2->clearModel(); + s2->setVTF(s1->getVTF()); + + // Main solver will be preconditioner in next cycle + std::swap(s1, s2); + } + + this->S2.setVTF(nullptr); + + return 0; + } + + //! \brief Load and initialize a simulator. + //! \param infile Input file to use + //! \param sim The simulator to initialize + //! \param setPetsc \e true to force a petsc solver + int initSystem(const char* infile, SIMKRefWrap& sim, bool setPetsc=true) + { + if (!sim.read(infile)) + return 1; + + if (setPetsc) + sim.opt.solver = SystemMatrix::PETSC; + else + oldSolver = sim.opt.solver; + + if (!sim.preprocess()) + return 2; + + sim.setQuadratureRule(sim.opt.nGauss[0],true); + if (setPetsc) { + sim.setLinSolParam("matrixfree", "1"); + std::stringstream str; + str << rtol; + sim.setLinSolParam("rtol", str.str(), &oldRtol); + } + + sim.setMode(SIM::STATIC); + return sim.initSystem(sim.opt.solver,1,1,0,true); + } + +protected: + //! \brief Parses a data section from an XML element. + bool parse(const TiXmlElement* elem) override + { + if (strcasecmp(elem->Value(),"krefinement")) + return this->SIMSolverStat>::parse(elem); + + if (utl::getAttribute(elem,"cycles",cycles)) + IFEM::cout << "\tDoing " << cycles << " k-refinement cycles." << std::endl; + if (utl::getAttribute(elem,"rtol",rtol)) + IFEM::cout << "\tK-solver tolerance " << rtol << std::endl; + + return true; + } + + int cycles = 1; //!< Number of k-refinement cycles + double rtol = 1e-6; //!< Relative tolerance in outer petsc solver + int oldSolver = -1; //!< Old matrix type before we overwrote with PETSc for k-refinement + std::string oldRtol; //!< Old relative tolerance before we overwrite with k-refinement tolerance + SIMKRefWrap& S2; //!< Reference to second simulator +}; + +#endif diff --git a/src/ASM/ASMs2D.h b/src/ASM/ASMs2D.h index 6781f9ab2..884e88da5 100644 --- a/src/ASM/ASMs2D.h +++ b/src/ASM/ASMs2D.h @@ -482,6 +482,12 @@ class ASMs2D : public ASMstruct, public ASM2D virtual bool evalSolution(Matrix& sField, const IntegrandBase& integrand, const RealArray* gpar, bool regular = true) const; + //! \brief Calculates parameter values for the Greville points. + //! \param[out] prm Parameter values in given direction for all points + //! \param[in] dir Parameter direction (0,1) + //! \param[in] basis Basis number (mixed) + bool getGrevilleParameters(RealArray& prm, int dir, int basis=1) const; + protected: // Internal utility methods @@ -519,12 +525,6 @@ class ASMs2D : public ASMstruct, public ASM2D const Vector& getGaussPointParameters(Matrix& uGP, int dir, int nGauss, const double* xi) const; - //! \brief Calculates parameter values for the Greville points. - //! \param[out] prm Parameter values in given direction for all points - //! \param[in] dir Parameter direction (0,1) - //! \param[in] basis Basis number (mixed) - bool getGrevilleParameters(RealArray& prm, int dir, int basis=1) const; - //! \brief Calculates parameter values for the Quasi-Interpolation points. //! \param[out] prm Parameter values in given direction for all points //! \param[in] dir Parameter direction (0,1) diff --git a/src/ASM/ASMs3D.h b/src/ASM/ASMs3D.h index 2d31d5eb6..4c736d8c1 100644 --- a/src/ASM/ASMs3D.h +++ b/src/ASM/ASMs3D.h @@ -544,6 +544,12 @@ class ASMs3D : public ASMstruct, public ASM3D virtual bool evalSolution(Matrix& sField, const IntegrandBase& integrand, const RealArray* gpar, bool regular = true) const; + //! \brief Calculates parameter values for the Greville points. + //! \param[out] prm Parameter values in given direction for all points + //! \param[in] dir Parameter direction (0,1,2) + //! \param[in] basis Basis number (mixed) + bool getGrevilleParameters(RealArray& prm, int dir, int basis=1) const; + protected: // Internal utility methods @@ -581,12 +587,6 @@ class ASMs3D : public ASMstruct, public ASM3D const Vector& getGaussPointParameters(Matrix& uGP, int dir, int nGauss, const double* xi) const; - //! \brief Calculates parameter values for the Greville points. - //! \param[out] prm Parameter values in given direction for all points - //! \param[in] dir Parameter direction (0,1,2) - //! \param[in] basis Basis number (mixed) - bool getGrevilleParameters(RealArray& prm, int dir, int basis=1) const; - //! \brief Calculates parameter values for the Quasi-Interpolation points. //! \param[out] prm Parameter values in given direction for all points //! \param[in] dir Parameter direction (0,1,2) diff --git a/src/ASM/SAMpatchPETSc.C b/src/ASM/SAMpatchPETSc.C index 7d9ad69fa..3da1dd13f 100644 --- a/src/ASM/SAMpatchPETSc.C +++ b/src/ASM/SAMpatchPETSc.C @@ -178,7 +178,7 @@ bool SAMpatchPETSc::expandSolution(const SystemVector& solVec, { PETScVector* Bptr = const_cast(dynamic_cast(&solVec)); if (!Bptr) - return false; + return this->SAMpatch::expandSolution(solVec, dofVec, scaleSD); #ifdef HAVE_MPI if (adm.isParallel()) { diff --git a/src/LinAlg/LinSolParams.C b/src/LinAlg/LinSolParams.C index 4c32f83e0..6ed1f0153 100644 --- a/src/LinAlg/LinSolParams.C +++ b/src/LinAlg/LinSolParams.C @@ -182,6 +182,8 @@ bool LinSolParams::read (const TiXmlElement* elem) addValue("dtol", value); else if ((value = utl::getValue(child,"maxits"))) addValue("maxits", value); + else if (!strcasecmp(child->Value(),"matrixfree")) + addValue("matrixfree", "1"); } if (parseblock == 0) diff --git a/src/LinAlg/PETScMatrix.C b/src/LinAlg/PETScMatrix.C index e9c421914..3ac0c9b5b 100644 --- a/src/LinAlg/PETScMatrix.C +++ b/src/LinAlg/PETScMatrix.C @@ -135,8 +135,8 @@ Real PETScVector::Linfnorm() const PETScMatrix::PETScMatrix (const ProcessAdm& padm, const LinSolParams& spar, LinAlg::LinearSystemType ltype) : - SparseMatrix(SUPERLU, 1), nsp(nullptr), adm(padm), solParams(spar, adm), - linsysType(ltype) + SparseMatrix(SUPERLU, 1), mxv(nullptr), mxvOwnMatrix(true), mfpc(nullptr), + nsp(nullptr), adm(padm), solParams(spar, adm), linsysType(ltype) { // Create matrix object, by default the matrix type is AIJ MatCreate(*adm.getCommunicator(),&pA); @@ -188,9 +188,12 @@ void PETScMatrix::initAssembly (const SAM& sam, bool delayLocking) return; // Get number of local equations in linear system - const PetscInt neq = adm.dd.getMaxEq()- adm.dd.getMinEq() + 1; + bool mxv = solParams.getIntValue("matrixfree") ? true : false; + if (mxv && !mxvOwnMatrix) + return; // Set correct number of rows and columns for matrix. + const PetscInt neq = adm.dd.getMaxEq() - adm.dd.getMinEq() + 1; MatSetSizes(pA,neq,neq,PETSC_DECIDE,PETSC_DECIDE); // Allocate sparsity pattern @@ -407,6 +410,9 @@ void PETScMatrix::initAssembly (const SAM& sam, bool delayLocking) bool PETScMatrix::beginAssembly() { + if (mxv && !mxvOwnMatrix) + return true; + if (matvec.empty()) { for (size_t j = 0; j < cols(); ++j) for (int i = IA[j]; i < IA[j+1]; ++i) @@ -432,6 +438,9 @@ bool PETScMatrix::beginAssembly() bool PETScMatrix::endAssembly() { + if (mxv && !mxvOwnMatrix) + return true; + // Finalizes parallel assembly process MatAssemblyEnd(pA,MAT_FINAL_ASSEMBLY); return true; @@ -441,6 +450,8 @@ bool PETScMatrix::endAssembly() void PETScMatrix::init () { SparseMatrix::init(); + if (mxv && !mxvOwnMatrix) + return; // Set all matrix elements to zero if (matvec.empty()) @@ -467,6 +478,9 @@ bool PETScMatrix::multiply (const SystemVector& B, SystemVector& C) const bool PETScMatrix::solve (SystemVector& B, bool newLHS, Real*) { + if (this->solveSparse) + return this->SparseMatrix::solve(B, newLHS, nullptr); + PETScVector* Bptr = dynamic_cast(&B); if (!Bptr) return false; @@ -474,13 +488,24 @@ bool PETScMatrix::solve (SystemVector& B, bool newLHS, Real*) if (A.empty()) return this->solveDirect(*Bptr); - Vec x; - VecDuplicate(Bptr->getVector(),&x); - VecCopy(Bptr->getVector(),x); + bool result; + if (extInitGuess) { + Vec x; + VecDuplicate(Bptr->getVector(),&x); + VecCopy(Bptr->getVector(), x); + VecCopy(extInitGuess->getVector(), Bptr->getVector()); + VecCopy(x, extInitGuess->getVector()); + VecDestroy(&x); + result = this->solve(extInitGuess->getVector(),Bptr->getVector(),newLHS,false); + } else { + Vec x; + VecDuplicate(Bptr->getVector(),&x); + VecCopy(Bptr->getVector(),x); - bool result = this->solve(x,Bptr->getVector(),newLHS, - solParams.getStringValue("type") != "preonly"); - VecDestroy(&x); + result = this->solve(x,Bptr->getVector(),newLHS, + solParams.getStringValue("type") != "preonly"); + VecDestroy(&x); + } return result; } @@ -488,6 +513,8 @@ bool PETScMatrix::solve (SystemVector& B, bool newLHS, Real*) bool PETScMatrix::solve (const SystemVector& b, SystemVector& x, bool newLHS) { + if (this->solveSparse) + return this->SparseMatrix::solve(b, x, newLHS); const PETScVector* Bptr = dynamic_cast(&b); if (!Bptr) return false; @@ -500,27 +527,63 @@ bool PETScMatrix::solve (const SystemVector& b, SystemVector& x, bool newLHS) } -bool PETScMatrix::solve (const Vec& b, Vec& x, bool newLHS, bool knoll) +void PETScMatrix::setMxV(PETScMxV* MxV, bool ownMatrix) { - // Reset linear solver - if (nLinSolves && solParams.getIntValue("reset_solves")) - if (nLinSolves%solParams.getIntValue("reset_solves") == 0) { - KSPDestroy(&ksp); - KSPCreate(*adm.getCommunicator(),&ksp); - setParams = true; - } + mxv = MxV; + mxvOwnMatrix = ownMatrix; + solveSparse = false; + forcedKSPType = "gmres"; +} + +void PETScMatrix::clearMxV() +{ + if (!mxv) + return; + + if (mxvOwnMatrix) { + MatDestroy(&pA); + MatDuplicate(matvec.front(), MAT_COPY_VALUES, &pA); + MatDestroy(&matvec.front()); + matvec.clear(); + } + + this->setMxV(nullptr); + this->setPC(nullptr); + setParams = true; + extInitGuess = nullptr; + forcedKSPType.clear(); +} + + +bool PETScMatrix::solve (const Vec& b, Vec& x, bool newLHS, bool knoll) +{ if (setParams) { + Mat P; + if (mxv) { + if (mxvOwnMatrix) { + matvec.resize(1); + MatDuplicate(pA, MAT_COPY_VALUES, &matvec.front()); + } + const PetscInt neq = adm.dd.getMaxEq()- adm.dd.getMinEq() + 1; + MatCreateShell(*adm.getCommunicator(), neq, neq, + PETSC_DETERMINE, PETSC_DETERMINE, mxv, &pA); + MatShellSetOperation(pA, MATOP_MULT, (void(*)(void))&PETScSIMMxV); + MatSetUp(pA); + if (!mfpc) + mxv->evalPC(P); + } #if PETSC_VERSION_MINOR < 5 - KSPSetOperators(ksp,pA,pA, newLHS ? SAME_NONZERO_PATTERN : SAME_PRECONDITIONER); + KSPSetOperators(ksp,pA,mxv ? P : pA, newLHS ? SAME_NONZERO_PATTERN : SAME_PRECONDITIONER); #else - KSPSetOperators(ksp,pA,pA); + KSPSetOperators(ksp,pA,mxv && !mfpc ? P : pA); KSPSetReusePreconditioner(ksp, newLHS ? PETSC_FALSE : PETSC_TRUE); #endif if (!setParameters()) return false; setParams = false; } + if (knoll) KSPSetInitialGuessKnoll(ksp,PETSC_TRUE); else @@ -682,7 +745,12 @@ bool PETScMatrix::setParameters(PETScMatrix* P, PETScVector* Pb) PC pc; KSPGetPC(ksp,&pc); - if (matvec.empty()) { + if (mfpc) { + PCSetType(pc,PCSHELL); + PCShellSetContext(pc, mfpc); + PCShellSetName(pc, mfpc->getName()); + PCShellSetApply(pc, PETScSIMPC); + } else if (matvec.empty()) { solParams.setupPC(pc, 0, "", std::set()); } else { if (matvec.size() > 4) { diff --git a/src/LinAlg/PETScMatrix.h b/src/LinAlg/PETScMatrix.h index 6300bfa9e..813c364b1 100644 --- a/src/LinAlg/PETScMatrix.h +++ b/src/LinAlg/PETScMatrix.h @@ -171,13 +171,31 @@ class PETScMatrix : public SparseMatrix //! \brief Get vector of block matrices. Used for tests only. const std::vector& getBlockMatrices() const { return matvec; } + //! \brief Remove pointers related to matrix-free operators. + void clearMxV(); + //! \brief Set the linear solver parameters (solver type, preconditioner, tolerances). //! \param[in] P Preconditioner matrix (ignored here) //! \param[in] Pb Preconditioner vector (ignored here) //! \return True on success virtual bool setParameters(PETScMatrix* P = nullptr, PETScVector* Pb = nullptr); + + //! \brief Set matrix-free operator. + void setMxV(PETScMxV* MxV, bool ownMatrix=false); + + //! \brief Set matrix-free preconditioner. + void setPC(PETScPC* pc) { mfpc = pc; } + + //! \brief Set externally provided initial guess for iterative solver. + void setInitialGuess(PETScVector* init) { extInitGuess = init; } + + //! \brief Enable/disable solution using sparse direct solver. + void setSolveSparse(bool direct) { solveSparse = direct; } + + //! \brief Returns a reference to our linear solver parameters. + PETScSolParams& getSolParams() { return solParams; } protected: - //! \brief Solve a linear system + //! \brief Solve a linear system. bool solve(const Vec& b, Vec& x, bool newLHS, bool knoll); //! \brief Solve system stored in the elem map. @@ -190,6 +208,9 @@ class PETScMatrix : public SparseMatrix PETScMatrix(const PETScMatrix& A) = delete; Mat pA; //!< The actual PETSc matrix + PETScMxV* mxv; //!< Matrix-free operator implementation + bool mxvOwnMatrix=true; //!< The matrix-free operator depends on the assembled matrix + PETScPC* mfpc; //!< Matrix-free preconditioner implementation KSP ksp; //!< Linear equation solver MatNullSpace* nsp; //!< Null-space of linear operator const ProcessAdm& adm; //!< Process administrator @@ -203,6 +224,8 @@ class PETScMatrix : public SparseMatrix LinAlg::LinearSystemType linsysType; //!< Linear system type IS glob2LocEq = nullptr; //!< Index set for global-to-local equations. std::vector matvec; //!< Blocks for block matrices. + PETScVector* extInitGuess = nullptr; //!< Externally provided initial guess + bool solveSparse = false; //!< Solve using sparse direct solver std::vector isvec; //!< Index sets for blocks. std::vector> glb2Blk; //!< Maps matrix entries in CSC order to block matrix entries. diff --git a/src/LinAlg/PETScSolParams.C b/src/LinAlg/PETScSolParams.C index 62dd00dfe..c1b1901b2 100644 --- a/src/LinAlg/PETScSolParams.C +++ b/src/LinAlg/PETScSolParams.C @@ -384,6 +384,7 @@ void PETScSolParams::setupAdditiveSchwarz(PC& pc, size_t block, } } + void PETScSolParams::setupSchurComplement(const std::vector& matvec) { PetscInt m1, n1; @@ -396,6 +397,7 @@ void PETScSolParams::setupSchurComplement(const std::vector& matvec) // TODO: non-SIMPLE schur preconditioners MatGetDiagonal(matvec[0],diagA00); + // TODO: Lumping VecReciprocal(diagA00); SPsetup = true; @@ -410,3 +412,25 @@ void PETScSolParams::setupSchurComplement(const std::vector& matvec) MatDestroy(&tmp); MatDestroy(&tmp2); } + + +extern "C" { + +PetscErrorCode PETScSIMMxV(Mat A, Vec x, Vec y) +{ + void* p; + MatShellGetContext(A, &p); + PETScMxV* sim = static_cast(p); + return sim->evalMxV(x,y) ? 0 : 1; +} + + +PetscErrorCode PETScSIMPC(PC pc, Vec x, Vec y) +{ + void* p; + PCShellGetContext(pc, &p); + PETScPC* sim = static_cast(p); + return sim->eval(x,y) ? 0 : 1; +} + +} diff --git a/src/LinAlg/PETScSolParams.h b/src/LinAlg/PETScSolParams.h index 2523727d7..867850c84 100644 --- a/src/LinAlg/PETScSolParams.h +++ b/src/LinAlg/PETScSolParams.h @@ -40,6 +40,41 @@ typedef std::vector ISMat; //!< Index set matrix enum SchurPrec { SIMPLE, MSIMPLER, PCD }; +/*! \brief Interface to use a simulator as a matrix-vector product. +*/ + +class PETScMxV { +public: + //! \brief Obtain the matrix to use as preconditioner. + virtual bool evalPC(Mat& P) = 0; + + //! \brief Evaluate the matrix-vector product y=A*x. + virtual bool evalMxV(Vec& x, Vec& y) = 0; +}; + + +/*! \brief Interface to use a simulator as a preconditioner. +*/ + +class PETScPC { +public: + //! \brief Evaluate y = P^-1*x + virtual bool eval(Vec& x, Vec& y) = 0; + + //! \brief Get the name of the preconditioner. + virtual const char* getName() const { return "SimPC"; } +}; + + +extern "C" { + //! \brief Evaluate the matrix-vector product y=A*x + PetscErrorCode PETScSIMMxV(Mat A, Vec x, Vec y); + + //! \brief Evaluate the preconditioner y=P^-1*x + PetscErrorCode PETScSIMPC(PC pc, Vec x, Vec y); +} + + /*! \brief Class for PETSc solver parameters. \details It contains information about solver method, preconditioner diff --git a/src/SIM/SIM2D.C b/src/SIM/SIM2D.C index e20d6bf84..09712d6a1 100644 --- a/src/SIM/SIM2D.C +++ b/src/SIM/SIM2D.C @@ -168,6 +168,7 @@ bool SIM2D::parseGeometryTag (const TiXmlElement* elem) int addu = 0, addv = 0; utl::getAttribute(elem,"u",addu); utl::getAttribute(elem,"v",addv); + addu += addRef, addv += addRef; for (int j : patches) { ASM2D* pch; if ((pch = dynamic_cast(this->getPatch(j,true)))) diff --git a/src/SIM/SIM3D.C b/src/SIM/SIM3D.C index ea8ca8072..13fb4df7f 100644 --- a/src/SIM/SIM3D.C +++ b/src/SIM/SIM3D.C @@ -148,6 +148,7 @@ bool SIM3D::parseGeometryTag (const TiXmlElement* elem) utl::getAttribute(elem,"u",addu); utl::getAttribute(elem,"v",addv); utl::getAttribute(elem,"w",addw); + addu += addRef, addv += addRef, addw += addRef; for (int j : patches) { ASM3D* pch; if ((pch = dynamic_cast(this->getPatch(j,true)))) diff --git a/src/SIM/SIMinput.h b/src/SIM/SIMinput.h index 5aea15ea8..9ab20562a 100644 --- a/src/SIM/SIMinput.h +++ b/src/SIM/SIMinput.h @@ -99,6 +99,9 @@ class SIMinput : public SIMbase bool setTracProperty(int code, Property::Type ptype, TractionFunc* field = nullptr); + //! \brief Increase additional refinement. + void increaseAdditionalRefinement() { ++addRef; } + private: //! \brief Parses a subelement of the \a geometry XML-tag. bool parseGeometryTag(const TiXmlElement* elem); @@ -242,6 +245,8 @@ class SIMinput : public SIMbase const InitialCondVec& info); protected: + int addRef = 0; //!< Additional refinement to perform + ModelGenerator* myGen; //!< Model generator TopologySet myEntitys; //!< Set of named topological entities