25#ifndef DUMUX_PARALLEL_AMGBACKEND_HH
26#define DUMUX_PARALLEL_AMGBACKEND_HH
30#include <dune/common/exceptions.hh>
31#include <dune/common/parallel/indexset.hh>
32#include <dune/common/parallel/mpicommunication.hh>
33#include <dune/common/parallel/mpihelper.hh>
34#include <dune/grid/common/capabilities.hh>
35#include <dune/istl/paamg/amg.hh>
36#include <dune/istl/paamg/pinfo.hh>
37#include <dune/istl/solvers.hh>
49template <
class LinearSolverTraits>
60 , isParallel_(
Dune::MPIHelper::getCollectiveCommunication().size() > 1)
63 DUNE_THROW(Dune::InvalidStateException,
"Using sequential constructor for parallel run. Use signature with gridView and dofMapper!");
65 checkAvailabilityOfDirectSolver_();
76 const typename LinearSolverTraits::DofMapper& dofMapper,
80 , isParallel_(
Dune::MPIHelper::getCollectiveCommunication().size() > 1)
85 phelper_ = std::make_unique<ParallelISTLHelper<LinearSolverTraits>>(gridView, dofMapper);
87 checkAvailabilityOfDirectSolver_();
97 template<
class Matrix,
class Vector>
98 bool solve(Matrix& A, Vector& x, Vector& b)
101 solveSequentialOrParallel_(A, x, b);
103 solveSequential_(A, x, b);
105 return result_.converged;
113 return "AMG-preconditioned BiCGSTAB solver";
119 const Dune::InverseOperatorResult&
result()
const
126 void checkAvailabilityOfDirectSolver_()
128#if !HAVE_SUPERLU && !HAVE_UMFPACK
129 std::cout <<
"\nAMGBiCGSTABBackend: No direct solver backend found. Using iterative solver as coarse grid solver.\n"
130 <<
"Note that dune-istl currently hard-codes a tolerance of 1e-2 for the iterative coarse grid solver.\n"
131 <<
"This may result in reduced accuracy or performance depending on your setup.\nConsider installing "
132 <<
"UMFPack (SuiteSparse) or SuperLU or apply the istl patch, see dumux/patches/README.md." << std::endl;
137 template<
class Matrix,
class Vector>
138 void solveSequentialOrParallel_(Matrix& A, Vector& x, Vector& b)
144 if (LinearSolverTraits::isNonOverlapping(phelper_->gridView()))
146 using PTraits =
typename LinearSolverTraits::template ParallelNonoverlapping<Matrix, Vector>;
147 solveParallel_<PTraits>(A, x, b);
151 using PTraits =
typename LinearSolverTraits::template ParallelOverlapping<Matrix, Vector>;
152 solveParallel_<PTraits>(A, x, b);
156 solveSequential_(A, x, b);
160 solveSequential_(A, x, b);
164 template<
class ParallelTraits,
class Matrix,
class Vector>
165 void solveParallel_(Matrix& A, Vector& x, Vector& b)
167 using Comm =
typename ParallelTraits::Comm;
168 using LinearOperator =
typename ParallelTraits::LinearOperator;
169 using ScalarProduct =
typename ParallelTraits::ScalarProduct;
171 std::shared_ptr<Comm> comm;
172 std::shared_ptr<LinearOperator> linearOperator;
173 std::shared_ptr<ScalarProduct> scalarProduct;
174 prepareLinearAlgebraParallel<LinearSolverTraits, ParallelTraits>(A, b, comm, linearOperator, scalarProduct, *phelper_);
176 using SeqSmoother = Dune::SeqSSOR<Matrix, Vector, Vector>;
177 using Smoother =
typename ParallelTraits::template Preconditioner<SeqSmoother>;
178 solveWithAmg_<Smoother>(A, x, b, linearOperator, comm, scalarProduct);
182 template<
class Matrix,
class Vector>
183 void solveSequential_(Matrix& A, Vector& x, Vector& b)
185 using Comm = Dune::Amg::SequentialInformation;
186 using Traits =
typename LinearSolverTraits::template Sequential<Matrix, Vector>;
187 using LinearOperator =
typename Traits::LinearOperator;
188 using ScalarProduct =
typename Traits::ScalarProduct;
190 auto comm = std::make_shared<Comm>();
191 auto linearOperator = std::make_shared<LinearOperator>(A);
192 auto scalarProduct = std::make_shared<ScalarProduct>();
194 using Smoother = Dune::SeqSSOR<Matrix, Vector, Vector>;
195 solveWithAmg_<Smoother>(A, x, b, linearOperator, comm, scalarProduct);
198 template<
class Smoother,
class Matrix,
class Vector,
class LinearOperator,
class Comm,
class ScalarProduct>
199 void solveWithAmg_(Matrix& A, Vector& x, Vector& b,
200 std::shared_ptr<LinearOperator>& linearOperator,
201 std::shared_ptr<Comm>& comm,
202 std::shared_ptr<ScalarProduct>& scalarProduct)
204 using SmootherArgs =
typename Dune::Amg::SmootherTraits<Smoother>::Arguments;
205 using Criterion = Dune::Amg::CoarsenCriterion<Dune::Amg::SymmetricCriterion<Matrix, Dune::Amg::FirstDiagonal>>;
209 Dune::Amg::Parameters params(15, 2000, 1.2, 1.6, Dune::Amg::atOnceAccu);
210 params.setDefaultValuesIsotropic(LinearSolverTraits::GridView::dimension);
212 Criterion criterion(params);
213 SmootherArgs smootherArgs;
214 smootherArgs.iterations = 1;
215 smootherArgs.relaxationFactor = 1;
217 using Amg = Dune::Amg::AMG<LinearOperator, Vector, Smoother, Comm>;
218 auto amg = std::make_shared<Amg>(*linearOperator, criterion, smootherArgs, *comm);
220 Dune::BiCGSTABSolver<Vector> solver(*linearOperator, *scalarProduct, *amg, this->
residReduction(), this->
maxIter(),
221 comm->communicator().rank() == 0 ? this->verbosity() : 0);
223 solver.apply(x, b, result_);
227 std::unique_ptr<ParallelISTLHelper<LinearSolverTraits>> phelper_;
229 Dune::InverseOperatorResult result_;
230 bool isParallel_ =
false;
Base class for linear solvers.
Provides a helper class for nonoverlapping decomposition.
static constexpr bool canCommunicate
Definition: gridcapabilities.hh:63
Definition: common/pdesolver.hh:36
A linear solver based on the ISTL AMG preconditioner and the ISTL BiCGSTAB solver.
Definition: amgbackend.hh:51
const Dune::InverseOperatorResult & result() const
The result containing the convergence history.
Definition: amgbackend.hh:119
std::string name() const
The name of the solver.
Definition: amgbackend.hh:111
AMGBiCGSTABBackend(const std::string ¶mGroup="")
Construct the backend for the sequential case only.
Definition: amgbackend.hh:58
AMGBiCGSTABBackend(const typename LinearSolverTraits::GridView &gridView, const typename LinearSolverTraits::DofMapper &dofMapper, const std::string ¶mGroup="")
Construct the backend for parallel or sequential runs.
Definition: amgbackend.hh:75
bool solve(Matrix &A, Vector &x, Vector &b)
Solve a linear system.
Definition: amgbackend.hh:98
Base class for linear solvers.
Definition: solver.hh:37
double residReduction() const
the linear solver residual reduction
Definition: solver.hh:99
int maxIter() const
the maximum number of linear solver iterations
Definition: solver.hh:91
const std::string & paramGroup() const
the parameter group for getting parameter from the parameter tree
Definition: solver.hh:79
int verbosity() const
the verbosity level
Definition: solver.hh:83