diff --git a/dumux/linear/pdesolver.hh b/dumux/linear/pdesolver.hh
index ce204474e949b65155ccd7906161384a6c048fc5..2303b8cbc50c07c68a42fe65f979641d5ca193ca 100644
--- a/dumux/linear/pdesolver.hh
+++ b/dumux/linear/pdesolver.hh
@@ -75,6 +75,7 @@ public:
                     const std::string& paramGroup = "")
     : ParentType(assembler, linearSolver)
     , paramGroup_(paramGroup)
+    , reuseMatrix_(false)
     {
         initParams_(paramGroup);
 
@@ -102,7 +103,10 @@ public:
 
         // linearize the problem at the current solution
         assembleTimer.start();
-        this->assembler().assembleJacobianAndResidual(uCurrentIter);
+        if (reuseMatrix_)
+            this->assembler().assembleResidual(uCurrentIter);
+        else
+            this->assembler().assembleJacobianAndResidual(uCurrentIter);
         assembleTimer.stop();
 
         ///////////////
@@ -194,6 +198,15 @@ public:
     const std::string& paramGroup() const
     { return paramGroup_; }
 
+    /*!
+     * \brief Set whether the matrix should be reused
+     * \note If this is set to true, the matrix will not be assembled. Make
+     *       sure there is an assembled matrix that can be resused before
+     *       setting this flag to true.
+     */
+    void reuseMatrix(bool reuse = true)
+    { reuseMatrix_ = reuse; }
+
 private:
 
     virtual bool solveLinearSystem_(SolutionVector& deltaU)
@@ -322,6 +335,9 @@ private:
 
     //! the parameter group for getting parameters from the parameter tree
     std::string paramGroup_;
+
+    //! check if the matrix is supposed to be reused
+    bool reuseMatrix_;
 };
 
 } // end namespace Dumux
diff --git a/test/porousmediumflow/tracer/constvel/main.cc b/test/porousmediumflow/tracer/constvel/main.cc
index ee423251b2be1afd26b48f93e196cd2fe7612022..ce969bab15b1d42aad0bffc13780234464402836 100644
--- a/test/porousmediumflow/tracer/constvel/main.cc
+++ b/test/porousmediumflow/tracer/constvel/main.cc
@@ -35,6 +35,7 @@
 #include <dumux/common/dumuxmessage.hh>
 
 #include <dumux/linear/seqsolverbackend.hh>
+#include <dumux/linear/pdesolver.hh>
 #include <dumux/assembly/fvassembler.hh>
 
 #include <dumux/io/vtkoutputmodule.hh>
@@ -107,15 +108,16 @@ int main(int argc, char** argv)
     //! the assembler with time loop for instationary problem
     using Assembler = FVAssembler<TypeTag, DiffMethod::analytic, IMPLICIT>;
     auto assembler = std::make_shared<Assembler>(problem, gridGeometry, gridVariables, timeLoop, xOld);
-    using JacobianMatrix = GetPropType<TypeTag, Properties::JacobianMatrix>;
-    auto A = std::make_shared<JacobianMatrix>();
-    auto r = std::make_shared<SolutionVector>();
-    assembler->setLinearSystem(A, r);
 
     //! the linear solver
     using LinearSolver = UMFPackBackend;
     auto linearSolver = std::make_shared<LinearSolver>();
 
+    //! pde solver (assemble, solve, update)
+    LinearPDESolver solver(assembler, linearSolver);
+    assembler->assembleJacobianAndResidual(x);
+    solver.reuseMatrix();
+
     //! intialize the vtk output module
     VtkOutputModule<GridVariables, SolutionVector> vtkWriter(*gridVariables, x, problem->name());
     using VelocityOutput = GetPropType<TypeTag, Properties::VelocityOutput>;
@@ -135,30 +137,8 @@ int main(int argc, char** argv)
     timeLoop->start();
     while (!timeLoop->finished())
     {
-        Dune::Timer assembleTimer;
-        assembler->assembleJacobianAndResidual(x);
-        assembleTimer.stop();
-
-        // solve the linear system A(xOld-xNew) = r
-        Dune::Timer solveTimer;
-        SolutionVector xDelta(x);
-        linearSolver->solve(*A, xDelta, *r);
-        solveTimer.stop();
-
-        // update solution and grid variables
-        Dune::Timer updateTimer;
-        x -= xDelta;
-        gridVariables->update(x);
-        updateTimer.stop();
-
-        // statistics
-        const auto elapsedTot = assembleTimer.elapsed() + solveTimer.elapsed() + updateTimer.elapsed();
-        if (mpiHelper.rank() == 0)
-            std::cout << "Assemble/solve/update time: "
-                      <<  assembleTimer.elapsed() << "(" << 100*assembleTimer.elapsed()/elapsedTot << "%)/"
-                      <<  solveTimer.elapsed() << "(" << 100*solveTimer.elapsed()/elapsedTot << "%)/"
-                      <<  updateTimer.elapsed() << "(" << 100*updateTimer.elapsed()/elapsedTot << "%)"
-                      <<  std::endl;
+        // assemble & solve & update
+        solver.solve(x);
 
         // make the new solution the old solution
         xOld = x;