Skip to content
Snippets Groups Projects
Commit e18fbe4d authored by Uli Sack's avatar Uli Sack Committed by usack
Browse files

numLevels_ is no longer set explicitely but inferred from the number of... 

numLevels_ is no longer set explicitely but inferred from the number of provided transfer operators.
members mat_, x_, rhs_ have been renamed to ~Hierarchy_ in order to have a clear distinction to the base class members x_,mat_,rhs_
This has been done in accordance with some tasks on the TODO and in cooperation with Carsten.
The changes in a little more detail are:
* constructors with argument numLevels have been deprecated
* resizing of various Hierarchy containers has been moved to preprocess()
* setNumberOfLevels() has been deprecated and does nothing now (s.above)
* setProblem with argument numLevels has been deprecated
* setProblem now only sets the base class members. the fine-level hierarchy entries are set in preprocess
* getSol() and getMatrix() are now implemented and used in/from the base class
* xHierarchy_ is now -same as matrixHierarchy_already was- a std_vector<shared_ptr<...> >
* implementation details of setting the smoothers have changed, but should not affect existing code

If this breaks your code please get back to me.
Uli

[[Imported from SVN: r7374]]
parent 558458da
Branches
Tags
No related merge requests found
Hide whitespace changes
Inline Side-by-side
dune/solvers/iterationsteps/multigridstep.cc
+ 96
− 38
View file @ e18fbe4d
......@@ -9,20 +9,20 @@
#include <dune/solvers/solvers/quadraticipopt.hh>
#endif
template <class MatrixType, class VectorType, class BitVectorType>
VectorType MultigridStep<MatrixType, VectorType, BitVectorType>::
getSol()
{
return this->x_[numLevels_-1];
}
template<class MatrixType, class VectorType, class BitVectorType>
inline
const MatrixType* MultigridStep<MatrixType, VectorType, BitVectorType>::
getMatrix()
{
return this->mat_.back().get();
}
//template <class MatrixType, class VectorType, class BitVectorType>
//VectorType MultigridStep<MatrixType, VectorType, BitVectorType>::
//getSol()
//{
// return *(this->x_);
//}
//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
......@@ -34,6 +34,23 @@ setMGType(int mu, int nu1, int nu2)
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()
{
......@@ -41,36 +58,76 @@ void MultigridStep<MatrixType, VectorType, BitVectorType>::preprocess()
std::cout << "Warning: Preprocess has already been called without calling iterate afterwards!" << std::endl;
preprocessCalled = true;
numLevels_ = mgTransfer_.size() + 1;
this->level_ = this->numLevels_-1;
// //////////////////////////////////////////////////////////
// Check that transfer operators have been supplied
// //////////////////////////////////////////////////////////
if (this->mgTransfer_.size()!=this->numLevels_-1)
DUNE_THROW(SolverError, "You requested " << this->numLevels_ << " levels "
<< "but you supplied " << this->mgTransfer_.size() << " transfer operators!");
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 numLevels_
// 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->mat_[i] = Dune::shared_ptr<MatrixType>(new MatrixType);
this->matrixHierarchy_[i] = Dune::shared_ptr<MatrixType>(new MatrixType);
this->xHierarchy_[i] = Dune::shared_ptr<VectorType>(new VectorType);
// Compute which entries are present in the (sparse) coarse grid stiffness
// matrices.
this->mgTransfer_[i]->galerkinRestrictSetOccupation(*this->mat_[i+1], *std::const_pointer_cast<MatrixType>(this->mat_[i]));
this->mgTransfer_[i]->galerkinRestrictSetOccupation(*this->matrixHierarchy_[i+1], *std::const_pointer_cast<MatrixType>(this->matrixHierarchy_[i]));
// setup matrix
this->mgTransfer_[i]->galerkinRestrict(*this->mat_[i+1], *std::const_pointer_cast<MatrixType>(this->mat_[i]));
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->x_[i],*this->mat_[i]);
GenericVector::resize(*(this->xHierarchy_[i]),*this->matrixHierarchy_[i]);
}
// /////////////////////////////////////////////////////
......@@ -78,7 +135,7 @@ void MultigridStep<MatrixType, VectorType, BitVectorType>::preprocess()
// /////////////////////////////////////////////////////
for (size_t i=0; i<ignoreNodesHierarchy_.size()-1; i++)
delete(ignoreNodesHierarchy_[i]);
if (this->ignoreNodes_!=0)
{
ignoreNodesHierarchy_[this->numLevels_-1] = const_cast<BitVectorType*>(this->ignoreNodes_);
......@@ -87,7 +144,7 @@ void MultigridStep<MatrixType, VectorType, BitVectorType>::preprocess()
ignoreNodesHierarchy_[i] = new BitVectorType();
this->mgTransfer_[i]->restrictToFathers(*ignoreNodesHierarchy_[i+1], *ignoreNodesHierarchy_[i]);
}
} else
} else
DUNE_THROW(SolverError, "We need a set of nodes to ignore");
// /////////////////////////////////////////////
......@@ -110,7 +167,7 @@ void MultigridStep<MatrixType, VectorType, BitVectorType>::preprocess()
typedef LinearIterationStep<MatrixType, VectorType> SmootherType;
assert(dynamic_cast<SmootherType*>(loopBaseSolver->iterationStep_));
dynamic_cast<SmootherType*>(loopBaseSolver->iterationStep_)->setProblem(*(this->mat_[0]), this->x_[0], this->rhs_[0]);
dynamic_cast<SmootherType*>(loopBaseSolver->iterationStep_)->setProblem(*(this->matrixHierarchy_[0]), *this->xHierarchy_[0], this->rhsHierarchy_[0]);
dynamic_cast<SmootherType*>(loopBaseSolver->iterationStep_)->ignoreNodes_ = ignoreNodesHierarchy_[0];
}
......@@ -119,7 +176,7 @@ void MultigridStep<MatrixType, VectorType, BitVectorType>::preprocess()
QuadraticIPOptSolver<MatrixType,VectorType>* ipoptBaseSolver = dynamic_cast<QuadraticIPOptSolver<MatrixType,VectorType>*> (this->basesolver_);
ipoptBaseSolver->setProblem(*(this->mat_[0]), this->x_[0], this->rhs_[0]);
ipoptBaseSolver->setProblem(*(this->matrixHierarchy_[0]), *this->xHierarchy_[0], this->rhsHierarchy_[0]);
}
#endif
else {
......@@ -129,7 +186,6 @@ else {
}
template<class MatrixType, class VectorType, class BitVectorType>
void MultigridStep<MatrixType, VectorType, BitVectorType>::nestedIteration()
{
......@@ -137,7 +193,7 @@ void MultigridStep<MatrixType, VectorType, BitVectorType>::nestedIteration()
{
iterate();
mgTransfer_[level_]->prolong(x_[level_], x_[level_+1]);
mgTransfer_[level_]->prolong(*xHierarchy_[level_], *xHierarchy_[level_+1]);
}
iterate();
......@@ -158,9 +214,9 @@ void MultigridStep<MatrixType, VectorType, BitVectorType>::iterate()
int& level = this->level_;
// Define references just for ease of notation
std::vector<Dune::shared_ptr<const MatrixType> > const &mat = this->mat_;
std::vector<VectorType>& x = this->x_;
std::vector<VectorType>& rhs = this->rhs_;
std::vector<Dune::shared_ptr<const MatrixType> > const &mat = this->matrixHierarchy_;
std::vector<Dune::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) {
......@@ -170,7 +226,7 @@ void MultigridStep<MatrixType, VectorType, BitVectorType>::iterate()
// Presmoothing
presmoother_[level]->setProblem(*(this->mat_[level]), x[level], rhs[level]);
presmoother_[level]->setProblem(*(this->matrixHierarchy_[level]), *x[level], rhs[level]);
presmoother_[level]->ignoreNodes_ = ignoreNodesHierarchy_[level];
for (int i=0; i<this->nu1_; i++)
......@@ -182,7 +238,7 @@ void MultigridStep<MatrixType, VectorType, BitVectorType>::iterate()
// compute residual
// fineResidual = rhs[level] - mat[level] * x[level];
VectorType fineResidual = rhs[level];
mat[level]->mmv(x[level], fineResidual);
mat[level]->mmv(*x[level], fineResidual);
// restrict residual
this->mgTransfer_[level-1]->restrict(fineResidual, rhs[level-1]);
......@@ -192,7 +248,7 @@ void MultigridStep<MatrixType, VectorType, BitVectorType>::iterate()
GenericVector::truncate(rhs[level-1], *ignoreNodesHierarchy_[level-1]);
// Choose all zeros as the initial correction
x[level-1] = 0;
*x[level-1] = 0;
// ///////////////////////////////////////
// Recursively solve the coarser system
......@@ -206,14 +262,16 @@ void MultigridStep<MatrixType, VectorType, BitVectorType>::iterate()
// add correction to the presmoothed solution
VectorType tmp;
this->mgTransfer_[level-1]->prolong(x[level-1], tmp);
x[level] += 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]->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();
}
dune/solvers/iterationsteps/multigridstep.hh
+ 90
− 53
View file @ e18fbe4d
......@@ -2,6 +2,7 @@
#define MULTIGRID_STEP_HH
#include <vector>
#include <map>
#include <dune/common/bitsetvector.hh>
#include <dune/solvers/transferoperators/multigridtransfer.hh>
......@@ -37,22 +38,41 @@
LinearIterationStep<MatrixType, VectorType>* preSmoother,
LinearIterationStep<MatrixType, VectorType>* postSmoother,
Solver* baseSolver,
const BitVectorType* ignoreNodes) :
const BitVectorType* ignoreNodes)
DUNE_DEPRECATED_MSG("The number of levels is no longer set explicitely, but instead inferred from the number of transfer operators.\n Just erase the number of levels from the argument list in your constructor call.") :
LinearIterationStep<MatrixType, VectorType>(mat, x, rhs),
presmoother_(numLevels), postsmoother_(numLevels),
mat_(numLevels),
matrixHierarchy_(numLevels),
ignoreNodesHierarchy_(numLevels),
x_(numLevels),
rhs_(numLevels),
xHierarchy_(numLevels),
rhsHierarchy_(numLevels),
preprocessCalled(false)
{
numLevels_ = numLevels;
level_ = numLevels-1;
mat_[level_] = &mat;
x_[level_] = x;
rhs_[level_] = rhs;
mu_ = mu;
nu1_ = nu1;
nu2_ = nu2;
setSmoother(preSmoother,postSmoother);
basesolver_ = baseSolver;
this->ignoreNodes_ = ignoreNodes;
}
MultigridStep(const MatrixType& mat,
VectorType& x,
VectorType& rhs,
int mu, int nu1, int nu2,
LinearIterationStep<MatrixType, VectorType>* preSmoother,
LinearIterationStep<MatrixType, VectorType>* postSmoother,
Solver* baseSolver,
const BitVectorType* ignoreNodes) :
LinearIterationStep<MatrixType, VectorType>(mat, x, rhs),
preprocessCalled(false)
{
mu_ = mu;
nu1_ = nu1;
nu2_ = nu2;
......@@ -64,29 +84,35 @@
this->ignoreNodes_ = ignoreNodes;
}
MultigridStep(const MatrixType& mat,
VectorType& x,
VectorType& rhs, int numLevels) :
VectorType& rhs, int numLevels)
DUNE_DEPRECATED_MSG("The number of levels is no longer set explicitely, but instead inferred from the number of transfer operators.\n Just erase the number of levels from the argument list in your constructor call.") :
LinearIterationStep<MatrixType, VectorType>(mat, x, rhs),
presmoother_(numLevels),postsmoother_(numLevels),
basesolver_(0),
mat_(numLevels),
matrixHierarchy_(numLevels),
ignoreNodesHierarchy_(numLevels),
x_(numLevels),
rhs_(numLevels),
xHierarchy_(numLevels),
rhsHierarchy_(numLevels),
preprocessCalled(false)
{
numLevels_ = numLevels;
level_ = numLevels-1;
mat_[level_] = Dune::stackobject_to_shared_ptr(mat);
x_[level_] = x;
rhs_[level_] = rhs;
}
MultigridStep(const MatrixType& mat,
VectorType& x,
VectorType& rhs) :
LinearIterationStep<MatrixType, VectorType>(mat, x, rhs),
basesolver_(0),
preprocessCalled(false)
{}
virtual ~MultigridStep() {
for (size_t i=0; i<mat_.size()-1; i++)
virtual ~MultigridStep()
{
for (size_t i=0; i<matrixHierarchy_.size()-1; i++)
delete(ignoreNodesHierarchy_[i]);
}
......@@ -94,12 +120,13 @@
VectorType& x,
VectorType& rhs,
int numLevels)
DUNE_DEPRECATED_MSG("Use setProblem(const MatrixType& mat,VectorType& x,const VectorType& rhs).\n The number of levels is no longer set explicitely, but instead inferred from the number of transfer operators.")
{
setNumberOfLevels(numLevels);
setProblem(mat, x, rhs);
}
virtual void setNumberOfLevels(int numLevels)
DUNE_DEPRECATED_MSG("The number of levels is no longer set explicitely, but instead inferred from the number of transfer operators.\n The functionality of this function is now performed in preprocess().")
{
numLevels_ = numLevels;
......@@ -107,31 +134,28 @@
if (ignoreNodesHierarchy_[i])
delete(ignoreNodesHierarchy_[i]);
mat_.resize(numLevels);
matrixHierarchy_.resize(numLevels);
ignoreNodesHierarchy_.resize(numLevels);
for (int i=0; i<int(mat_.size())-1; i++)
for (int i=0; i<int(matrixHierarchy_.size())-1; i++)
{
mat_[i].reset();
matrixHierarchy_[i].reset();
ignoreNodesHierarchy_[i] = NULL;
}
presmoother_.resize(numLevels);
postsmoother_.resize(numLevels);
x_.resize(numLevels);
rhs_.resize(numLevels);
xHierarchy_.resize(numLevels);
rhsHierarchy_.resize(numLevels);
}
virtual void setProblem(const MatrixType& mat,
VectorType& x,
virtual void setProblem(const MatrixType& mat,
VectorType& x,
const VectorType& rhs)
{
level_ = numLevels_-1;
mat_[level_] = Dune::stackobject_to_shared_ptr(mat);
x_[level_] = x;
rhs_[level_] = rhs;
LinearIterationStep<MatrixType, VectorType, BitVectorType>::setProblem(mat,x,rhs);
// Preprocess must be called again, to create the new matrix hierarchy
preprocessCalled = false;
......@@ -139,7 +163,9 @@
void setRhs(const VectorType& rhs)
{
rhs_[numLevels_-1] = rhs;
level_ = numLevels_-1;
this->rhs_ = &rhs;
rhsHierarchy_.back() = rhs;
}
template <class DerivedTransferHierarchy>
......@@ -158,9 +184,9 @@
virtual void postprocess();
virtual VectorType getSol();
// virtual VectorType getSol();
virtual const MatrixType* getMatrix();
// virtual const MatrixType* getMatrix();
virtual int level() const {return level_;}
......@@ -173,36 +199,41 @@
virtual void setMGType(int mu, int nu1, int nu2);
/** \brief Set the smoother iteration step */
virtual void setSmoother(LinearIterationStep<MatrixType, VectorType>* smoother) {
for (size_t i=0; i<presmoother_.size(); i++)
presmoother_[i] = postsmoother_[i] = Dune::stackobject_to_shared_ptr(*smoother);
virtual void setSmoother(LinearIterationStep<MatrixType, VectorType>* smoother)
{
presmootherDefault_ = postsmootherDefault_ = Dune::stackobject_to_shared_ptr(*smoother);
levelWiseSmoothers_.clear();
}
/** \brief Set the smoother iteration step from a smart pointer*/
virtual void setSmoother(Dune::shared_ptr<LinearIterationStep<MatrixType, VectorType> > smoother) {
for (size_t i=0; i<presmoother_.size(); i++)
presmoother_[i] = postsmoother_[i] = smoother;
virtual void setSmoother(Dune::shared_ptr<LinearIterationStep<MatrixType, VectorType> > smoother)
{
presmootherDefault_ = postsmootherDefault_ = smoother;
levelWiseSmoothers_.clear();
}
/** \brief Set pre- and post smoothers individually */
virtual void setSmoother(LinearIterationStep<MatrixType, VectorType>* preSmoother,
LinearIterationStep<MatrixType, VectorType>* postSmoother) {
for (size_t i=0; i<presmoother_.size(); i++) {
presmoother_[i] = Dune::stackobject_to_shared_ptr(*preSmoother);
postsmoother_[i] = Dune::stackobject_to_shared_ptr(*postSmoother);
}
LinearIterationStep<MatrixType, VectorType>* postSmoother)
{
presmootherDefault_ = Dune::stackobject_to_shared_ptr(*preSmoother);
postsmootherDefault_ = Dune::stackobject_to_shared_ptr(*postSmoother);
levelWiseSmoothers_.clear();
}
/** \brief Set the smoother iteration step for a particular level */
virtual void setSmoother(LinearIterationStep<MatrixType, VectorType>* smoother,
int level) {
presmoother_[level] = postsmoother_[level] = Dune::stackobject_to_shared_ptr(*smoother);
virtual void setSmoother(LinearIterationStep<MatrixType, VectorType>* smoother, std::size_t level)
{
levelWiseSmoothers_[level] = Dune::stackobject_to_shared_ptr(*smoother);
}
/** \brief Set the smoother iteration step for a particular level, from a smart pointer */
virtual void setSmoother(Dune::shared_ptr<LinearIterationStep<MatrixType, VectorType> > smoother,
int level) {
presmoother_[level] = postsmoother_[level] = smoother;
virtual void setSmoother(Dune::shared_ptr<LinearIterationStep<MatrixType, VectorType> > smoother, std::size_t level)
{
levelWiseSmoothers_[level] = smoother;
}
protected:
......@@ -232,23 +263,29 @@
//! The total number of levels of the underlying grid
int numLevels_;
public:
//! The linear operators on each level
std::vector<Dune::shared_ptr<const MatrixType> > mat_;
std::vector<Dune::shared_ptr<const MatrixType> > matrixHierarchy_;
protected:
//! Flags specifying the dirichlet nodes on each level
std::vector<BitVectorType*> ignoreNodesHierarchy_;
public:
std::vector<VectorType> x_;
std::vector<Dune::shared_ptr<VectorType> > xHierarchy_;
std::vector<VectorType> rhs_;
std::vector<VectorType> rhsHierarchy_;
//protected:
std::vector<MultigridTransfer<VectorType, BitVectorType, MatrixType>* > mgTransfer_;
protected:
Dune::shared_ptr<LinearIterationStep<MatrixType, VectorType> > presmootherDefault_;
Dune::shared_ptr<LinearIterationStep<MatrixType, VectorType> > postsmootherDefault_;
typedef std::map<std::size_t, Dune::shared_ptr<LinearIterationStep<MatrixType, VectorType> > > SmootherCache;
SmootherCache levelWiseSmoothers_;
bool preprocessCalled;
};
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment