levelpatchpreconditioner.hh

#ifndef LEVEL_PATCH_PRECONDITIONER_HH
#define LEVEL_PATCH_PRECONDITIONER_HH

#include <string>

#include <dune/common/timer.hh>
#include <dune/common/fvector.hh>
#include <dune/common/bitsetvector.hh>

#include <dune/solvers/iterationsteps/lineariterationstep.hh>

#include <dune/faultnetworks/assemblers/globalfaultassembler.hh>
#include <dune/faultnetworks/patchfactories/supportpatchfactory.hh>
#include <dune/faultnetworks/localproblem.hh>
#include <dune/faultnetworks/levelinterfacenetwork.hh>
#include <dune/faultnetworks/utils/debugutils.hh>

#include <dune/fufem/boundarypatch.hh>
#include <dune/fufem/functiontools/boundarydofs.hh>


template <class BasisType, class LocalAssembler, class LocalInterfaceAssembler, class MatrixType, class VectorType>
class LevelPatchPreconditioner : public LinearIterationStep<MatrixType, VectorType> {

public:
    enum Mode {additive, multiplicative};
    enum BoundaryMode {homogeneous, fromIterate};

private:
    typedef typename BasisType::GridView GridView;
    typedef typename GridView::Grid GridType;

    const LevelInterfaceNetwork<GridView>& levelInterfaceNetwork_;
    const BasisType& patchLevelBasis_;
    const LocalAssembler& localAssembler_;
    const std::vector<std::shared_ptr<LocalInterfaceAssembler>>& localInterfaceAssemblers_;
    const Mode mode_;

    const GridType& grid_;
    const int level_;
    const BasisType basis_;

    size_t patchDepth_;
    BoundaryMode boundaryMode_;

    MatrixType matrix_;
    std::vector<std::shared_ptr<OscLocalProblem<MatrixType, VectorType>> > localProblems_;

public:

    // for each coarse patch given by patchLevelBasis: set up local problem, compute local correction
    LevelPatchPreconditioner(const LevelInterfaceNetwork<GridView>& levelInterfaceNetwork,
                             const BasisType& patchLevelBasis,
                             const LocalAssembler& localAssembler,
                             const std::vector<std::shared_ptr<LocalInterfaceAssembler>>& localInterfaceAssemblers,
                             const Mode mode = LevelPatchPreconditioner::Mode::additive) :
          levelInterfaceNetwork_(levelInterfaceNetwork),
          patchLevelBasis_(patchLevelBasis),
          localAssembler_(localAssembler),
          localInterfaceAssemblers_(localInterfaceAssemblers),
          mode_(mode),
          grid_(levelInterfaceNetwork_.grid()),
          level_(levelInterfaceNetwork_.level()),
          basis_(levelInterfaceNetwork_)
    {
        setPatchDepth();
        setBoundaryMode();

        assert(localInterfaceAssemblers_.size() == levelInterfaceNetwork_.size());
    }

    void build() {
        // assemble stiffness matrix for entire level including all faults
        GlobalFaultAssembler<BasisType, BasisType> globalFaultAssembler(basis_, basis_, levelInterfaceNetwork_);
        globalFaultAssembler.assembleOperator(localAssembler_, localInterfaceAssemblers_, matrix_);

        // set boundary conditions
        Dune::BitSetVector<1> globalBoundaryDofs;
        BoundaryPatch<GridView> boundaryPatch(levelInterfaceNetwork_.levelGridView(), true);
        constructBoundaryDofs(boundaryPatch, basis_, globalBoundaryDofs);

        typedef typename MatrixType::row_type RowType;
        typedef typename RowType::ConstIterator ColumnIterator;

        for(size_t i=0; i<globalBoundaryDofs.size(); i++) {
            if(!globalBoundaryDofs[i][0])
                continue;

            RowType& row = matrix_[i];

            ColumnIterator cIt    = row.begin();
            ColumnIterator cEndIt = row.end();

            for(; cIt!=cEndIt; ++cIt) {
                row[cIt.index()] = 0;
            }
            row[i] = 1;
        }

        // init vertexInElements
        const int dim = GridType::dimension;
        typedef typename GridType::template Codim<0>::Entity EntityType;
        std::vector<std::vector<EntityType>> vertexInElements;

        const GridView& patchLevelGridView = patchLevelBasis_.getGridView();
        vertexInElements.resize(patchLevelGridView.size(dim));
        for (size_t i=0; i<vertexInElements.size(); i++) {
            vertexInElements[i].resize(0);
        }

        typedef typename BasisType::LocalFiniteElement FE;
        typedef  typename GridView::template Codim <0>::Iterator  ElementLevelIterator;
        ElementLevelIterator endElemIt = patchLevelGridView.template end <0>();
        for (ElementLevelIterator  it = patchLevelGridView. template begin <0>(); it!=endElemIt; ++it) {

            // compute coarseGrid vertexInElements
            const FE& coarseFE = patchLevelBasis_.getLocalFiniteElement(*it);

            for (size_t i=0; i<coarseFE.localBasis().size(); ++i) {
                int dofIndex = patchLevelBasis_.indexInGridView(*it, i);
                vertexInElements[dofIndex].push_back(*it);
            }
        }

        localProblems_.resize(vertexInElements.size());

        std::cout << std::endl;
        std::cout << "---------------------------------------------" << std::endl;
        std::cout << "Initializing local fine grid corrections! Level: " << level_ << std::endl;

        // init local fine level corrections
        Dune::Timer timer;
        timer.reset();
        timer.start();

        const int patchLevel = patchLevelBasis_.faultNetwork().level();
        for (size_t i=0; i<vertexInElements.size(); i++) {
            //std::cout << i << std::endl;
            //std::cout << "---------------" << std::endl;

            SupportPatchFactory<BasisType> patchFactory(patchLevelBasis_, basis_, patchLevel, level_, vertexInElements, i, patchDepth_);
            std::vector<int>& localToGlobal = patchFactory.getLocalToGlobal();
            Dune::BitSetVector<1>& boundaryDofs = patchFactory.getBoundaryDofs();
            VectorType rhs;
            rhs.resize(matrix_.N());
            rhs = 0;

            //print(localToGlobal, "localToGlobal: ");
            //print(boundaryDofs, "boundaryDofs: ");

            localProblems_[i] = std::make_shared<OscLocalProblem<MatrixType, VectorType>>(matrix_, rhs, localToGlobal, boundaryDofs);

            /*
            if ((i+1) % 10 == 0) {
                std::cout << '\r' << std::floor((i+1.0)/patchLevelBasis_.size()*100) << " % done. Elapsed time: " << timer.elapsed() << " seconds. Predicted total setup time: " << timer.elapsed()/(i+1.0)*patchLevelBasis_.size() << " seconds." << std::flush;
            }*/
        }

        std::cout << std::endl;
        std::cout << "Total setup time: " << timer.elapsed() << " seconds." << std::endl;
        std::cout << "---------------------------------------------" << std::endl;
        timer.stop();
    }

    void setPatchDepth(const size_t patchDepth = 0) {
        patchDepth_ = patchDepth;
    }

    void setBoundaryMode(const BoundaryMode boundaryMode = LevelPatchPreconditioner::BoundaryMode::homogeneous) {
        boundaryMode_ = boundaryMode;
    }

    virtual void setProblem(const MatrixType& mat, VectorType& x, const VectorType& rhs) {
        this->x_ = &x;
        this->rhs_ = &rhs;
        this->mat_ = Dune::stackobject_to_shared_ptr(mat);

        for (size_t i=0; i<localProblems_.size(); i++) {
            if (boundaryMode_ == BoundaryMode::homogeneous)
                localProblems_[i]->updateRhs(rhs);
            else
                localProblems_[i]->updateRhsAndBoundary(rhs, x);
        }
    }

    virtual void iterate() {
        if (mode_ == additive)
            iterateAdd();
        else
            iterateMult();
    }

    void iterateAdd() {
        *(this->x_) = 0;	

        VectorType it, x;
        for (size_t i=0; i<localProblems_.size(); i++) {
            localProblems_[i]->solve(it);
            localProblems_[i]->prolong(it, x);

            /*if (i==5) {
                writeToVTK(basis_, x, "/storage/mi/podlesny/data/faultnetworks/iterates/", "exactvertexdata_patchDepth_"+std::to_string(patchDepth_));
            }*/

            *(this->x_) += x;
        }
    }

    void iterateMult() {
        *(this->x_) = 0;

        VectorType it, x;
        for (size_t i=0; i<localProblems_.size(); i++) {
            VectorType updatedRhs(*(this->rhs_));
            matrix_.mmv(*(this->x_), updatedRhs);

            //step(i);
            //print(updatedRhs, "localRhs: ");
            //writeToVTK(basis_, updatedRhs, "/storage/mi/podlesny/data/faultnetworks/rhs/local", "exactvertexdata_step_"+std::to_string(i));

            if (boundaryMode_ == BoundaryMode::homogeneous)
                localProblems_[i]->updateRhs(updatedRhs);
            else
                localProblems_[i]->updateRhsAndBoundary(updatedRhs, *(this->x_));

            localProblems_[i]->solve(it);
            localProblems_[i]->prolong(it, x);


            //print(it, "it: ");
            //print(x, "x: ");

            //writeToVTK(basis_, x, "/storage/mi/podlesny/data/faultnetworks/sol/local", "exactvertexdata_step_"+std::to_string(i));

            /*if (i==5) {
                writeToVTK(basis_, x, "/storage/mi/podlesny/data/faultnetworks/iterates/", "exactvertexdata_patchDepth_"+std::to_string(patchDepth_));
            }*/

            *(this->x_) += x;
        }
    }

    const BasisType& basis() const {
        return basis_;
    }

    const GridView& gridView() const {
        return basis_.getGridView();
    }

    const LevelInterfaceNetwork<GridView>& levelInterfaceNetwork() const {
        return levelInterfaceNetwork_;
    }

    size_t size() const {
        return localProblems_.size();
    }
};

#endif