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>
42#include <opm/simulators/linalg/gpuistl_hip/SolverAdapter.hpp>
44#include <opm/simulators/linalg/gpuistl/SolverAdapter.hpp>
51 template <
class Operator>
55 const std::function<VectorType()>& weightsCalculator,
56 std::size_t pressureIndex)
58 init(op, Dune::Amg::SequentialInformation(), prm, weightsCalculator,
63 template <
class Operator>
69 const std::function<VectorType()>& weightsCalculator,
70 std::size_t pressureIndex)
72 init(op, comm, prm, weightsCalculator, pressureIndex);
75 template <
class Operator>
78 apply(VectorType& x, VectorType& rhs, Dune::InverseOperatorResult& res)
81 recreateDirectSolver();
83 linsolver_->apply(x, rhs, res);
86 template <
class Operator>
88 FlexibleSolver<Operator>::
89 apply(VectorType& x, VectorType& rhs,
double reduction, Dune::InverseOperatorResult& res)
92 recreateDirectSolver();
94 linsolver_->apply(x, rhs, reduction, res);
98 template <
class Operator>
103 return *preconditioner_;
106 template <
class Operator>
107 Dune::SolverCategory::Category
111 return linearoperator_for_solver_->category();
115 template <
class Operator>
116 template <
class Comm>
118 FlexibleSolver<Operator>::
119 initOpPrecSp(Operator& op,
121 const std::function<VectorType()> weightsCalculator,
123 std::size_t pressureIndex)
126 linearoperator_for_solver_ = &op;
127 auto child = prm.get_child_optional(
"preconditioner");
129 child ? *child :
Opm::PropertyTree(),
133 scalarproduct_ = Dune::createScalarProduct<VectorType, Comm>(comm, op.category());
136 template <
class Operator>
138 FlexibleSolver<Operator>::
139 initOpPrecSp(Operator& op,
141 const std::function<VectorType()> weightsCalculator,
142 const Dune::Amg::SequentialInformation&,
143 std::size_t pressureIndex)
146 linearoperator_for_solver_ = &op;
147 auto child = prm.get_child_optional(
"preconditioner");
149 child ? *child :
Opm::PropertyTree(),
152 scalarproduct_ = std::make_shared<Dune::SeqScalarProduct<VectorType>>();
155 template <
class Operator>
156 template <
class Comm>
158 FlexibleSolver<Operator>::
161 const bool is_iorank = comm.communicator().rank() == 0;
162 const double tol = prm.get<
double>(
"tol", 1e-2);
163 const int maxiter = prm.get<
int>(
"maxiter", 200);
164 const int verbosity = is_iorank ? prm.get<
int>(
"verbosity", 0) : 0;
165 const std::string solver_type = prm.get<std::string>(
"solver",
"bicgstab");
166 if (solver_type ==
"bicgstab") {
167 linsolver_ = std::make_shared<Dune::BiCGSTABSolver<VectorType>>(*linearoperator_for_solver_,
173 }
else if (solver_type ==
"loopsolver") {
174 linsolver_ = std::make_shared<Dune::LoopSolver<VectorType>>(*linearoperator_for_solver_,
180 }
else if (solver_type ==
"gmres") {
181 int restart = prm.get<
int>(
"restart", 15);
182 linsolver_ = std::make_shared<Dune::RestartedGMResSolver<VectorType>>(*linearoperator_for_solver_,
189 }
else if (solver_type ==
"flexgmres") {
190 int restart = prm.get<
int>(
"restart", 15);
191 linsolver_ = std::make_shared<Dune::RestartedFlexibleGMResSolver<VectorType>>(*linearoperator_for_solver_,
198#if HAVE_SUITESPARSE_UMFPACK
199 }
else if (solver_type ==
"umfpack") {
200 if constexpr (std::is_same_v<typename VectorType::field_type,float>) {
201 OPM_THROW(std::invalid_argument,
"UMFPack cannot be used with floats");
203 using MatrixType = std::remove_const_t<std::remove_reference_t<
decltype(linearoperator_for_solver_->getmat())>>;
204 linsolver_ = std::make_shared<Dune::UMFPack<MatrixType>>(linearoperator_for_solver_->getmat(), verbosity,
false);
205 direct_solver_ =
true;
209 }
else if (solver_type ==
"gpubicgstab") {
211 *linearoperator_for_solver_,
220 OPM_THROW(std::invalid_argument,
221 "Properties: Solver " + solver_type +
" not known.");
231 template <
class Operator>
233 FlexibleSolver<Operator>::
234 recreateDirectSolver()
236#if HAVE_SUITESPARSE_UMFPACK
237 if constexpr (std::is_same_v<typename VectorType::field_type, float>) {
238 OPM_THROW(std::invalid_argument,
"UMFPack cannot be used with floats");
240 using MatrixType = std::remove_const_t<std::remove_reference_t<
decltype(linearoperator_for_solver_->getmat())>>;
241 linsolver_ = std::make_shared<Dune::UMFPack<MatrixType>>(linearoperator_for_solver_->getmat(), 0,
false);
244 OPM_THROW(std::logic_error,
"Direct solver specified, but the FlexibleSolver class was not compiled with SuiteSparse support.");
251 template <
class Operator>
252 template <
class Comm>
254 FlexibleSolver<Operator>::
258 const std::function<VectorType()> weightsCalculator,
259 std::size_t pressureIndex)
261 initOpPrecSp(op, prm, weightsCalculator, comm, pressureIndex);
262 initSolver(prm, comm);
271template<
class Scalar,
int N>
272using BV = Dune::BlockVector<Dune::FieldVector<Scalar, N>>;
273template<
class Scalar,
int N>
274using OBM = Dune::BCRSMatrix<Opm::MatrixBlock<Scalar, N, N>>;
277template<
class Scalar,
int N>
278using SeqOpM = Dune::MatrixAdapter<OBM<Scalar,N>, BV<Scalar,N>, BV<Scalar,N>>;
279template<
class Scalar,
int N>
285using Comm = Dune::OwnerOverlapCopyCommunication<int, int>;
286template<
class Scalar,
int N>
288template<
class Scalar,
int N>
290template<
class Scalar,
int N>
291using ParOpD = Dune::OverlappingSchwarzOperator<OBM<Scalar,N>, BV<Scalar,N>, BV<Scalar,N>, Comm>;
297#define INSTANTIATE_FLEXIBLESOLVER_OP(...) \
298 template class Dune::FlexibleSolver<__VA_ARGS__>; \
299 template Dune::FlexibleSolver<__VA_ARGS__>:: \
300 FlexibleSolver(__VA_ARGS__& op, \
302 const Opm::PropertyTree& prm, \
303 const std::function<typename __VA_ARGS__::domain_type()>& weightsCalculator, \
304 std::size_t pressureIndex);
306#define INSTANTIATE_FLEXIBLESOLVER(T,N) \
307 INSTANTIATE_FLEXIBLESOLVER_OP(SeqOpM<T,N>); \
308 INSTANTIATE_FLEXIBLESOLVER_OP(SeqOpW<T,N>); \
309 INSTANTIATE_FLEXIBLESOLVER_OP(ParOpM<T,N>); \
310 INSTANTIATE_FLEXIBLESOLVER_OP(ParOpW<T,N>); \
311 INSTANTIATE_FLEXIBLESOLVER_OP(ParOpD<T,N>);
315#define INSTANTIATE_FLEXIBLESOLVER_OP(...) \
316 template class Dune::FlexibleSolver<__VA_ARGS__>;
318#define INSTANTIATE_FLEXIBLESOLVER(T,N) \
319 INSTANTIATE_FLEXIBLESOLVER_OP(SeqOpM<T,N>); \
320 INSTANTIATE_FLEXIBLESOLVER_OP(SeqOpW<T,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:53
AbstractPrecondType & preconditioner()
Access the contained preconditioner.
Definition FlexibleSolver_impl.hpp:101
Dune linear operator that assumes ghost rows are ordered after interior rows.
Definition WellOperators.hpp:376
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:739
Definition PropertyTree.hpp:37
Adapter to combine a matrix and another linear operator into a combined linear operator.
Definition WellOperators.hpp:273
Adapter to combine a matrix and another linear operator into a combined linear operator.
Definition WellOperators.hpp:179
Wraps a CUDA solver to work with CPU data.
Definition SolverAdapter.hpp:52
This file contains a set of helper functions used by VFPProd / VFPInj.
Definition blackoilboundaryratevector.hh:37