mmgstep.cc

// -*- tab-width: 8; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=8 sw=4 sts=4:
#include <dune/solvers/transferoperators/truncatedmgtransfer.hh>
#include <dune/solvers/iterationsteps/projectedblockgsstep.hh>

#include <dune/solvers/computeenergy.hh>


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

template <class MatrixType, class VectorType>
void MonotoneMGStep<MatrixType, VectorType>::
preprocess()
{
    // Unset the recompute bitfields, so we compute the full stiffness matrix hierarchy at the beginning
    for (size_t i=0; i<this->mgTransfer_.size(); i++) {
        std::dynamic_pointer_cast<TruncatedMGTransfer<VectorType> >(this->mgTransfer_[i])->setRecomputeBitField(nullptr);
        std::dynamic_pointer_cast<TruncatedMGTransfer<VectorType> >(this->mgTransfer_[i])->setCriticalBitField(nullptr);
    }

    // Call preprocess of the base class
    MultigridStep<MatrixType,VectorType>::preprocess();

    // Then setup the obstacle flags which specify the subset of entries to be reassembled at each iteration
    hasObstacleHierarchy_.resize(this->numLevels());
    hasObstacleHierarchy_.back() = hasObstacle_;

    for (int i=this->mgTransfer_.size()-1; i>=0; i--) {
        hasObstacleHierarchy_[i] = std::make_shared<Dune::BitSetVector<dim> >();
        this->mgTransfer_[i]->restrict(*hasObstacleHierarchy_[i+1], *hasObstacleHierarchy_[i]);
        for (size_t j=0; j<hasObstacleHierarchy_[i]->size(); j++)
            if ((*this->ignoreNodesHierarchy_[i])[j].any())
                (*hasObstacleHierarchy_[i])[j][0]=false;
    }

    recompute_.resize(this->mgTransfer_.size());
    for (size_t i=0; i<this->mgTransfer_.size(); i++)
    {
        recompute_[i].resize(hasObstacleHierarchy_[i]->size());
        recompute_[i].unsetAll();
        std::dynamic_pointer_cast<TruncatedMGTransfer<VectorType> >(this->mgTransfer_[i])->setRecomputeBitField(&recompute_[i]);
    }

    oldCritical_.resize(this->numLevels());
    for (size_t i=0; i<this->numLevels(); i++)
    {
      oldCritical_[i].resize(hasObstacleHierarchy_[i]->size());
      oldCritical_[i].unsetAll();
    }

    // /////////////////////////////////////////////
    //   Set up base solver
    // /////////////////////////////////////////////

    obstacleHierarchy_.resize(this->numLevels());
    obstacleHierarchy_.back() = obstacles_;

    for (size_t i=0; i<obstacleHierarchy_.size()-1; i++)
        obstacleHierarchy_[i] = std::make_shared<ObstacleVectorType>();

    if (typeid(*this->basesolver_) == typeid(::LoopSolver<VectorType>)) {

        ::LoopSolver<VectorType>* loopBaseSolver = dynamic_cast< ::LoopSolver<VectorType>* > (this->basesolver_.get());

        typedef ProjectedBlockGSStep<MatrixType, VectorType> SmootherType;

        dynamic_cast<SmootherType*>(&loopBaseSolver->getIterationStep())->hasObstacle_ = hasObstacleHierarchy_[0].get();
        dynamic_cast<SmootherType*>(&loopBaseSolver->getIterationStep())->obstacles_   = obstacleHierarchy_[0].get();

#if HAVE_IPOPT
    } else if (typeid(*this->basesolver_) == typeid(QuadraticIPOptSolver<MatrixType,VectorType>)) {

        QuadraticIPOptSolver<MatrixType,VectorType>* ipoptBaseSolver = dynamic_cast<QuadraticIPOptSolver<MatrixType,VectorType>*> (this->basesolver_.get());

        ipoptBaseSolver->setObstacles(*obstacleHierarchy_[0]);
#endif
    } else {
        DUNE_THROW(SolverError, "You can't use " << typeid(*this->basesolver_).name()
                   << " as a base solver for obstacle problems!");
    }
}

template<class MatrixType, class VectorType>
void MonotoneMGStep<MatrixType, VectorType>::nestedIteration()
{
    for (int i=this->numLevels()-1; i>0; i--) {

        Dune::BitSetVector<dim> dummy(this->rhsHierarchy_[i].size(), false);

        // /////////////////////////////
        //   Restrict obstacles
        // //////////////////////////////
        obstacleRestrictor_->restrict(*obstacleHierarchy_[i], *obstacleHierarchy_[i-1],
                                    *hasObstacleHierarchy_[i], *hasObstacleHierarchy_[i-1],
                                    *this->mgTransfer_[i-1], dummy);

        // Restrict right hand side
        this->mgTransfer_[i-1]->restrict(this->rhsHierarchy_[i], this->rhsHierarchy_[i-1]);
    }

    for (this->level_ = 0; this->level_<(int) this->numLevels()-1; this->level_++) {

        // If we start from an infeasible configuration, the restricted
        // obstacles may be inconsistent.  We do an ad hoc correction here.
        // The true obstacles on the maxlevel are not touched.
        for (size_t i=0; i<obstacleHierarchy_[this->level_]->size(); i++) {
            BoxConstraint<field_type,dim>& cO = (*obstacleHierarchy_[this->level_])[i];
            for (int j=0; j<dim; j++)
                if (cO.lower(j) > cO.upper(j))
                    cO.lower(j) = cO.upper(j) = (cO.lower(j) + cO.upper(j)) / 2;
        }

        std::cout << "Nested iteration on level " << this->level_ << std::endl;
        iterate();
        iterate();
        //std::cout << this->xHierarchy_[this->level_] << std::endl;
        this->mgTransfer_[this->level_]->prolong(*this->xHierarchy_[this->level_], *this->xHierarchy_[this->level_+1]);

    }

}


template <class MatrixType, class VectorType>
void MonotoneMGStep<MatrixType, VectorType>::iterate()
{
    this->preprocessCalled = false;

    int& level = this->level_;

    // Define references just for ease of notation
    std::vector<std::shared_ptr<const MatrixType> >& mat = this->matrixHierarchy_;
    std::vector<std::shared_ptr<VectorType> >& x   = this->xHierarchy_;
    std::vector<VectorType>& rhs = this->rhsHierarchy_;
    auto& obstacles = obstacleHierarchy_;

    Dune::BitSetVector<dim> critical(x[level]->size(), false);

    // Solve directly if we're looking at the coarse problem
    if (level == 0) {

        this->basesolver_->solve();

        // Determine critical nodes  (only for debugging)
        const double eps = 1e-12;
        for (size_t i=0; i<obstacles[level]->size(); i++) {

            for (int j=0; j<dim; j++) {

                if ((*obstacles[level])[i].lower(j) >= (*x[level])[i][j] - eps
                    || (*obstacles[level])[i].upper(j) <= (*x[level])[i][j] + eps) {
                    critical[i][j] = true;
                }

            }
        }

    } else {

        // Presmoothing
        ProjectedBlockGSStep<MatrixType,VectorType>* presmoother = dynamic_cast<ProjectedBlockGSStep<MatrixType, VectorType>*>(this->presmoother_[level].get());
        assert(presmoother);
        presmoother->setProblem(*(mat[level]), *x[level], rhs[level]);
        presmoother->obstacles_ = obstacles[level].get();
        presmoother->hasObstacle_ = hasObstacleHierarchy_[level].get();
        presmoother->ignoreNodes_ = this->ignoreNodesHierarchy_[level];

        for (int i=0; i<this->nu1_; i++)
            presmoother->iterate();

        // First, a backup of the obstacles for postsmoothing
        std::vector<BoxConstraint<field_type,dim> > obstacleBackup = *obstacles[level];

        // Compute defect obstacles
        for (size_t i=0; i<obstacles[level]->size(); i++)
            (*obstacles[level])[i] -= (*x[level])[i];

        // ///////////////////////
        //    Truncation
        // ///////////////////////

        // Determine critical nodes
        const double eps = 1e-12;
        for (size_t i=0; i<obstacles[level]->size(); i++) {

            for (int j=0; j<dim; j++) {

                if ((*obstacles[level])[i].lower(j) >= -eps
                    || (*obstacles[level])[i].upper(j) <= eps) {
                    critical[i][j] = true;
                }

            }
        }

        ///////////////////////////////////////////////////////////////////////////////
        // Compute the part of the coarse grid matrix that needs to be recomputed.
        // There are two reasons why a matrix entry may need to be recomputed:
        // 1) A corresponding fine grid vertex switched from critical to non-critical or vice versa
        // 2) A corresponding fine grid matrix entry got recomputed
        ///////////////////////////////////////////////////////////////////////////////

        Dune::BitSetVector<dim> changed(critical.size());
        for (size_t i=0; i<changed.size(); i++)
          for (int j=0; j<dim; j++)
            changed[i][j] = (critical[i][j] != oldCritical_[level][i][j]);

        if (level < (int) this->numLevels()-1 )
          for (size_t i=0; i<changed.size(); i++)
            for (int j=0; j<dim; j++)
              changed[i][j] = (changed[i][j] || recompute_[level][i][j]);

        std::dynamic_pointer_cast<TruncatedMGTransfer<VectorType> >(this->mgTransfer_[level-1])->restrict(changed, recompute_[level-1]);

        oldCritical_[level] = critical;

        // Set bitfield of nodes that will be truncated
        std::dynamic_pointer_cast<TruncatedMGTransfer<VectorType> >(this->mgTransfer_[level-1])->setCriticalBitField(&critical);

        // Restrict stiffness matrix
        this->mgTransfer_[level-1]->galerkinRestrict(*(mat[level]), *(const_cast<MatrixType*>(mat[level-1].get())));

        // Restrict obstacles
        obstacleRestrictor_->restrict(*obstacles[level], *obstacles[level-1],
                                    *hasObstacleHierarchy_[level], *hasObstacleHierarchy_[level-1],
                                    *this->mgTransfer_[level-1],
                                    critical);

        // //////////////////////////////////////////////////////////////////////
        // Restriction:  fineResiduum = 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]);

        // Choose all zeros as the initial correction
        *x[level-1] = 0;

        // ///////////////////////////////////////
        // Recursively solve the coarser system
        // ///////////////////////////////////////
        level--;
        for (int i=0; i<this->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;

        // Remove pointer to the temporary critical bitfield
        // this avoids a memory problem when the same mmg step is reused
        std::dynamic_pointer_cast<TruncatedMGTransfer<VectorType> >(this->mgTransfer_[level-1])->setCriticalBitField(nullptr);

        // restore the true (non-defect) obstacles
        *obstacles[level] = obstacleBackup;

#ifndef NDEBUG
        // Debug: is the current iterate really admissible?
        for (size_t i=0; i<obstacles[level]->size(); i++)
            for (int j=0; j<VectorType::block_type::dimension; j++)
                if (((*x[level])[i][j] - (*obstacles[level])[i].lower(j)<-1e-14
                    || (*x[level])[i][j] - (*obstacles[level])[i].upper(j) >1e-14 )
                    && (!(*this->ignoreNodesHierarchy_[level])[i][j])) {

                    std::cout << "Obstacle disregarded!\n";
                    std::cout << (*x[level])[i] << std::endl << (*obstacles[level])[i] << std::endl;
                    std::cout << "level: " << level << "   index: " << i << "   komponent: " << j << std::endl;
                    std::cout << "is " << ((critical[i][j]) ? "" : "not") << "critical" << std::endl;
                }
#endif

        // Postsmoothing
        ProjectedBlockGSStep<MatrixType,VectorType>* postsmoother = dynamic_cast<ProjectedBlockGSStep<MatrixType, VectorType>*>(this->postsmoother_[level].get());
        assert(postsmoother);
        postsmoother->setProblem(*(mat[level]), *x[level], rhs[level]);
        postsmoother->obstacles_ = obstacles[level].get();
        postsmoother->hasObstacle_ = hasObstacleHierarchy_[level].get();
        postsmoother->ignoreNodes_ = this->ignoreNodesHierarchy_[level];

        for (int i=0; i<this->nu2_; i++)
            postsmoother->iterate();

    }

    // ////////////////////////////////////////////////////////////////////
    //   Track the number of critical nodes found during this iteration
    // ////////////////////////////////////////////////////////////////////
    if (level==(int) this->numLevels()-1 && this->verbosity_==NumProc::FULL) {

        std::cout << critical.count() << " critical nodes found on level " << level;

        int changes = 0;
        for (unsigned int i=0; i<oldCritical.size(); i++)
            for (int j=0; j<dim; j++)
                if (oldCritical[i][j] !=critical[i][j])
                    changes++;

        std::cout << ", and " << changes << " changes." << std::endl;
        oldCritical = critical;
    }

    // Debug: output energy
    if (level==(int) this->numLevels()-1 && this->verbosity_==NumProc::FULL) {
        std::streamsize const oldPrecision = std::cout.precision();
        std::cout << "Total energy: "
                  << std::setprecision(10)
                  << computeEnergy(*mat[level], *x[level], rhs[level]) << std::endl
                  << std::setprecision(oldPrecision);
    }

}