From 3871114492c874129e9996865a6bc091df18c7a1 Mon Sep 17 00:00:00 2001 From: Timo Koch Date: Sun, 3 May 2020 21:10:04 +0200 Subject: [PATCH 01/32] [linear] Add a general istl solver backend not based on the factory --- dumux/linear/istlsolvers.hh | 295 ++++++++++++++++++++++++++++++++++++ 1 file changed, 295 insertions(+) create mode 100644 dumux/linear/istlsolvers.hh diff --git a/dumux/linear/istlsolvers.hh b/dumux/linear/istlsolvers.hh new file mode 100644 index 0000000000..90f271f898 --- /dev/null +++ b/dumux/linear/istlsolvers.hh @@ -0,0 +1,295 @@ +// -*- 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 . * + *****************************************************************************/ +/*! + * \file + * \ingroup Linear + * \brief Linear solvers from dune-istl + */ +#ifndef DUMUX_LINEAR_ISTL_SOLVERS_HH +#define DUMUX_LINEAR_ISTL_SOLVERS_HH + +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include + +namespace Dumux::Detail { + +/*! + * \ingroup Linear + * \brief Returns the block level for the preconditioner for a given matrix + * + * \tparam M The matrix. + */ +template +constexpr std::size_t preconditionerBlockLevel() noexcept +{ + return isMultiTypeBlockMatrix::value ? 2 : 1; +} + +template +Dune::SolverCategory::Category solverCategory(const typename LinearSolverTraits::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 { + +/*! + * \ingroup Linear + * \brief Base class for linear solvers + */ +template +class IstlLinearSolver +{ + using XVector = typename InverseOperator::domain_type; + using BVector = typename InverseOperator::range_type; + using Scalar = typename InverseOperator::real_type; +#if HAVE_MPI + using Comm = Dune::OwnerOverlapCopyCommunication, int>; + using ScalarProduct = Dune::ScalarProduct; +#endif + +public: + /*! + * \brief Constructor for sequential solvers + */ + IstlLinearSolver(const std::string& paramGroup = "") + { + if (Dune::MPIHelper::getCollectiveCommunication().size() > 1) + DUNE_THROW(Dune::InvalidStateException, "Using sequential constructor for parallel run. Use signature with gridView and dofMapper!"); + + initializeParameters_(paramGroup); + solverCategory_ = Dune::SolverCategory::sequential; + scalarProduct_ = std::make_shared(); + } + + /*! + * \brief Constructor for parallel and sequential solvers + */ + IstlLinearSolver(const typename LinearSolverTraits::GridView& gridView, + const typename LinearSolverTraits::DofMapper& dofMapper, + const std::string& paramGroup = "") + { + initializeParameters_(paramGroup); +#if HAVE_MPI + solverCategory_ = Detail::solverCategory(gridView); + if (solverCategory_ != Dune::SolverCategory::sequential) + { + parallelHelper_ = std::make_unique>(gridView, dofMapper); + communication_ = std::make_shared(gridView.comm(), solverCategory_); + scalarProduct_ = Dune::createScalarProduct(*communication_, solverCategory_); + parallelHelper_->createParallelIndexSet(*communication_); + } + else + scalarProduct_ = std::make_shared(); +#else + solverCategory_ = Dune::SolverCategory::sequential; +#endif + } + +#if HAVE_MPI + /*! + * \brief Constructor with custom scalar product and communication + */ + IstlLinearSolver(std::shared_ptr communication, + std::shared_ptr scalarProduct, + const typename LinearSolverTraits::GridView& gridView, + const typename LinearSolverTraits::DofMapper& dofMapper, + const std::string& paramGroup = "") + { + initializeParameters_(paramGroup); + solverCategory_ = Detail::solverCategory(gridView); + scalarProduct_ = scalarProduct; + communication_ = communication; + if (solverCategory_ != Dune::SolverCategory::sequential) + { + parallelHelper_ = std::make_unique>(gridView, dofMapper); + parallelHelper_->createParallelIndexSet(communication); + } + } +#endif + + bool solve(Matrix& A, XVector& x, BVector& b) + { +#if HAVE_MPI + return solveSequentialOrParallel_(A, x, b); +#else + return solveSequential_(A, x, b); +#endif + } + + Scalar norm(const XVector& x) const + { +#if HAVE_MPI + 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 vectorHelper(parallelHelper_->gridView(), parallelHelper_->dofMapper()); + vectorHelper.makeNonOverlappingConsistent(y); + return scalarProduct_->norm(y); + } + else +#endif + { + return scalarProduct_->norm(x); + } + } + + const Dune::InverseOperatorResult& result() const + { + return result_; + } + + const std::string& name() const + { + return name_; + } + + void setResidualReduction(double residReduction) + { params_["reduction"] = std::to_string(residReduction); } + +private: + //! reset the linear solver factory + void initializeParameters_(const std::string& paramGroup) + { + params_ = LinearSolverParameters::createParameterTree(paramGroup); + } + + bool solveSequential_(Matrix& A, XVector& x, BVector& b) + { + // construct solver from linear operator + using SequentialTraits = typename LinearSolverTraits::template Sequential; + auto linearOperator = std::make_shared(A); + auto solver = constructPreconditionedSolver_(linearOperator); + + // solve linear system + solver.apply(x, b, result_); + return result_.converged; + } + +#if HAVE_MPI + bool solveSequentialOrParallel_(Matrix& A, XVector& x, BVector& 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 (isMultiTypeBlockMatrix::value) + return solveSequential_(A, x, b); + + else + { + switch (solverCategory_) + { + case Dune::SolverCategory::sequential: + return solveSequential_(A, x, b); + case Dune::SolverCategory::nonoverlapping: + using NOTraits = typename LinearSolverTraits::template ParallelNonoverlapping; + return solveParallel_(A, x, b); + case Dune::SolverCategory::overlapping: + using OTraits = typename LinearSolverTraits::template ParallelOverlapping; + return solveParallel_(A, x, b); + default: DUNE_THROW(Dune::InvalidStateException, "Unknown solver category"); + } + } + } + + template + bool solveParallel_(Matrix& A, XVector& x, BVector& b) + { +#if DUNE_VERSION_GT(DUNE_ISTL,2,7) + // make linear algebra consistent + prepareLinearAlgebraParallel(A, b, *parallelHelper_); + + // construct solver from linear operator + auto linearOperator = std::make_shared(A, *communication_); + auto solver = constructPreconditionedSolver_(linearOperator); + + // solve linear system + solver.apply(x, b, result_); + return result_.converged; +#else + DUNE_THROW(Dune::NotImplemented, "Parallel solvers only available for dune-istl > 2.7"); +#endif + } +#endif // HAVE_MPI + + template + InverseOperator constructPreconditionedSolver_(std::shared_ptr& op) + { + const auto& params = params_.sub("preconditioner"); + using Prec = Dune::Preconditioner; + std::shared_ptr prec = std::make_shared(op, params); + +#if HAVE_MPI && DUNE_VERSION_GT(DUNE_ISTL,2,7) + if (prec->category() != op->category() && prec->category() == Dune::SolverCategory::sequential) + prec = Dune::wrapPreconditioner4Parallel(prec, op); +#endif + return {op, scalarProduct_, prec, params_}; + } + +#if HAVE_MPI + std::unique_ptr> parallelHelper_; + std::shared_ptr communication_; +#endif + Dune::SolverCategory::Category solverCategory_; + std::shared_ptr scalarProduct_; + + Dune::InverseOperatorResult result_; + Dune::ParameterTree params_; + std::string name_; +}; + +template +using ILUBiCGSTABIstlSolver = IstlLinearSolver, + Dune::SeqILU()> + >; + +} // end namespace Dumux + +#endif -- GitLab From 8ae3fdd18708abcd73e16841de8a102dd2247de1 Mon Sep 17 00:00:00 2001 From: Timo Koch Date: Sun, 3 May 2020 21:10:48 +0200 Subject: [PATCH 02/32] [test][swe] Use ILUBiCGSTABIstlSolver in dambreak test --- test/freeflow/shallowwater/dambreak/main.cc | 15 ++++----------- 1 file changed, 4 insertions(+), 11 deletions(-) diff --git a/test/freeflow/shallowwater/dambreak/main.cc b/test/freeflow/shallowwater/dambreak/main.cc index 3d9aa15439..792ddac8cb 100644 --- a/test/freeflow/shallowwater/dambreak/main.cc +++ b/test/freeflow/shallowwater/dambreak/main.cc @@ -41,12 +41,7 @@ #include #include - -#if DUNE_VERSION_GTE(DUNE_ISTL,2,8) -#include -#else -#include -#endif +#include #include @@ -129,11 +124,9 @@ int main(int argc, char** argv) auto assembler = std::make_shared(problem, gridGeometry, gridVariables, timeLoop, xOld); // the linear solver -#if DUNE_VERSION_GTE(DUNE_ISTL,2,8) - using LinearSolver = IstlSolverFactoryBackend>; -#else - using LinearSolver = AMGBiCGSTABBackend>; -#endif + constexpr auto blockSize = std::decay_t::dimension; + using BlockType = Dune::FieldVector; + using LinearSolver = ILUBiCGSTABIstlSolver, Assembler::JacobianMatrix, Dune::BlockVector>; auto linearSolver = std::make_shared(leafGridView, gridGeometry->dofMapper()); // the non-linear solver -- GitLab From 7b9284376274676c6d0961eb0303f2b6d7455432 Mon Sep 17 00:00:00 2001 From: Kilian Weishaupt Date: Tue, 8 Sep 2020 10:19:57 +0200 Subject: [PATCH 03/32] [assembler] Add deprecation warning to residualNorm() --- dumux/assembly/fvassembler.hh | 1 + dumux/multidomain/fvassembler.hh | 1 + 2 files changed, 2 insertions(+) diff --git a/dumux/assembly/fvassembler.hh b/dumux/assembly/fvassembler.hh index e1a6614341..acc41140b2 100644 --- a/dumux/assembly/fvassembler.hh +++ b/dumux/assembly/fvassembler.hh @@ -202,6 +202,7 @@ public: } //! compute the residual and return it's vector norm + [[deprecated("Use norm(curSol) provided by the linear solver class instead. Will be deleted after 3.3")]] Scalar residualNorm(const SolutionVector& curSol) const { ResidualType residual(numDofs()); diff --git a/dumux/multidomain/fvassembler.hh b/dumux/multidomain/fvassembler.hh index f1e6bd9718..1725b0daf9 100644 --- a/dumux/multidomain/fvassembler.hh +++ b/dumux/multidomain/fvassembler.hh @@ -235,6 +235,7 @@ public: } //! compute the residual and return it's vector norm + [[deprecated("Use norm(curSol) provided by the linear solver class instead. Will be deleted after 3.3")]] Scalar residualNorm(const SolutionVector& curSol) { ResidualType residual; -- GitLab From e022365a5576e7e5d65ef485c27380b416107fcb Mon Sep 17 00:00:00 2001 From: Bernd Flemisch Date: Fri, 4 Dec 2020 08:55:12 +0100 Subject: [PATCH 04/32] [linear] add linear algebra traits, use them for the istl solver --- dumux/linear/algebratraits.hh | 46 +++++++++++++++++++++ dumux/linear/istlsolvers.hh | 11 +++-- test/freeflow/shallowwater/dambreak/main.cc | 6 +-- 3 files changed, 56 insertions(+), 7 deletions(-) create mode 100644 dumux/linear/algebratraits.hh diff --git a/dumux/linear/algebratraits.hh b/dumux/linear/algebratraits.hh new file mode 100644 index 0000000000..51715a4b33 --- /dev/null +++ b/dumux/linear/algebratraits.hh @@ -0,0 +1,46 @@ +// -*- 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 . * + *****************************************************************************/ +/*! + * \file + * \ingroup Linear + * \brief Define traits for linear algebra. + */ +#ifndef DUMUX_LINEAR_ALGEBRA_TRAITS_HH +#define DUMUX_LINEAR_ALGEBRA_TRAITS_HH + + +namespace Dumux { + +template +struct LinearAlgebraTraits +{ + using Matrix = M; + using Vector = V; +}; + +template +struct LinearAlgebraTraitsFromAssembler +{ + using Matrix = typename Assembler::JacobianMatrix; + using Vector = typename Assembler::ResidualType; +}; + +} // end namespace Dumux + +#endif diff --git a/dumux/linear/istlsolvers.hh b/dumux/linear/istlsolvers.hh index 90f271f898..bc86493fa0 100644 --- a/dumux/linear/istlsolvers.hh +++ b/dumux/linear/istlsolvers.hh @@ -284,10 +284,13 @@ private: std::string name_; }; -template -using ILUBiCGSTABIstlSolver = IstlLinearSolver, - Dune::SeqILU()> +template +using ILUBiCGSTABIstlSolver = IstlLinearSolver, + Dune::SeqILU()> >; } // end namespace Dumux diff --git a/test/freeflow/shallowwater/dambreak/main.cc b/test/freeflow/shallowwater/dambreak/main.cc index 792ddac8cb..641d034881 100644 --- a/test/freeflow/shallowwater/dambreak/main.cc +++ b/test/freeflow/shallowwater/dambreak/main.cc @@ -41,6 +41,7 @@ #include #include +#include #include #include @@ -124,9 +125,8 @@ int main(int argc, char** argv) auto assembler = std::make_shared(problem, gridGeometry, gridVariables, timeLoop, xOld); // the linear solver - constexpr auto blockSize = std::decay_t::dimension; - using BlockType = Dune::FieldVector; - using LinearSolver = ILUBiCGSTABIstlSolver, Assembler::JacobianMatrix, Dune::BlockVector>; + using LinearSolver = ILUBiCGSTABIstlSolver, + LinearAlgebraTraitsFromAssembler>; auto linearSolver = std::make_shared(leafGridView, gridGeometry->dofMapper()); // the non-linear solver -- GitLab From 4c41292a16fe3073f4639559ac70ac02d5ba6dfd Mon Sep 17 00:00:00 2001 From: Bernd Flemisch Date: Fri, 4 Dec 2020 09:44:13 +0100 Subject: [PATCH 05/32] [linear] make linear algebra traits work for residual vectors with state Apply the logic from the Newton solver. Test in porousmediumflow/2p2c/implicit/injection. --- dumux/linear/algebratraits.hh | 23 +++++++++++++++++++- dumux/linear/istlsolvers.hh | 1 + dumux/nonlinear/newtonsolver.hh | 1 + test/porousmediumflow/2p2c/injection/main.cc | 6 +++-- 4 files changed, 28 insertions(+), 3 deletions(-) diff --git a/dumux/linear/algebratraits.hh b/dumux/linear/algebratraits.hh index 51715a4b33..74240b6021 100644 --- a/dumux/linear/algebratraits.hh +++ b/dumux/linear/algebratraits.hh @@ -24,9 +24,24 @@ #ifndef DUMUX_LINEAR_ALGEBRA_TRAITS_HH #define DUMUX_LINEAR_ALGEBRA_TRAITS_HH +#include +#include + +#include +#include namespace Dumux { +namespace Detail { + +template>::value, int> = 0> +constexpr std::size_t blockSize() { return 1; } + +template>::value, int> = 0> +constexpr std::size_t blockSize() { return std::decay_t::size(); } + +} // end namespace Detail + template struct LinearAlgebraTraits { @@ -37,8 +52,14 @@ struct LinearAlgebraTraits template struct LinearAlgebraTraitsFromAssembler { +private: + using VectorPossiblyWithState = typename Assembler::ResidualType; + using Scalar = std::decay_t()[0][0])>; + static constexpr auto blockSize = Detail::blockSize()[0])>(); + using BlockType = Dune::FieldVector; +public: + using Vector = Dune::BlockVector; using Matrix = typename Assembler::JacobianMatrix; - using Vector = typename Assembler::ResidualType; }; } // end namespace Dumux diff --git a/dumux/linear/istlsolvers.hh b/dumux/linear/istlsolvers.hh index bc86493fa0..1e95b0bf40 100644 --- a/dumux/linear/istlsolvers.hh +++ b/dumux/linear/istlsolvers.hh @@ -26,6 +26,7 @@ #include +#include #include #include #include diff --git a/dumux/nonlinear/newtonsolver.hh b/dumux/nonlinear/newtonsolver.hh index 0701f314f9..d78141353f 100644 --- a/dumux/nonlinear/newtonsolver.hh +++ b/dumux/nonlinear/newtonsolver.hh @@ -53,6 +53,7 @@ #include #include #include +#include #include #include "newtonconvergencewriter.hh" diff --git a/test/porousmediumflow/2p2c/injection/main.cc b/test/porousmediumflow/2p2c/injection/main.cc index 8ba188a6a9..e8b8ae2dfe 100644 --- a/test/porousmediumflow/2p2c/injection/main.cc +++ b/test/porousmediumflow/2p2c/injection/main.cc @@ -32,7 +32,8 @@ #include #include -#include +#include +#include #include #include @@ -128,7 +129,8 @@ int main(int argc, char** argv) auto assembler = std::make_shared(problem, gridGeometry, gridVariables, timeLoop, xOld); // the linear solver - using LinearSolver = AMGBiCGSTABBackend>; + using LinearSolver = ILUBiCGSTABIstlSolver, + LinearAlgebraTraitsFromAssembler>; auto linearSolver = std::make_shared(leafGridView, gridGeometry->dofMapper()); // the non-linear solver -- GitLab From 1722fad1396b1451dd4154ba2f04d0a37ea8e04b Mon Sep 17 00:00:00 2001 From: Bernd Flemisch Date: Mon, 7 Dec 2020 11:51:32 +0100 Subject: [PATCH 06/32] [linear] use algebra traits for ISTL solver --- dumux/linear/istlsolvers.hh | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/dumux/linear/istlsolvers.hh b/dumux/linear/istlsolvers.hh index 1e95b0bf40..8169b4025f 100644 --- a/dumux/linear/istlsolvers.hh +++ b/dumux/linear/istlsolvers.hh @@ -80,11 +80,11 @@ namespace Dumux { * \ingroup Linear * \brief Base class for linear solvers */ -template +template class IstlLinearSolver { + using Matrix = typename LinearAlgebraTraits::Matrix; using XVector = typename InverseOperator::domain_type; using BVector = typename InverseOperator::range_type; using Scalar = typename InverseOperator::real_type; @@ -286,7 +286,7 @@ private: }; template -using ILUBiCGSTABIstlSolver = IstlLinearSolver, Dune::SeqILU Date: Wed, 9 Dec 2020 16:04:38 +0100 Subject: [PATCH 07/32] [linear] make the ISTL solver backend work for multitype --- dumux/linear/algebratraits.hh | 17 ++- dumux/linear/istlsolvers.hh | 133 +++++++++++++----- .../linearsolveracceptsmultitypematrix.hh | 9 ++ dumux/linear/linearsolverparameters.hh | 13 +- dumux/linear/linearsolvertraits.hh | 8 ++ dumux/nonlinear/newtonsolver.hh | 4 +- test/multidomain/poromechanics/el1p/main.cc | 8 +- 7 files changed, 152 insertions(+), 40 deletions(-) diff --git a/dumux/linear/algebratraits.hh b/dumux/linear/algebratraits.hh index 74240b6021..a11399ae59 100644 --- a/dumux/linear/algebratraits.hh +++ b/dumux/linear/algebratraits.hh @@ -30,6 +30,8 @@ #include #include +#include + namespace Dumux { namespace Detail { @@ -49,8 +51,8 @@ struct LinearAlgebraTraits using Vector = V; }; -template -struct LinearAlgebraTraitsFromAssembler +template +struct LATraitsFromAssemblerImpl { private: using VectorPossiblyWithState = typename Assembler::ResidualType; @@ -62,6 +64,17 @@ public: using Matrix = typename Assembler::JacobianMatrix; }; +template +struct LATraitsFromAssemblerImpl +{ + using Vector = typename Assembler::ResidualType; + using Matrix = typename Assembler::JacobianMatrix; +}; + +template +using LinearAlgebraTraitsFromAssembler = LATraitsFromAssemblerImpl::value>; + } // end namespace Dumux #endif diff --git a/dumux/linear/istlsolvers.hh b/dumux/linear/istlsolvers.hh index 8169b4025f..f37b2d2da0 100644 --- a/dumux/linear/istlsolvers.hh +++ b/dumux/linear/istlsolvers.hh @@ -38,7 +38,10 @@ #include #include +#include #include +#include +#include #include namespace Dumux::Detail { @@ -55,8 +58,8 @@ constexpr std::size_t preconditionerBlockLevel() noexcept return isMultiTypeBlockMatrix::value ? 2 : 1; } -template -Dune::SolverCategory::Category solverCategory(const typename LinearSolverTraits::GridView& gridView) +template +Dune::SolverCategory::Category solverCategory(const GridView& gridView) { if constexpr (LinearSolverTraits::canCommunicate) { @@ -85,12 +88,15 @@ template, int>; - using ScalarProduct = Dune::ScalarProduct; + using ScalarProduct = Dune::ScalarProduct; + using ParallelHelper = std::conditional_t::value, + double, + ParallelISTLHelper>; #endif public: @@ -110,8 +116,9 @@ public: /*! * \brief Constructor for parallel and sequential solvers */ - IstlLinearSolver(const typename LinearSolverTraits::GridView& gridView, - const typename LinearSolverTraits::DofMapper& dofMapper, + template + IstlLinearSolver(const GridView& gridView, + const DofMapper& dofMapper, const std::string& paramGroup = "") { initializeParameters_(paramGroup); @@ -135,10 +142,11 @@ public: /*! * \brief Constructor with custom scalar product and communication */ + template IstlLinearSolver(std::shared_ptr communication, std::shared_ptr scalarProduct, - const typename LinearSolverTraits::GridView& gridView, - const typename LinearSolverTraits::DofMapper& dofMapper, + const GridView& gridView, + const DofMapper& dofMapper, const std::string& paramGroup = "") { initializeParameters_(paramGroup); @@ -165,20 +173,27 @@ public: Scalar norm(const XVector& x) const { #if HAVE_MPI - if (solverCategory_ == Dune::SolverCategory::nonoverlapping) + if constexpr (LinearSolverTraits::canCommunicate) { - 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 vectorHelper(parallelHelper_->gridView(), parallelHelper_->dofMapper()); - vectorHelper.makeNonOverlappingConsistent(y); - return scalarProduct_->norm(y); + 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 vectorHelper(parallelHelper_->gridView(), parallelHelper_->dofMapper()); + vectorHelper.makeNonOverlappingConsistent(y); + return scalarProduct_->norm(y); + } } - else #endif + if constexpr (!preconditionerAcceptsMultiTypeMatrix::value + && isMultiTypeBlockVector::value) { - return scalarProduct_->norm(x); + auto y = VectorConverter::multiTypeToBlockVector(x); + return scalarProduct_->norm(y); } + else + return scalarProduct_->norm(x); } const Dune::InverseOperatorResult& result() const @@ -203,13 +218,46 @@ private: bool solveSequential_(Matrix& A, XVector& x, BVector& b) { - // construct solver from linear operator - using SequentialTraits = typename LinearSolverTraits::template Sequential; - auto linearOperator = std::make_shared(A); - auto solver = constructPreconditionedSolver_(linearOperator); + if constexpr (preconditionerAcceptsMultiTypeMatrix::value + || !isMultiTypeBlockMatrix::value) + { + // construct solver from linear operator + using SequentialTraits = typename LinearSolverTraits::template Sequential; + auto linearOperator = std::make_shared(A); + auto solver = constructPreconditionedSolver_(linearOperator); + + // solve linear system + solver.apply(x, b, result_); + } + else + { + // create the bcrs matrix the IterativeSolver backend can handle + auto M = MatrixConverter::multiTypeToBCRSMatrix(A); + + // get the new matrix sizes + const std::size_t numRows = M.N(); + assert(numRows == M.M()); + + // create the vector the IterativeSolver backend can handle + auto bTmp = VectorConverter::multiTypeToBlockVector(b); + assert(bTmp.size() == numRows); + + // create a blockvector to which the linear solver writes the solution + using VectorBlock = typename Dune::FieldVector; + using BlockVector = typename Dune::BlockVector; + BlockVector y(numRows); + + auto linearOperator = std::make_shared>(M); + auto solver = constructPreconditionedSolver_(linearOperator); + + // solve linear system + solver.apply(y, bTmp, result_); + + // copy back the result y into x + if(result_.converged) + VectorConverter::retrieveValues(x, y); + } - // solve linear system - solver.apply(x, b, result_); return result_.converged; } @@ -221,7 +269,6 @@ private: // TODO: This can be adapted once the situation in Dune ISTL changes. if constexpr (isMultiTypeBlockMatrix::value) return solveSequential_(A, x, b); - else { switch (solverCategory_) @@ -274,7 +321,7 @@ private: } #if HAVE_MPI - std::unique_ptr> parallelHelper_; + std::unique_ptr parallelHelper_; std::shared_ptr communication_; #endif Dune::SolverCategory::Category solverCategory_; @@ -285,14 +332,34 @@ private: std::string name_; }; +template +struct ILUBiCGSTABIstlSolverImpl; + +template +struct ILUBiCGSTABIstlSolverImpl +{ + using type = IstlLinearSolver, + Dune::SeqILU + >; +}; + +template +struct ILUBiCGSTABIstlSolverImpl +{ + using type = IstlLinearSolver::multiTypeToBlockVector(std::declval()))>, + Dune::SeqILU::multiTypeToBCRSMatrix(std::declval())), + decltype(VectorConverter::multiTypeToBlockVector(std::declval())), + decltype(VectorConverter::multiTypeToBlockVector(std::declval()))> + >; +}; + template -using ILUBiCGSTABIstlSolver = IstlLinearSolver, - Dune::SeqILU()> - >; +using ILUBiCGSTABIstlSolver = typename ILUBiCGSTABIstlSolverImpl::value>::type; } // end namespace Dumux diff --git a/dumux/linear/linearsolveracceptsmultitypematrix.hh b/dumux/linear/linearsolveracceptsmultitypematrix.hh index de471c0639..57128c6cfc 100644 --- a/dumux/linear/linearsolveracceptsmultitypematrix.hh +++ b/dumux/linear/linearsolveracceptsmultitypematrix.hh @@ -24,7 +24,10 @@ #ifndef DUMUX_LINEAR_SOLVER_ACCEPTS_MULTITYPEMATRIX_HH #define DUMUX_LINEAR_SOLVER_ACCEPTS_MULTITYPEMATRIX_HH +#include + #include +#include namespace Dumux { @@ -64,6 +67,12 @@ template<> struct linearSolverAcceptsMultiTypeMatrix : public std::false_type {}; #endif // HAVE_UMFPACK +template +struct preconditionerAcceptsMultiTypeMatrix : public std::true_type {}; + +template +struct preconditionerAcceptsMultiTypeMatrix> : public std::false_type {}; + } // end namespace Dumux #endif diff --git a/dumux/linear/linearsolverparameters.hh b/dumux/linear/linearsolverparameters.hh index b222c9dc71..1a230607e8 100644 --- a/dumux/linear/linearsolverparameters.hh +++ b/dumux/linear/linearsolverparameters.hh @@ -34,6 +34,14 @@ namespace Dumux { +namespace Detail { +template +using GVDetector = typename T::GridView; + +template +using hasGridView = Dune::Std::is_detected; +} + template class LinearSolverParameters { @@ -62,7 +70,10 @@ public: params["preconditioner.relaxation"] = "1.0"; params["preconditioner.verbosity"] = "0"; params["preconditioner.defaultAggregationSizeMode"] = "isotropic"; - params["preconditioner.defaultAggregationDimension"] = std::to_string(LinearSolverTraits::GridView::dimension); + if constexpr (Detail::hasGridView::value) + params["preconditioner.defaultAggregationDimension"] = std::to_string(LinearSolverTraits::GridView::dimension); + else + params["preconditioner.defaultAggregationDimension"] = 3; params["preconditioner.maxLevel"] = "100"; params["ParameterGroup"] = paramGroup; params["preconditioner.ParameterGroup"] = paramGroup; diff --git a/dumux/linear/linearsolvertraits.hh b/dumux/linear/linearsolvertraits.hh index 146e63b40c..4e8556f167 100644 --- a/dumux/linear/linearsolvertraits.hh +++ b/dumux/linear/linearsolvertraits.hh @@ -57,6 +57,14 @@ struct SequentialSolverTraits using Preconditioner = SeqPreconditioner; }; +struct SeqLinearSolverTraits +{ + template + using Sequential = SequentialSolverTraits; + + static constexpr bool canCommunicate = false; +}; + #if HAVE_MPI template struct NonoverlappingSolverTraits diff --git a/dumux/nonlinear/newtonsolver.hh b/dumux/nonlinear/newtonsolver.hh index d78141353f..534473e4ec 100644 --- a/dumux/nonlinear/newtonsolver.hh +++ b/dumux/nonlinear/newtonsolver.hh @@ -1223,14 +1223,14 @@ private: assert(this->checkSizesOfSubMatrices(A) && "Sub-blocks of MultiTypeBlockMatrix have wrong sizes!"); // create the bcrs matrix the IterativeSolver backend can handle - const auto M = MatrixConverter::multiTypeToBCRSMatrix(A); + auto M = MatrixConverter::multiTypeToBCRSMatrix(A); // get the new matrix sizes const std::size_t numRows = M.N(); assert(numRows == M.M()); // create the vector the IterativeSolver backend can handle - const auto bTmp = VectorConverter::multiTypeToBlockVector(b); + auto bTmp = VectorConverter::multiTypeToBlockVector(b); assert(bTmp.size() == numRows); // create a blockvector to which the linear solver writes the solution diff --git a/test/multidomain/poromechanics/el1p/main.cc b/test/multidomain/poromechanics/el1p/main.cc index 0c8fd89651..a5ee0334c7 100644 --- a/test/multidomain/poromechanics/el1p/main.cc +++ b/test/multidomain/poromechanics/el1p/main.cc @@ -31,7 +31,10 @@ #include #include -#include +#include +#include +#include + #include #include #include @@ -141,7 +144,8 @@ int main(int argc, char** argv) couplingManager); // the linear solver - using LinearSolver = ILU0BiCGSTABBackend; + using LinearSolver = ILUBiCGSTABIstlSolver>; auto linearSolver = std::make_shared(); // the non-linear solver -- GitLab From 00cf670c40a77cfb3a86c0cb0446a93d5791c8ee Mon Sep 17 00:00:00 2001 From: Bernd Flemisch Date: Wed, 16 Dec 2020 16:03:22 +0100 Subject: [PATCH 08/32] [linear] Specialize parallel helper for sequential settings --- dumux/linear/istlsolvers.hh | 4 +--- dumux/linear/parallelhelpers.hh | 18 ++++++++++++++---- 2 files changed, 15 insertions(+), 7 deletions(-) diff --git a/dumux/linear/istlsolvers.hh b/dumux/linear/istlsolvers.hh index f37b2d2da0..31786c5892 100644 --- a/dumux/linear/istlsolvers.hh +++ b/dumux/linear/istlsolvers.hh @@ -94,9 +94,7 @@ class IstlLinearSolver #if HAVE_MPI using Comm = Dune::OwnerOverlapCopyCommunication, int>; using ScalarProduct = Dune::ScalarProduct; - using ParallelHelper = std::conditional_t::value, - double, - ParallelISTLHelper>; + using ParallelHelper = ParallelISTLHelper; #endif public: diff --git a/dumux/linear/parallelhelpers.hh b/dumux/linear/parallelhelpers.hh index f8cf755489..26ba484f8b 100644 --- a/dumux/linear/parallelhelpers.hh +++ b/dumux/linear/parallelhelpers.hh @@ -41,9 +41,11 @@ namespace Dumux { * \brief A parallel helper class providing a nonoverlapping * decomposition of all degrees of freedom */ -// operator that resets result to zero at constrained DOFS +template +class ParallelISTLHelperImpl; + template -class ParallelISTLHelper +class ParallelISTLHelperImpl { using GridView = typename LinearSolverTraits::GridView; using DofMapper = typename LinearSolverTraits::DofMapper; @@ -276,7 +278,7 @@ class ParallelISTLHelper public: - ParallelISTLHelper(const GridView& gridView, const DofMapper& mapper) + ParallelISTLHelperImpl(const GridView& gridView, const DofMapper& mapper) : gridView_(gridView), mapper_(mapper) { if constexpr (Detail::canCommunicate) @@ -423,7 +425,15 @@ private: std::vector isOwned_; //!< vector to identify unique decomposition std::vector isGhost_; //!< vector to identify ghost dofs -}; // class ParallelISTLHelper +}; // class ParallelISTLHelperImpl + +template +class ParallelISTLHelperImpl +{}; + +template +using ParallelISTLHelper = ParallelISTLHelperImpl; template class ParallelVectorHelper -- GitLab From 4c3b6e02ef6ebacbdee71de2f30e6b8099a0be88 Mon Sep 17 00:00:00 2001 From: Bernd Flemisch Date: Wed, 16 Dec 2020 16:16:06 +0100 Subject: [PATCH 09/32] [linear] improve hasGridView construct --- dumux/linear/linearsolverparameters.hh | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/dumux/linear/linearsolverparameters.hh b/dumux/linear/linearsolverparameters.hh index 1a230607e8..0a6c9979ce 100644 --- a/dumux/linear/linearsolverparameters.hh +++ b/dumux/linear/linearsolverparameters.hh @@ -39,7 +39,7 @@ template using GVDetector = typename T::GridView; template -using hasGridView = Dune::Std::is_detected; +constexpr bool hasGridView = Dune::Std::is_detected::value; } template @@ -70,7 +70,7 @@ public: params["preconditioner.relaxation"] = "1.0"; params["preconditioner.verbosity"] = "0"; params["preconditioner.defaultAggregationSizeMode"] = "isotropic"; - if constexpr (Detail::hasGridView::value) + if constexpr (Detail::hasGridView) params["preconditioner.defaultAggregationDimension"] = std::to_string(LinearSolverTraits::GridView::dimension); else params["preconditioner.defaultAggregationDimension"] = 3; -- GitLab From 1e32631a3bcf065adcb49c1a4be5ba01c93356c5 Mon Sep 17 00:00:00 2001 From: Bernd Flemisch Date: Wed, 16 Dec 2020 16:49:37 +0100 Subject: [PATCH 10/32] [linear] simplify branching between single and multitype matrices/vectors --- dumux/linear/algebratraits.hh | 5 ++++- dumux/linear/istlsolvers.hh | 41 ++++++++++++++++++----------------- 2 files changed, 25 insertions(+), 21 deletions(-) diff --git a/dumux/linear/algebratraits.hh b/dumux/linear/algebratraits.hh index a11399ae59..f7fe9930d0 100644 --- a/dumux/linear/algebratraits.hh +++ b/dumux/linear/algebratraits.hh @@ -30,7 +30,8 @@ #include #include -#include +#include +#include namespace Dumux { @@ -69,6 +70,8 @@ struct LATraitsFromAssemblerImpl { using Vector = typename Assembler::ResidualType; using Matrix = typename Assembler::JacobianMatrix; + using SingleTypeVector = decltype(VectorConverter::multiTypeToBlockVector(std::declval())); + using SingleTypeMatrix = decltype(MatrixConverter::multiTypeToBCRSMatrix(std::declval())); }; template diff --git a/dumux/linear/istlsolvers.hh b/dumux/linear/istlsolvers.hh index 31786c5892..6da8ee8712 100644 --- a/dumux/linear/istlsolvers.hh +++ b/dumux/linear/istlsolvers.hh @@ -91,6 +91,9 @@ class IstlLinearSolver using XVector = typename LinearAlgebraTraits::Vector; using BVector = typename LinearAlgebraTraits::Vector; using Scalar = typename InverseOperator::real_type; + static constexpr bool convertMatrixAndVector + = !preconditionerAcceptsMultiTypeMatrix::value + && isMultiTypeBlockVector::value; #if HAVE_MPI using Comm = Dune::OwnerOverlapCopyCommunication, int>; using ScalarProduct = Dune::ScalarProduct; @@ -184,8 +187,7 @@ public: } } #endif - if constexpr (!preconditionerAcceptsMultiTypeMatrix::value - && isMultiTypeBlockVector::value) + if constexpr (convertMatrixAndVector) { auto y = VectorConverter::multiTypeToBlockVector(x); return scalarProduct_->norm(y); @@ -216,18 +218,7 @@ private: bool solveSequential_(Matrix& A, XVector& x, BVector& b) { - if constexpr (preconditionerAcceptsMultiTypeMatrix::value - || !isMultiTypeBlockMatrix::value) - { - // construct solver from linear operator - using SequentialTraits = typename LinearSolverTraits::template Sequential; - auto linearOperator = std::make_shared(A); - auto solver = constructPreconditionedSolver_(linearOperator); - - // solve linear system - solver.apply(x, b, result_); - } - else + if constexpr (convertMatrixAndVector) { // create the bcrs matrix the IterativeSolver backend can handle auto M = MatrixConverter::multiTypeToBCRSMatrix(A); @@ -255,6 +246,16 @@ private: if(result_.converged) VectorConverter::retrieveValues(x, y); } + else + { + // construct solver from linear operator + using SequentialTraits = typename LinearSolverTraits::template Sequential; + auto linearOperator = std::make_shared(A); + auto solver = constructPreconditionedSolver_(linearOperator); + + // solve linear system + solver.apply(x, b, result_); + } return result_.converged; } @@ -347,12 +348,12 @@ struct ILUBiCGSTABIstlSolverImpl template struct ILUBiCGSTABIstlSolverImpl { - using type = IstlLinearSolver::multiTypeToBlockVector(std::declval()))>, - Dune::SeqILU::multiTypeToBCRSMatrix(std::declval())), - decltype(VectorConverter::multiTypeToBlockVector(std::declval())), - decltype(VectorConverter::multiTypeToBlockVector(std::declval()))> - >; + using type = IstlLinearSolver, + Dune::SeqILU + >; }; template -- GitLab From ba2f0ecd45d3ffd1f2f2f1c869b36765b9fb2b38 Mon Sep 17 00:00:00 2001 From: Bernd Flemisch Date: Mon, 21 Dec 2020 08:48:13 +0100 Subject: [PATCH 11/32] [linear] remove circular dependency --- dumux/linear/linearsolveracceptsmultitypematrix.hh | 1 - 1 file changed, 1 deletion(-) diff --git a/dumux/linear/linearsolveracceptsmultitypematrix.hh b/dumux/linear/linearsolveracceptsmultitypematrix.hh index 57128c6cfc..abfd44b623 100644 --- a/dumux/linear/linearsolveracceptsmultitypematrix.hh +++ b/dumux/linear/linearsolveracceptsmultitypematrix.hh @@ -27,7 +27,6 @@ #include #include -#include namespace Dumux { -- GitLab From 1dbe1fc363114feb9c401a46c4717adb4ac59ba9 Mon Sep 17 00:00:00 2001 From: Bernd Flemisch Date: Mon, 21 Dec 2020 08:48:41 +0100 Subject: [PATCH 12/32] [linear] add norm to base class --- dumux/linear/solver.hh | 12 ++++++++++++ 1 file changed, 12 insertions(+) diff --git a/dumux/linear/solver.hh b/dumux/linear/solver.hh index 0e73ac6b21..75f64be3f3 100644 --- a/dumux/linear/solver.hh +++ b/dumux/linear/solver.hh @@ -25,6 +25,7 @@ #define DUMUX_LINEAR_SOLVER_HH #include +#include #include namespace Dumux { @@ -71,6 +72,17 @@ public: DUNE_THROW(Dune::NotImplemented, "Linear solver doesn't implement a solve method!"); } + template + auto norm(const Vector& x) const + { + if (Dune::MPIHelper::getCollectiveCommunication().size() > 1) + DUNE_THROW(Dune::NotImplemented, "norm in parallel"); + + static const Dune::SeqScalarProduct sp; + + return sp.norm(x); + } + //! the name of the linear solver std::string name() const { return "unknown solver"; } -- GitLab From f3b17c2909657019ad30d4613d01ff62fd1afe2c Mon Sep 17 00:00:00 2001 From: Bernd Flemisch Date: Mon, 21 Dec 2020 17:03:57 +0100 Subject: [PATCH 13/32] [linear] adapt the AMG backend to the new structure --- dumux/linear/amgbackend.hh | 252 +++++++++++++++++++- dumux/linear/istlsolvers.hh | 18 +- dumux/linear/solvercategory.hh | 45 ++++ test/porousmediumflow/2p2c/waterair/main.cc | 4 +- 4 files changed, 289 insertions(+), 30 deletions(-) create mode 100644 dumux/linear/solvercategory.hh diff --git a/dumux/linear/amgbackend.hh b/dumux/linear/amgbackend.hh index 468d4f7dd9..14285f3c9d 100644 --- a/dumux/linear/amgbackend.hh +++ b/dumux/linear/amgbackend.hh @@ -38,16 +38,14 @@ #include #include +#include namespace Dumux { -/*! - * \ingroup Linear - * \brief A linear solver based on the ISTL AMG preconditioner - * and the ISTL BiCGSTAB solver. - */ +namespace Detail { + template -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 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 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>(gridView, dofMapper); + if constexpr (LinearSolverTraits::canCommunicate) + { + if (isParallel_) + phelper_ = std::make_unique>(gridView, dofMapper); + } #endif checkAvailabilityOfDirectSolver_(); } @@ -121,6 +124,9 @@ public: return result_; } + template + 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 NewAMGBiCGSTABBackend : public LinearSolver +{ + using Matrix = typename LinearAlgebraTraits::Matrix; + using Vector = typename LinearAlgebraTraits::Vector; + using ScalarProduct = Dune::ScalarProduct; +#if HAVE_MPI + using Comm = Dune::OwnerOverlapCopyCommunication, 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(); + } + + /*! + * \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(gridView); + if (solverCategory_ != Dune::SolverCategory::sequential) + { + phelper_ = std::make_unique>(gridView, dofMapper); + comm_ = std::make_shared(gridView.comm(), solverCategory_); + scalarProduct_ = Dune::createScalarProduct(*comm_, solverCategory_); + phelper_->createParallelIndexSet(*comm_); + } + else + scalarProduct_ = std::make_shared(); +#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 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; + solveParallel_(A, x, b); + } + else + { + using PTraits = typename LinearSolverTraits::template ParallelOverlapping; + solveParallel_(A, x, b); + } + } + else + solveSequential_(A, x, b); + } + else + { + solveSequential_(A, x, b); + } + } + + template + void solveParallel_(Matrix& A, Vector& x, Vector& b) + { + prepareLinearAlgebraParallel(A, b, *phelper_); + auto linearOperator = std::make_shared(A, *comm_); + + using SeqSmoother = Dune::SeqSSOR; + using Smoother = typename ParallelTraits::template Preconditioner; + solveWithAmg_(A, x, b, linearOperator, comm_); + } +#endif // HAVE_MPI + + void solveSequential_(Matrix& A, Vector& x, Vector& b) + { + using SequentialTraits = typename LinearSolverTraits::template Sequential; + auto linearOperator = std::make_shared(A); + + auto comm = std::make_shared(); + + using Smoother = Dune::SeqSSOR; + solveWithAmg_(A, x, b, linearOperator, comm); + } + + template + void solveWithAmg_(Matrix& A, Vector& x, Vector& b, + std::shared_ptr& linearOperator, + std::shared_ptr& comm) + { + using SmootherArgs = typename Dune::Amg::SmootherTraits::Arguments; + using Criterion = Dune::Amg::CoarsenCriterion>; + + //! \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; + auto amg = std::make_shared(*linearOperator, criterion, smootherArgs, *comm); + + Dune::BiCGSTABSolver 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> phelper_; + std::shared_ptr comm_; +#endif + Dune::SolverCategory::Category solverCategory_; + std::shared_ptr scalarProduct_; + Dune::InverseOperatorResult result_; + bool isParallel_ = false; +}; + +template +struct AMGImplHelper +{ + template + using type = NewAMGBiCGSTABBackend; +}; + +template<> +struct AMGImplHelper<1> +{ + template + using type = OldAMGBiCGSTABBackend; +}; + +} // end namespace Detail + +/*! + * \ingroup Linear + * \brief A linear solver based on the ISTL AMG preconditioner + * and the ISTL BiCGSTAB solver. + */ +template +using AMGBiCGSTABBackend = typename Detail::AMGImplHelper::template type; + } // end namespace Dumux #endif diff --git a/dumux/linear/istlsolvers.hh b/dumux/linear/istlsolvers.hh index 6da8ee8712..15c6375418 100644 --- a/dumux/linear/istlsolvers.hh +++ b/dumux/linear/istlsolvers.hh @@ -43,6 +43,7 @@ #include #include #include +#include namespace Dumux::Detail { @@ -58,23 +59,6 @@ constexpr std::size_t preconditionerBlockLevel() noexcept return isMultiTypeBlockMatrix::value ? 2 : 1; } -template -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 { diff --git a/dumux/linear/solvercategory.hh b/dumux/linear/solvercategory.hh new file mode 100644 index 0000000000..725cdc0139 --- /dev/null +++ b/dumux/linear/solvercategory.hh @@ -0,0 +1,45 @@ +// -*- 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 . * + *****************************************************************************/ +#ifndef DUMUX_SOLVERCATEGORY_HH +#define DUMUX_SOLVERCATEGORY_HH + +#include + +namespace Dumux::Detail { + +template +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 diff --git a/test/porousmediumflow/2p2c/waterair/main.cc b/test/porousmediumflow/2p2c/waterair/main.cc index 6faa549539..74830f3ed8 100644 --- a/test/porousmediumflow/2p2c/waterair/main.cc +++ b/test/porousmediumflow/2p2c/waterair/main.cc @@ -34,6 +34,7 @@ #include #include +#include #include #include @@ -111,7 +112,8 @@ int main(int argc, char** argv) auto assembler = std::make_shared(problem, gridGeometry, gridVariables, timeLoop, xOld); // the linear solver - using LinearSolver = AMGBiCGSTABBackend>; + using LinearSolver = AMGBiCGSTABBackend, + LinearAlgebraTraitsFromAssembler>; auto linearSolver = std::make_shared(leafGridView, gridGeometry->dofMapper()); // the non-linear solver -- GitLab From 1f6033479c8347935fda20b0df25cba3b48fa66f Mon Sep 17 00:00:00 2001 From: Bernd Flemisch Date: Mon, 21 Dec 2020 21:23:34 +0100 Subject: [PATCH 14/32] [linear] adapt ISTL solver factory backend to new structure --- dumux/linear/istlsolverfactorybackend.hh | 268 +++++++++++++++++- .../navierstokes/sincos/CMakeLists.txt | 2 +- test/freeflow/navierstokes/sincos/main.cc | 1 + test/porousmediumflow/richards/lens/main.cc | 8 +- 4 files changed, 267 insertions(+), 12 deletions(-) diff --git a/dumux/linear/istlsolverfactorybackend.hh b/dumux/linear/istlsolverfactorybackend.hh index 9c90f0f83e..f1d87478f9 100644 --- a/dumux/linear/istlsolverfactorybackend.hh +++ b/dumux/linear/istlsolverfactorybackend.hh @@ -47,6 +47,7 @@ #include #include #include +#include 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 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 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 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 NewIstlSolverFactoryBackend : public LinearSolver +{ + using Matrix = typename LinearAlgebraTraits::Matrix; + using Vector = typename LinearAlgebraTraits::Vector; + using ScalarProduct = Dune::ScalarProduct; +#if HAVE_MPI + using Comm = Dune::OwnerOverlapCopyCommunication, 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(); + 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(gridView); + if (solverCategory_ != Dune::SolverCategory::sequential) + { + parallelHelper_ = std::make_unique>(gridView, dofMapper); + comm_ = std::make_shared(gridView.comm(), solverCategory_); + scalarProduct_ = Dune::createScalarProduct(*comm_, solverCategory_); + parallelHelper_->createParallelIndexSet(*comm_); + } + else + scalarProduct_ = std::make_shared(); +#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 vectorHelper(parallelHelper_->gridView(), parallelHelper_->dofMapper()); + vectorHelper.makeNonOverlappingConsistent(y); + return scalarProduct_->norm(y); + } + } +#endif + return scalarProduct_->norm(x); + } + +private: + + void initializeParameters_() + { + params_ = LinearSolverParameters::createParameterTree(paramGroup_); + checkMandatoryParameters_(); + name_ = params_.get("preconditioner.type") + "-preconditioned " + params_.get("type"); + if (params_.get("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::value) + { + if (isParallel_) + { + if (LinearSolverTraits::isNonOverlapping(parallelHelper_->gridView())) + { + using PTraits = typename LinearSolverTraits::template ParallelNonoverlapping; + solveParallel_(A, x, b); + } + else + { + using PTraits = typename LinearSolverTraits::template ParallelOverlapping; + solveParallel_(A, x, b); + } + } + else + solveSequential_(A, x, b); + } + else + { + solveSequential_(A, x, b); + } + } + + template + 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(); + + prepareLinearAlgebraParallel(A, b, *parallelHelper_); + auto linearOperator = std::make_shared(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; + using LinearOperator = typename Traits::LinearOperator; + auto linearOperator = std::make_shared(A); + + if (firstCall_) + initSolverFactories(); + + // construct solver + auto solver = getSolverFromFactory_(linearOperator); + + // solve linear system + solver->apply(x, b, result_); + } + + template + auto getSolverFromFactory_(std::shared_ptr& 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> parallelHelper_; + std::shared_ptr comm_; +#endif + bool isParallel_ = false; + bool firstCall_; + + std::shared_ptr scalarProduct_; + Dune::SolverCategory::Category solverCategory_; + Dune::InverseOperatorResult result_; + Dune::ParameterTree params_; + std::string name_; +}; + +template +struct BackendImplHelper +{ + template + using type = NewIstlSolverFactoryBackend; +}; + +template<> +struct BackendImplHelper<1> +{ + template + using type = OldIstlSolverFactoryBackend; +}; + +} // 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 +using IstlSolverFactoryBackend = typename Detail::BackendImplHelper::template type; + } // end namespace Dumux #endif diff --git a/test/freeflow/navierstokes/sincos/CMakeLists.txt b/test/freeflow/navierstokes/sincos/CMakeLists.txt index 8f0d7d3d5b..2ac2f497a6 100644 --- a/test/freeflow/navierstokes/sincos/CMakeLists.txt +++ b/test/freeflow/navierstokes/sincos/CMakeLists.txt @@ -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> + COMPILE_DEFINITIONS LINEARSOLVER=IstlSolverFactoryBackend,LinearAlgebraTraitsFromAssembler> 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 diff --git a/test/freeflow/navierstokes/sincos/main.cc b/test/freeflow/navierstokes/sincos/main.cc index 563a6536b7..b5e0e5ad9a 100644 --- a/test/freeflow/navierstokes/sincos/main.cc +++ b/test/freeflow/navierstokes/sincos/main.cc @@ -42,6 +42,7 @@ #include #include +#include #include #include diff --git a/test/porousmediumflow/richards/lens/main.cc b/test/porousmediumflow/richards/lens/main.cc index cdf033b772..57fa47a8e0 100644 --- a/test/porousmediumflow/richards/lens/main.cc +++ b/test/porousmediumflow/richards/lens/main.cc @@ -38,6 +38,7 @@ #include #include #include +#include #if DUNE_VERSION_GTE(DUNE_ISTL,2,8) #include @@ -152,11 +153,8 @@ int main(int argc, char** argv) auto assembler = std::make_shared(problem, gridGeometry, gridVariables, timeLoop, xOld); // the linear solver -#if DUNE_VERSION_GTE(DUNE_ISTL,2,8) - using LinearSolver = IstlSolverFactoryBackend>; -#else - using LinearSolver = AMGBiCGSTABBackend>; -#endif + using LinearSolver = IstlSolverFactoryBackend, + LinearAlgebraTraitsFromAssembler>; auto linearSolver = std::make_shared(leafGridView, gridGeometry->dofMapper()); -- GitLab From 7e01b9f944597cfda126e6e8ba5ff54dffff5a02 Mon Sep 17 00:00:00 2001 From: Bernd Flemisch Date: Mon, 21 Dec 2020 21:28:58 +0100 Subject: [PATCH 15/32] Apply 3 suggestion(s) to 2 file(s) --- dumux/linear/parallelhelpers.hh | 8 ++------ dumux/linear/solver.hh | 4 +--- 2 files changed, 3 insertions(+), 9 deletions(-) diff --git a/dumux/linear/parallelhelpers.hh b/dumux/linear/parallelhelpers.hh index 26ba484f8b..50fa016ea2 100644 --- a/dumux/linear/parallelhelpers.hh +++ b/dumux/linear/parallelhelpers.hh @@ -41,8 +41,8 @@ namespace Dumux { * \brief A parallel helper class providing a nonoverlapping * decomposition of all degrees of freedom */ -template -class ParallelISTLHelperImpl; +template +class ParallelISTLHelperImpl {}; template class ParallelISTLHelperImpl @@ -427,10 +427,6 @@ private: }; // class ParallelISTLHelperImpl -template -class ParallelISTLHelperImpl -{}; - template using ParallelISTLHelper = ParallelISTLHelperImpl; diff --git a/dumux/linear/solver.hh b/dumux/linear/solver.hh index 75f64be3f3..1ed4ebaef5 100644 --- a/dumux/linear/solver.hh +++ b/dumux/linear/solver.hh @@ -78,9 +78,7 @@ public: if (Dune::MPIHelper::getCollectiveCommunication().size() > 1) DUNE_THROW(Dune::NotImplemented, "norm in parallel"); - static const Dune::SeqScalarProduct sp; - - return sp.norm(x); + return Dune::SeqScalarProduct().norm(x); } //! the name of the linear solver -- GitLab From a53201bc184b7b964ad36527c5472d72cb78f049 Mon Sep 17 00:00:00 2001 From: Bernd Flemisch Date: Mon, 21 Dec 2020 23:21:13 +0100 Subject: [PATCH 16/32] Apply 2 suggestion(s) to 1 file(s) --- dumux/nonlinear/newtonsolver.hh | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/dumux/nonlinear/newtonsolver.hh b/dumux/nonlinear/newtonsolver.hh index 534473e4ec..d78141353f 100644 --- a/dumux/nonlinear/newtonsolver.hh +++ b/dumux/nonlinear/newtonsolver.hh @@ -1223,14 +1223,14 @@ private: assert(this->checkSizesOfSubMatrices(A) && "Sub-blocks of MultiTypeBlockMatrix have wrong sizes!"); // create the bcrs matrix the IterativeSolver backend can handle - auto M = MatrixConverter::multiTypeToBCRSMatrix(A); + const auto M = MatrixConverter::multiTypeToBCRSMatrix(A); // get the new matrix sizes const std::size_t numRows = M.N(); assert(numRows == M.M()); // create the vector the IterativeSolver backend can handle - auto bTmp = VectorConverter::multiTypeToBlockVector(b); + const auto bTmp = VectorConverter::multiTypeToBlockVector(b); assert(bTmp.size() == numRows); // create a blockvector to which the linear solver writes the solution -- GitLab From 22934c795f850cb2dc6e5ea7db5e1f40006e9052 Mon Sep 17 00:00:00 2001 From: Bernd Flemisch Date: Mon, 21 Dec 2020 23:22:58 +0100 Subject: [PATCH 17/32] [linear] don't call makeNonOverlappingConsistent for multitype --- dumux/linear/istlsolverfactorybackend.hh | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/dumux/linear/istlsolverfactorybackend.hh b/dumux/linear/istlsolverfactorybackend.hh index f1d87478f9..43536d7230 100644 --- a/dumux/linear/istlsolverfactorybackend.hh +++ b/dumux/linear/istlsolverfactorybackend.hh @@ -44,6 +44,7 @@ #include #include +#include #include #include #include @@ -404,7 +405,7 @@ public: Scalar norm(const Vector& x) const { #if HAVE_MPI - if constexpr (LinearSolverTraits::canCommunicate) + if constexpr (LinearSolverTraits::canCommunicate && !isMultiTypeBlockVector::value) { if (solverCategory_ == Dune::SolverCategory::nonoverlapping) { -- GitLab From 6aaf48c085550d7e309afa41d728b25265012fd4 Mon Sep 17 00:00:00 2001 From: Timo Koch Date: Thu, 31 Dec 2020 19:12:30 +0100 Subject: [PATCH 18/32] [latraits] Always provide single block type alrebra types --- dumux/linear/algebratraits.hh | 42 ++++++++++++++++++++++------------- 1 file changed, 27 insertions(+), 15 deletions(-) diff --git a/dumux/linear/algebratraits.hh b/dumux/linear/algebratraits.hh index f7fe9930d0..7c88f99031 100644 --- a/dumux/linear/algebratraits.hh +++ b/dumux/linear/algebratraits.hh @@ -33,9 +33,7 @@ #include #include -namespace Dumux { - -namespace Detail { +namespace Dumux::Detail { template>::value, int> = 0> constexpr std::size_t blockSize() { return 1; } @@ -43,16 +41,7 @@ constexpr std::size_t blockSize() { return 1; } template>::value, int> = 0> constexpr std::size_t blockSize() { return std::decay_t::size(); } -} // end namespace Detail - -template -struct LinearAlgebraTraits -{ - using Matrix = M; - using Vector = V; -}; - -template +template struct LATraitsFromAssemblerImpl { private: @@ -63,6 +52,8 @@ private: public: using Vector = Dune::BlockVector; using Matrix = typename Assembler::JacobianMatrix; + using SingleTypeVector = Vector; + using SingleTypeMatrix = Matrix; }; template @@ -74,9 +65,30 @@ struct LATraitsFromAssemblerImpl using SingleTypeMatrix = decltype(MatrixConverter::multiTypeToBCRSMatrix(std::declval())); }; +} // end namespace Dumux::Detail + +namespace Dumux { + +/* + * \ingroup Linear + * \brief Traits providing linear algebra types (vector, matrix) + */ +template +struct LinearAlgebraTraits +{ + using Matrix = M; + using Vector = V; + using SingleTypeMatrix = STM; + using SingleTypeVector = STV; +}; + +/* + * \ingroup Linear + * \brief Helper to extract linear algebra types from an assembler + */ template -using LinearAlgebraTraitsFromAssembler = LATraitsFromAssemblerImpl::value>; +using LinearAlgebraTraitsFromAssembler + = Detail::LATraitsFromAssemblerImpl::value>; } // end namespace Dumux -- GitLab From 7493e99e3f52dcef01a1fdef93a90ac5f3c09cf7 Mon Sep 17 00:00:00 2001 From: Timo Koch Date: Thu, 31 Dec 2020 19:21:22 +0100 Subject: [PATCH 19/32] [istlsolvers] Simplify solver alias --- dumux/linear/istlsolvers.hh | 38 ++++++++++++------------------------- 1 file changed, 12 insertions(+), 26 deletions(-) diff --git a/dumux/linear/istlsolvers.hh b/dumux/linear/istlsolvers.hh index 15c6375418..f81f7926de 100644 --- a/dumux/linear/istlsolvers.hh +++ b/dumux/linear/istlsolvers.hh @@ -61,6 +61,7 @@ constexpr std::size_t preconditionerBlockLevel() noexcept } // end namespace Dumux::Detail + namespace Dumux { /*! @@ -315,34 +316,19 @@ private: std::string name_; }; -template -struct ILUBiCGSTABIstlSolverImpl; - -template -struct ILUBiCGSTABIstlSolverImpl -{ - using type = IstlLinearSolver, - Dune::SeqILU - >; -}; - -template -struct ILUBiCGSTABIstlSolverImpl -{ - using type = IstlLinearSolver, - Dune::SeqILU - >; -}; +/*! + * \ingroup Linear + * \brief An ILU preconditioned BiCGSTAB solver using dune-istl + */ template -using ILUBiCGSTABIstlSolver = typename ILUBiCGSTABIstlSolverImpl::value>::type; +using ILUBiCGSTABIstlSolver = + IstlLinearSolver, + Dune::SeqILU + >; } // end namespace Dumux -- GitLab From 9c9e9fd43ad48d2f2c4ca299110338e3dbc9b807 Mon Sep 17 00:00:00 2001 From: Timo Koch Date: Thu, 31 Dec 2020 21:48:22 +0100 Subject: [PATCH 20/32] [temp][istlsolver] Use preconditioner factory instead of specifying preconditioner directly --- dumux/linear/istlsolvers.hh | 44 ++++++++++++++++++++++++++++++++----- 1 file changed, 38 insertions(+), 6 deletions(-) diff --git a/dumux/linear/istlsolvers.hh b/dumux/linear/istlsolvers.hh index f81f7926de..c7cb0c656e 100644 --- a/dumux/linear/istlsolvers.hh +++ b/dumux/linear/istlsolvers.hh @@ -59,6 +59,39 @@ constexpr std::size_t preconditionerBlockLevel() noexcept return isMultiTypeBlockMatrix::value ? 2 : 1; } +template class Preconditioner, int blockLevel = 1> +class IstlDefaultBlockLevelPreconditionerFactory +{ +public: + template + auto operator() (TL typeList, const M& matrix, const Dune::ParameterTree& config) + { + using Matrix = typename Dune::TypeListElement<0, decltype(typeList)>::type; + using Domain = typename Dune::TypeListElement<1, decltype(typeList)>::type; + using Range = typename Dune::TypeListElement<2, decltype(typeList)>::type; + std::shared_ptr> preconditioner + = std::make_shared>(matrix, config); + return preconditioner; + } +}; + +template class Preconditioner> +class IstlDefaultPreconditionerFactory +{ + template + auto operator() (TL typeList, const M& matrix, const Dune::ParameterTree& config) + { + using Matrix = typename Dune::TypeListElement<0, decltype(typeList)>::type; + using Domain = typename Dune::TypeListElement<1, decltype(typeList)>::type; + using Range = typename Dune::TypeListElement<2, decltype(typeList)>::type; + std::shared_ptr> preconditioner + = std::make_shared>(matrix, config); + return preconditioner; + } +}; + +using IstlAmgPreconditionerFactory = Dune::AMGCreator; + } // end namespace Dumux::Detail @@ -69,7 +102,7 @@ namespace Dumux { * \brief Base class for linear solvers */ template + class InverseOperator, class PreconditionerFactory> class IstlLinearSolver { using Matrix = typename LinearAlgebraTraits::Matrix; @@ -77,7 +110,7 @@ class IstlLinearSolver using BVector = typename LinearAlgebraTraits::Vector; using Scalar = typename InverseOperator::real_type; static constexpr bool convertMatrixAndVector - = !preconditionerAcceptsMultiTypeMatrix::value + = !preconditionerAcceptsMultiTypeMatrix::value && isMultiTypeBlockVector::value; #if HAVE_MPI using Comm = Dune::OwnerOverlapCopyCommunication, int>; @@ -295,7 +328,8 @@ private: { const auto& params = params_.sub("preconditioner"); using Prec = Dune::Preconditioner; - std::shared_ptr prec = std::make_shared(op, params); + using TL = Dune::TypeList; + std::shared_ptr prec = PreconditionerFactory{}(TL{}, op, params); #if HAVE_MPI && DUNE_VERSION_GT(DUNE_ISTL,2,7) if (prec->category() != op->category() && prec->category() == Dune::SolverCategory::sequential) @@ -325,9 +359,7 @@ template using ILUBiCGSTABIstlSolver = IstlLinearSolver, - Dune::SeqILU + Detail::IstlDefaultBlockLevelPreconditionerFactory >; } // end namespace Dumux -- GitLab From deebda56ab98769bf3c0998627ccddbc13d00cd0 Mon Sep 17 00:00:00 2001 From: Timo Koch Date: Fri, 1 Jan 2021 22:18:34 +0100 Subject: [PATCH 21/32] [amg] Use istl solver class to implement amg backend --- dumux/linear/amgbackend.hh | 212 +++---------------------------------- 1 file changed, 13 insertions(+), 199 deletions(-) diff --git a/dumux/linear/amgbackend.hh b/dumux/linear/amgbackend.hh index 14285f3c9d..398bef7f54 100644 --- a/dumux/linear/amgbackend.hh +++ b/dumux/linear/amgbackend.hh @@ -39,6 +39,7 @@ #include #include #include +#include namespace Dumux { @@ -236,209 +237,22 @@ private: bool isParallel_ = false; }; -template -class NewAMGBiCGSTABBackend : public LinearSolver -{ - using Matrix = typename LinearAlgebraTraits::Matrix; - using Vector = typename LinearAlgebraTraits::Vector; - using ScalarProduct = Dune::ScalarProduct; -#if HAVE_MPI - using Comm = Dune::OwnerOverlapCopyCommunication, 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(); - } - - /*! - * \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(gridView); - if (solverCategory_ != Dune::SolverCategory::sequential) - { - phelper_ = std::make_unique>(gridView, dofMapper); - comm_ = std::make_shared(gridView.comm(), solverCategory_); - scalarProduct_ = Dune::createScalarProduct(*comm_, solverCategory_); - phelper_->createParallelIndexSet(*comm_); - } - else - scalarProduct_ = std::make_shared(); -#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 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; - solveParallel_(A, x, b); - } - else - { - using PTraits = typename LinearSolverTraits::template ParallelOverlapping; - solveParallel_(A, x, b); - } - } - else - solveSequential_(A, x, b); - } - else - { - solveSequential_(A, x, b); - } - } - - template - void solveParallel_(Matrix& A, Vector& x, Vector& b) - { - prepareLinearAlgebraParallel(A, b, *phelper_); - auto linearOperator = std::make_shared(A, *comm_); - - using SeqSmoother = Dune::SeqSSOR; - using Smoother = typename ParallelTraits::template Preconditioner; - solveWithAmg_(A, x, b, linearOperator, comm_); - } -#endif // HAVE_MPI - - void solveSequential_(Matrix& A, Vector& x, Vector& b) - { - using SequentialTraits = typename LinearSolverTraits::template Sequential; - auto linearOperator = std::make_shared(A); - - auto comm = std::make_shared(); - - using Smoother = Dune::SeqSSOR; - solveWithAmg_(A, x, b, linearOperator, comm); - } - - template - void solveWithAmg_(Matrix& A, Vector& x, Vector& b, - std::shared_ptr& linearOperator, - std::shared_ptr& comm) - { - using SmootherArgs = typename Dune::Amg::SmootherTraits::Arguments; - using Criterion = Dune::Amg::CoarsenCriterion>; - - //! \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; - auto amg = std::make_shared(*linearOperator, criterion, smootherArgs, *comm); - - Dune::BiCGSTABSolver 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> phelper_; - std::shared_ptr comm_; -#endif - Dune::SolverCategory::Category solverCategory_; - std::shared_ptr scalarProduct_; - Dune::InverseOperatorResult result_; - bool isParallel_ = false; -}; +/*! + * \ingroup Linear + * \brief An AMG preconditioned BiCGSTAB solver using dune-istl + */ +template +using NewAMGBiCGSTABBackend = + IstlLinearSolver, + Detail::IstlAmgPreconditionerFactory + >; template struct AMGImplHelper { - template - using type = NewAMGBiCGSTABBackend; + template + using type = NewAMGBiCGSTABBackend; }; template<> -- GitLab From a10dc4ca8287ddc1de0274e30bac061d09468383 Mon Sep 17 00:00:00 2001 From: Timo Koch Date: Sat, 2 Jan 2021 13:21:41 +0100 Subject: [PATCH 22/32] [istlsolvers] Simplify new AMG solver impl by renaming class to AMGBiCGSTABIstlSolver --- dumux/linear/amgbackend.hh | 55 ++++----------------- dumux/linear/istlsolvers.hh | 15 +++++- test/porousmediumflow/2p2c/waterair/main.cc | 6 +-- 3 files changed, 26 insertions(+), 50 deletions(-) diff --git a/dumux/linear/amgbackend.hh b/dumux/linear/amgbackend.hh index 398bef7f54..66ccb7a1dd 100644 --- a/dumux/linear/amgbackend.hh +++ b/dumux/linear/amgbackend.hh @@ -25,6 +25,8 @@ #ifndef DUMUX_PARALLEL_AMGBACKEND_HH #define DUMUX_PARALLEL_AMGBACKEND_HH +#warning "This header is deprecated. Use the AMG solver from dumux/linear/istlsolvers.hh" + #include #include @@ -39,14 +41,12 @@ #include #include #include -#include namespace Dumux { -namespace Detail { - +// OLD AMG Backend. Use the new one from dumux/linear/istlsolvers.hh template -class OldAMGBiCGSTABBackend : public LinearSolver +class AMGBiCGSTABBackend : public LinearSolver { public: /*! @@ -54,8 +54,8 @@ public: * * \param paramGroup the parameter group for parameter lookup */ - [[deprecated("Use new AMGBiCGSTABBackend with 2nd template parameter.")]] - OldAMGBiCGSTABBackend(const std::string& paramGroup = "") + [[deprecated("Use new AMGBiCGSTABIstlSolver with 2nd template parameter from dumux/linear/istlsolvers.hh.")]] + AMGBiCGSTABBackend(const std::string& paramGroup = "") : LinearSolver(paramGroup) , isParallel_(Dune::MPIHelper::getCommunication().size() > 1) { @@ -72,10 +72,10 @@ public: * \param dofMapper an index mapper for dof entities * \param paramGroup the parameter group for parameter lookup */ - [[deprecated("Use new AMGBiCGSTABBackend with 2nd template parameter.")]] - OldAMGBiCGSTABBackend(const typename LinearSolverTraits::GridView& gridView, - const typename LinearSolverTraits::DofMapper& dofMapper, - const std::string& paramGroup = "") + [[deprecated("Use new AMGBiCGSTABIstlSolver with 2nd template parameter from dumux/linear/istlsolvers.hh.")]] + AMGBiCGSTABBackend(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) @@ -237,41 +237,6 @@ private: bool isParallel_ = false; }; -/*! - * \ingroup Linear - * \brief An AMG preconditioned BiCGSTAB solver using dune-istl - */ -template -using NewAMGBiCGSTABBackend = - IstlLinearSolver, - Detail::IstlAmgPreconditionerFactory - >; - -template -struct AMGImplHelper -{ - template - using type = NewAMGBiCGSTABBackend; -}; - -template<> -struct AMGImplHelper<1> -{ - template - using type = OldAMGBiCGSTABBackend; -}; - -} // end namespace Detail - -/*! - * \ingroup Linear - * \brief A linear solver based on the ISTL AMG preconditioner - * and the ISTL BiCGSTAB solver. - */ -template -using AMGBiCGSTABBackend = typename Detail::AMGImplHelper::template type; - } // end namespace Dumux #endif diff --git a/dumux/linear/istlsolvers.hh b/dumux/linear/istlsolvers.hh index c7cb0c656e..4af95121d3 100644 --- a/dumux/linear/istlsolvers.hh +++ b/dumux/linear/istlsolvers.hh @@ -31,11 +31,11 @@ #include #include #include -#include -#include #include #include #include +#include +#include #include #include @@ -362,6 +362,17 @@ using ILUBiCGSTABIstlSolver = Detail::IstlDefaultBlockLevelPreconditionerFactory >; +/*! + * \ingroup Linear + * \brief An AMG preconditioned BiCGSTAB solver using dune-istl + */ +template +using AMGBiCGSTABIstlSolver = + IstlIterativeLinearSolver, + Detail::IstlAmgPreconditionerFactory + >; + } // end namespace Dumux #endif diff --git a/test/porousmediumflow/2p2c/waterair/main.cc b/test/porousmediumflow/2p2c/waterair/main.cc index 74830f3ed8..836e959f13 100644 --- a/test/porousmediumflow/2p2c/waterair/main.cc +++ b/test/porousmediumflow/2p2c/waterair/main.cc @@ -32,7 +32,7 @@ #include #include -#include +#include #include #include #include @@ -112,8 +112,8 @@ int main(int argc, char** argv) auto assembler = std::make_shared(problem, gridGeometry, gridVariables, timeLoop, xOld); // the linear solver - using LinearSolver = AMGBiCGSTABBackend, - LinearAlgebraTraitsFromAssembler>; + using LinearSolver = AMGBiCGSTABIstlSolver, + LinearAlgebraTraitsFromAssembler>; auto linearSolver = std::make_shared(leafGridView, gridGeometry->dofMapper()); // the non-linear solver -- GitLab From e26f416ab5ef5436675f564cd4a714d5fb37b806 Mon Sep 17 00:00:00 2001 From: Timo Koch Date: Sat, 2 Jan 2021 13:23:24 +0100 Subject: [PATCH 23/32] [istlsolvers] Specify acceptsMultiType logic explicitly for the specfic solvers --- dumux/linear/istlsolvers.hh | 44 ++++++++++++------- .../linearsolveracceptsmultitypematrix.hh | 8 ---- 2 files changed, 28 insertions(+), 24 deletions(-) diff --git a/dumux/linear/istlsolvers.hh b/dumux/linear/istlsolvers.hh index 4af95121d3..83aee155e4 100644 --- a/dumux/linear/istlsolvers.hh +++ b/dumux/linear/istlsolvers.hh @@ -40,7 +40,6 @@ #include #include #include -#include #include #include #include @@ -99,19 +98,19 @@ namespace Dumux { /*! * \ingroup Linear - * \brief Base class for linear solvers + * \brief Standard dune-istl iterative linear solvers */ template -class IstlLinearSolver + class InverseOperator, class PreconditionerFactory, + bool convertMultiTypeLATypes = false> +class IstlIterativeLinearSolver { using Matrix = typename LinearAlgebraTraits::Matrix; using XVector = typename LinearAlgebraTraits::Vector; using BVector = typename LinearAlgebraTraits::Vector; using Scalar = typename InverseOperator::real_type; - static constexpr bool convertMatrixAndVector - = !preconditionerAcceptsMultiTypeMatrix::value - && isMultiTypeBlockVector::value; + static constexpr bool convertMultiTypeVectorAndMatrix + = convertMultiTypeLATypes && isMultiTypeBlockVector::value; #if HAVE_MPI using Comm = Dune::OwnerOverlapCopyCommunication, int>; using ScalarProduct = Dune::ScalarProduct; @@ -122,7 +121,7 @@ public: /*! * \brief Constructor for sequential solvers */ - IstlLinearSolver(const std::string& paramGroup = "") + IstlIterativeLinearSolver(const std::string& paramGroup = "") { if (Dune::MPIHelper::getCollectiveCommunication().size() > 1) DUNE_THROW(Dune::InvalidStateException, "Using sequential constructor for parallel run. Use signature with gridView and dofMapper!"); @@ -136,9 +135,9 @@ public: * \brief Constructor for parallel and sequential solvers */ template - IstlLinearSolver(const GridView& gridView, - const DofMapper& dofMapper, - const std::string& paramGroup = "") + IstlIterativeLinearSolver(const GridView& gridView, + const DofMapper& dofMapper, + const std::string& paramGroup = "") { initializeParameters_(paramGroup); #if HAVE_MPI @@ -162,7 +161,7 @@ public: * \brief Constructor with custom scalar product and communication */ template - IstlLinearSolver(std::shared_ptr communication, + IstlIterativeLinearSolver(std::shared_ptr communication, std::shared_ptr scalarProduct, const GridView& gridView, const DofMapper& dofMapper, @@ -205,7 +204,7 @@ public: } } #endif - if constexpr (convertMatrixAndVector) + if constexpr (convertMultiTypeVectorAndMatrix) { auto y = VectorConverter::multiTypeToBlockVector(x); return scalarProduct_->norm(y); @@ -236,7 +235,7 @@ private: bool solveSequential_(Matrix& A, XVector& x, BVector& b) { - if constexpr (convertMatrixAndVector) + if constexpr (convertMultiTypeVectorAndMatrix) { // create the bcrs matrix the IterativeSolver backend can handle auto M = MatrixConverter::multiTypeToBCRSMatrix(A); @@ -354,12 +353,25 @@ private: /*! * \ingroup Linear * \brief An ILU preconditioned BiCGSTAB solver using dune-istl + * + * Solver: The BiCGSTAB (stabilized biconjugate gradients method) solver has + * faster and smoother convergence than the original BiCG. It can be applied to + * nonsymmetric matrices.\n + * See: Van der Vorst, H. A. (1992). "Bi-CGSTAB: A Fast and Smoothly Converging + * Variant of Bi-CG for the Solution of Nonsymmetric Linear Systems". + * SIAM J. Sci. and Stat. Comput. 13 (2): 631–644. doi:10.1137/0913035. + * + * Preconditioner: ILU(n) incomplete LU factorization. The order n indicates + * fill-in. It can be damped by the relaxation parameter + * LinearSolver.PreconditionerRelaxation.\n + * See: Golub, G. H., and Van Loan, C. F. (2012). Matrix computations. JHU Press. */ template using ILUBiCGSTABIstlSolver = - IstlLinearSolver, - Detail::IstlDefaultBlockLevelPreconditionerFactory + Detail::IstlDefaultBlockLevelPreconditionerFactory, + /*accepts multi-type istl types?*/ false >; /*! diff --git a/dumux/linear/linearsolveracceptsmultitypematrix.hh b/dumux/linear/linearsolveracceptsmultitypematrix.hh index abfd44b623..de471c0639 100644 --- a/dumux/linear/linearsolveracceptsmultitypematrix.hh +++ b/dumux/linear/linearsolveracceptsmultitypematrix.hh @@ -24,8 +24,6 @@ #ifndef DUMUX_LINEAR_SOLVER_ACCEPTS_MULTITYPEMATRIX_HH #define DUMUX_LINEAR_SOLVER_ACCEPTS_MULTITYPEMATRIX_HH -#include - #include namespace Dumux { @@ -66,12 +64,6 @@ template<> struct linearSolverAcceptsMultiTypeMatrix : public std::false_type {}; #endif // HAVE_UMFPACK -template -struct preconditionerAcceptsMultiTypeMatrix : public std::true_type {}; - -template -struct preconditionerAcceptsMultiTypeMatrix> : public std::false_type {}; - } // end namespace Dumux #endif -- GitLab From 1e1794edb3c6df349e449f3ec301cad4f67319ca Mon Sep 17 00:00:00 2001 From: Timo Koch Date: Sat, 2 Jan 2021 13:23:40 +0100 Subject: [PATCH 24/32] [istlsolvers] Add new direct solvers accepting multitypematrices and state vectors --- dumux/linear/istlsolvers.hh | 152 ++++++++++++++++++++++++++++++++++++ 1 file changed, 152 insertions(+) diff --git a/dumux/linear/istlsolvers.hh b/dumux/linear/istlsolvers.hh index 83aee155e4..2db84bcbba 100644 --- a/dumux/linear/istlsolvers.hh +++ b/dumux/linear/istlsolvers.hh @@ -39,10 +39,12 @@ #include #include +#include #include #include #include #include +#include namespace Dumux::Detail { @@ -385,6 +387,156 @@ using AMGBiCGSTABIstlSolver = Detail::IstlAmgPreconditionerFactory >; + +/*! + * \ingroup Linear + * \brief Direct dune-istl linear solvers + */ +template class Solver> +class DirectIstlSolver : public LinearSolver +{ +public: + using LinearSolver::LinearSolver; + + template + bool solve(const Matrix& A, Vector& x, const Vector& b) + { + // support dune-istl multi-type block vector/matrix + if constexpr (isMultiTypeBlockVector()) + { + auto [AA, xx, bb] = convertIstlMultiTypeToBCRSSystem_(A, x, b); + bool converged = solve_(AA, xx, bb); + if (converged) + VectorConverter::retrieveValues(x, xx); + return converged; + } + + /* + Copy vectors into standard istl block vector. + This is necessary for all model _not_ using a FieldVector + Could be avoided for vectors that already have the right type using another if constexpr branch + but it shouldn't impact performance too much + */ + else + { + constexpr auto blockSize = Detail::blockSize(); + using Scalar = std::decay_t; + using BlockType = Dune::FieldVector; + Dune::BlockVector xTmp; xTmp.resize(b.size()); + Dune::BlockVector bTmp(xTmp); + + bTmp = b; + bool converged = solve_(A, xTmp, bTmp); + if (converged) + x = xTmp; + return converged; + } + } + + std::string name() const + { + return "Direct solver"; + } + + const Dune::InverseOperatorResult& result() const + { + return result_; + } + +private: + Dune::InverseOperatorResult result_; + + template + bool solve_(const Matrix& A, Vector& x, const Vector& b) + { + static_assert(isBCRSMatrix::value, "Direct solver only works with BCRS matrices!"); + using BlockType = typename Matrix::block_type; + static_assert(BlockType::rows == BlockType::cols, "Matrix block must be quadratic!"); + constexpr auto blockSize = BlockType::rows; + + Solver solver(A, this->verbosity() > 0); + + Vector bTmp(b); + solver.apply(x, bTmp, result_); + + int size = x.size(); + for (int i = 0; i < size; i++) + { + for (int j = 0; j < blockSize; j++) + { + using std::isnan; + using std::isinf; + if (isnan(x[i][j]) || isinf(x[i][j])) + { + result_.converged = false; + break; + } + } + } + + return result_.converged; + } + + template + auto convertIstlMultiTypeToBCRSSystem_(const Matrix& A, Vector& x, const Vector& b) + { + const auto AA = MatrixConverter::multiTypeToBCRSMatrix(A); + + // get the new matrix sizes + const std::size_t numRows = AA.N(); + assert(numRows == AA.M()); + + // create the vector the IterativeSolver backend can handle + const auto bb = VectorConverter::multiTypeToBlockVector(b); + assert(bb.size() == numRows); + + // create a blockvector to which the linear solver writes the solution + using VectorBlock = typename Dune::FieldVector; + using BlockVector = typename Dune::BlockVector; + BlockVector xx(numRows); + + return std::make_tuple(std::move(AA), std::move(xx), std::move(bb)); + } +}; + +} // end namespace Dumux + +#if HAVE_SUPERLU +#include + +namespace Dumux { + +/*! + * \ingroup Linear + * \brief Direct linear solver using the SuperLU library. + * + * See: Li, X. S. (2005). "An overview of SuperLU: Algorithms, implementation, + * and user interface." ACM Transactions on Mathematical Software (TOMS) 31(3): 302-325. + * http://crd-legacy.lbl.gov/~xiaoye/SuperLU/ + */ +using SuperLUIstlSolver = DirectIstlSolver; + } // end namespace Dumux +#endif // HAVE_SUPERLU + +#if HAVE_UMFPACK +#include + +namespace Dumux { + +/*! + * \ingroup Linear + * \brief Direct linear solver using the UMFPack library. + * + * See: Davis, Timothy A. (2004). "Algorithm 832". ACM Transactions on + * Mathematical Software 30 (2): 196–199. doi:10.1145/992200.992206. + * http://faculty.cse.tamu.edu/davis/suitesparse.html + */ +using UMFPackIstlSolver = DirectIstlSolver; + +} // end namespace Dumux + +#endif // HAVE_UMFPACK + #endif -- GitLab From 6c63bb6dd7bfa2c227ff8343d668419ee1357347 Mon Sep 17 00:00:00 2001 From: Timo Koch Date: Sat, 2 Jan 2021 13:23:46 +0100 Subject: [PATCH 25/32] [test][f] Use new UMFPack solver in donea test --- test/freeflow/navierstokes/donea/main.cc | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/test/freeflow/navierstokes/donea/main.cc b/test/freeflow/navierstokes/donea/main.cc index 7986ac5af9..80bd1a3561 100644 --- a/test/freeflow/navierstokes/donea/main.cc +++ b/test/freeflow/navierstokes/donea/main.cc @@ -36,6 +36,7 @@ #include #include #include +#include #include #include @@ -137,7 +138,7 @@ int main(int argc, char** argv) couplingManager); // the linear solver - using LinearSolver = Dumux::UMFPackBackend; + using LinearSolver = UMFPackIstlSolver; auto linearSolver = std::make_shared(); // the non-linear solver -- GitLab From 986f3c7c0675f2cbd57304ab9008a63fead8c26e Mon Sep 17 00:00:00 2001 From: Timo Koch Date: Sat, 2 Jan 2021 13:31:42 +0100 Subject: [PATCH 26/32] [linear] Deprecate obsolete traits header --- dumux/linear/linearsolveracceptsmultitypematrix.hh | 2 ++ 1 file changed, 2 insertions(+) diff --git a/dumux/linear/linearsolveracceptsmultitypematrix.hh b/dumux/linear/linearsolveracceptsmultitypematrix.hh index de471c0639..0c3e25b258 100644 --- a/dumux/linear/linearsolveracceptsmultitypematrix.hh +++ b/dumux/linear/linearsolveracceptsmultitypematrix.hh @@ -24,6 +24,8 @@ #ifndef DUMUX_LINEAR_SOLVER_ACCEPTS_MULTITYPEMATRIX_HH #define DUMUX_LINEAR_SOLVER_ACCEPTS_MULTITYPEMATRIX_HH +#warning "This header is deprecated and will be removed after release 3.4." + #include namespace Dumux { -- GitLab From 98331c559681bc2dec66f9170380cdb0ee2dc45b Mon Sep 17 00:00:00 2001 From: Timo Koch Date: Sat, 2 Jan 2021 13:52:06 +0100 Subject: [PATCH 27/32] [newton] Deprecate multitype conversion in Newton solver --- dumux/nonlinear/newtonsolver.hh | 1 + 1 file changed, 1 insertion(+) diff --git a/dumux/nonlinear/newtonsolver.hh b/dumux/nonlinear/newtonsolver.hh index d78141353f..ad831e67fc 100644 --- a/dumux/nonlinear/newtonsolver.hh +++ b/dumux/nonlinear/newtonsolver.hh @@ -1213,6 +1213,7 @@ private: * */ template + [[deprecated("After 3.4 Newton will no longer support conversion of multitype matrices for solvers that don't support this feature!")]] typename std::enable_if_t() && isMultiTypeBlockVector(), bool> solveLinearSystemImpl_(LinearSolver& ls, -- GitLab From 0d214d220d643ffd50aa7dfcf93f593411e3cbff Mon Sep 17 00:00:00 2001 From: Timo Koch Date: Sat, 2 Jan 2021 17:20:44 +0100 Subject: [PATCH 28/32] [linear] Silence deprecation warning by using forward declarations --- .../linear/linearsolveracceptsmultitypematrix.hh | 16 ++++++++++++++-- 1 file changed, 14 insertions(+), 2 deletions(-) diff --git a/dumux/linear/linearsolveracceptsmultitypematrix.hh b/dumux/linear/linearsolveracceptsmultitypematrix.hh index 0c3e25b258..6fc59812a3 100644 --- a/dumux/linear/linearsolveracceptsmultitypematrix.hh +++ b/dumux/linear/linearsolveracceptsmultitypematrix.hh @@ -26,8 +26,6 @@ #warning "This header is deprecated and will be removed after release 3.4." -#include - namespace Dumux { //! The default @@ -38,30 +36,44 @@ struct linearSolverAcceptsMultiTypeMatrix : public std::true_type {}; //! Those are all with ILU preconditioner that doesn't support the additional block level //! And the direct solvers that have BCRS Matrix hardcoded +class ILUnBiCGSTABBackend; + template<> struct linearSolverAcceptsMultiTypeMatrix : public std::false_type {}; +class ILUnCGBackend; + template<> struct linearSolverAcceptsMultiTypeMatrix : public std::false_type {}; +class ILU0BiCGSTABBackend; + template<> struct linearSolverAcceptsMultiTypeMatrix : public std::false_type {}; +class ILU0CGBackend; + template<> struct linearSolverAcceptsMultiTypeMatrix : public std::false_type {}; +class ILU0RestartedGMResBackend; + template<> struct linearSolverAcceptsMultiTypeMatrix : public std::false_type {}; +class ILUnRestartedGMResBackend; + template<> struct linearSolverAcceptsMultiTypeMatrix : public std::false_type {}; #if HAVE_SUPERLU +class SuperLUBackend; template<> struct linearSolverAcceptsMultiTypeMatrix : public std::false_type {}; #endif // HAVE_SUPERLU #if HAVE_UMFPACK +class UMFPackBackend; template<> struct linearSolverAcceptsMultiTypeMatrix : public std::false_type {}; #endif // HAVE_UMFPACK -- GitLab From 97419c6f14a7b4990dca5a9dfa1f1f15483e9970 Mon Sep 17 00:00:00 2001 From: Timo Koch Date: Sat, 2 Jan 2021 17:21:17 +0100 Subject: [PATCH 29/32] [linear][istl] Add AMG-CG solver --- dumux/linear/istlsolvers.hh | 11 +++++++++++ test/porousmediumflow/1p/incompressible/main.cc | 7 +++++-- 2 files changed, 16 insertions(+), 2 deletions(-) diff --git a/dumux/linear/istlsolvers.hh b/dumux/linear/istlsolvers.hh index 2db84bcbba..b2c13d1a6d 100644 --- a/dumux/linear/istlsolvers.hh +++ b/dumux/linear/istlsolvers.hh @@ -387,6 +387,17 @@ using AMGBiCGSTABIstlSolver = Detail::IstlAmgPreconditionerFactory >; +/*! + * \ingroup Linear + * \brief An AMG preconditioned CG solver using dune-istl + */ +template +using AMGCGIstlSolver = + IstlIterativeLinearSolver, + Detail::IstlAmgPreconditionerFactory + >; + /*! * \ingroup Linear diff --git a/test/porousmediumflow/1p/incompressible/main.cc b/test/porousmediumflow/1p/incompressible/main.cc index 2bf4beeae8..379c40b504 100644 --- a/test/porousmediumflow/1p/incompressible/main.cc +++ b/test/porousmediumflow/1p/incompressible/main.cc @@ -32,7 +32,9 @@ #include #include -#include +#include +#include +#include #include #include @@ -143,7 +145,8 @@ int main(int argc, char** argv) using Assembler = FVAssembler; auto assembler = std::make_shared(problem, gridGeometry, gridVariables); - using LinearSolver = SSORCGBackend; + using LinearSolver = AMGCGIstlSolver, + LinearAlgebraTraitsFromAssembler>; auto linearSolver = std::make_shared(); // solver the linear problem -- GitLab From e526764bce22f8327ba2dfb702a4ebd00c9b273a Mon Sep 17 00:00:00 2001 From: Bernd Flemisch Date: Fri, 28 Jan 2022 09:10:53 +0100 Subject: [PATCH 30/32] [nonlinear] fix newton solver --- dumux/nonlinear/newtonsolver.hh | 6 ------ 1 file changed, 6 deletions(-) diff --git a/dumux/nonlinear/newtonsolver.hh b/dumux/nonlinear/newtonsolver.hh index ad831e67fc..5f1745d9ca 100644 --- a/dumux/nonlinear/newtonsolver.hh +++ b/dumux/nonlinear/newtonsolver.hh @@ -182,12 +182,6 @@ void assign(To& to, const From& from) DUNE_THROW(Dune::Exception, "Values are not assignable to each other!"); } -template>::value, int> = 0> -constexpr std::size_t blockSize() { return 1; } - -template>::value, int> = 0> -constexpr std::size_t blockSize() { return std::decay_t::size(); } - template struct BlockTypeHelper { using type = std::decay_t()[0])>; }; -- GitLab From c634d4f31a7cea5c2f3efdf0d5375997abbe932d Mon Sep 17 00:00:00 2001 From: Bernd Flemisch Date: Fri, 28 Jan 2022 09:11:38 +0100 Subject: [PATCH 31/32] [test][navierstokes][donea] use new istl solver factory backend --- test/freeflow/navierstokes/donea/main_momentum.cc | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/test/freeflow/navierstokes/donea/main_momentum.cc b/test/freeflow/navierstokes/donea/main_momentum.cc index c9248abf19..14c01c0b47 100644 --- a/test/freeflow/navierstokes/donea/main_momentum.cc +++ b/test/freeflow/navierstokes/donea/main_momentum.cc @@ -42,6 +42,7 @@ #include #include +#include #include #include "properties_momentum.hh" @@ -153,7 +154,8 @@ int main(int argc, char** argv) using Assembler = FVAssembler; auto assembler = std::make_shared(problem, gridGeometry, gridVariables); - using LinearSolver = IstlSolverFactoryBackend>; + using LinearSolver = IstlSolverFactoryBackend, + LinearAlgebraTraitsFromAssembler>; const auto dofMapper = LinearSolverTraits::DofMapper(gridGeometry->gridView()); auto linearSolver = std::make_shared(gridGeometry->gridView(), dofMapper); -- GitLab From c8646266cac0d955f233d5d48365ceedec4c9ee8 Mon Sep 17 00:00:00 2001 From: Bernd Flemisch Date: Fri, 28 Jan 2022 09:38:32 +0100 Subject: [PATCH 32/32] [linear][istl] fix acceptance of multitype block vectors/matrices --- dumux/linear/istlsolvers.hh | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/dumux/linear/istlsolvers.hh b/dumux/linear/istlsolvers.hh index b2c13d1a6d..e1781b934b 100644 --- a/dumux/linear/istlsolvers.hh +++ b/dumux/linear/istlsolvers.hh @@ -373,7 +373,7 @@ using ILUBiCGSTABIstlSolver = IstlIterativeLinearSolver, Detail::IstlDefaultBlockLevelPreconditionerFactory, - /*accepts multi-type istl types?*/ false + /*accepts multi-type istl types?*/ true >; /*! -- GitLab