CPRPreconditioner.hpp

Go to the documentation of this file.

/*
  Copyright 2014 SINTEF ICT, Applied Mathematics.
  Copyright 2014 IRIS AS.
  Copyright 2015 Dr. Blatt - HPC-Simulation-Software & Services
  Copyright 2015 NTNU
  Copyright 2015 Statoil AS

  This file is part of the Open Porous Media project (OPM).

  OPM is free software: you can redistribute it and/or modify
  it under the terms of the GNU General Public License as published by
  the Free Software Foundation, either version 3 of the License, or
  (at your option) any later version.

  OPM is distributed in the hope that it will be useful,
  but WITHOUT ANY WARRANTY; without even the implied warranty of
  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  GNU General Public License for more details.

  You should have received a copy of the GNU General Public License
  along with OPM.  If not, see <http://www.gnu.org/licenses/>.
*/

#ifndef OPM_CPRPRECONDITIONER_HEADER_INCLUDED
#define OPM_CPRPRECONDITIONER_HEADER_INCLUDED

#include <memory>
#include <type_traits>

#include <opm/common/utility/platform_dependent/disable_warnings.h>
#include <opm/core/utility/parameters/ParameterGroup.hpp>

#include <dune/istl/bvector.hh>
#include <dune/istl/bcrsmatrix.hh>
#include <dune/istl/operators.hh>
#include <dune/istl/io.hh>
#include <dune/istl/owneroverlapcopy.hh>
#include <dune/istl/preconditioners.hh>
#include <dune/istl/schwarz.hh>
#include <dune/istl/solvers.hh>
#include <dune/istl/paamg/amg.hh>
#include <dune/istl/paamg/kamg.hh>
#include <dune/istl/paamg/pinfo.hh>

#include <opm/common/utility/platform_dependent/reenable_warnings.h>

#include <opm/common/ErrorMacros.hpp>
#include <opm/common/Exceptions.hpp>

#include <opm/autodiff/AdditionalObjectDeleter.hpp>
namespace Opm
{
namespace
{
template<class M, class X, class Y, class P>
struct CPRSelector
{
    typedef Dune::Amg::SequentialInformation ParallelInformation;
    typedef Dune::MatrixAdapter<M, X, Y> Operator;
    typedef Dune::SeqILU0<M,X,X> EllipticPreconditioner;
    typedef std::unique_ptr<EllipticPreconditioner> EllipticPreconditionerPointer;

    typedef EllipticPreconditioner Smoother;
    typedef Dune::Amg::AMG<Operator, X, Smoother, ParallelInformation> AMG;

    static Operator* makeOperator(const M& m, const P&)
    {
        return new Operator(m);
    }
};

#if HAVE_MPI
template<class M, class X, class Y, class I1, class I2>
struct CPRSelector<M,X,Y,Dune::OwnerOverlapCopyCommunication<I1,I2> >
{
    typedef Dune::OwnerOverlapCopyCommunication<I1,I2> ParallelInformation;
    typedef Dune::OverlappingSchwarzOperator<M,X,X,ParallelInformation> Operator;
    typedef Dune::BlockPreconditioner<X, X, ParallelInformation, Dune::SeqILU0<M,X,X> >
    EllipticPreconditioner;
    typedef std::unique_ptr<EllipticPreconditioner,
                            AdditionalObjectDeleter<Dune::SeqILU0<M,X,X> > >
    EllipticPreconditionerPointer;

    typedef EllipticPreconditioner Smoother;
    typedef Dune::Amg::AMG<Operator, X, Smoother, ParallelInformation> AMG;

    static Operator* makeOperator(const M& m, const ParallelInformation& p)
    {
        return new Operator(m, p);
    }
};

template<class ILU, class I1, class I2>
AdditionalObjectDeleter<ILU>
createParallelDeleter(ILU& ilu, const Dune::OwnerOverlapCopyCommunication<I1,I2>& p)
    {
        (void) p;
        return AdditionalObjectDeleter<ILU>(ilu);
    }
#endif

template<class M, class X>
std::shared_ptr<Dune::SeqILU0<M,X,X> >
createILU0Ptr(const M& A, double relax, const Dune::Amg::SequentialInformation&)
{
    return std::shared_ptr<Dune::SeqILU0<M,X,X> >(new Dune::SeqILU0<M,X,X>( A, relax) );
}
template<class M, class X>
std::shared_ptr<Dune::SeqILUn<M,X,X> >
createILUnPtr(const M& A, int ilu_n, double relax, const Dune::Amg::SequentialInformation&)
{
    return std::shared_ptr<Dune::SeqILUn<M,X,X> >(new Dune::SeqILUn<M,X,X>( A, ilu_n, relax) );
}

#if HAVE_MPI
template<class ILU, class I1, class I2>
struct SelectParallelILUSharedPtr
{
    typedef std::shared_ptr<
        Dune::BlockPreconditioner<
            typename ILU::range_type,
            typename ILU::domain_type,
            Dune::OwnerOverlapCopyCommunication<I1,I2>,
            ILU
            >
        > type;
};

template<class M, class X, class I1, class I2>
typename SelectParallelILUSharedPtr<Dune::SeqILU0<M,X,X>, I1, I2>::type
createILU0Ptr(const M& A, double relax,
              const Dune::OwnerOverlapCopyCommunication<I1,I2>& comm)
{
    typedef Dune::BlockPreconditioner<
        X,
        X,
        Dune::OwnerOverlapCopyCommunication<I1,I2>,
        Dune::SeqILU0<M,X,X>
        > PointerType;
    Dune::SeqILU0<M,X,X>* ilu = nullptr;
    int ilu_setup_successful = 1;
    std::string message;
    try {
        ilu = new Dune::SeqILU0<M,X,X>(A, relax);
    }
    catch ( Dune::MatrixBlockError error )
    {
                message = error.what();
                std::cerr<<"Exception occured on process " <<
                    comm.communicator().rank() << " during " <<
                    "setup of ILU0 preconditioner with message: " <<
                    message<<std::endl;
                ilu_setup_successful = 0;
    }
    // Check whether there was a problem on some process
    if ( comm.communicator().min(ilu_setup_successful) == 0 )
    {
        if ( ilu ) // not null if constructor succeeded
        {
            // prevent memory leak
            delete ilu;
        }
        throw Dune::MatrixBlockError();
    }
    return typename SelectParallelILUSharedPtr<Dune::SeqILU0<M,X,X>, I1, I2>
        ::type ( new PointerType(*ilu, comm), createParallelDeleter(*ilu, comm));
}

template<class M, class X, class I1, class I2>
typename SelectParallelILUSharedPtr<Dune::SeqILUn<M,X,X>, I1, I2>::type
createILUnPtr(const M& A, int ilu_n, double relax,
              const Dune::OwnerOverlapCopyCommunication<I1,I2>& comm)
{
    typedef Dune::BlockPreconditioner<
        X,
        X,
        Dune::OwnerOverlapCopyCommunication<I1,I2>,
        Dune::SeqILUn<M,X,X>
        > PointerType;
    Dune::SeqILUn<M,X,X>* ilu = new Dune::SeqILUn<M,X,X>( A, ilu_n, relax);

    return typename SelectParallelILUSharedPtr<Dune::SeqILUn<M,X,X>, I1, I2>::type
        (new PointerType(*ilu, comm),createParallelDeleter(*ilu, comm));
}
#endif

template<class M, class X=typename M::range_type>
std::unique_ptr<Dune::SeqILU0<M,X,X> >
createEllipticPreconditionerPointer(const M& Ae, double relax,
                                    const Dune::Amg::SequentialInformation&)
{
    return std::unique_ptr<Dune::SeqILU0<M,X,X> >(new Dune::SeqILU0<M,X,X>(Ae, relax));
}

#if HAVE_MPI
template<class M, class X=typename M::range_type, class I1, class I2>
typename CPRSelector<M,X,X,Dune::OwnerOverlapCopyCommunication<I1,I2> >
::EllipticPreconditionerPointer
createEllipticPreconditionerPointer(const M& Ae, double relax,
                                    const Dune::OwnerOverlapCopyCommunication<I1,I2>& comm)
{
    typedef Dune::BlockPreconditioner<X, X,
                                      Dune::OwnerOverlapCopyCommunication<I1,I2>,
                                      Dune::SeqILU0<M,X,X> >
    ParallelPreconditioner;

    Dune::SeqILU0<M,X,X>* ilu=new Dune::SeqILU0<M,X,X>(Ae, relax);
    typedef typename CPRSelector<M,X,X,Dune::OwnerOverlapCopyCommunication<I1,I2> >
        ::EllipticPreconditionerPointer EllipticPreconditionerPointer;
    return EllipticPreconditionerPointer(new ParallelPreconditioner(*ilu, comm),
                                         createParallelDeleter(*ilu, comm));
}
#endif
} // end namespace

    struct CPRParameter
    {
        double cpr_relax_;
        double cpr_solver_tol_;
        int cpr_ilu_n_;
        int cpr_max_ell_iter_;
        bool cpr_use_amg_;
        bool cpr_use_bicgstab_;
        bool cpr_solver_verbose_;

        CPRParameter() { reset(); }

        CPRParameter( const parameter::ParameterGroup& param)
        {
            // reset values to default
            reset();

            cpr_relax_          = param.getDefault("cpr_relax", cpr_relax_);
            cpr_solver_tol_     = param.getDefault("cpr_solver_tol", cpr_solver_tol_);
            cpr_ilu_n_          = param.getDefault("cpr_ilu_n", cpr_ilu_n_);
            cpr_max_ell_iter_   = param.getDefault("cpr_max_elliptic_iter",cpr_max_ell_iter_);
            cpr_use_amg_        = param.getDefault("cpr_use_amg", cpr_use_amg_);
            cpr_use_bicgstab_   = param.getDefault("cpr_use_bicgstab", cpr_use_bicgstab_);
            cpr_solver_verbose_ = param.getDefault("cpr_solver_verbose", cpr_solver_verbose_);
        }

        void reset()
        {
            cpr_relax_          = 1.0;
            cpr_solver_tol_     = 1e-2;
            cpr_ilu_n_          = 0;
            cpr_max_ell_iter_   = 25;
            cpr_use_amg_        = true;
            cpr_use_bicgstab_   = true;
            cpr_solver_verbose_ = false;
        }
    };


    template<class M, class X, class Y,
             class P=Dune::Amg::SequentialInformation>
    class CPRPreconditioner : public Dune::Preconditioner<X,Y>
    {
        // prohibit copying for now
        CPRPreconditioner( const CPRPreconditioner& );

    public:
        typedef P ParallelInformation;
        typedef typename Dune::remove_const<M>::type matrix_type;
        typedef X domain_type;
        typedef Y range_type;
        typedef typename X::field_type field_type;

        // define the category
        enum {
            category = std::is_same<P,Dune::Amg::SequentialInformation>::value?
            Dune::SolverCategory::sequential:Dune::SolverCategory::overlapping
        };

        typedef typename CPRSelector<M,X,X,P>::Operator Operator;

        typedef Dune::Preconditioner<X,X> WholeSystemPreconditioner;

        typedef typename CPRSelector<M,X,X,P>::EllipticPreconditionerPointer
        EllipticPreconditionerPointer;

        typedef typename CPRSelector<M,X,X,P>::AMG  AMG;

        CPRPreconditioner (const CPRParameter& param, const M& A, const M& Ae,
                           const ParallelInformation& comm=ParallelInformation(),
                           const ParallelInformation& commAe=ParallelInformation())
            : param_( param ),
              A_(A),
              Ae_(Ae),
              de_( Ae_.N() ),
              ve_( Ae_.M() ),
              dmodified_( A_.N() ),
              opAe_(CPRSelector<M,X,Y,P>::makeOperator(Ae_, commAe)),
              precond_(), // ilu0 preconditioner for elliptic system
              amg_(),     // amg  preconditioner for elliptic system
              pre_(), // copy A will be made be the preconditioner
              vilu_( A_.N() ),
              comm_(comm),
              commAe_(commAe)
        {
            // create appropriate preconditioner for elliptic system
            createEllipticPreconditioner( param_.cpr_use_amg_, commAe_ );

            // create the preconditioner for the whole system.
            if( param_.cpr_ilu_n_ == 0 ) {
                pre_ = createILU0Ptr<M,X>( A_, param_.cpr_relax_, comm_ );
            }
            else {
                pre_ = createILUnPtr<M,X>( A_, param_.cpr_ilu_n_, param_.cpr_relax_, comm_ );
            }
        }

        virtual void pre (X& /*x*/, Y& /*b*/)
        {
        }

        virtual void apply (X& v, const Y& d)
        {
            // Extract part of d corresponding to elliptic part.
            // Note: Assumes that the elliptic part comes first.
            std::copy_n(d.begin(), de_.size(), de_.begin());

            // Solve elliptic part, extend solution to full.
            // reset result
            ve_ = 0;
            solveElliptic( ve_, de_ );

            //reset return value
            v = 0.0;
            // Again assuming that the elliptic part comes first.
            std::copy(ve_.begin(), ve_.end(), v.begin());

            // Subtract elliptic residual from initial residual.
            // dmodified = d - A * vfull
            dmodified_ = d;
            A_.mmv(v, dmodified_);
            // A is not parallel, do communication manually.
            comm_.copyOwnerToAll(dmodified_, dmodified_);

            // Apply Preconditioner for whole system (relax will be applied already)
            pre_->apply( vilu_, dmodified_);

            // don't apply relaxation if relax_ == 1
            if( std::abs( param_.cpr_relax_ - 1.0 ) < 1e-12 ) {
                v += vilu_;
            }
            else {
                v *= param_.cpr_relax_;
                v += vilu_;
            }
        }

        virtual void post (X& /*x*/)
        {
        }

     protected:
        void solveElliptic(Y& x, Y& de)
        {
            // Linear solver parameters
            const double tolerance = param_.cpr_solver_tol_;
            const int maxit        = param_.cpr_max_ell_iter_;
            const int verbosity    = ( param_.cpr_solver_verbose_ &&
                                       comm_.communicator().rank()==0 ) ? 1 : 0;

            // operator result containing iterations etc.
            Dune::InverseOperatorResult result;

            // the scalar product chooser
            typedef Dune::ScalarProductChooser<X,ParallelInformation,category>
                ScalarProductChooser;
            // the scalar product.
            std::unique_ptr<typename ScalarProductChooser::ScalarProduct>
                sp(ScalarProductChooser::construct(commAe_));

            if( amg_ )
            {
                // Solve system with AMG
                if( param_.cpr_use_bicgstab_ ) {
                    Dune::BiCGSTABSolver<X> linsolve(*opAe_, *sp, (*amg_), tolerance, maxit, verbosity);
                    linsolve.apply(x, de, result);
                }
                else {
                    Dune::CGSolver<X> linsolve(*opAe_, *sp, (*amg_), tolerance, maxit, verbosity);
                    linsolve.apply(x, de, result);
                }
            }
            else
            {
                assert( precond_ );
                // Solve system with ILU-0
                if( param_.cpr_use_bicgstab_ ) {
                    Dune::BiCGSTABSolver<X> linsolve(*opAe_, *sp, (*precond_), tolerance, maxit, verbosity);
                    linsolve.apply(x, de, result);
                }
                else {
                    Dune::CGSolver<X> linsolve(*opAe_, *sp, (*precond_), tolerance, maxit, verbosity);
                    linsolve.apply(x, de, result);
                }

            }

            if (!result.converged) {
                OPM_THROW(LinearSolverProblem, "CPRPreconditioner failed to solve elliptic subsystem.");
            }
        }

        const CPRParameter& param_;

        const matrix_type& A_;
        const matrix_type& Ae_;

        Y de_, ve_, dmodified_;

        std::unique_ptr<Operator> opAe_;

        EllipticPreconditionerPointer precond_;
        std::unique_ptr< AMG > amg_;

        std::shared_ptr< WholeSystemPreconditioner > pre_;

        Y vilu_;

        const P& comm_;
        const P& commAe_;
     protected:
        void createEllipticPreconditioner( const bool amg, const P& comm )
        {
            if( amg )
            {
              // The coupling metric used in the AMG
              typedef Dune::Amg::FirstDiagonal CouplingMetric;

              // The coupling criterion used in the AMG
              typedef Dune::Amg::SymmetricCriterion<M, CouplingMetric> CritBase;

              // The coarsening criterion used in the AMG
              typedef Dune::Amg::CoarsenCriterion<CritBase> Criterion;

              // TODO: revise choice of parameters
              int coarsenTarget=1200;
              Criterion criterion(15,coarsenTarget);
              criterion.setDebugLevel( 0 ); // no debug information, 1 for printing hierarchy information
              criterion.setDefaultValuesIsotropic(2);
              criterion.setNoPostSmoothSteps( 1 );
              criterion.setNoPreSmoothSteps( 1 );

              // for DUNE 2.2 we also need to pass the smoother args
              typedef typename AMG::Smoother Smoother;
              typedef typename Dune::Amg::SmootherTraits<Smoother>::Arguments  SmootherArgs;
              SmootherArgs  smootherArgs;
              smootherArgs.iterations = 1;
              smootherArgs.relaxationFactor = param_.cpr_relax_;

              amg_ = std::unique_ptr< AMG > (new AMG(*opAe_, criterion, smootherArgs, comm ));
            }
            else
            {
                precond_ = createEllipticPreconditionerPointer<M,X>( Ae_, param_.cpr_relax_, comm);
            }
       }
    };


} // namespace Opm

#endif // OPM_CPRPRECONDITIONER_HEADER_INCLUDED