Skip to content

Commit

Permalink
added: support for K-refinement cycles
Browse files Browse the repository at this point in the history
  • Loading branch information
@akva2
akva2 committed Aug 28, 2018
1 parent f815348 commit 5650bc9
Show file tree
Hide file tree
Showing 14 changed files with 648 additions and 34 deletions.
228 changes: 228 additions & 0 deletions Apps/Common/SIMKCyclePC.h
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,228 @@
//==============================================================================
//!
//! \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 "ASMs2D.h"
#include "ASMs3D.h"
#include "SAM.h"
#include "PETScSolParams.h"
#include <GoTools/geometry/SplineSurface.h>
#include <GoTools/trivariate/SplineVolume.h>


/*!
\brief Template for K-cycle preconditioning.
\details Currently only single patch models.
*/

template<class Sim>
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<const ASMs2D*>(fSim.getPatch(1));
const ASMs2D* pch = dynamic_cast<const ASMs2D*>(S1.getPatch(1));
std::array<RealArray,2> 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<Go::BasisDerivsSf> 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<const ASMs3D*>(fSim.getPatch(1));
const ASMs3D* pch = dynamic_cast<const ASMs3D*>(S1.getPatch(1));
std::array<RealArray,3> 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<Go::BasisDerivs> 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<PetscInt> 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 <class V1, class V2>
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
};

#endif

#endif
69 changes: 69 additions & 0 deletions Apps/Common/SIMMxV.h
View file
Original file line number Diff line number Diff line change
@@ -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 Sim>
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
Loading

0 comments on commit 5650bc9

Please sign in to comment.