Code owners
Assign users and groups as approvers for specific file changes. Learn more.
multigridstep.cc 9.45 KiB
// -*- tab-width: 8; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=8 sw=4 sts=4:
#include <typeinfo>
#include <dune/solvers/transferoperators/multigridtransfer.hh>
#include <dune/solvers/solvers/loopsolver.hh>
#include <dune/solvers/solvers/linearsolver.hh>
#include <dune/solvers/common/genericvectortools.hh>
//template<class MatrixType, class VectorType, class BitVectorType>
//inline
//const MatrixType* MultigridStep<MatrixType, VectorType, BitVectorType>::
//getMatrix()
//{
// return this->mat_;
//}
template<class MatrixType, class VectorType, class BitVectorType>
inline
void MultigridStep<MatrixType, VectorType, BitVectorType>::
setMGType(int mu, int nu1, int nu2)
{
mu_ = mu;
nu1_ = nu1;
nu2_ = nu2;
}
/** \brief Sets up the MultigridStep for computation
*
* Resizes all hierarchy containers to appropriate size
* Sets up:
* \li matrix hierarchy
* \li ignore nodes hierarchy
* \li the base solver
*
* \b NOTE: All data, i.e.
* \li transfer operators
* \li fine level matrix, solution vector and right hand side
* \li smoothers
* have to be set prior to the call to preprocess().
* If any is changed subsequently, preprocess() has to be called again before iterate().
* There is one exception. The right hand side may be changed by setRhs() without a subsequent call to preprocess().
* This is in order to be able to compute solutions to multiple rhs w/o having to reassemble the same matrix hierarchy time and again.
*/
template<class MatrixType, class VectorType, class BitVectorType>
void MultigridStep<MatrixType, VectorType, BitVectorType>::preprocess()
{
if(preprocessCalled)
std::cout << "Warning: Preprocess has already been called without calling iterate afterwards!" << std::endl;
preprocessCalled = true;
this->level_ = numLevels()-1;
// //////////////////////////////////////////////////////////
// Check that transfer operators have been supplied
// //////////////////////////////////////////////////////////
for (size_t i=0; i<this->mgTransfer_.size(); i++)
if (this->mgTransfer_[i] == NULL)
DUNE_THROW(SolverError, "You have not supplied a multigrid restriction operator "
"for level " << i);
// //////////////////////////////////////////////////////////
// Resize hierarchy containers to current number of levels
// and set fine level data where needed.
// //////////////////////////////////////////////////////////
for (int i=0; i<int(ignoreNodesHierarchy_.size())-1; i++)
if (ignoreNodesHierarchy_[i])
delete(ignoreNodesHierarchy_[i]);
xHierarchy_.resize(numLevels());
rhsHierarchy_.resize(numLevels());
matrixHierarchy_.resize(numLevels());
ignoreNodesHierarchy_.resize(numLevels());
for (int i=0; i<int(matrixHierarchy_.size())-1; i++)
{
matrixHierarchy_[i].reset();
xHierarchy_[i].reset();
ignoreNodesHierarchy_[i] = NULL;
}
matrixHierarchy_.back() = this->mat_;
xHierarchy_.back() = Dune::stackobject_to_shared_ptr(*(this->x_));
rhsHierarchy_.back() = *(this->rhs_);
presmoother_.resize(numLevels());
postsmoother_.resize(numLevels());
typename SmootherCache::iterator end=levelWiseSmoothers_.end();
for (size_t i=0; i<numLevels(); ++i)
{
typename SmootherCache::iterator levelSmoother = levelWiseSmoothers_.find(i);
if (levelSmoother != end)
presmoother_[i] = postsmoother_[i] = levelSmoother->second;
else
{
presmoother_[i] = presmootherDefault_;
postsmoother_[i] = postsmootherDefault_;
}
}
// /////////////////////////////////////////////////////
// Assemble the complete hierarchy of matrices
// /////////////////////////////////////////////////////
for (int i=this->numLevels()-2; i>=0; i--)
{
this->matrixHierarchy_[i] = std::shared_ptr<MatrixType>(new MatrixType);
this->xHierarchy_[i] = std::shared_ptr<VectorType>(new VectorType);
// Compute which entries are present in the (sparse) coarse grid stiffness
// matrices.
this->mgTransfer_[i]->galerkinRestrictSetOccupation(*this->matrixHierarchy_[i+1], *std::const_pointer_cast<MatrixType>(this->matrixHierarchy_[i]));
// setup matrix
this->mgTransfer_[i]->galerkinRestrict(*this->matrixHierarchy_[i+1], *std::const_pointer_cast<MatrixType>(this->matrixHierarchy_[i]));
// Set solution vector sizes for the lower levels
GenericVector::resize(*(this->xHierarchy_[i]),*this->matrixHierarchy_[i]);
}
// /////////////////////////////////////////////////////
// Setup dirichlet bitfields
// /////////////////////////////////////////////////////
if (this->ignoreNodes_!=0)
{
ignoreNodesHierarchy_[numLevels()-1] = const_cast<BitVectorType*>(this->ignoreNodes_);
for (int i=numLevels()-2; i>=0; --i)
{
ignoreNodesHierarchy_[i] = new BitVectorType();
this->mgTransfer_[i]->restrictToFathers(*ignoreNodesHierarchy_[i+1], *ignoreNodesHierarchy_[i]);
}
} else
DUNE_THROW(SolverError, "We need a set of nodes to ignore");
// /////////////////////////////////////////////
// Set up base solver
// /////////////////////////////////////////////
using SmootherType = LinearIterationStep<MatrixType, VectorType, BitVectorType>;
using LinearSolverType = LinearSolver<MatrixType, VectorType>;
if (basesolver_)
{
// If the base solver can ignore dofs give it the ignoreNodes field
if (dynamic_cast<CanIgnore<BitVectorType>*>(this->basesolver_.get()))
dynamic_cast<CanIgnore<BitVectorType>*>(this->basesolver_.get())->setIgnore(*ignoreNodesHierarchy_[0]);
typedef ::LoopSolver<VectorType> DuneSolversLoopSolver;
if (dynamic_cast<DuneSolversLoopSolver*>(this->basesolver_.get())) {
DuneSolversLoopSolver* loopBaseSolver = dynamic_cast<DuneSolversLoopSolver*> (this->basesolver_.get());
assert(dynamic_cast<SmootherType*>(&loopBaseSolver->getIterationStep()));
dynamic_cast<SmootherType*>(&loopBaseSolver->getIterationStep())->setProblem(*(this->matrixHierarchy_[0]), *this->xHierarchy_[0], this->rhsHierarchy_[0]);
dynamic_cast<SmootherType*>(&loopBaseSolver->getIterationStep())->setIgnore(*ignoreNodesHierarchy_[0]);
}
else if (dynamic_cast<LinearSolverType*>(this->basesolver_.get())) {
LinearSolverType* linearBaseSolver = dynamic_cast<LinearSolverType*> (this->basesolver_.get());
linearBaseSolver->setProblem(*(this->matrixHierarchy_[0]), *this->xHierarchy_[0], this->rhsHierarchy_[0]);
}
else {
DUNE_THROW(SolverError, "You can't use " << typeid(*this->basesolver_).name()
<< " as a base solver in a MultigridStep!");
}
}
}
template<class MatrixType, class VectorType, class BitVectorType>
void MultigridStep<MatrixType, VectorType, BitVectorType>::nestedIteration()
{
for (level_ = 0; level_<(int) numLevels(); level_++)
{
iterate();
mgTransfer_[level_]->prolong(*xHierarchy_[level_], *xHierarchy_[level_+1]);
}
iterate();
}
template<class MatrixType, class VectorType, class BitVectorType>
void MultigridStep<MatrixType, VectorType, BitVectorType>::postprocess()
{
}
template<class MatrixType, class VectorType, class BitVectorType>
void MultigridStep<MatrixType, VectorType, BitVectorType>::iterate()
{
preprocessCalled = false;
int& level = this->level_;
// Define references just for ease of notation
std::vector<std::shared_ptr<const MatrixType> > const &mat = this->matrixHierarchy_;
std::vector<std::shared_ptr<VectorType> >& x = this->xHierarchy_;
std::vector<VectorType>& rhs = this->rhsHierarchy_;
// Solve directly if we're looking at the coarse problem
if ((level == 0) and (basesolver_)) {
basesolver_->solve();
return; }
// Presmoothing
presmoother_[level]->setProblem(*(this->matrixHierarchy_[level]), *x[level], rhs[level]);
presmoother_[level]->ignoreNodes_ = ignoreNodesHierarchy_[level];
for (int i=0; i<this->nu1_; i++)
presmoother_[level]->iterate();
if (level>0)
{
// /////////////////////////
// Restriction
// compute residual
// fineResidual = rhs[level] - mat[level] * x[level];
VectorType fineResidual = rhs[level];
mat[level]->mmv(*x[level], fineResidual);
// restrict residual
this->mgTransfer_[level-1]->restrict(fineResidual, rhs[level-1]);
// Set Dirichlet values.
GenericVector::truncate(rhs[level-1], *ignoreNodesHierarchy_[level-1]);
// Choose all zeros as the initial correction
*x[level-1] = 0;
// ///////////////////////////////////////
// Recursively solve the coarser system
level--;
for (int i=0; i<mu_; i++)
iterate();
level++;
// ////////////////////////////////////////
// Prolong
// add correction to the presmoothed solution
VectorType tmp;
this->mgTransfer_[level-1]->prolong(*x[level-1], tmp);
*x[level] += tmp;
}
// Postsmoothing
postsmoother_[level]->setProblem(*(mat[level]), *x[level], rhs[level]);
postsmoother_[level]->ignoreNodes_ = ignoreNodesHierarchy_[level];
for (int i=0; i<this->nu2_; i++)
postsmoother_[level]->iterate();
// *(this->x_) = xHierarchy_.back();
}