Commit deebda56 authored by Timo Koch's avatar Timo Koch Committed by Bernd Flemisch
Browse files

[amg] Use istl solver class to implement amg backend

parent 9c9e9fd4
......@@ -39,6 +39,7 @@
#include <dumux/linear/solver.hh>
#include <dumux/linear/parallelhelpers.hh>
#include <dumux/linear/solvercategory.hh>
#include <dumux/linear/istlsolvers.hh>
namespace Dumux {
......@@ -236,209 +237,22 @@ private:
bool isParallel_ = false;
};
template <class LinearSolverTraits, class LinearAlgebraTraits>
class NewAMGBiCGSTABBackend : public LinearSolver
{
using Matrix = typename LinearAlgebraTraits::Matrix;
using Vector = typename LinearAlgebraTraits::Vector;
using ScalarProduct = Dune::ScalarProduct<Vector>;
#if HAVE_MPI
using Comm = Dune::OwnerOverlapCopyCommunication<Dune::bigunsignedint<96>, int>;
#endif
public:
/*!
* \brief Construct the backend for the sequential case only
*
* \param paramGroup the parameter group for parameter lookup
*/
NewAMGBiCGSTABBackend(const std::string& paramGroup = "")
: LinearSolver(paramGroup)
, isParallel_(Dune::MPIHelper::getCollectiveCommunication().size() > 1)
{
if (isParallel_)
DUNE_THROW(Dune::InvalidStateException, "Using sequential constructor for parallel run. Use signature with gridView and dofMapper!");
solverCategory_ = Dune::SolverCategory::sequential;
scalarProduct_ = std::make_shared<ScalarProduct>();
}
/*!
* \brief Construct the backend for parallel or sequential runs
*
* \param gridView the grid view on which we are performing the multi-grid
* \param dofMapper an index mapper for dof entities
* \param paramGroup the parameter group for parameter lookup
*/
NewAMGBiCGSTABBackend(const typename LinearSolverTraits::GridView& gridView,
const typename LinearSolverTraits::DofMapper& dofMapper,
const std::string& paramGroup = "")
: LinearSolver(paramGroup)
#if HAVE_MPI
, isParallel_(Dune::MPIHelper::getCollectiveCommunication().size() > 1)
#endif
{
#if HAVE_MPI
solverCategory_ = Detail::solverCategory<LinearSolverTraits>(gridView);
if (solverCategory_ != Dune::SolverCategory::sequential)
{
phelper_ = std::make_unique<ParallelISTLHelper<LinearSolverTraits>>(gridView, dofMapper);
comm_ = std::make_shared<Comm>(gridView.comm(), solverCategory_);
scalarProduct_ = Dune::createScalarProduct<Vector>(*comm_, solverCategory_);
phelper_->createParallelIndexSet(*comm_);
}
else
scalarProduct_ = std::make_shared<ScalarProduct>();
#else
solverCategory_ = Dune::SolverCategory::sequential;
#endif
}
/*!
* \brief Solve a linear system.
*
* \param A the matrix
* \param x the seeked solution vector, containing the initial solution upon entry
* \param b the right hand side vector
*/
bool solve(Matrix& A, Vector& x, Vector& b)
{
#if HAVE_MPI
solveSequentialOrParallel_(A, x, b);
#else
solveSequential_(A, x, b);
#endif
return result_.converged;
}
/*!
* \brief The name of the solver
*/
std::string name() const
{
return "AMG-preconditioned BiCGSTAB solver";
}
/*!
* \brief The result containing the convergence history.
*/
const Dune::InverseOperatorResult& result() const
{
return result_;
}
Scalar norm(const Vector& x) const
{
#if HAVE_MPI
if constexpr (LinearSolverTraits::canCommunicate)
{
if (solverCategory_ == Dune::SolverCategory::nonoverlapping)
{
auto y(x); // make a copy because the vector needs to be made consistent
using GV = typename LinearSolverTraits::GridView;
using DM = typename LinearSolverTraits::DofMapper;
ParallelVectorHelper<GV, DM, LinearSolverTraits::dofCodim> vectorHelper(phelper_->gridView(), phelper_->dofMapper());
vectorHelper.makeNonOverlappingConsistent(y);
return scalarProduct_->norm(y);
}
}
#endif
return scalarProduct_->norm(x);
}
private:
#if HAVE_MPI
void solveSequentialOrParallel_(Matrix& A, Vector& x, Vector& b)
{
if constexpr (LinearSolverTraits::canCommunicate)
{
if (isParallel_)
{
if (LinearSolverTraits::isNonOverlapping(phelper_->gridView()))
{
using PTraits = typename LinearSolverTraits::template ParallelNonoverlapping<Matrix, Vector>;
solveParallel_<PTraits>(A, x, b);
}
else
{
using PTraits = typename LinearSolverTraits::template ParallelOverlapping<Matrix, Vector>;
solveParallel_<PTraits>(A, x, b);
}
}
else
solveSequential_(A, x, b);
}
else
{
solveSequential_(A, x, b);
}
}
template<class ParallelTraits>
void solveParallel_(Matrix& A, Vector& x, Vector& b)
{
prepareLinearAlgebraParallel<LinearSolverTraits, ParallelTraits>(A, b, *phelper_);
auto linearOperator = std::make_shared<typename ParallelTraits::LinearOperator>(A, *comm_);
using SeqSmoother = Dune::SeqSSOR<Matrix, Vector, Vector>;
using Smoother = typename ParallelTraits::template Preconditioner<SeqSmoother>;
solveWithAmg_<Smoother>(A, x, b, linearOperator, comm_);
}
#endif // HAVE_MPI
void solveSequential_(Matrix& A, Vector& x, Vector& b)
{
using SequentialTraits = typename LinearSolverTraits::template Sequential<Matrix, Vector>;
auto linearOperator = std::make_shared<typename SequentialTraits::LinearOperator>(A);
auto comm = std::make_shared<Dune::Amg::SequentialInformation>();
using Smoother = Dune::SeqSSOR<Matrix, Vector, Vector>;
solveWithAmg_<Smoother>(A, x, b, linearOperator, comm);
}
template<class Smoother, class LinearOperator, class Comm>
void solveWithAmg_(Matrix& A, Vector& x, Vector& b,
std::shared_ptr<LinearOperator>& linearOperator,
std::shared_ptr<Comm>& comm)
{
using SmootherArgs = typename Dune::Amg::SmootherTraits<Smoother>::Arguments;
using Criterion = Dune::Amg::CoarsenCriterion<Dune::Amg::SymmetricCriterion<Matrix, Dune::Amg::FirstDiagonal>>;
//! \todo Check whether the default accumulation mode atOnceAccu is needed.
//! \todo make parameters changeable at runtime from input file / parameter tree
Dune::Amg::Parameters params(15, 2000, 1.2, 1.6, Dune::Amg::atOnceAccu);
params.setDefaultValuesIsotropic(LinearSolverTraits::GridView::dimension);
params.setDebugLevel(this->verbosity());
Criterion criterion(params);
SmootherArgs smootherArgs;
smootherArgs.iterations = 1;
smootherArgs.relaxationFactor = 1;
using Amg = Dune::Amg::AMG<LinearOperator, Vector, Smoother, Comm>;
auto amg = std::make_shared<Amg>(*linearOperator, criterion, smootherArgs, *comm);
Dune::BiCGSTABSolver<Vector> solver(*linearOperator, *scalarProduct_, *amg, this->residReduction(), this->maxIter(),
comm->communicator().rank() == 0 ? this->verbosity() : 0);
solver.apply(x, b, result_);
}
#if HAVE_MPI
std::unique_ptr<ParallelISTLHelper<LinearSolverTraits>> phelper_;
std::shared_ptr<Comm> comm_;
#endif
Dune::SolverCategory::Category solverCategory_;
std::shared_ptr<ScalarProduct> scalarProduct_;
Dune::InverseOperatorResult result_;
bool isParallel_ = false;
};
/*!
* \ingroup Linear
* \brief An AMG preconditioned BiCGSTAB solver using dune-istl
*/
template<class LSTraits, class LATraits>
using NewAMGBiCGSTABBackend =
IstlLinearSolver<LSTraits, LATraits,
Dune::BiCGSTABSolver<typename LATraits::SingleTypeVector>,
Detail::IstlAmgPreconditionerFactory
>;
template<int i>
struct AMGImplHelper
{
template<class ...T>
using type = NewAMGBiCGSTABBackend<T...>;
template<class LSTraits, class LATraits>
using type = NewAMGBiCGSTABBackend<LSTraits, LATraits>;
};
template<>
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment