multigridstep.cc

#include <typeinfo>

#include <dune/solvers/transferoperators/multigridtransfer.hh>
#include <dune/solvers/solvers/loopsolver.hh>
#include <dune/solvers/common/genericvectortools.hh>
#include "blockgsstep.hh"

#ifdef HAVE_IPOPT
#include <dune/solvers/solvers/quadraticipopt.hh>
#endif

//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
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;

    numLevels_ = mgTransfer_.size() + 1;

    this->level_ = this->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 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->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->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
    // /////////////////////////////////////////////////////
    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_);
        for (int i=this->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
    // /////////////////////////////////////////////

    if (this->basesolver_ == NULL)
        DUNE_THROW(SolverError, "You have not provided a base solver!");

    // If the base solver can ignore dofs give it the ignoreNodes field
    if (dynamic_cast<CanIgnore<BitVectorType>*>(this->basesolver_))
        dynamic_cast<CanIgnore<BitVectorType>*>(this->basesolver_)->ignoreNodes_ = ignoreNodesHierarchy_[0];

    typedef ::LoopSolver<VectorType> DuneSolversLoopSolver;

    if (typeid(*this->basesolver_) == typeid(DuneSolversLoopSolver)) {

        DuneSolversLoopSolver* loopBaseSolver = dynamic_cast<DuneSolversLoopSolver*> (this->basesolver_);
     
        typedef LinearIterationStep<MatrixType, VectorType> SmootherType;
        assert(dynamic_cast<SmootherType*>(loopBaseSolver->iterationStep_));

        dynamic_cast<SmootherType*>(loopBaseSolver->iterationStep_)->setProblem(*(this->matrixHierarchy_[0]), *this->xHierarchy_[0], this->rhsHierarchy_[0]);
        dynamic_cast<SmootherType*>(loopBaseSolver->iterationStep_)->ignoreNodes_ = ignoreNodesHierarchy_[0];
        
    }
#ifdef HAVE_IPOPT
    else if (typeid(*this->basesolver_) == typeid(QuadraticIPOptSolver<MatrixType,VectorType>)) {

        QuadraticIPOptSolver<MatrixType,VectorType>* ipoptBaseSolver = dynamic_cast<QuadraticIPOptSolver<MatrixType,VectorType>*> (this->basesolver_);
     
        ipoptBaseSolver->setProblem(*(this->matrixHierarchy_[0]), *this->xHierarchy_[0], this->rhsHierarchy_[0]);
    } 
#endif
else {
        DUNE_THROW(SolverError, "You can't use " << typeid(*this->basesolver_).name()
                   << " as a base solver for contact problems!");
    }
}


template<class MatrixType, class VectorType, class BitVectorType>
void MultigridStep<MatrixType, VectorType, BitVectorType>::nestedIteration()
{
    for (level_ = 0; level_<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<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) {
        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();

    // /////////////////////////
    // 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();

}