21#ifndef OPM_FLEXIBLE_SOLVER_IMPL_HEADER_INCLUDED
22#define OPM_FLEXIBLE_SOLVER_IMPL_HEADER_INCLUDED
24#include <opm/common/ErrorMacros.hpp>
26#include <opm/simulators/linalg/matrixblock.hh>
27#include <opm/simulators/linalg/ilufirstelement.hh>
28#include <opm/simulators/linalg/FlexibleSolver.hpp>
29#include <opm/simulators/linalg/PreconditionerFactory.hpp>
30#include <opm/simulators/linalg/PropertyTree.hpp>
31#include <opm/simulators/linalg/WellOperators.hpp>
33#include <dune/common/fmatrix.hh>
34#include <dune/istl/bcrsmatrix.hh>
35#include <dune/istl/solvers.hh>
36#include <dune/istl/umfpack.hh>
37#include <dune/istl/owneroverlapcopy.hh>
38#include <dune/istl/paamg/pinfo.hh>
43 template <
class Operator>
47 const std::function<VectorType()>& weightsCalculator,
48 std::size_t pressureIndex)
50 init(op, Dune::Amg::SequentialInformation(), prm, weightsCalculator,
55 template <
class Operator>
61 const std::function<VectorType()>& weightsCalculator,
62 std::size_t pressureIndex)
64 init(op, comm, prm, weightsCalculator, pressureIndex);
67 template <
class Operator>
70 apply(VectorType& x, VectorType& rhs, Dune::InverseOperatorResult& res)
72 linsolver_->apply(x, rhs, res);
75 template <
class Operator>
77 FlexibleSolver<Operator>::
78 apply(VectorType& x, VectorType& rhs,
double reduction, Dune::InverseOperatorResult& res)
80 linsolver_->apply(x, rhs, reduction, res);
84 template <
class Operator>
89 return *preconditioner_;
92 template <
class Operator>
93 Dune::SolverCategory::Category
97 return linearoperator_for_solver_->category();
101 template <
class Operator>
102 template <
class Comm>
104 FlexibleSolver<Operator>::
105 initOpPrecSp(Operator& op,
107 const std::function<VectorType()> weightsCalculator,
109 std::size_t pressureIndex)
112 linearoperator_for_solver_ = &op;
113 auto child = prm.get_child_optional(
"preconditioner");
119 scalarproduct_ = Dune::createScalarProduct<VectorType, Comm>(comm, op.category());
122 template <
class Operator>
124 FlexibleSolver<Operator>::
125 initOpPrecSp(Operator& op,
127 const std::function<VectorType()> weightsCalculator,
128 const Dune::Amg::SequentialInformation&,
129 std::size_t pressureIndex)
132 linearoperator_for_solver_ = &op;
133 auto child = prm.get_child_optional(
"preconditioner");
138 scalarproduct_ = std::make_shared<Dune::SeqScalarProduct<VectorType>>();
141 template <
class Operator>
143 FlexibleSolver<Operator>::
146 const double tol = prm.get<
double>(
"tol", 1e-2);
147 const int maxiter = prm.get<
int>(
"maxiter", 200);
148 const int verbosity = is_iorank ? prm.get<
int>(
"verbosity", 0) : 0;
149 const std::string solver_type = prm.get<std::string>(
"solver",
"bicgstab");
150 if (solver_type ==
"bicgstab") {
151 linsolver_.reset(
new Dune::BiCGSTABSolver<VectorType>(*linearoperator_for_solver_,
157 }
else if (solver_type ==
"loopsolver") {
158 linsolver_.reset(
new Dune::LoopSolver<VectorType>(*linearoperator_for_solver_,
164 }
else if (solver_type ==
"gmres") {
165 int restart = prm.get<
int>(
"restart", 15);
166 linsolver_.reset(
new Dune::RestartedGMResSolver<VectorType>(*linearoperator_for_solver_,
173#if HAVE_SUITESPARSE_UMFPACK
174 }
else if (solver_type ==
"umfpack") {
176 using MatrixType = std::remove_const_t<std::remove_reference_t<
decltype(linearoperator_for_solver_->getmat())>>;
180 OPM_THROW(std::invalid_argument,
"Properties: Solver " << solver_type <<
" not known.");
187 template <
class Operator>
188 template <
class Comm>
190 FlexibleSolver<Operator>::
194 const std::function<VectorType()> weightsCalculator,
195 std::size_t pressureIndex)
197 initOpPrecSp(op, prm, weightsCalculator, comm, pressureIndex);
198 initSolver(prm, comm.communicator().rank() == 0);
208using BV = Dune::BlockVector<Dune::FieldVector<double, N>>;
210using OBM = Dune::BCRSMatrix<Opm::MatrixBlock<double, N, N>>;
214using SeqOpM = Dune::MatrixAdapter<OBM<N>, BV<N>, BV<N>>;
221using Comm = Dune::OwnerOverlapCopyCommunication<int, int>;
223using ParOpM = Dune::OverlappingSchwarzOperator<OBM<N>, BV<N>, BV<N>, Comm>;
231#define INSTANTIATE_FLEXIBLESOLVER_OP(Operator) \
232template class Dune::FlexibleSolver<Operator>; \
233template Dune::FlexibleSolver<Operator>::FlexibleSolver(Operator& op, \
235 const Opm::PropertyTree& prm, \
236 const std::function<typename Operator::domain_type()>& weightsCalculator, \
237 std::size_t pressureIndex);
238#define INSTANTIATE_FLEXIBLESOLVER(N) \
239INSTANTIATE_FLEXIBLESOLVER_OP(SeqOpM<N>); \
240INSTANTIATE_FLEXIBLESOLVER_OP(SeqOpW<N>); \
241INSTANTIATE_FLEXIBLESOLVER_OP(ParOpM<N>); \
242INSTANTIATE_FLEXIBLESOLVER_OP(ParOpW<N>);
246#define INSTANTIATE_FLEXIBLESOLVER_OP(Operator) \
247template class Dune::FlexibleSolver<Operator>;
249#define INSTANTIATE_FLEXIBLESOLVER(N) \
250INSTANTIATE_FLEXIBLESOLVER_OP(SeqOpM<N>); \
251INSTANTIATE_FLEXIBLESOLVER_OP(SeqOpW<N>);
A solver class that encapsulates all needed objects for a linear solver (operator,...
Definition: FlexibleSolver.hpp:45
FlexibleSolver(Operator &op, const Opm::PropertyTree &prm, const std::function< VectorType()> &weightsCalculator, std::size_t pressureIndex)
Create a sequential solver.
Definition: FlexibleSolver_impl.hpp:45
AbstractPrecondType & preconditioner()
Access the contained preconditioner.
Definition: FlexibleSolver_impl.hpp:87
Definition: MSWellHelpers.hpp:29
static PrecPtr create(const Operator &op, const PropertyTree &prm, const std::function< Vector()> &weightsCalculator={}, std::size_t pressureIndex=std::numeric_limits< std::size_t >::max())
Create a new serial preconditioner and return a pointer to it.
Definition: PreconditionerFactory_impl.hpp:500
Definition: PropertyTree.hpp:37
Adapter to combine a matrix and another linear operator into a combined linear operator.
Definition: WellOperators.hpp:209
Adapter to combine a matrix and another linear operator into a combined linear operator.
Definition: WellOperators.hpp:123