coarse_prec.hh

#ifndef COARSE_PREC_HH
#define COARSE_PREC_HH
 
#if HAVE_ARPACKPP
 
#include <dune/pdelab/boilerplate/pdelab.hh>
 
#include <dune/common/timer.hh>
 
#include "coarsespace.hh"
#include "subdomainprojectedcoarsespace.hh"
 
 
#include <dune/common/deprecated.hh>
#include <dune/common/parallel/mpihelper.hh>
 
#include <dune/istl/owneroverlapcopy.hh>
#include <dune/istl/solvercategory.hh>
#include <dune/istl/operators.hh>
#include <dune/istl/solvers.hh>
#include <dune/istl/solver.hh>
#include <dune/istl/solvertype.hh>
#include <dune/istl/preconditioners.hh>
#include <dune/istl/preconditioner.hh>
#include <dune/istl/scalarproducts.hh>
#include <dune/istl/paamg/amg.hh>
#include <dune/istl/paamg/pinfo.hh>
#include <dune/istl/io.hh>
#include <dune/istl/superlu.hh>
#include <dune/istl/umfpack.hh>
 
#include <dune/istl/istlexception.hh>
#include <dune/istl/matrixutils.hh>
#include <dune/istl/gsetc.hh>
#include <dune/istl/ilu.hh>
 
#include <dune/pdelab/constraints/common/constraints.hh>
#include <dune/pdelab/gridfunctionspace/genericdatahandle.hh>
#include <dune/pdelab/backend/solver.hh>
#include <dune/pdelab/backend/istl/vector.hh>
#include <dune/pdelab/backend/istl/bcrsmatrix.hh>
 
 
namespace Dune {
    namespace Geneo {
 
        template<class Coarse, class COARSE_V, class COARSE_W>
            class CoarseTransfer{
                public:
                    CoarseTransfer(std::shared_ptr<Coarse>& coarse, std::shared_ptr<Coarse>& coarse_small, int ranks)
                        : coarse_(coarse),
                        coarse_small_(coarse_small),
                        ranks_(ranks)
                {}
 
 
                    /*
                     *  Global Mappings
                     */
                    /*
                     * Restriction mapping from coarse to coarse_small
                     */
                    COARSE_V restrict(const COARSE_V& v){
                        //create restricted v
                        int cnt_small=0, cnt=0;
                        COARSE_V v0(coarse_small_->basis_size(),coarse_small_->basis_size());
                        for (int rank = 0; rank < ranks_; rank++){
                            for (int i = 0; i< coarse_->get_local_basis_sizes(rank); i++){
                                if(i < coarse_small_->get_local_basis_sizes(rank)){
                                    v0[cnt_small] = v[cnt];
                                    cnt_small++;
                                }
                                cnt++;
                            }
                        }
                        return v0;
                    }
 
                    /*
                     * Prolongation mapping from coarse_small to coarse
                     */
                    COARSE_W prolongate(COARSE_W& w, const COARSE_W& v){
                        //create prolongated w
                        COARSE_W w0(coarse_->basis_size(),coarse_->basis_size());
                        int cnt_small = 0, cnt = 0;
                        for (int rank = 0; rank < ranks_; rank++){
                            for (int i = 0; i< coarse_->get_local_basis_sizes(rank); i++){
                                w0[cnt] = v[cnt];
                                if(i < coarse_small_->get_local_basis_sizes(rank)){
                                    w0[cnt] = w[cnt_small];
                                    cnt_small++;
                                }
                                cnt++;
                            }
                        }
                        return w0;
                    }
                    /*
                     *  Local Mappings
                     */
                    /*
                     * Restriction mapping from coarse to coarse_small
                     */
                    COARSE_V restrict_local(const COARSE_V& v, int offset, int size){
                        //create restricted v
                        COARSE_V v0(size,size);
                        for (int i = 0; i< size; i++){
                            v0[i] = v[offset+i];
                        }
                        return v0;
                    }
 
                    /*
                     * Prolongation mapping from coarse_small to coarse
                     */
                    COARSE_W prolongate_local(COARSE_W& w, int offset, int size){
                        //create prolongated w
                        COARSE_W w0(coarse_->basis_size(),coarse_->basis_size());
                        int cnt_small = 0, cnt = 0;
                        for (int i = 0; i< coarse_->basis_size(); i++){
                            w0[i] = 0;
                            if(i >= offset && i < offset + size)
                                w0[i] = w[i-offset];
                        }
                        return w0;
                    }
 
 
                private:
                    std::shared_ptr<Coarse> coarse_;
                    std::shared_ptr<Coarse> coarse_small_;
                    int ranks_;
            };
 
 
        template<class GFS, class COARSE_M, class COARSE_V, class COARSE_W, class M, class X, class Y>
            class CoarsePrec : public Preconditioner<COARSE_V,COARSE_W> {
                public:
                    typedef Dune::PDELab::Backend::Native<COARSE_M> matrix_type;
 
                    // define the category
                    enum {
                        category=SolverCategory::sequential
                    };
 
                    CoarsePrec (const GFS& gfs, const COARSE_M& A, std::shared_ptr<CoarseSpace<M,X> > coarse_space, std::shared_ptr<CoarseSpace<M,X>> coarse_space_small) 
                        :   my_rank_(gfs.gridView().comm().rank()),
                        ranks_(gfs.gridView().comm().size()),
                        gfs_(gfs),
                        coarse_transfer(coarse_space_,coarse_space_small_,ranks_),
                        CA(Dune::PDELab::Backend::native(A)),
                        coarse_space_(coarse_space), 
                        coarse_space_small_(coarse_space_small),
                        solverc(CA),
                        local_matrix(ranks_), solver_r(ranks_)
                {
                    for (int rank = 0; rank < ranks_; rank++){
                        //local size and offset
                        int local_size = coarse_space_->get_local_basis_sizes(rank);
                        auto offset = coarse_space_->basis_array_offset(rank);
 
                        //Set up local matrix
                        local_matrix[rank] = std::make_shared<COARSE_M>(setupLocalMatrix(offset, local_size));
 
                        solver_r[rank] = std::make_shared<Dune::UMFPack<COARSE_M>>(*local_matrix[rank]);
                    }
                }
 
                    virtual void pre (COARSE_V& x, COARSE_W& b)
                    {
                        DUNE_UNUSED_PARAMETER(x);
                        DUNE_UNUSED_PARAMETER(b);
                    }
 
                    virtual void apply (COARSE_V& v, const COARSE_W& d)
                    {
                        Dune::InverseOperatorResult result;
                        //Loop over subspaces and invert local matrix
                        if(configuration.get<bool>("prec_subdomain",true)!=true){
                            for (int rank = 0; rank < ranks_; rank++){
                                //local size and offset
                                int local_size = coarse_space_->get_local_basis_sizes(rank);
                                auto offset = coarse_space_->basis_array_offset(rank);
 
                                //restrict the right hand side to the local subspace
                                auto btilde = coarse_transfer.restrict_local(d,offset,local_size);
 
                                //solve the linear system
                                COARSE_V vl(local_size,local_size);
                                solver_r[rank]->apply(vl,btilde,result);
 
                                //prolongate the solution and add to v
                                auto local_correction = coarse_transfer.prolongate_local(vl,offset,local_size);
                                v += local_correction;
                            }
                        }
                        // Solve coarse system
                        auto btilde = coarse_transfer.restrict(d);
                        COARSE_V v0(coarse_space_small_->basis_size(),coarse_space_small_->basis_size());
                        solverc.apply(v0,btilde,result);
                        auto coarse_correction = coarse_transfer.prolongate(v0,d);
                        v = coarse_correction;
 
                    }
                private:
                    COARSE_M setupLocalMatrix(int offset, int local_size){
                        //Set up local matrix for rank
                        COARSE_M local_matrix(local_size,local_size,COARSE_M::row_wise );
                        auto A = *coarse_space_->get_coarse_system();
                        auto setup_row = local_matrix.createbegin();
                        for (int idx = 0; idx < local_size; idx++){
                            for (int j = 0; j < local_size; j++)
                                setup_row.insert(j);
                            ++setup_row;
                            for (int j = 0; j < local_size; j++){
                                if(A.exists(offset+idx,offset+j))
                                    local_matrix[idx][j] = A[offset+idx][offset+j];
                                else
                                    local_matrix[idx][j] = 0;
                            }
                        }
                        return local_matrix;
                    }
 
                    virtual void post (COARSE_V& x)
                    {
                        DUNE_UNUSED_PARAMETER(x);
                    }
 
                private:
                    const GFS& gfs_;
                    int my_rank_;
                    int ranks_;
                    matrix_type CA;
                    std::shared_ptr<CoarseSpace<M,X> > coarse_space_;
                    std::shared_ptr<CoarseSpace<M,X>> coarse_space_small_;
                    Dune::UMFPack<matrix_type> solverc;
                    CoarseTransfer< CoarseSpace<M,X>,COARSE_V,COARSE_W > coarse_transfer;
                    std::vector<std::shared_ptr<COARSE_M> > local_matrix;
                    std::vector<std::shared_ptr<Dune::UMFPack<COARSE_M>> > solver_r;
            };
 
 
    }//end namespace Dune
 
}
#endif //have arpack
 
#endif