Go to the documentation of this file.
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>
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>
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>
81 recreateDirectSolver();
83 linsolver_->apply(x, rhs, res);
86 template < class Operator>
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>
121 const std::function<VectorType()> weightsCalculator,
123 std::size_t pressureIndex)
126 linearoperator_for_solver_ = &op;
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;
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>
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>
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>);
Dune::OverlappingSchwarzOperator< OBM< Scalar, N >, BV< Scalar, N >, BV< Scalar, N >, Comm > ParOpD Definition: FlexibleSolver_impl.hpp:291
Dune::BCRSMatrix< Opm::MatrixBlock< Scalar, N, N > > OBM Definition: FlexibleSolver_impl.hpp:274
Dune::BlockVector< Dune::FieldVector< Scalar, N > > BV Definition: FlexibleSolver_impl.hpp:272
Dune::OwnerOverlapCopyCommunication< int, int > Comm Definition: FlexibleSolver_impl.hpp:285
Dune::MatrixAdapter< OBM< Scalar, N >, BV< Scalar, N >, BV< Scalar, N > > SeqOpM Definition: FlexibleSolver_impl.hpp:278
Definition: FlexibleSolver.hpp:45
virtual Dune::SolverCategory::Category category() const override Definition: FlexibleSolver_impl.hpp:109
virtual void apply(VectorType &x, VectorType &rhs, Dune::InverseOperatorResult &res) override Definition: FlexibleSolver_impl.hpp:78
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:53
AbstractPrecondType & preconditioner() Access the contained preconditioner. Definition: FlexibleSolver_impl.hpp:101
Dune linear operator that assumes ghost rows are ordered after interior rows. Avoids some computation... Definition: WellOperators.hpp:340
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:735
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:237
Adapter to combine a matrix and another linear operator into a combined linear operator. Definition: WellOperators.hpp:143
Wraps a CUDA solver to work with CPU data. Definition: SolverAdapter.hpp:52
Definition: fvbaseprimaryvariables.hh:141
Dune::InverseOperatorResult InverseOperatorResult Definition: BdaBridge.hpp:32
|