21#ifndef OPM_FLEXIBLE_SOLVER_IMPL_HEADER_INCLUDED
22#define OPM_FLEXIBLE_SOLVER_IMPL_HEADER_INCLUDED
24#include <opm/common/ErrorMacros.hpp>
25#include <opm/common/TimingMacros.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>
41#include <opm/simulators/linalg/cuistl/SolverAdapter.hpp>
47 template <
class Operator>
51 const std::function<VectorType()>& weightsCalculator,
52 std::size_t pressureIndex)
54 init(op, Dune::Amg::SequentialInformation(), prm, weightsCalculator,
59 template <
class Operator>
65 const std::function<VectorType()>& weightsCalculator,
66 std::size_t pressureIndex)
68 init(op, comm, prm, weightsCalculator, pressureIndex);
71 template <
class Operator>
74 apply(VectorType& x, VectorType& rhs, Dune::InverseOperatorResult& res)
77 recreateDirectSolver();
79 linsolver_->apply(x, rhs, res);
82 template <
class Operator>
84 FlexibleSolver<Operator>::
85 apply(VectorType& x, VectorType& rhs,
double reduction, Dune::InverseOperatorResult& res)
88 recreateDirectSolver();
90 linsolver_->apply(x, rhs, reduction, res);
94 template <
class Operator>
99 return *preconditioner_;
102 template <
class Operator>
103 Dune::SolverCategory::Category
107 return linearoperator_for_solver_->category();
111 template <
class Operator>
112 template <
class Comm>
114 FlexibleSolver<Operator>::
115 initOpPrecSp(Operator& op,
117 const std::function<VectorType()> weightsCalculator,
119 std::size_t pressureIndex)
122 linearoperator_for_solver_ = &op;
123 auto child = prm.get_child_optional(
"preconditioner");
125 child ? *child :
Opm::PropertyTree(),
129 scalarproduct_ = Dune::createScalarProduct<VectorType, Comm>(comm, op.category());
132 template <
class Operator>
134 FlexibleSolver<Operator>::
135 initOpPrecSp(Operator& op,
137 const std::function<VectorType()> weightsCalculator,
138 const Dune::Amg::SequentialInformation&,
139 std::size_t pressureIndex)
142 linearoperator_for_solver_ = &op;
143 auto child = prm.get_child_optional(
"preconditioner");
145 child ? *child :
Opm::PropertyTree(),
148 scalarproduct_ = std::make_shared<Dune::SeqScalarProduct<VectorType>>();
151 template <
class Operator>
153 FlexibleSolver<Operator>::
156 const double tol = prm.get<
double>(
"tol", 1e-2);
157 const int maxiter = prm.get<
int>(
"maxiter", 200);
158 const int verbosity = is_iorank ? prm.get<
int>(
"verbosity", 0) : 0;
159 const std::string solver_type = prm.get<std::string>(
"solver",
"bicgstab");
160 if (solver_type ==
"bicgstab") {
161 linsolver_ = std::make_shared<Dune::BiCGSTABSolver<VectorType>>(*linearoperator_for_solver_,
167 }
else if (solver_type ==
"loopsolver") {
168 linsolver_ = std::make_shared<Dune::LoopSolver<VectorType>>(*linearoperator_for_solver_,
174 }
else if (solver_type ==
"gmres") {
175 int restart = prm.get<
int>(
"restart", 15);
176 linsolver_ = std::make_shared<Dune::RestartedGMResSolver<VectorType>>(*linearoperator_for_solver_,
183 }
else if (solver_type ==
"flexgmres") {
184 int restart = prm.get<
int>(
"restart", 15);
185 linsolver_ = std::make_shared<Dune::RestartedFlexibleGMResSolver<VectorType>>(*linearoperator_for_solver_,
192#if HAVE_SUITESPARSE_UMFPACK
193 }
else if (solver_type ==
"umfpack") {
194 using MatrixType = std::remove_const_t<std::remove_reference_t<
decltype(linearoperator_for_solver_->getmat())>>;
195 linsolver_ = std::make_shared<Dune::UMFPack<MatrixType>>(linearoperator_for_solver_->getmat(), verbosity,
false);
196 direct_solver_ =
true;
199 }
else if (solver_type ==
"cubicgstab") {
201 *linearoperator_for_solver_,
209 OPM_THROW(std::invalid_argument,
210 "Properties: Solver " + solver_type +
" not known.");
220 template <
class Operator>
222 FlexibleSolver<Operator>::
223 recreateDirectSolver()
225#if HAVE_SUITESPARSE_UMFPACK
226 using MatrixType = std::remove_const_t<std::remove_reference_t<
decltype(linearoperator_for_solver_->getmat())>>;
227 linsolver_ = std::make_shared<Dune::UMFPack<MatrixType>>(linearoperator_for_solver_->getmat(), 0,
false);
229 OPM_THROW(std::logic_error,
"Direct solver specified, but the FlexibleSolver class was not compiled with SuiteSparse support.");
236 template <
class Operator>
237 template <
class Comm>
239 FlexibleSolver<Operator>::
243 const std::function<VectorType()> weightsCalculator,
244 std::size_t pressureIndex)
246 initOpPrecSp(op, prm, weightsCalculator, comm, pressureIndex);
247 initSolver(prm, comm.communicator().rank() == 0);
257using BV = Dune::BlockVector<Dune::FieldVector<double, N>>;
259using OBM = Dune::BCRSMatrix<Opm::MatrixBlock<double, N, N>>;
263using SeqOpM = Dune::MatrixAdapter<OBM<N>, BV<N>, BV<N>>;
270using Comm = Dune::OwnerOverlapCopyCommunication<int, int>;
272using ParOpM = Dune::OverlappingSchwarzOperator<OBM<N>, BV<N>, BV<N>, Comm>;
280#define INSTANTIATE_FLEXIBLESOLVER_OP(Operator) \
281template class Dune::FlexibleSolver<Operator>; \
282template Dune::FlexibleSolver<Operator>::FlexibleSolver(Operator& op, \
284 const Opm::PropertyTree& prm, \
285 const std::function<typename Operator::domain_type()>& weightsCalculator, \
286 std::size_t pressureIndex);
287#define INSTANTIATE_FLEXIBLESOLVER(N) \
288INSTANTIATE_FLEXIBLESOLVER_OP(SeqOpM<N>); \
289INSTANTIATE_FLEXIBLESOLVER_OP(SeqOpW<N>); \
290INSTANTIATE_FLEXIBLESOLVER_OP(ParOpM<N>); \
291INSTANTIATE_FLEXIBLESOLVER_OP(ParOpW<N>);
295#define INSTANTIATE_FLEXIBLESOLVER_OP(Operator) \
296template class Dune::FlexibleSolver<Operator>;
298#define INSTANTIATE_FLEXIBLESOLVER(N) \
299INSTANTIATE_FLEXIBLESOLVER_OP(SeqOpM<N>); \
300INSTANTIATE_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:49
AbstractPrecondType & preconditioner()
Access the contained preconditioner.
Definition FlexibleSolver_impl.hpp:97
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:702
Definition PropertyTree.hpp:37
Adapter to combine a matrix and another linear operator into a combined linear operator.
Definition WellOperators.hpp:224
Adapter to combine a matrix and another linear operator into a combined linear operator.
Definition WellOperators.hpp:134
Wraps a CUDA solver to work with CPU data.
Definition SolverAdapter.hpp:47
This file contains a set of helper functions used by VFPProd / VFPInj.
Definition BlackoilPhases.hpp:27