diff --git a/examples/1protationsymmetry/doc/main.md b/examples/1protationsymmetry/doc/main.md
index 25e76e0a13756d432ed6bae036baaa8f6fc73f1e..ec84262262520d2bf0b81d710e45a3d8de7323c8 100644
--- a/examples/1protationsymmetry/doc/main.md
+++ b/examples/1protationsymmetry/doc/main.md
@@ -37,8 +37,10 @@ and compute the convergence rates.
 #include <dumux/common/parameters.hh> // for getParam
 #include <dumux/common/integrate.hh>  // for integrateL2Error
 
-#include <dumux/linear/seqsolverbackend.hh> // for ILU0BiCGSTABBackend
-#include <dumux/linear/pdesolver.hh>        // for LinearPDESolver
+#include <dumux/linear/istlsolvers.hh>
+#include <dumux/linear/linearsolvertraits.hh>
+#include <dumux/linear/linearalgebratraits.hh>
+#include <dumux/linear/pdesolver.hh>
 #include <dumux/assembly/fvassembler.hh>
 #include <dumux/assembly/diffmethod.hh>
 
@@ -145,8 +147,9 @@ solution and stores the result therein.
     using Assembler = FVAssembler<TypeTag, DiffMethod::analytic>;
     auto assembler = std::make_shared<Assembler>(problem, gridGeometry, gridVariables);
 
-    using LinearSolver = ILU0BiCGSTABBackend;
-    auto linearSolver = std::make_shared<LinearSolver>();
+    using LinearSolver = ILUBiCGSTABIstlSolver<LinearSolverTraits<GridGeometry>,
+                                               LinearAlgebraTraitsFromAssembler<Assembler>>;
+    auto linearSolver = std::make_shared<LinearSolver>(gridGeometry->gridView(), gridGeometry->dofMapper());
     LinearPDESolver<Assembler, LinearSolver> solver(assembler,  linearSolver);
     solver.setVerbose(false); // suppress output during solve()
 ```
diff --git a/examples/1protationsymmetry/main.cc b/examples/1protationsymmetry/main.cc
index 3dfbf9c23b55886360968186af13c4874c4a1cee..7e24bebd63e8aba4bd7eebe2baee614bea694e34 100644
--- a/examples/1protationsymmetry/main.cc
+++ b/examples/1protationsymmetry/main.cc
@@ -34,8 +34,10 @@
 #include <dumux/common/parameters.hh> // for getParam
 #include <dumux/common/integrate.hh>  // for integrateL2Error
 
-#include <dumux/linear/seqsolverbackend.hh> // for ILU0BiCGSTABBackend
-#include <dumux/linear/pdesolver.hh>        // for LinearPDESolver
+#include <dumux/linear/istlsolvers.hh>
+#include <dumux/linear/linearsolvertraits.hh>
+#include <dumux/linear/linearalgebratraits.hh>
+#include <dumux/linear/pdesolver.hh>
 #include <dumux/assembly/fvassembler.hh>
 #include <dumux/assembly/diffmethod.hh>
 
@@ -134,8 +136,9 @@ int main(int argc, char** argv) try
     using Assembler = FVAssembler<TypeTag, DiffMethod::analytic>;
     auto assembler = std::make_shared<Assembler>(problem, gridGeometry, gridVariables);
 
-    using LinearSolver = ILU0BiCGSTABBackend;
-    auto linearSolver = std::make_shared<LinearSolver>();
+    using LinearSolver = ILUBiCGSTABIstlSolver<LinearSolverTraits<GridGeometry>,
+                                               LinearAlgebraTraitsFromAssembler<Assembler>>;
+    auto linearSolver = std::make_shared<LinearSolver>(gridGeometry->gridView(), gridGeometry->dofMapper());
     LinearPDESolver<Assembler, LinearSolver> solver(assembler,  linearSolver);
     solver.setVerbose(false); // suppress output during solve()
     // [[/codeblock]]
diff --git a/examples/1ptracer/doc/main.md b/examples/1ptracer/doc/main.md
index a011f31d69168a6d66d204d932a1796e039b299b..a66dcf25fccc7c8d143f703ad111117c6511b6e4 100644
--- a/examples/1ptracer/doc/main.md
+++ b/examples/1ptracer/doc/main.md
@@ -45,7 +45,9 @@ The following files contains the available linear solver backends and the assemb
 systems arising from finite volume discretizations (box-scheme, tpfa-approximation, mpfa-approximation).
 
 ```cpp
-#include <dumux/linear/seqsolverbackend.hh>
+#include <dumux/linear/istlsolvers.hh>
+#include <dumux/linear/linearsolvertraits.hh>
+#include <dumux/linear/linearalgebratraits.hh>
 #include <dumux/assembly/fvassembler.hh>
 #include <dumux/assembly/diffmethod.hh> // analytic or numeric differentiation
 ```
@@ -158,7 +160,7 @@ The grid variables are used store variables (primary and secondary variables) on
 ```
 
 We now instantiate the assembler class, assemble the linear system and solve it with the linear
-solver ILUnBiCGSTABBackend (a bi-conjugate gradient solver preconditioned by an incomplete LU-factorization preconditioner).
+solver AMGCGIstlSolver (a conjugate gradient solver preconditioned by an algebraic multigrid).
 Besides that, the time needed for assembly and solve is measured and printed.
 
 ```cpp
@@ -173,10 +175,10 @@ Besides that, the time needed for assembly and solve is measured and printed.
 
     (*r) *= -1.0; // We want to solve `Ax = -r`.
 
-    using LinearSolver = ILUnBiCGSTABBackend;
+    using OnePLinearSolver = AMGCGIstlSolver<LinearSolverTraits<GridGeometry>, LinearAlgebraTraitsFromAssembler<OnePAssembler>>;
     Dune::Timer solverTimer; std::cout << "Solving linear system ..." << std::flush;
-    auto linearSolver = std::make_shared<LinearSolver>();
-    linearSolver->solve(*A, p, *r);
+    auto onePLinearSolver = std::make_shared<OnePLinearSolver>(gridGeometry->gridView(), gridGeometry->dofMapper());
+    onePLinearSolver->solve(*A, p, *r);
     solverTimer.stop(); std::cout << " took " << solverTimer.elapsed() << " seconds." << std::endl;
 
     Dune::Timer updateTimer; std::cout << "Updating variables ..." << std::flush;
@@ -199,7 +201,7 @@ problem defined in `problem_1p.hh`. Let us now write this solution to a VTK file
 
     // print overall CPU time required for assembling and solving the 1p problem.
     timer.stop();
-    const auto& comm = Dune::MPIHelper::getCommunication();
+    const auto& comm = leafGridView.comm();
     std::cout << "Simulation took " << timer.elapsed() << " seconds on "
               << comm.size() << " processes.\n"
               << "The cumulative CPU time was " << timer.elapsed()*comm.size() << " seconds.\n";
@@ -291,11 +293,15 @@ Moreover, we define 10 check points in the time loop at which we will write the
 
 We create and initialize the assembler with a time loop for the transient problem.
 Within the time loop, we will use this assembler in each time step to assemble the linear system.
+As solver we use the ILUBiCGSTABIstlSolver (a bi-conjugate gradient solver preconditioned by an incomplete LU-factorization).
 
 ```cpp
     using TracerAssembler = FVAssembler<TracerTypeTag, DiffMethod::analytic, /*implicit=*/false>;
     auto assembler = std::make_shared<TracerAssembler>(tracerProblem, gridGeometry, gridVariables, timeLoop, xOld);
     assembler->setLinearSystem(A, r);
+
+    using TracerLinearSolver = ILUBiCGSTABIstlSolver<LinearSolverTraits<GridGeometry>, LinearAlgebraTraitsFromAssembler<TracerAssembler>>;
+    auto tracerLinearSolver = std::make_shared<TracerLinearSolver>(gridGeometry->gridView(), gridGeometry->dofMapper());
 ```
 
 The following lines of code initialize the vtk output module, add the velocity output facility
@@ -334,7 +340,7 @@ and the time step sizes used is printed to the terminal.
         // We solve the linear system `A(xOld-xNew) = r`.
         Dune::Timer solveTimer;
         SolutionVector xDelta(x);
-        linearSolver->solve(*A, xDelta, *r);
+        tracerLinearSolver->solve(*A, xDelta, *r);
         solveTimer.stop();
 
         // update the solution vector and the grid variables.
diff --git a/examples/1ptracer/main.cc b/examples/1ptracer/main.cc
index 31cbca6cdee009341580854930b43efe11b0f6de..94cf78169e99afde1b4164286c091e15e7784268 100644
--- a/examples/1ptracer/main.cc
+++ b/examples/1ptracer/main.cc
@@ -39,7 +39,9 @@
 
 // The following files contains the available linear solver backends and the assembler for the linear
 // systems arising from finite volume discretizations (box-scheme, tpfa-approximation, mpfa-approximation).
-#include <dumux/linear/seqsolverbackend.hh>
+#include <dumux/linear/istlsolvers.hh>
+#include <dumux/linear/linearsolvertraits.hh>
+#include <dumux/linear/linearalgebratraits.hh>
 #include <dumux/assembly/fvassembler.hh>
 #include <dumux/assembly/diffmethod.hh> // analytic or numeric differentiation
 
@@ -124,7 +126,7 @@ int main(int argc, char** argv) try
     onePGridVariables->init(p);
 
     // We now instantiate the assembler class, assemble the linear system and solve it with the linear
-    // solver ILUnBiCGSTABBackend (a bi-conjugate gradient solver preconditioned by an incomplete LU-factorization preconditioner).
+    // solver AMGCGIstlSolver (a conjugate gradient solver preconditioned by an algebraic multigrid).
     // Besides that, the time needed for assembly and solve is measured and printed.
     using OnePAssembler = FVAssembler<OnePTypeTag, DiffMethod::analytic>;
     auto assemblerOneP = std::make_shared<OnePAssembler>(problemOneP, gridGeometry, onePGridVariables);
@@ -137,10 +139,10 @@ int main(int argc, char** argv) try
 
     (*r) *= -1.0; // We want to solve `Ax = -r`.
 
-    using LinearSolver = ILUnBiCGSTABBackend;
+    using OnePLinearSolver = AMGCGIstlSolver<LinearSolverTraits<GridGeometry>, LinearAlgebraTraitsFromAssembler<OnePAssembler>>;
     Dune::Timer solverTimer; std::cout << "Solving linear system ..." << std::flush;
-    auto linearSolver = std::make_shared<LinearSolver>();
-    linearSolver->solve(*A, p, *r);
+    auto onePLinearSolver = std::make_shared<OnePLinearSolver>(gridGeometry->gridView(), gridGeometry->dofMapper());
+    onePLinearSolver->solve(*A, p, *r);
     solverTimer.stop(); std::cout << " took " << solverTimer.elapsed() << " seconds." << std::endl;
 
     Dune::Timer updateTimer; std::cout << "Updating variables ..." << std::flush;
@@ -161,7 +163,7 @@ int main(int argc, char** argv) try
 
     // print overall CPU time required for assembling and solving the 1p problem.
     timer.stop();
-    const auto& comm = Dune::MPIHelper::getCommunication();
+    const auto& comm = leafGridView.comm();
     std::cout << "Simulation took " << timer.elapsed() << " seconds on "
               << comm.size() << " processes.\n"
               << "The cumulative CPU time was " << timer.elapsed()*comm.size() << " seconds.\n";
@@ -245,10 +247,14 @@ int main(int argc, char** argv) try
 
     // We create and initialize the assembler with a time loop for the transient problem.
     // Within the time loop, we will use this assembler in each time step to assemble the linear system.
+    // As solver we use the ILUBiCGSTABIstlSolver (a bi-conjugate gradient solver preconditioned by an incomplete LU-factorization).
     using TracerAssembler = FVAssembler<TracerTypeTag, DiffMethod::analytic, /*implicit=*/false>;
     auto assembler = std::make_shared<TracerAssembler>(tracerProblem, gridGeometry, gridVariables, timeLoop, xOld);
     assembler->setLinearSystem(A, r);
 
+    using TracerLinearSolver = ILUBiCGSTABIstlSolver<LinearSolverTraits<GridGeometry>, LinearAlgebraTraitsFromAssembler<TracerAssembler>>;
+    auto tracerLinearSolver = std::make_shared<TracerLinearSolver>(gridGeometry->gridView(), gridGeometry->dofMapper());
+
     // The following lines of code initialize the vtk output module, add the velocity output facility
     // and write out the initial solution. At each checkpoint, we will use the output module to write
     // the solution of a time step into a corresponding vtk file.
@@ -283,7 +289,7 @@ int main(int argc, char** argv) try
         // We solve the linear system `A(xOld-xNew) = r`.
         Dune::Timer solveTimer;
         SolutionVector xDelta(x);
-        linearSolver->solve(*A, xDelta, *r);
+        tracerLinearSolver->solve(*A, xDelta, *r);
         solveTimer.stop();
 
         // update the solution vector and the grid variables.
diff --git a/examples/embedded_network_1d3d/bloodflow.hh b/examples/embedded_network_1d3d/bloodflow.hh
index f0c821b16deaf8ff08a94b82477152d93a2ff23a..2c051b86b81d3d512cdeb8e6ac9fc4525d52980e 100644
--- a/examples/embedded_network_1d3d/bloodflow.hh
+++ b/examples/embedded_network_1d3d/bloodflow.hh
@@ -42,7 +42,9 @@
 #include <dumux/material/components/constant.hh>
 #include <dumux/material/fluidsystems/1pliquid.hh>
 
-#include <dumux/linear/seqsolverbackend.hh>
+#include <dumux/linear/istlsolvers.hh>
+#include <dumux/linear/linearsolvertraits.hh>
+#include <dumux/linear/linearalgebratraits.hh>
 #include <dumux/assembly/fvassembler.hh>
 #include <dumux/nonlinear/newtonsolver.hh>
 
@@ -462,7 +464,7 @@ std::vector<double> computeBloodVolumeFluxes(const GridGeometry& gridGeometry, c
 
     using Assembler = FVAssembler<TypeTag, DiffMethod::analytic>;
     auto assembler = std::make_shared<Assembler>(problem, gridGeometry, gridVariables);
-    using LinearSolver = UMFPackBackend;
+    using LinearSolver = UMFPackIstlSolver<SeqLinearSolverTraits, LinearAlgebraTraitsFromAssembler<Assembler>>;
     auto linearSolver = std::make_shared<LinearSolver>();
     using NewtonSolver = Dumux::NewtonSolver<Assembler, LinearSolver>;
     NewtonSolver nonLinearSolver(assembler, linearSolver);
diff --git a/examples/porenetwork_upscaling/doc/main.md b/examples/porenetwork_upscaling/doc/main.md
index d484fe132dd19b6ea5958004d9f1a94ba9f29d8f..0d8505c2ce9dddf736cc67cc7c6db7d4c263d941 100644
--- a/examples/porenetwork_upscaling/doc/main.md
+++ b/examples/porenetwork_upscaling/doc/main.md
@@ -39,8 +39,10 @@ Pore-Network-Model to evaluate the upscaled Darcy permeability of a given networ
 #include <dumux/common/parameters.hh> // for getParam
 #include <dumux/common/initialize.hh>
 
-#include <dumux/linear/seqsolverbackend.hh> // for ILU0BiCGSTABBackend
-#include <dumux/linear/pdesolver.hh>        // for LinearPDESolver
+#include <dumux/linear/istlsolvers.hh>
+#include <dumux/linear/linearsolvertraits.hh>
+#include <dumux/linear/linearalgebratraits.hh>
+#include <dumux/linear/pdesolver.hh>
 #include <dumux/nonlinear/newtonsolver.hh>
 #include <dumux/assembly/fvassembler.hh>
 
@@ -114,7 +116,7 @@ solver classes to assemble and solve the non-linear system.
     using Assembler = FVAssembler<TypeTag, DiffMethod::numeric>;
     auto assembler = std::make_shared<Assembler>(problem, gridGeometry, gridVariables);
 
-    using LinearSolver = UMFPackBackend;
+    using LinearSolver = Dumux::UMFPackIstlSolver<SeqLinearSolverTraits, LinearAlgebraTraitsFromAssembler<Assembler>>;
     auto linearSolver = std::make_shared<LinearSolver>();
 
     using NewtonSolver = NewtonSolver<Assembler, LinearSolver>;
diff --git a/examples/porenetwork_upscaling/main.cc b/examples/porenetwork_upscaling/main.cc
index 540109b475ddf409ebde8a375a9d49f5cb7d993e..75d4b3d331d0278fa036eceedf3d499a4a7891de 100644
--- a/examples/porenetwork_upscaling/main.cc
+++ b/examples/porenetwork_upscaling/main.cc
@@ -36,8 +36,10 @@
 #include <dumux/common/parameters.hh> // for getParam
 #include <dumux/common/initialize.hh>
 
-#include <dumux/linear/seqsolverbackend.hh> // for ILU0BiCGSTABBackend
-#include <dumux/linear/pdesolver.hh>        // for LinearPDESolver
+#include <dumux/linear/istlsolvers.hh>
+#include <dumux/linear/linearsolvertraits.hh>
+#include <dumux/linear/linearalgebratraits.hh>
+#include <dumux/linear/pdesolver.hh>
 #include <dumux/nonlinear/newtonsolver.hh>
 #include <dumux/assembly/fvassembler.hh>
 
@@ -103,7 +105,7 @@ void runExample()
     using Assembler = FVAssembler<TypeTag, DiffMethod::numeric>;
     auto assembler = std::make_shared<Assembler>(problem, gridGeometry, gridVariables);
 
-    using LinearSolver = UMFPackBackend;
+    using LinearSolver = Dumux::UMFPackIstlSolver<SeqLinearSolverTraits, LinearAlgebraTraitsFromAssembler<Assembler>>;
     auto linearSolver = std::make_shared<LinearSolver>();
 
     using NewtonSolver = NewtonSolver<Assembler, LinearSolver>;