21 #ifndef OPM_FLEXIBLE_SOLVER_IMPL_HEADER_INCLUDED
22 #define OPM_FLEXIBLE_SOLVER_IMPL_HEADER_INCLUDED
24 #include <opm/simulators/linalg/matrixblock.hh>
25 #include <opm/simulators/linalg/ilufirstelement.hh>
26 #include <opm/simulators/linalg/FlexibleSolver.hpp>
27 #include <opm/simulators/linalg/PreconditionerFactory.hpp>
29 #include <dune/common/fmatrix.hh>
30 #include <dune/istl/bcrsmatrix.hh>
31 #include <dune/istl/solvers.hh>
32 #include <dune/istl/umfpack.hh>
33 #include <dune/istl/owneroverlapcopy.hh>
34 #include <dune/istl/paamg/pinfo.hh>
39 template <
class MatrixType,
class VectorType>
43 const std::function<VectorType()>& weightsCalculator,
44 std::size_t pressureIndex)
46 init(op, Dune::Amg::SequentialInformation(), prm, weightsCalculator,
51 template <
class MatrixType,
class VectorType>
57 const std::function<VectorType()>& weightsCalculator,
58 std::size_t pressureIndex)
60 init(op, comm, prm, weightsCalculator, pressureIndex);
63 template <
class MatrixType,
class VectorType>
66 apply(VectorType& x, VectorType& rhs, Dune::InverseOperatorResult& res)
68 linsolver_->apply(x, rhs, res);
71 template <
class MatrixType,
class VectorType>
73 FlexibleSolver<MatrixType, VectorType>::
74 apply(VectorType& x, VectorType& rhs,
double reduction, Dune::InverseOperatorResult& res)
76 linsolver_->apply(x, rhs, reduction, res);
80 template <
class MatrixType,
class VectorType>
85 return *preconditioner_;
88 template <
class MatrixType,
class VectorType>
89 Dune::SolverCategory::Category
93 return linearoperator_for_solver_->category();
97 template <
class MatrixType,
class VectorType>
100 FlexibleSolver<MatrixType, VectorType>::
101 initOpPrecSp(AbstractOperatorType& op,
103 const std::function<VectorType()> weightsCalculator,
105 std::size_t pressureIndex)
108 using ParOperatorType = Dune::OverlappingSchwarzOperator<MatrixType, VectorType, VectorType, Comm>;
109 linearoperator_for_solver_ = &op;
110 auto op_prec = std::make_shared<ParOperatorType>(op.getmat(), comm);
111 auto child = prm.get_child_optional(
"preconditioner");
117 scalarproduct_ = Dune::createScalarProduct<VectorType, Comm>(comm, op.category());
118 linearoperator_for_precond_ = op_prec;
121 template <
class MatrixType,
class VectorType>
123 FlexibleSolver<MatrixType, VectorType>::
124 initOpPrecSp(AbstractOperatorType& op,
126 const std::function<VectorType()> weightsCalculator,
127 const Dune::Amg::SequentialInformation&,
128 std::size_t pressureIndex)
131 using SeqOperatorType = Dune::MatrixAdapter<MatrixType, VectorType, VectorType>;
132 linearoperator_for_solver_ = &op;
133 auto op_prec = std::make_shared<SeqOperatorType>(op.getmat());
134 auto child = prm.get_child_optional(
"preconditioner");
139 scalarproduct_ = std::make_shared<Dune::SeqScalarProduct<VectorType>>();
140 linearoperator_for_precond_ = op_prec;
143 template <
class MatrixType,
class VectorType>
145 FlexibleSolver<MatrixType, VectorType>::
148 const double tol = prm.get<
double>(
"tol", 1e-2);
149 const int maxiter = prm.get<
int>(
"maxiter", 200);
150 const int verbosity = is_iorank ? prm.get<
int>(
"verbosity", 0) : 0;
151 const std::string solver_type = prm.get<std::string>(
"solver",
"bicgstab");
152 if (solver_type ==
"bicgstab") {
153 linsolver_.reset(
new Dune::BiCGSTABSolver<VectorType>(*linearoperator_for_solver_,
159 }
else if (solver_type ==
"loopsolver") {
160 linsolver_.reset(
new Dune::LoopSolver<VectorType>(*linearoperator_for_solver_,
166 }
else if (solver_type ==
"gmres") {
167 int restart = prm.get<
int>(
"restart", 15);
168 linsolver_.reset(
new Dune::RestartedGMResSolver<VectorType>(*linearoperator_for_solver_,
175 #if HAVE_SUITESPARSE_UMFPACK
176 }
else if (solver_type ==
"umfpack") {
178 linsolver_.reset(
new Dune::UMFPack<MatrixType>(linearoperator_for_solver_->getmat(), verbosity, dummy));
181 OPM_THROW(std::invalid_argument,
"Properties: Solver " << solver_type <<
" not known.");
188 template <
class MatrixType,
class VectorType>
189 template <
class Comm>
191 FlexibleSolver<MatrixType, VectorType>::
192 init(AbstractOperatorType& op,
195 const std::function<VectorType()> weightsCalculator,
196 std::size_t pressureIndex)
198 initOpPrecSp(op, prm, weightsCalculator, comm, pressureIndex);
199 initSolver(prm, comm.communicator().rank() == 0);
208 using BV = Dune::BlockVector<Dune::FieldVector<double, N>>;
210 using BM = Dune::BCRSMatrix<Dune::FieldMatrix<double, N, N>>;
212 using OBM = Dune::BCRSMatrix<Opm::MatrixBlock<double, N, N>>;
216 using Comm = Dune::OwnerOverlapCopyCommunication<int, int>;
222 #define INSTANTIATE_FLEXIBLESOLVER(N) \
223 template class Dune::FlexibleSolver<BM<N>, BV<N>>; \
224 template class Dune::FlexibleSolver<OBM<N>, BV<N>>; \
225 template Dune::FlexibleSolver<BM<N>, BV<N>>::FlexibleSolver(AbstractOperatorType& op, \
227 const Opm::PropertyTree& prm, \
228 const std::function<BV<N>()>& weightsCalculator, \
230 template Dune::FlexibleSolver<OBM<N>, BV<N>>::FlexibleSolver(AbstractOperatorType& op, \
232 const Opm::PropertyTree& prm, \
233 const std::function<BV<N>()>& weightsCalculator, \
238 #define INSTANTIATE_FLEXIBLESOLVER(N) \
239 template class Dune::FlexibleSolver<BM<N>, BV<N>>; \
240 template class Dune::FlexibleSolver<OBM<N>, BV<N>>;
A solver class that encapsulates all needed objects for a linear solver (operator,...
Definition: FlexibleSolver.hpp:40
Dune::AssembledLinearOperator< MatrixType, VectorType, VectorType > AbstractOperatorType
Base class type of the operator passed to the solver.
Definition: FlexibleSolver.hpp:46
FlexibleSolver(AbstractOperatorType &op, const Opm::PropertyTree &prm, const std::function< VectorType()> &weightsCalculator, std::size_t pressureIndex)
Create a sequential solver.
Definition: FlexibleSolver_impl.hpp:41
AbstractPrecondType & preconditioner()
Access the contained preconditioner.
Definition: FlexibleSolver_impl.hpp:83
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.hpp:71
Definition: PropertyTree.hpp:37