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>
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>
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>
77 recreateDirectSolver();
79 linsolver_->apply(x, rhs, res);
82 template <
class Operator>
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>
117 const std::function<VectorType()> weightsCalculator,
119 std::size_t pressureIndex)
122 linearoperator_for_solver_ = &op;
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;
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>>;
270using Comm = Dune::OwnerOverlapCopyCommunication<int, int>;
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>);
Dune::OverlappingSchwarzOperator< OBM< N >, BV< N >, BV< N >, Comm > ParOpM
Definition: FlexibleSolver_impl.hpp:272
Dune::BlockVector< Dune::FieldVector< double, N > > BV
Definition: FlexibleSolver_impl.hpp:257
Dune::OwnerOverlapCopyCommunication< int, int > Comm
Definition: FlexibleSolver_impl.hpp:270
Dune::BCRSMatrix< Opm::MatrixBlock< double, N, N > > OBM
Definition: FlexibleSolver_impl.hpp:259
Dune::MatrixAdapter< OBM< N >, BV< N >, BV< N > > SeqOpM
Definition: FlexibleSolver_impl.hpp:263
Definition: FlexibleSolver.hpp:45
virtual Dune::SolverCategory::Category category() const override
Definition: FlexibleSolver_impl.hpp:105
virtual void apply(VectorType &x, VectorType &rhs, Dune::InverseOperatorResult &res) override
Definition: FlexibleSolver_impl.hpp:74
typename Operator::domain_type VectorType
Definition: FlexibleSolver.hpp:47
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())
Definition: PreconditionerFactory_impl.hpp:702
Definition: PropertyTree.hpp:37
T get(const std::string &key) const
std::optional< PropertyTree > get_child_optional(const std::string &key) const
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:51
Definition: SupportsFaceTag.hpp:27
Dune::InverseOperatorResult InverseOperatorResult
Definition: BdaBridge.hpp:32