TpsaPreconditioner.hpp

Go to the documentation of this file.

// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=4 sw=4 sts=4:
/*
  Copyright 2025 NORCE 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_TPSA_PRECONDITIONER_HPP
#define OPM_TPSA_PRECONDITIONER_HPP
 
#include <dune/istl/operators.hh>
#include <dune/istl/owneroverlapcopy.hh>
#include <dune/istl/schwarz.hh>
#include <dune/istl/paamg/pinfo.hh>
 
#include <opm/simulators/linalg/FlexibleSolver.hpp>
#include <opm/simulators/linalg/PreconditionerWithUpdate.hpp>
#include <opm/simulators/linalg/PropertyTree.hpp>
#include <opm/simulators/linalg/tpsa/TpsaTypes.hpp>
 
#include <functional>
#include <memory>
#include <type_traits>
 
namespace Opm
{
 
template <typename Scalar>
using SeqDispDispOperatorT = Dune::MatrixAdapter<Linear::DispDispMatrix00T<Scalar>,
                                                 Linear::DispVector0T<Scalar>,
                                                 Linear::DispVector0T<Scalar> >;
template <typename Scalar>
using SeqRotRotOperatorT = Dune::MatrixAdapter<Linear::RotRotMatrixT<Scalar>,
                                               Linear::RotVectorT<Scalar>,
                                               Linear::RotVectorT<Scalar> >;
template <typename Scalar>
using SeqSPresSPresOperatorT = Dune::MatrixAdapter<Linear::SPresSPresMatrixT<Scalar>,
                                                   Linear::SPresVectorT<Scalar>,
                                                   Linear::SPresVectorT<Scalar> >;
 
#if HAVE_MPI
using TpsaParComm = Dune::OwnerOverlapCopyCommunication<int, int>;
 
template <typename Scalar>
using ParDispDispOperatorT = Dune::OverlappingSchwarzOperator<Linear::DispDispMatrix00T<Scalar>,
                                                              Linear::DispVector0T<Scalar>,
                                                              Linear::DispVector0T<Scalar>,
                                                              TpsaParComm>;
template <typename Scalar>
using ParRotRotOperatorT = Dune::OverlappingSchwarzOperator<Linear::RotRotMatrixT<Scalar>,
                                                            Linear::RotVectorT<Scalar>,
                                                            Linear::RotVectorT<Scalar>,
                                                            TpsaParComm>;
template <typename Scalar>
using ParSPresSPresOperatorT = Dune::OverlappingSchwarzOperator<Linear::SPresSPresMatrixT<Scalar>,
                                                                Linear::SPresVectorT<Scalar>,
                                                                Linear::SPresVectorT<Scalar>,
                                                                TpsaParComm>;
#endif
 
template <class Scalar, class DispOp, class RotOp, class SPresOp,
          class Comm = Dune::Amg::SequentialInformation>
class TpsaPreconditioner
    : public Dune::PreconditionerWithUpdate<Linear::TpsaMultiVector<Scalar>,
                                            Linear::TpsaMultiVector<Scalar> >
{
    using MultiVector = Linear::TpsaMultiVector<Scalar>;
    using DispSolver = Dune::FlexibleSolver<DispOp>;
    using RotSolver = Dune::FlexibleSolver<RotOp>;
    using SPresSolver = Dune::FlexibleSolver<SPresOp>;
 
public:
    static constexpr bool isParallel = !std::is_same_v<Comm, Dune::Amg::SequentialInformation>;
 
    static constexpr auto _0 = Dune::Indices::_0;
    static constexpr auto _1 = Dune::Indices::_1;
    static constexpr auto _2 = Dune::Indices::_2;
    static constexpr auto _3 = Dune::Indices::_3;
    static constexpr auto _4 = Dune::Indices::_4;
 
    static constexpr std::size_t pressureIdx = 0;
 
    TpsaPreconditioner(const Linear::TpsaMatrixView<Scalar>& S, const PropertyTree& prm)
        requires (!isParallel);
 
    TpsaPreconditioner(const Linear::TpsaMatrixView<Scalar>& S,
                       const PropertyTree& prm,
                       const Comm& comm)
        requires (isParallel);
 
    void pre(MultiVector&, MultiVector&) override;
 
    void post(MultiVector&) override;
 
    Dune::SolverCategory::Category category() const override;
 
    void update() override;
 
    bool hasPerfectUpdate() const override;
 
    void apply(MultiVector& v, const MultiVector& d) override;
 
private:
    void initSubSolvers_(const PropertyTree& prm);
 
    const Linear::TpsaMatrixView<Scalar>& S_;
    const Comm* comm_{nullptr};
 
    // Operators of the diagonal blocks, referenced by the solvers below.
    std::unique_ptr<DispOp> dispOp0_;
    std::unique_ptr<DispOp> dispOp1_;
    std::unique_ptr<DispOp> dispOp2_;
    std::unique_ptr<RotOp> rotOp_;
    std::unique_ptr<SPresOp> sPresOp_;
 
    // Solvers of the diagonal blocks.
    std::unique_ptr<DispSolver> dispSolver0_;
    std::unique_ptr<DispSolver> dispSolver1_;
    std::unique_ptr<DispSolver> dispSolver2_;
    std::unique_ptr<RotSolver> rotSolver_;
    std::unique_ptr<SPresSolver> sPresSolver_;
};
 
} // namespace Opm
 
#endif // OPM_TPSA_PRECONDITIONER_HPP