diff --git a/dumux/linear/parallelmatrixadapter.hh b/dumux/linear/parallelmatrixadapter.hh
new file mode 100644
index 0000000000000000000000000000000000000000..4394c29ec28ef3cb931638cc0d63d570ada8d48e
--- /dev/null
+++ b/dumux/linear/parallelmatrixadapter.hh
@@ -0,0 +1,118 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*****************************************************************************
+ * See the file COPYING for full copying permissions. *
+ * *
+ * This program is free software: you can redistribute it and/or modify *
+ * it under the terms of the GNU General Public License as published by *
+ * the Free Software Foundation, either version 3 of the License, or *
+ * (at your option) any later version. *
+ * *
+ * This program is distributed in the hope that it will be useful, *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of *
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
+ * GNU General Public License for more details. *
+ * *
+ * You should have received a copy of the GNU General Public License *
+ * along with this program. If not, see <http://www.gnu.org/licenses/>. *
+ *****************************************************************************/
+/*!
+ * \file
+ * \ingroup Linear
+ * \brief A parallel version of a linear operator
+ */
+#ifndef DUMUX_LINEAR_PARALLEL_MATRIX_ADAPTER_HH
+#define DUMUX_LINEAR_PARALLEL_MATRIX_ADAPTER_HH
+
+#include <dune/common/hybridutilities.hh>
+#include <dune/istl/operators.hh>
+#include <dune/istl/bcrsmatrix.hh>
+#include <dumux/parallel/parallel_for.hh>
+
+namespace Dumux {
+
+/*!
+ * \brief Adapter to turn a multi-type matrix into a thread-parallel linear operator.
+ * Adapts a matrix to the assembled linear operator interface
+ */
+template<class M, class X, class Y>
+class ParallelMultiTypeMatrixAdapter : public Dune::MatrixAdapter<M,X,Y>
+{
+ using ParentType = Dune::MatrixAdapter<M,X,Y>;
+public:
+ //! export types
+ typedef M matrix_type;
+ typedef X domain_type;
+ typedef Y range_type;
+ typedef typename X::field_type field_type;
+
+ //! constructor: just store a reference to a matrix
+ explicit ParallelMultiTypeMatrixAdapter (const M& A)
+ : ParentType(A) {}
+
+ //! constructor: store an std::shared_ptr to a matrix
+ explicit ParallelMultiTypeMatrixAdapter (std::shared_ptr<const M> A)
+ : ParentType(A) {}
+
+ //! apply operator to x: \f$ y = A(x) \f$
+ void apply (const X& x, Y& y) const override
+ {
+ y = 0.0;
+
+ auto& A = this->getmat();
+ using namespace Dune::Hybrid;
+ forEach(integralRange(Dune::Hybrid::size(y)), [&](auto&& i)
+ {
+ forEach(integralRange(Dune::Hybrid::size(x)), [&](auto&& j)
+ {
+ const auto& mat = A[i][j];
+ const auto& xx = x[j];
+ auto& yy = y[i];
+
+ Dumux::parallelFor(mat.N(), [&](const std::size_t ii)
+ {
+ const auto& row = mat[ii];
+ const auto endj = row.end();
+ for (auto jj=row.begin(); jj!=endj; ++jj)
+ {
+ const auto& xj = Dune::Impl::asVector(xx[jj.index()]);
+ auto&& yi = Dune::Impl::asVector(yy[ii]);
+ Dune::Impl::asMatrix(*jj).umv(xj, yi);
+ }
+ });
+ });
+ });
+ }
+
+ //! apply operator to x, scale and add: \f$ y = y + \alpha A(x) \f$
+ void applyscaleadd (field_type alpha, const X& x, Y& y) const override
+ {
+ auto& A = this->getmat();
+ using namespace Dune::Hybrid;
+ forEach(integralRange(Dune::Hybrid::size(y)), [&](auto&& i)
+ {
+ forEach(integralRange(Dune::Hybrid::size(x)), [&](auto&& j)
+ {
+ const auto& mat = A[i][j];
+ const auto& xx = x[j];
+ auto& yy = y[i];
+
+ Dumux::parallelFor(mat.N(), [&](const std::size_t ii)
+ {
+ const auto& row = mat[ii];
+ const auto endj = row.end();
+ for (auto jj=row.begin(); jj!=endj; ++jj)
+ {
+ const auto& xj = Dune::Impl::asVector(xx[jj.index()]);
+ auto&& yi = Dune::Impl::asVector(yy[ii]);
+ Dune::Impl::asMatrix(*jj).usmv(alpha, xj, yi);
+ }
+ });
+ });
+ });
+ }
+};
+
+} // end namespace Dumux
+
+#endif
diff --git a/dumux/linear/seqsolverbackend.hh b/dumux/linear/seqsolverbackend.hh
index 6612e9daf7e1e4bb55cb86d0ffeae6013cb688f5..6827728cf1ce4a2df3aa3dc8ecf2bf7e34724ae1 100644
--- a/dumux/linear/seqsolverbackend.hh
+++ b/dumux/linear/seqsolverbackend.hh
@@ -42,6 +42,7 @@
#include <dumux/linear/solver.hh>
#include <dumux/linear/preconditioners.hh>
#include <dumux/linear/linearsolverparameters.hh>
+#include <dumux/linear/parallelmatrixadapter.hh>
namespace Dumux {
@@ -1040,7 +1041,7 @@ public:
bool solve(const Matrix& M, Vector& x, const Vector& b)
{
BlockDiagILU0Preconditioner<Matrix, Vector, Vector> preconditioner(M);
- Dune::MatrixAdapter<Matrix, Vector, Vector> op(M);
+ Dumux::ParallelMultiTypeMatrixAdapter<Matrix, Vector, Vector> op(M);
Dune::BiCGSTABSolver<Vector> solver(op, preconditioner, this->residReduction(),
this->maxIter(), this->verbosity());
auto bTmp(b);
@@ -1078,7 +1079,7 @@ public:
bool solve(const Matrix& M, Vector& x, const Vector& b)
{
BlockDiagILU0Preconditioner<Matrix, Vector, Vector> preconditioner(M);
- Dune::MatrixAdapter<Matrix, Vector, Vector> op(M);
+ Dumux::ParallelMultiTypeMatrixAdapter<Matrix, Vector, Vector> op(M);
static const int restartGMRes = getParamFromGroup<int>(this->paramGroup(), "LinearSolver.GMResRestart");
Dune::RestartedGMResSolver<Vector> solver(op, preconditioner, this->residReduction(), restartGMRes,
this->maxIter(), this->verbosity());
@@ -1284,7 +1285,7 @@ public:
BlockDiagAMGPreconditioner<Matrix, Vector, Vector> preconditioner(linearOperator, criterion, smootherArgs);
- Dune::MatrixAdapter<Matrix, Vector, Vector> op(m);
+ Dumux::ParallelMultiTypeMatrixAdapter<Matrix, Vector, Vector> op(m);
Dune::BiCGSTABSolver<Vector> solver(op, preconditioner, this->residReduction(),
this->maxIter(), this->verbosity());
auto bTmp(b);
diff --git a/examples/embedded_network_1d3d/solver.hh b/examples/embedded_network_1d3d/solver.hh
index cc957c09cc4764f1ce21a8e9717e4a1a59e58048..3b0e6193ac14537912b41e85c925722ffd493294 100644
--- a/examples/embedded_network_1d3d/solver.hh
+++ b/examples/embedded_network_1d3d/solver.hh
@@ -37,9 +37,9 @@
#include <dumux/common/typetraits/matrix.hh>
#include <dumux/common/typetraits/utility.hh>
#include <dumux/linear/solver.hh>
-#include <dumux/linear/amgbackend.hh>
#include <dumux/linear/preconditioners.hh>
#include <dumux/linear/linearsolverparameters.hh>
+#include <dumux/linear/parallelmatrixadapter.hh>
namespace Dumux::Example {
@@ -134,7 +134,7 @@ public:
void setup(const MatrixType& m)
{
preconditioner_ = std::make_unique<BlockDiagILU0Preconditioner<MatrixType, VectorType, VectorType>>(m);
- linearOperator_ = std::make_shared<Dune::MatrixAdapter<MatrixType, VectorType, VectorType>>(m);
+ linearOperator_ = std::make_shared<Dumux::ParallelMultiTypeMatrixAdapter<MatrixType, VectorType, VectorType>>(m);
solver_ = std::make_unique<Dune::BiCGSTABSolver<VectorType>>(*linearOperator_, *preconditioner_, this->residReduction(), this->maxIter(), this->verbosity());
isSetup_ = true;