diff --git a/python/dune/perftool/pdelab/driver.py b/python/dune/perftool/pdelab/driver.py
index a3d7e88aa86c4c3f6900c43918b27eb32be0ec20..6bed1ab630fe7e54fed481d504a45af85d78f28a 100644
--- a/python/dune/perftool/pdelab/driver.py
+++ b/python/dune/perftool/pdelab/driver.py
@@ -1387,6 +1387,17 @@ def time_loop():
     elif not linear and not matrix_free:
         assert False
 
+    print_residual()
+    print_matrix()
+    if get_option("exact_solution_expression"):
+        from ufl import MixedElement, VectorElement, TensorElement
+        if isinstance(_driver_data['form'].coefficients()[0].ufl_element(), (MixedElement, VectorElement, TensorElement)):
+            raise NotImplementedError("Comparing to exact solution is olny implemented for scalar elements.")
+
+        if get_option("compare_dofs"):
+            compare_dofs()
+        if get_option("compare_l2errorsquared"):
+            compare_L2_squared()
 
 def generate_driver(formdatas, data):
     # The driver module uses a global dictionary for storing necessary data
diff --git a/test/heatequation/CMakeLists.txt b/test/heatequation/CMakeLists.txt
index 3ac7152738434df1d396fa582d18ee8f7963bd8f..543dd2d08e6b2c862c2f990b1cc992bb559a0c26 100644
--- a/test/heatequation/CMakeLists.txt
+++ b/test/heatequation/CMakeLists.txt
@@ -1,6 +1,6 @@
 add_generated_executable(UFLFILE heatequation.ufl
                          TARGET heatequation_implicit
-                         FORM_COMPILER_ARGS
+                         FORM_COMPILER_ARGS --exact-solution-expression g --compare-l2errorsquared 1e-7
                          )
 
 
@@ -10,7 +10,7 @@ dune_add_system_test(TARGET heatequation_implicit
 
 add_generated_executable(UFLFILE heatequation.ufl
                          TARGET heatequation_explicit
-                         FORM_COMPILER_ARGS --explicit-time-stepping
+                         FORM_COMPILER_ARGS --explicit-time-stepping --exact-solution-expression g --compare-l2errorsquared 1e-7
                          )
 
 
@@ -18,14 +18,10 @@ dune_add_system_test(TARGET heatequation_explicit
                      INIFILE heatequation_explicit.mini
                      )
 
+# Hand written reference solutions
 add_executable(heatequation_implicit_ref heatequation_implicit_ref.cc)
 dune_symlink_to_source_files(FILES heatequation_implicit_ref.ini)
 set_target_properties(heatequation_implicit_ref PROPERTIES EXCLUDE_FROM_ALL 1)
 
 add_executable(heatequation_explicit_ref heatequation_explicit_ref.cc)
 set_target_properties(heatequation_explicit_ref PROPERTIES EXCLUDE_FROM_ALL 1)
-
-# # Add the reference code with the PDELab localoperator that produced
-# # the reference vtk file
-# add_executable(heatequation_ref heatequation_main.cc)
-# set_target_properties(heatequation_ref PROPERTIES EXCLUDE_FROM_ALL 1)
diff --git a/test/heatequation/heatequation.ufl b/test/heatequation/heatequation.ufl
index 5495df8a009101d17991e69627f18eedbf1fffc0..b6d3f65c7e107483029acd67226bd8e8ec6ab455 100644
--- a/test/heatequation/heatequation.ufl
+++ b/test/heatequation/heatequation.ufl
@@ -1,3 +1,4 @@
+f = Expression("Dune::FieldVector<double,2> c(0.5); c-= x; return 4.*(1.-c.two_norm2())*std::exp(-1.*c.two_norm2());")
 g = Expression("Dune::FieldVector<double,2> c(0.5); c-= x; return std::exp(-1.*c.two_norm2());")
 
 V = FiniteElement("CG", "triangle", 1, dirichlet_expression=g)
@@ -5,6 +6,6 @@ u = TrialFunction(V)
 v = TestFunction(V)
 
 mass = (u*v)*dx
-poisson = (inner(grad(u), grad(v)))*dx
+poisson = (inner(grad(u), grad(v)) - f*v)*dx
 
 forms = [mass, poisson]
\ No newline at end of file
diff --git a/test/heatequation/heatequation_explicit_ref.cc b/test/heatequation/heatequation_explicit_ref.cc
index 28470723359409f3f055c74a17fb8fa79f698eaf..9a07f04958156b7a93d1f17106695c6c81d76312 100644
--- a/test/heatequation/heatequation_explicit_ref.cc
+++ b/test/heatequation/heatequation_explicit_ref.cc
@@ -53,9 +53,13 @@ public:
 
   //! source term
   typename Traits::RangeFieldType
-  f (const typename Traits::ElementType& e, const typename Traits::DomainType& x) const
+  f (const typename Traits::ElementType& e, const typename Traits::DomainType& xlocal) const
   {
-    return 0.0;
+    typename Traits::DomainType x = e.geometry().global(xlocal);
+    Dune::FieldVector<double,2> c(0.5);
+    c-= x;
+    return 4.*(1.-c.two_norm2())*std::exp(-1.*c.two_norm2());
+    // return 0.0;
   }
 
   //! boundary condition type function
diff --git a/test/heatequation/problem.hh b/test/heatequation/problem.hh
index 37c93f23fa916467d0cffacfbf090431b64478c7..8c1a13c6d53b2c822b1a4c276f23c4f47a80f0d3 100644
--- a/test/heatequation/problem.hh
+++ b/test/heatequation/problem.hh
@@ -29,9 +29,9 @@ public:
   template<typename E, typename X>
   Number f (const E& e, const X& local) const
   {
-    return 0.0;
-    // auto x = e.geometry().global(local);
-    // Dune::FieldVector<double,2> c(0.5); c-= x; return 4.*(1.-c.two_norm2())*std::exp(-1.*c.two_norm2());
+    // return 0.0;
+    auto x = e.geometry().global(local);
+    Dune::FieldVector<double,2> c(0.5); c-= x; return 4.*(1.-c.two_norm2())*std::exp(-1.*c.two_norm2());
   }
 
   //! boundary condition type function (true = Dirichlet)