From eb38c414ab3e17b33bfc5ba685c9a5a925f4b2d0 Mon Sep 17 00:00:00 2001
From: Kilian Weishaupt <kilian.weishaupt@iws.uni-stuttgart.de>
Date: Tue, 17 Mar 2020 12:15:28 +0100
Subject: [PATCH] [istlsolverfactorybackend] Introduce init-function for
dumux-specific preconditioners
* check matrix type, chose appropriate function (dune or dumux)
* MultiTypeBlock matrix only works sequnetially so far
---
dumux/linear/istlsolverfactorybackend.hh | 56 ++++++++++++++++++++++--
1 file changed, 53 insertions(+), 3 deletions(-)
diff --git a/dumux/linear/istlsolverfactorybackend.hh b/dumux/linear/istlsolverfactorybackend.hh
index 776ba71309..2de5db9694 100644
--- a/dumux/linear/istlsolverfactorybackend.hh
+++ b/dumux/linear/istlsolverfactorybackend.hh
@@ -37,16 +37,63 @@
// preconditioners
#include <dune/istl/preconditioners.hh>
#include <dune/istl/paamg/amg.hh>
+#include "preconditioners.hh"
// solvers
#include <dune/istl/solvers.hh>
#include <dune/istl/solverfactory.hh>
+#include <dumux/common/typetraits/matrix.hh>
#include <dumux/linear/solver.hh>
#include <dumux/linear/parallelhelpers.hh>
namespace Dumux {
+// initSolverFactoriesForMultiTypeBlockMatrix differs in different compilation units,
+// so we have it in an anonymous namespace
+namespace {
+
+/*!
+ * \brief Initializes the direct solvers, preconditioners and iterative solvers in
+ * the factories with the corresponding Matrix and Vector types.
+ * \note We currently consider only direct solvers and preconditioners provided by
+ * Dumux which, unlike the ones implemented in Dune, also support MultiTypeBlockMatrices.
+ * \note This function could be removed once Dune::initSolverFactories supports MultiTypeBlockMatrices.
+ * \tparam LinearOperator the linear operator type
+ */
+template<class LinearOperator>
+int initSolverFactoriesForMultiTypeBlockMatrix()
+{
+ using M = typename LinearOperator::matrix_type;
+ using X = typename LinearOperator::range_type;
+ using Y = typename LinearOperator::domain_type;
+ using TL = Dune::TypeList<M,X,Y>;
+ auto& dsfac = Dune::DirectSolverFactory<M,X,Y>::instance();
+ Dune::addRegistryToFactory<TL>(dsfac, Dumux::MultiTypeBlockMatrixDirectSolverTag{});
+ auto& pfac = Dune::PreconditionerFactory<LinearOperator,X,Y>::instance();
+ Dune::addRegistryToFactory<TL>(pfac, Dumux::MultiTypeBlockMatrixPreconditionerTag{});
+ using TLS = Dune::TypeList<X,Y>;
+ auto& isfac = Dune::IterativeSolverFactory<X,Y>::instance();
+ return Dune::addRegistryToFactory<TLS>(isfac, Dune::IterativeSolverTag{});
+}
+} // end namespace
+
+/*!
+ * \brief Initialize the solver factories for regular matrices or MultiTypeBlockMatrices
+ * \tparam Matrix the matrix
+ * \tparam LinearOperator the linear operator
+ *
+ * \note This function could be removed once Dune::initSolverFactories supports MultiTypeBlockMatrices.
+ */
+template<class Matrix, class LinearOperator>
+void initSolverFactories()
+{
+ if constexpr (isMultiTypeBlockMatrix<Matrix>::value)
+ initSolverFactoriesForMultiTypeBlockMatrix<LinearOperator>();
+ else
+ Dune::initSolverFactories<LinearOperator>();
+}
+
/*!
* \ingroup Linear
* \brief A linear solver using the dune-istl solver factory
@@ -148,7 +195,10 @@ private:
template<class Matrix, class Vector>
void solveSequentialOrParallel_(Matrix& A, Vector& x, Vector& b)
{
- if constexpr (LinearSolverTraits::canCommunicate)
+ // Dune::MultiTypeBlockMatrix does not provide a ColIterator which is needed by Dune::NonoverlappingSchwarzOperator.
+ // We therefore can only solve these types of systems sequentially.
+ // TODO: This can be adapted once the situation in Dune ISTL changes.
+ if constexpr (LinearSolverTraits::canCommunicate && !isMultiTypeBlockMatrix<Matrix>::value)
{
if (isParallel_)
{
@@ -181,7 +231,7 @@ private:
if (firstCall_)
{
- Dune::initSolverFactories<LinearOperator>();
+ initSolverFactories<Matrix, LinearOperator>();
parallelHelper_->initGhostsAndOwners();
}
@@ -207,7 +257,7 @@ private:
auto linearOperator = std::make_shared<LinearOperator>(A);
if (firstCall_)
- Dune::initSolverFactories<LinearOperator>();
+ initSolverFactories<Matrix, LinearOperator>();
// construct solver
auto solver = getSolverFromFactory_(linearOperator);
--
GitLab