PressureSolverPolicy.hpp

Go to the documentation of this file.

/*
  Copyright 2019 SINTEF Digital, Mathematics and Cybernetics.
 
  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_PRESSURE_SOLVER_POLICY_HEADER_INCLUDED
#define OPM_PRESSURE_SOLVER_POLICY_HEADER_INCLUDED
 
#include <opm/simulators/linalg/PressureTransferPolicy.hpp>
#include <opm/simulators/linalg/PropertyTree.hpp>
 
#include <dune/istl/solver.hh>
#include <dune/istl/owneroverlapcopy.hh>
 
namespace Dune
{
namespace Amg
{
    template <class OperatorType, class Solver, class LevelTransferPolicy>
    class PressureSolverPolicy
    {
    public:
        using Operator = OperatorType;
        explicit PressureSolverPolicy(const Opm::PropertyTree& prm)
            : prm_(prm)
        {
        }
 
    private:
        using X = typename Operator::range_type;
        struct PressureInverseOperator : public Dune::InverseOperator<X, X>
        {
#if HAVE_MPI
            template <typename GlobalIndex, typename LocalIndex>
            PressureInverseOperator(Operator& op,
                                    const Opm::PropertyTree& prm,
                                    const Dune::OwnerOverlapCopyCommunication<GlobalIndex, LocalIndex>& comm)
                : linsolver_()
            {
                assert(op.category() == Dune::SolverCategory::overlapping);
                // Assuming that we do not use Cpr as Pressure solver and use hard
                // coded pressure index that might be wrong but should be unused.
                linsolver_ = std::make_unique<Solver>(op, comm, prm, std::function<X()>(),
                                                      /* pressureIndex = */ 1);
            }
#endif // HAVE_MPI
 
            PressureInverseOperator(Operator& op,
                                    const Opm::PropertyTree& prm,
                                    const SequentialInformation&)
                : linsolver_()
            {
                assert(op.category() != Dune::SolverCategory::overlapping);
                // Assuming that we do not use Cpr as Pressure solver and use hard
                // coded pressure index that might be wrong but should be unused.
                linsolver_ = std::make_unique<Solver>(op, prm, std::function<X()>(),
                                                      /* pressureIndex = */ 1);
            }
 
 
            Dune::SolverCategory::Category category() const override
            {
                return linsolver_->category();
            }
 
            void apply(X& x, X& b, double reduction, Dune::InverseOperatorResult& res) override
            {
                linsolver_->apply(x, b, reduction, res);
            }
 
            void apply(X& x, X& b, Dune::InverseOperatorResult& res) override
            {
                linsolver_->apply(x, b, res);
            }
 
            void updatePreconditioner()
            {
                linsolver_->preconditioner().update();
            }
 
        private:
            std::unique_ptr<Solver> linsolver_;
        };
 
    public:
        using CoarseLevelSolver = PressureInverseOperator;
 
        template <class LTP>
        void setCoarseOperator(LTP& transferPolicy)
        {
            coarseOperator_ = transferPolicy.getCoarseLevelOperator();
        }
        template <class LTP>
        CoarseLevelSolver* createCoarseLevelSolver(LTP& transferPolicy)
        {
            coarseOperator_ = transferPolicy.getCoarseLevelOperator();
            auto& tp = dynamic_cast<LevelTransferPolicy&>(transferPolicy); // TODO: make this unnecessary.
            PressureInverseOperator* inv
                = new PressureInverseOperator(*coarseOperator_, prm_, tp.getCoarseLevelCommunication());
            return inv;
        }
 
    private:
        std::shared_ptr<Operator> coarseOperator_;
        Opm::PropertyTree prm_;
    };
} // namespace Amg
} // namespace Dune
 
 
#endif // OPM_PRESSURE_SOLVER_POLICY_HEADER_INCLUDED