geneo.hh

#ifndef GENEO_SOLVER_HH
#define GENEO_SOLVER_HH
 
#include "conbase_fork.hh"
#include "cg_fork.hh"
#include "neumann_boundary_condition.hh"
 
#include "two_level_schwarz.hh"
 
#include "subdomainprojectedcoarsespace.hh"
 
#include "partitionofunity.hh"
#include "localoperator_ovlp_region.hh"
 
#include "geneobasis.hh"
#include "liptonbabuskabasis.hh"
#include "partitionofunitybasis.hh"
#include "zembasis.hh"
 
namespace Dune{
    namespace Geneo{
 
        template<class V, class GFS, class GO, class LOP, class CON, class MBE>
            class CG_GenEO  :
                public PDELab::OVLPScalarProductImplementation<GFS>,
                public PDELab::LinearResultStorage 
        {
 
                public:
 
                    typedef double RF;
 
                    CG_GenEO(V& x0_, const GFS& gfs_, LOP& lop_, CON& constraints_, MBE& mbe_, Composites::SolverInfo & si_, Dune::MPIHelper& helper_, int cells_) :
                        PDELab::OVLPScalarProductImplementation<GFS>(gfs_),
                        x0(x0_),
                        gfs(gfs_),
                        lop(lop_),
                        constraints(constraints_),
                        mbe(mbe_),
                        solver_info(si_),
                        helper(helper_),
                        cells(cells_) { }
 
                    void inline apply(){
 
                        typedef typename GFS::template ConstraintsContainer<double>::Type C;
                        C cg;           cg.clear();     
                        C cg_ext;       cg_ext.clear();
 
                        Dune::PDELab::constraints(constraints,gfs,cg);
 
                        std::vector<double> times;
                        Dune::Timer watch;
                        watch.reset();
 
#if HAVE_SUITESPARSE
                        //Neumann boudary conditions between processors
                        Dune::PDELab::constraints_exterior(constraints,gfs,cg_ext);
                        auto cc_bnd_neu_int_dir = C();
                        Dune::PDELab::PureNeumannBoundaryCondition pnbc;
                        cc_bnd_neu_int_dir.clear();
                        Dune::PDELab::constraints(pnbc,gfs,cc_bnd_neu_int_dir);
 
                        //Set up Grid operator
                        typedef typename GO::Jacobian J;
 
                        //Matrix with "correct" boundary conditions
                        GO go(gfs,cg,gfs,cg,lop,mbe);
                        J AF(go);
                        AF = 0.0;
                        go.jacobian(x0,AF);
 
                        //Matrix with pure neumann boundary conditions
                        GO go_neu(gfs,cg_ext,gfs,cg_ext,lop,mbe);
                        J AF_neu(go_neu);
                        AF_neu = 0.0;
                        go_neu.jacobian(x0,AF_neu);
 
                        //Create local operator on overlap
                        typedef Dune::PDELab::LocalOperatorOvlpRegion<LOP, GFS> LOP_ovlp;
                        LOP_ovlp lop_ovlp(lop,gfs);
                        typedef Dune::PDELab::GridOperator<GFS,GFS,LOP_ovlp,MBE,RF,RF,RF,C,C> GO_ovlp;
                        GO_ovlp go_ovlp(gfs,cg_ext,gfs,cg_ext,lop_ovlp,mbe);
                        typedef typename GO_ovlp::Jacobian J2;
                        J2 AF_ovlp(go_ovlp);
                        AF_ovlp = 0.0;
                        go_ovlp.jacobian(x0,AF_ovlp);
 
                        //Set up solution vector and some necessary operators
                        typedef Dune::PDELab::OverlappingOperator<C,J,V,V> POP;
                        POP popf(cg,AF);
                        typedef Dune::PDELab::istl::ParallelHelper<GFS> PIH;
                        PIH pihf(gfs);
                        typedef Dune::PDELab::OverlappingScalarProduct<GFS,V> OSP;
                        OSP ospf(gfs,pihf);
 
                        //Preconditioner
 
                        auto gv = gfs.gridView();
 
                        double eigenvalue_threshold = solver_info.eigenvalue_threshold*(double)solver_info.overlap/(cells+solver_info.overlap); //eigenvalue threshhold
                        if(helper.rank() == 0) std::cout << "Eigenvalue threshhold: " << eigenvalue_threshold << std::endl;
                        int verb = 10;
                        //if (gfs.gridView().comm().rank()==0) verb=solver_info.verb;
                        typedef Dune::PDELab::LocalFunctionSpace<GFS, Dune::PDELab::AnySpaceTag> LFSU;
                        typedef typename LFSU::template Child<0>::Type LFS;
                        LFSU lfsu(gfs);
                        LFS lfs = lfsu.template child<0>();
 
                        //file << "Time for geneo setup:" << timer.elapsed() << "\n";
                        //times.push_back(watch.elapsed());
                        //watch.reset();
 
                        std::shared_ptr<V> part_unity;
                        if (solver_info.widlund_part_unity == false)
                          part_unity = sarkisPartitionOfUnity<3,V>(gfs, lfs, cc_bnd_neu_int_dir);
                        else
                            part_unity = standardPartitionOfUnity<3,V>(gfs, lfs, cc_bnd_neu_int_dir);
 
                        //times.push_back(watch.elapsed());
                        //watch.reset();
 
                        V dummy(gfs, 1);
                        Dune::PDELab::set_constrained_dofs(cg_ext,0.0,dummy); // Zero on subdomain boundary
 
                        int nev = solver_info.nev;
 
                        std::shared_ptr<SubdomainBasis<V> > subdomain_basis;
                        //SubdomainBasis<V> subdomain_basis2;
                        if (nev > 0)
                            subdomain_basis = std::make_shared<GenEOBasis<GFS,J,V,V,Composites::SolverInfo,3> >(gfs, AF_neu, AF_ovlp, eigenvalue_threshold, *part_unity, solver_info);//,subdomain_basis2);
                        else if (nev == 0){
                            subdomain_basis = std::make_shared<ZEMBasis<GFS,LFS,V,3,3> >(gfs, lfs, *part_unity);
                            //subdomain_basis2 = *subdomain_basis;
                        }
                        else{
                            subdomain_basis = std::make_shared<PartitionOfUnityBasis<V> >(*part_unity);
                            //subdomain_basis2 = *subdomain_basis;
                        }
 
                        auto partunityspace = std::make_shared<SubdomainProjectedCoarseSpace<GFS,J,V,V,3> >(gfs, AF_neu, subdomain_basis, verb);
 
                        //typedef TwoLevelOverlappingAdditiveSchwarz<GFS,J,V,V> PREC_PCG;
                        typedef TwoLevelOverlappingAdditiveSchwarzDC<GFS,J,V,V> PREC_DC;
                        //std::shared_ptr<PREC_PCG> prec;
                        std::shared_ptr<PREC_DC> prec;
 
                        // std::shared_ptr<SubdomainProjectedCoarseSpace<GFS,J,V,V,3>> partunityspace2;
                        //partunityspace2 = std::make_shared<SubdomainProjectedCoarseSpace<GFS,J,V,V,3> >(gfs, AF_neu, std::make_shared<SubdomainBasis<V>>(subdomain_basis2), verb);
                        prec = std::make_shared<PREC_DC>(gfs, AF, partunityspace,solver_info.coarseSpaceActive);
 
                        MPI_Barrier(gfs.gridView().comm());
                        if(helper.rank()==0) std::cout << "Geneo setup time " << watch.elapsed() << std::endl;
                        times.push_back(watch.elapsed());
                        watch.reset();
 
                        // set up and assemble right hand side w.r.t. l(v)-a(u_g,v)
                        V d(gfs,0.0);
                        go.residual(x0,d);
 
                        // now solve defect equation A*v = d
                        V v(gfs,0.0);
                        //Solve using CG
 
                        std::shared_ptr<CGSolverFork<V>> solver;
                        solver = std::make_shared<CGSolverFork<V> >(popf,ospf,*prec,solver_info.KrylovTol,solver_info.MaxIt,solver_info.verb,true);
                        Dune::InverseOperatorResult result;
                        solver->apply(v,d,result);
 
                        times.push_back(watch.elapsed());
                        watch.reset();
                        x0 -= v;
 
                        solver_info.setTimes(times);
                        solver_info.recordResult(result);
 
#else
                        std::cout << "Solver not available . . . Please install UMFPACK as part of SuiteSparse" << std::endl;
                        return;
#endif
                    }
 
                    //Templated apply for Newton solver
                    template<class M, class V2, class W>
                        void apply(M& A, V2& z, W& r, typename Dune::template FieldTraits<typename V2::ElementType >::real_type reduction)
                        {
                            this->apply(z,r,reduction);
                        }
 
                private:
                    V & x0;
                    const GFS& gfs;
                    LOP& lop;
                    const CON& constraints;
                    const MBE& mbe;
                    Composites::SolverInfo& solver_info;
                    Dune::MPIHelper& helper;
                    int cells;
            };
    }
}
 
#endif