Commit f3b17c29 authored by Bernd Flemisch's avatar Bernd Flemisch
Browse files

[linear] adapt the AMG backend to the new structure

parent 1dbe1fc3
......@@ -38,16 +38,14 @@
#include <dumux/linear/solver.hh>
#include <dumux/linear/parallelhelpers.hh>
#include <dumux/linear/solvercategory.hh>
namespace Dumux {
/*!
* \ingroup Linear
* \brief A linear solver based on the ISTL AMG preconditioner
* and the ISTL BiCGSTAB solver.
*/
namespace Detail {
template <class LinearSolverTraits>
class AMGBiCGSTABBackend : public LinearSolver
class OldAMGBiCGSTABBackend : public LinearSolver
{
public:
/*!
......@@ -55,7 +53,8 @@ public:
*
* \param paramGroup the parameter group for parameter lookup
*/
AMGBiCGSTABBackend(const std::string& paramGroup = "")
[[deprecated("Use new AMGBiCGSTABBackend<LinearSolverTraits, LinearAlgebraTraits> with 2nd template parameter.")]]
OldAMGBiCGSTABBackend(const std::string& paramGroup = "")
: LinearSolver(paramGroup)
, isParallel_(Dune::MPIHelper::getCommunication().size() > 1)
{
......@@ -72,17 +71,21 @@ public:
* \param dofMapper an index mapper for dof entities
* \param paramGroup the parameter group for parameter lookup
*/
AMGBiCGSTABBackend(const typename LinearSolverTraits::GridView& gridView,
const typename LinearSolverTraits::DofMapper& dofMapper,
const std::string& paramGroup = "")
[[deprecated("Use new AMGBiCGSTABBackend<LinearSolverTraits, LinearAlgebraTraits> with 2nd template parameter.")]]
OldAMGBiCGSTABBackend(const typename LinearSolverTraits::GridView& gridView,
const typename LinearSolverTraits::DofMapper& dofMapper,
const std::string& paramGroup = "")
: LinearSolver(paramGroup)
#if HAVE_MPI
, isParallel_(Dune::MPIHelper::getCommunication().size() > 1)
#endif
{
#if HAVE_MPI
if (isParallel_)
phelper_ = std::make_unique<ParallelISTLHelper<LinearSolverTraits>>(gridView, dofMapper);
if constexpr (LinearSolverTraits::canCommunicate)
{
if (isParallel_)
phelper_ = std::make_unique<ParallelISTLHelper<LinearSolverTraits>>(gridView, dofMapper);
}
#endif
checkAvailabilityOfDirectSolver_();
}
......@@ -121,6 +124,9 @@ public:
return result_;
}
template<class Vector>
Scalar norm(const Vector& x) const = delete;
private:
//! see https://gitlab.dune-project.org/core/dune-istl/-/issues/62
void checkAvailabilityOfDirectSolver_()
......@@ -230,6 +236,228 @@ 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;
};
template<int i>
struct AMGImplHelper
{
template<class ...T>
using type = NewAMGBiCGSTABBackend<T...>;
};
template<>
struct AMGImplHelper<1>
{
template<class T>
using type = OldAMGBiCGSTABBackend<T>;
};
} // end namespace Detail
/*!
* \ingroup Linear
* \brief A linear solver based on the ISTL AMG preconditioner
* and the ISTL BiCGSTAB solver.
*/
template<class ...T>
using AMGBiCGSTABBackend = typename Detail::AMGImplHelper<sizeof...(T)>::template type<T...>;
} // end namespace Dumux
#endif
......@@ -43,6 +43,7 @@
#include <dumux/linear/linearsolveracceptsmultitypematrix.hh>
#include <dumux/linear/matrixconverter.hh>
#include <dumux/linear/parallelhelpers.hh>
#include <dumux/linear/solvercategory.hh>
namespace Dumux::Detail {
......@@ -58,23 +59,6 @@ constexpr std::size_t preconditionerBlockLevel() noexcept
return isMultiTypeBlockMatrix<M>::value ? 2 : 1;
}
template<class LinearSolverTraits, class GridView>
Dune::SolverCategory::Category solverCategory(const GridView& gridView)
{
if constexpr (LinearSolverTraits::canCommunicate)
{
if (gridView.comm().size() <= 1)
return Dune::SolverCategory::sequential;
if (LinearSolverTraits::isNonOverlapping(gridView))
return Dune::SolverCategory::nonoverlapping;
else
return Dune::SolverCategory::overlapping;
}
else
return Dune::SolverCategory::sequential;
}
} // end namespace Dumux::Detail
namespace Dumux {
......
// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=4 sw=4 sts=4:
/*****************************************************************************
* See the file COPYING for full copying permissions. *
* *
* This program is free software: you can redistribute it and/or modify *
* it under the terms of the GNU General Public License as published by *
* the Free Software Foundation, either version 3 of the License, or *
* (at your option) any later version. *
* *
* This program is distributed in the hope that it will be useful, *
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
* GNU General Public License for more details. *
* *
* You should have received a copy of the GNU General Public License *
* along with this program. If not, see <http://www.gnu.org/licenses/>. *
*****************************************************************************/
#ifndef DUMUX_SOLVERCATEGORY_HH
#define DUMUX_SOLVERCATEGORY_HH
#include <dune/istl/solvers.hh>
namespace Dumux::Detail {
template<class LinearSolverTraits, class GridView>
Dune::SolverCategory::Category solverCategory(const GridView& gridView)
{
if constexpr (LinearSolverTraits::canCommunicate)
{
if (gridView.comm().size() <= 1)
return Dune::SolverCategory::sequential;
if (LinearSolverTraits::isNonOverlapping(gridView))
return Dune::SolverCategory::nonoverlapping;
else
return Dune::SolverCategory::overlapping;
}
else
return Dune::SolverCategory::sequential;
}
} // end namespace Dumux::Detail
#endif
......@@ -34,6 +34,7 @@
#include <dumux/linear/amgbackend.hh>
#include <dumux/linear/linearsolvertraits.hh>
#include <dumux/linear/algebratraits.hh>
#include <dumux/nonlinear/newtonsolver.hh>
#include <dumux/assembly/fvassembler.hh>
......@@ -111,7 +112,8 @@ int main(int argc, char** argv)
auto assembler = std::make_shared<Assembler>(problem, gridGeometry, gridVariables, timeLoop, xOld);
// the linear solver
using LinearSolver = AMGBiCGSTABBackend<LinearSolverTraits<GridGeometry>>;
using LinearSolver = AMGBiCGSTABBackend<LinearSolverTraits<GridGeometry>,
LinearAlgebraTraitsFromAssembler<Assembler>>;
auto linearSolver = std::make_shared<LinearSolver>(leafGridView, gridGeometry->dofMapper());
// the non-linear solver
......
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