geneo.hh

// -*- tab-width: 4; indent-tabs-mode: nil -*-
// vi: set et ts=4 sw=4 sts=4:
#ifndef GENEO_SOLVER_HH
#define GENEO_SOLVER_HH
 
#include <dune/pdelab/backend/istl/geneo/geneo.hh>
#include "zembasis.hh"
 
namespace Dune{
    namespace Geneo{
 
        template<class V, class GFS, class GFS_EXT, class GO, class GO_EXT, 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_, const GFS_EXT gfs_ext_, LOP& lop_, CON& constraints_, MBE& mbe_, Composites::SolverInfo & si_, Dune::MPIHelper& helper_, int cells_) :
                PDELab::OVLPScalarProductImplementation<GFS>(gfs_),
                x0(x0_),
                gfs(gfs_),
                gfs_ext(gfs_ext_),
                lop(lop_),
                constraints(constraints_),
                mbe(mbe_),
                solver_info(si_),
                helper(helper_),
                cells(cells_) { }
 
            inline void apply()
            {
#if HAVE_SUITESPARSE
                typedef typename GFS::template ConstraintsContainer<RF>::Type CC;
                auto cg = CC();
                typedef typename GFS_EXT::template ConstraintsContainer<RF>::Type CC_EXTERIOR;
                auto cg_ext = CC_EXTERIOR();
                auto cc_bnd_neu_int_dir = CC();
 
                // assemble constraints
                Dune::PDELab::constraints(constraints,gfs,cg);
                Dune::PDELab::constraints(constraints,gfs_ext,cg_ext);
 
                Dune::PDELab::NoDirichletConstraintsParameters pnbc;
                Dune::PDELab::constraints(pnbc,gfs,cc_bnd_neu_int_dir);
 
                std::vector<double> times;
                Dune::Timer watch;
                watch.reset();
 
                // also create a wrapper using the same data based on the GFS without processor boundary constraints
                typedef Dune::PDELab::Backend::Vector<GFS_EXT, RF> V_EXTERIOR;
                V_EXTERIOR x0_ext(gfs_ext, Dune::PDELab::Backend::unattached_container());
                x0_ext.attach(x0.storage());
 
                //Set up Grid operator
                typedef typename GO::Jacobian J;
                typedef typename GO_EXT::Jacobian J_ext;
 
                //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_EXT go_neu(gfs_ext,cg_ext,gfs_ext,cg_ext,lop,mbe);
                J_ext AF_neu(go_neu);
                AF_neu = 0.0;
                go_neu.jacobian(x0_ext,AF_neu);
 
                //Create local operator on overlap
                typedef Dune::PDELab::LocalOperatorOvlpRegion<LOP, GFS_EXT> LOP_ovlp;
                LOP_ovlp lop_ovlp(lop,gfs_ext);
                typedef Dune::PDELab::GridOperator<GFS_EXT,GFS_EXT,LOP_ovlp,MBE,RF,RF,RF,CC,CC> GO_ovlp;
                GO_ovlp go_ovlp(gfs_ext,cg_ext,gfs_ext,cg_ext,lop_ovlp,mbe);
                J_ext AF_ovlp(go_ovlp);
                AF_ovlp = 0.0;
                go_ovlp.jacobian(x0_ext,AF_ovlp);
 
                //Set up solution vector and some necessary operators
                typedef Dune::PDELab::OverlappingOperator<CC,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);
 
                // 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){
                    if(eigenvalue_threshold > 0) std::cout << "Eigenvalue threshhold: " << eigenvalue_threshold << std::endl;
                    else std::cout << "Using " << solver_info.nev << " EV in each subdomain" << std::endl;
                }
                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>();
 
                std::shared_ptr<V> part_unity;
                part_unity = std::make_shared<V>(standardPartitionOfUnity<V>(gfs, cc_bnd_neu_int_dir));
 
                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<Dune::PDELab::SubdomainBasis<V> > subdomain_basis;
                if (nev > 0)
                    subdomain_basis = std::make_shared<Dune::PDELab::GenEOBasis<GFS,J_ext,V,3> >(gfs, AF_neu, AF_ovlp, eigenvalue_threshold, *part_unity, solver_info.nev, solver_info.nev_arpack, 0.001, false, solver_info.verb);//,subdomain_basis2);
                else if (nev == 0){
                    subdomain_basis = std::make_shared<Dune::PDELab::ZEMBasis<GFS,LFS,V,3,3> >(gfs, lfs, *part_unity);
                }
                else{
                    subdomain_basis = std::make_shared<Dune::PDELab::SubdomainBasis<V> >(*part_unity);
                }
 
                auto partunityspace = std::make_shared<Dune::PDELab::SubdomainProjectedCoarseSpace<GFS,J_ext,V,PIH> >(gfs, AF_neu, subdomain_basis, pihf);
                typedef Dune::PDELab::ISTL::TwoLevelOverlappingAdditiveSchwarz<GFS,J,V,V> PREC_PCG;
                std::shared_ptr<PREC_PCG> prec;
                prec = std::make_shared<PREC_PCG>(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);
                V v(gfs,0.0);
 
                //Solve using CG
                std::shared_ptr<CGSolver<V>> solver;
                solver = std::make_shared<CGSolver<V> >(popf,ospf,*prec,solver_info.KrylovTol,solver_info.MaxIt,(helper.rank()==0) ? solver_info.verb : 0,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;
            const GFS_EXT& gfs_ext;
            LOP& lop;
            const CON& constraints;
            const MBE& mbe;
            Composites::SolverInfo& solver_info;
            Dune::MPIHelper& helper;
            int cells;
        };
    }
}
 
#endif