14#ifndef DUMUX_LINEAR_ISTL_SOLVERFACTORYBACKEND_HH
15#define DUMUX_LINEAR_ISTL_SOLVERFACTORYBACKEND_HH
19#include <dune/common/parallel/mpihelper.hh>
20#include <dune/common/parametertree.hh>
25#include <dune/istl/preconditioners.hh>
26#include <dune/istl/paamg/amg.hh>
30#include <dune/istl/solvers.hh>
31#include <dune/istl/solverfactory.hh>
54template<
class LinearOperator>
55int initSolverFactoriesForMultiTypeBlockMatrix()
57 using M =
typename LinearOperator::matrix_type;
58 using X =
typename LinearOperator::range_type;
59 using Y =
typename LinearOperator::domain_type;
60 using TL = Dune::TypeList<M,X,Y>;
61 auto& dsfac = Dune::DirectSolverFactory<M,X,Y>::instance();
62 Dune::addRegistryToFactory<TL>(dsfac, Dumux::MultiTypeBlockMatrixDirectSolverTag{});
63 auto& pfac = Dune::PreconditionerFactory<LinearOperator,X,Y>::instance();
64 Dune::addRegistryToFactory<TL>(pfac, Dumux::MultiTypeBlockMatrixPreconditionerTag{});
65 using TLS = Dune::TypeList<X,Y>;
66 auto& isfac = Dune::IterativeSolverFactory<X,Y>::instance();
67 return Dune::addRegistryToFactory<TLS>(isfac, Dune::IterativeSolverTag{});
78template<
class Matrix,
class LinearOperator>
82 initSolverFactoriesForMultiTypeBlockMatrix<LinearOperator>();
84 Dune::initSolverFactories<LinearOperator>();
93template <
class LinearSolverTraits,
class LinearAlgebraTraits>
98 using ScalarProduct = Dune::ScalarProduct<Vector>;
100 using Comm = Dune::OwnerOverlapCopyCommunication<Dune::bigunsignedint<96>,
int>;
111 , isParallel_(
Dune::MPIHelper::getCommunication().size() > 1)
114 DUNE_THROW(Dune::InvalidStateException,
"Using sequential constructor for parallel run. Use signature with gridView and dofMapper!");
117 initializeParameters_();
118 scalarProduct_ = std::make_shared<ScalarProduct>();
119 solverCategory_ = Dune::SolverCategory::sequential;
130 const typename LinearSolverTraits::DofMapper& dofMapper,
134 , isParallel_(gridView.comm().size() > 1)
138 initializeParameters_();
140 solverCategory_ = Detail::solverCategory<LinearSolverTraits>(gridView);
141 if (solverCategory_ != Dune::SolverCategory::sequential)
143 parallelHelper_ = std::make_unique<ParallelISTLHelper<LinearSolverTraits>>(gridView, dofMapper);
144 comm_ = std::make_shared<Comm>(gridView.comm(), solverCategory_);
145 scalarProduct_ = Dune::createScalarProduct<Vector>(*comm_, solverCategory_);
146 parallelHelper_->createParallelIndexSet(*comm_);
149 scalarProduct_ = std::make_shared<ScalarProduct>();
151 solverCategory_ = Dune::SolverCategory::sequential;
152 scalarProduct_ = std::make_shared<ScalarProduct>();
163 bool solve(Matrix& A, Vector& x, Vector& b)
166 solveSequentialOrParallel_(A, x, b);
168 solveSequential_(A, x, b);
171 return result_.converged;
174 const Dune::InverseOperatorResult&
result()
const
179 const std::string&
name()
const
189 if (solverCategory_ == Dune::SolverCategory::nonoverlapping)
192 using GV =
typename LinearSolverTraits::GridView;
193 using DM =
typename LinearSolverTraits::DofMapper;
196 return scalarProduct_->norm(y);
200 return scalarProduct_->norm(x);
205 void initializeParameters_()
208 checkMandatoryParameters_();
209 name_ = params_.get<std::string>(
"preconditioner.type") +
"-preconditioned " + params_.get<std::string>(
"type");
210 if (params_.get<
int>(
"verbose", 0) > 0)
211 std::cout <<
"Initialized linear solver of type: " << name_ << std::endl;
214 void checkMandatoryParameters_()
216 if (!params_.hasKey(
"type"))
217 DUNE_THROW(Dune::InvalidStateException,
"Solver factory needs \"LinearSolver.Type\" parameter to select the solver");
219 if (!params_.hasKey(
"preconditioner.type"))
220 DUNE_THROW(Dune::InvalidStateException,
"Solver factory needs \"LinearSolver.Preconditioner.Type\" parameter to select the preconditioner");
224 void solveSequentialOrParallel_(Matrix& A, Vector& x, Vector& b)
233 if (LinearSolverTraits::isNonOverlapping(parallelHelper_->gridView()))
235 using PTraits =
typename LinearSolverTraits::template ParallelNonoverlapping<Matrix, Vector>;
236 solveParallel_<PTraits>(A, x, b);
240 using PTraits =
typename LinearSolverTraits::template ParallelOverlapping<Matrix, Vector>;
241 solveParallel_<PTraits>(A, x, b);
245 solveSequential_(A, x, b);
249 solveSequential_(A, x, b);
253 template<
class ParallelTraits>
254 void solveParallel_(Matrix& A, Vector& x, Vector& b)
256 using LinearOperator =
typename ParallelTraits::LinearOperator;
259 initSolverFactories<Matrix, LinearOperator>();
261 prepareLinearAlgebraParallel<LinearSolverTraits, ParallelTraits>(A, b, *parallelHelper_);
262 auto linearOperator = std::make_shared<LinearOperator>(A, *comm_);
265 apply_(linearOperator, x, b);
269 void solveSequential_(Matrix& A, Vector& x, Vector& b)
272 using Traits =
typename LinearSolverTraits::template Sequential<Matrix, Vector>;
273 using LinearOperator =
typename Traits::LinearOperator;
274 auto linearOperator = std::make_shared<LinearOperator>(A);
277 initSolverFactories<Matrix, LinearOperator>();
280 apply_(linearOperator, x, b);
283 template<
class LinearOperator>
284 auto getSolverFromFactory_(std::shared_ptr<LinearOperator>& fop)
287 return Dune::getSolverFromFactory(fop, params_);
288 }
catch(Dune::Exception& e) {
289 std::cerr <<
"Dune::Exception during solver construction: " << e.what() << std::endl;
290 return std::decay_t<decltype(Dune::getSolverFromFactory(fop, params_))>();
294 template<
class LinearOperator>
295 void apply_(std::shared_ptr<LinearOperator>& fop, Vector& x, Vector& b)
297 auto solver = getSolverFromFactory_(fop);
301 bool success =
static_cast<bool>(solver);
303 int successRemote = success;
305 successRemote = comm_->communicator().min(success);
308 DUNE_THROW(Dune::Exception,
"Could not create ISTL solver");
309 else if (!successRemote)
310 DUNE_THROW(Dune::Exception,
"Could not create ISTL solver on remote process");
313 DUNE_THROW(Dune::Exception,
"Could not create ISTL solver");
317 solver->apply(x, b, result_);
320 const std::string paramGroup_;
322 std::unique_ptr<ParallelISTLHelper<LinearSolverTraits>> parallelHelper_;
323 std::shared_ptr<Comm> comm_;
325 bool isParallel_ =
false;
328 std::shared_ptr<ScalarProduct> scalarProduct_;
329 Dune::SolverCategory::Category solverCategory_;
330 Dune::InverseOperatorResult result_;
331 Dune::ParameterTree params_;
A linear solver using the dune-istl solver factory to choose the solver and preconditioner at runtime...
Definition: istlsolverfactorybackend.hh:95
const std::string & name() const
Definition: istlsolverfactorybackend.hh:179
IstlSolverFactoryBackend(const std::string ¶mGroup="")
Construct the backend for the sequential case only.
Definition: istlsolverfactorybackend.hh:109
IstlSolverFactoryBackend(const typename LinearSolverTraits::GridView &gridView, const typename LinearSolverTraits::DofMapper &dofMapper, const std::string ¶mGroup="")
Construct the backend for parallel or sequential runs.
Definition: istlsolverfactorybackend.hh:129
bool solve(Matrix &A, Vector &x, Vector &b)
Solve a linear system.
Definition: istlsolverfactorybackend.hh:163
Scalar norm(const Vector &x) const
Definition: istlsolverfactorybackend.hh:184
const Dune::InverseOperatorResult & result() const
Definition: istlsolverfactorybackend.hh:174
Base class for linear solvers.
Definition: solver.hh:27
const std::string & paramGroup() const
the parameter group for getting parameter from the parameter tree
Definition: solver.hh:78
double Scalar
Definition: solver.hh:31
static Dune::ParameterTree createParameterTree(const std::string ¶mGroup="")
Create a tree containing parameters required for the linear solvers and precondioners of the Dune IST...
Definition: linearsolverparameters.hh:47
Definition: parallelhelpers.hh:476
void makeNonOverlappingConsistent(Dune::BlockVector< Block, Alloc > &v) const
Make a vector consistent for non-overlapping domain decomposition methods.
Definition: parallelhelpers.hh:484
The specialized Dumux macro and tag for the ISTL registry to choose the solver and preconditioner at ...
Generates a parameter tree required for the linear solvers and precondioners of the Dune ISTL.
Type traits to be used with matrix types.
static constexpr bool canCommunicate
Definition: gridcapabilities.hh:51
void initSolverFactories()
Initialize the solver factories for regular matrices or MultiTypeBlockMatrices.
Definition: istlsolverfactorybackend.hh:79
Definition: common/pdesolver.hh:24
Provides a helper class for nonoverlapping decomposition.
Dumux preconditioners for iterative solvers.
Base class for linear solvers.
V Vector
Definition: linearalgebratraits.hh:53
M Matrix
Definition: linearalgebratraits.hh:52
Helper type to determine whether a given type is a Dune::MultiTypeBlockMatrix.
Definition: matrix.hh:37
Helper type to determine whether a given type is a Dune::MultiTypeBlockVector.
Definition: vector.hh:22
Type traits to be used with vector types.