Commit 1f603347 authored by Bernd Flemisch's avatar Bernd Flemisch
Browse files

[linear] adapt ISTL solver factory backend to new structure

parent f3b17c29
......@@ -47,6 +47,7 @@
#include <dumux/linear/solver.hh>
#include <dumux/linear/parallelhelpers.hh>
#include <dumux/linear/istlsolverregistry.hh>
#include <dumux/linear/solvercategory.hh>
namespace Dumux {
......@@ -112,10 +113,11 @@ void initSolverFactories()
* \brief A linear solver using the dune-istl solver factory
* to choose the solver and preconditioner at runtime.
* \note the solvers are configured via the input file
* \note requires Dune version 2.7.1 or newer and 2.8 for parallel solvers
*/
namespace Detail {
template <class LinearSolverTraits>
class IstlSolverFactoryBackend : public LinearSolver
class OldIstlSolverFactoryBackend : public LinearSolver
{
public:
......@@ -124,7 +126,8 @@ public:
*
* \param paramGroup the parameter group for parameter lookup
*/
IstlSolverFactoryBackend(const std::string& paramGroup = "")
[[deprecated("Use new IstlSolverFactoryBackend<LinearSolverTraits, LinearAlgebraTraits> with 2nd template parameter.")]]
OldIstlSolverFactoryBackend(const std::string& paramGroup = "")
: paramGroup_(paramGroup)
, isParallel_(Dune::MPIHelper::getCommunication().size() > 1)
{
......@@ -142,9 +145,10 @@ public:
* \param dofMapper an index mapper for dof entities
* \param paramGroup the parameter group for parameter lookup
*/
IstlSolverFactoryBackend(const typename LinearSolverTraits::GridView& gridView,
const typename LinearSolverTraits::DofMapper& dofMapper,
const std::string& paramGroup = "")
[[deprecated("Use new IstlSolverFactoryBackend<LinearSolverTraits, LinearAlgebraTraits> with 2nd template parameter.")]]
OldIstlSolverFactoryBackend(const typename LinearSolverTraits::GridView& gridView,
const typename LinearSolverTraits::DofMapper& dofMapper,
const std::string& paramGroup = "")
: paramGroup_(paramGroup)
#if HAVE_MPI
, isParallel_(Dune::MPIHelper::getCommunication().size() > 1)
......@@ -307,6 +311,258 @@ private:
std::string name_;
};
template <class LinearSolverTraits, class LinearAlgebraTraits>
class NewIstlSolverFactoryBackend : 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
*/
NewIstlSolverFactoryBackend(const std::string& paramGroup = "")
: paramGroup_(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!");
firstCall_ = true;
initializeParameters_();
scalarProduct_ = std::make_shared<ScalarProduct>();
solverCategory_ = Dune::SolverCategory::sequential;
}
/*!
* \brief Construct the backend for parallel or sequential runs
*
* \param gridView the grid view for parallel communication via the grid
* \param dofMapper an index mapper for dof entities
* \param paramGroup the parameter group for parameter lookup
*/
NewIstlSolverFactoryBackend(const typename LinearSolverTraits::GridView& gridView,
const typename LinearSolverTraits::DofMapper& dofMapper,
const std::string& paramGroup = "")
: paramGroup_(paramGroup)
#if HAVE_MPI
, isParallel_(Dune::MPIHelper::getCollectiveCommunication().size() > 1)
#endif
{
firstCall_ = true;
initializeParameters_();
#if HAVE_MPI
solverCategory_ = Detail::solverCategory<LinearSolverTraits>(gridView);
if (solverCategory_ != Dune::SolverCategory::sequential)
{
parallelHelper_ = std::make_unique<ParallelISTLHelper<LinearSolverTraits>>(gridView, dofMapper);
comm_ = std::make_shared<Comm>(gridView.comm(), solverCategory_);
scalarProduct_ = Dune::createScalarProduct<Vector>(*comm_, solverCategory_);
parallelHelper_->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
firstCall_ = false;
return result_.converged;
}
const Dune::InverseOperatorResult& result() const
{
return result_;
}
const std::string& name() const
{
return name_;
}
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(parallelHelper_->gridView(), parallelHelper_->dofMapper());
vectorHelper.makeNonOverlappingConsistent(y);
return scalarProduct_->norm(y);
}
}
#endif
return scalarProduct_->norm(x);
}
private:
void initializeParameters_()
{
params_ = LinearSolverParameters<LinearSolverTraits>::createParameterTree(paramGroup_);
checkMandatoryParameters_();
name_ = params_.get<std::string>("preconditioner.type") + "-preconditioned " + params_.get<std::string>("type");
if (params_.get<int>("verbose", 0) > 0)
std::cout << "Initialized linear solver of type: " << name_ << std::endl;
}
void checkMandatoryParameters_()
{
if (!params_.hasKey("type"))
DUNE_THROW(Dune::InvalidStateException, "Solver factory needs \"LinearSolver.Type\" parameter to select the solver");
if (!params_.hasKey("preconditioner.type"))
DUNE_THROW(Dune::InvalidStateException, "Solver factory needs \"LinearSolver.Preconditioner.Type\" parameter to select the preconditioner");
}
#if HAVE_MPI
void solveSequentialOrParallel_(Matrix& A, Vector& x, Vector& b)
{
// Dune::MultiTypeBlockMatrix does not provide a ColIterator which is needed by Dune::NonoverlappingSchwarzOperator.
// We therefore can only solve these types of systems sequentially.
// TODO: This can be adapted once the situation in Dune ISTL changes.
if constexpr (LinearSolverTraits::canCommunicate && !isMultiTypeBlockMatrix<Matrix>::value)
{
if (isParallel_)
{
if (LinearSolverTraits::isNonOverlapping(parallelHelper_->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)
{
#if DUNE_VERSION_GT_REV(DUNE_ISTL,2,7,0)
using LinearOperator = typename ParallelTraits::LinearOperator;
if (firstCall_)
initSolverFactories<Matrix, LinearOperator>();
prepareLinearAlgebraParallel<LinearSolverTraits, ParallelTraits>(A, b, *parallelHelper_);
auto linearOperator = std::make_shared<LinearOperator>(A, *comm_);
// construct solver
auto solver = getSolverFromFactory_(linearOperator);
// solve linear system
solver->apply(x, b, result_);
#else
DUNE_THROW(Dune::NotImplemented, "Parallel solvers only available for dune-istl > 2.7.0");
#endif
}
#endif // HAVE_MPI
void solveSequential_(Matrix& A, Vector& x, Vector& b)
{
// construct linear operator
using Traits = typename LinearSolverTraits::template Sequential<Matrix, Vector>;
using LinearOperator = typename Traits::LinearOperator;
auto linearOperator = std::make_shared<LinearOperator>(A);
if (firstCall_)
initSolverFactories<Matrix, LinearOperator>();
// construct solver
auto solver = getSolverFromFactory_(linearOperator);
// solve linear system
solver->apply(x, b, result_);
}
template<class LinearOperator>
auto getSolverFromFactory_(std::shared_ptr<LinearOperator>& fop)
{
try { return Dune::getSolverFromFactory(fop, params_); }
catch(Dune::Exception& e)
{
std::cerr << "Could not create solver with factory" << std::endl;
std::cerr << e.what() << std::endl;
throw e;
}
}
const std::string paramGroup_;
#if HAVE_MPI
std::unique_ptr<ParallelISTLHelper<LinearSolverTraits>> parallelHelper_;
std::shared_ptr<Comm> comm_;
#endif
bool isParallel_ = false;
bool firstCall_;
std::shared_ptr<ScalarProduct> scalarProduct_;
Dune::SolverCategory::Category solverCategory_;
Dune::InverseOperatorResult result_;
Dune::ParameterTree params_;
std::string name_;
};
template<int i>
struct BackendImplHelper
{
template<class ...T>
using type = NewIstlSolverFactoryBackend<T...>;
};
template<>
struct BackendImplHelper<1>
{
template<class T>
using type = OldIstlSolverFactoryBackend<T>;
};
} // end namespace Detail
/*!
* \ingroup Linear
* \brief A linear solver using the dune-istl solver factory
* to choose the solver and preconditioner at runtime.
* \note the solvers are configured via the input file
* \note requires Dune version 2.7.1 or newer
*/
template<class ...T>
using IstlSolverFactoryBackend = typename Detail::BackendImplHelper<sizeof...(T)>::template type<T...>;
} // end namespace Dumux
#endif
......@@ -38,7 +38,7 @@ dumux_add_test(NAME test_ff_navierstokes_sincos_uzawapreconditioner_factory
LABELS freeflow
TIMEOUT 5000
CMAKE_GUARD "( HAVE_UMFPACK AND ( ( DUNE_ISTL_VERSION VERSION_GREATER 2.7 ) OR ( DUNE_ISTL_VERSION VERSION_EQUAL 2.7 ) ) )"
COMPILE_DEFINITIONS LINEARSOLVER=IstlSolverFactoryBackend<LinearSolverTraits<MassGridGeometry>>
COMPILE_DEFINITIONS LINEARSOLVER=IstlSolverFactoryBackend<LinearSolverTraits<MassGridGeometry>,LinearAlgebraTraitsFromAssembler<Assembler>>
COMMAND ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py
CMD_ARGS --script fuzzy
--files ${CMAKE_SOURCE_DIR}/test/references/test_ff_navierstokes_sincos_instationary-reference.vtu
......
......@@ -42,6 +42,7 @@
#include <dumux/linear/seqsolverbackend.hh>
#include <dumux/linear/linearsolvertraits.hh>
#include <dumux/linear/algebratraits.hh>
#include <dumux/linear/istlsolverfactorybackend.hh>
#include <dumux/multidomain/fvassembler.hh>
......
......@@ -38,6 +38,7 @@
#include <dumux/common/dumuxmessage.hh>
#include <dumux/common/initialize.hh>
#include <dumux/linear/linearsolvertraits.hh>
#include <dumux/linear/algebratraits.hh>
#if DUNE_VERSION_GTE(DUNE_ISTL,2,8)
#include <dumux/linear/istlsolverfactorybackend.hh>
......@@ -152,11 +153,8 @@ int main(int argc, char** argv)
auto assembler = std::make_shared<Assembler>(problem, gridGeometry, gridVariables, timeLoop, xOld);
// the linear solver
#if DUNE_VERSION_GTE(DUNE_ISTL,2,8)
using LinearSolver = IstlSolverFactoryBackend<LinearSolverTraits<GridGeometry>>;
#else
using LinearSolver = AMGBiCGSTABBackend<LinearSolverTraits<GridGeometry>>;
#endif
using LinearSolver = IstlSolverFactoryBackend<LinearSolverTraits<GridGeometry>,
LinearAlgebraTraitsFromAssembler<Assembler>>;
auto linearSolver = std::make_shared<LinearSolver>(leafGridView, gridGeometry->dofMapper());
......
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