diff --git a/applications/poisson_dg/donkey.sbatch b/applications/poisson_dg/donkey.sbatch
index 50a3570163421c471b282ef5a36cfde07ccf0f65..abf8b1aa2795fcb894d0e7dccb5bf9f12d33e0ae 100755
--- a/applications/poisson_dg/donkey.sbatch
+++ b/applications/poisson_dg/donkey.sbatch
@@ -31,6 +31,8 @@ srun $SRUNOPT ./app_poisson_dg_deg5_opcount app_poisson_dg_3d_deg5_opcount.ini
 srun $SRUNOPT ./app_poisson_dg_deg6_opcount app_poisson_dg_3d_deg6_opcount.ini
 srun $SRUNOPT ./app_poisson_dg_deg7_opcount app_poisson_dg_3d_deg7_opcount.ini
 srun $SRUNOPT ./app_poisson_dg_deg8_opcount app_poisson_dg_3d_deg8_opcount.ini
+srun $SRUNOPT ./app_poisson_dg_deg9_opcount app_poisson_dg_3d_deg9_opcount.ini
+srun $SRUNOPT ./app_poisson_dg_deg10_opcount app_poisson_dg_3d_deg10_opcount.ini
 
 # Run the timing executables
 COUNT=0
@@ -42,5 +44,7 @@ while [ $COUNT -lt 2 ]; do
     srun $SRUNOPT ./app_poisson_dg_deg6_nonopcount app_poisson_dg_3d_deg6_nonopcount.ini
     srun $SRUNOPT ./app_poisson_dg_deg7_nonopcount app_poisson_dg_3d_deg7_nonopcount.ini
     srun $SRUNOPT ./app_poisson_dg_deg8_nonopcount app_poisson_dg_3d_deg8_nonopcount.ini
+    srun $SRUNOPT ./app_poisson_dg_deg9_nonopcount app_poisson_dg_3d_deg9_nonopcount.ini
+    srun $SRUNOPT ./app_poisson_dg_deg10_nonopcount app_poisson_dg_3d_deg10_nonopcount.ini
     COUNT=$((COUNT + 1))
 done
diff --git a/applications/poisson_dg/poisson_dg.mini b/applications/poisson_dg/poisson_dg.mini
index 26c9b94ea8ae55f39ca50ca18df2716fc79409a6..a8aac4c2a2cef789f5380010c3cfd1fa6111500a 100644
--- a/applications/poisson_dg/poisson_dg.mini
+++ b/applications/poisson_dg/poisson_dg.mini
@@ -6,7 +6,7 @@ opcount_suffix = opcount, nonopcount | expand opcount
 # Calculate the size of the grid to equlibritate it to 100 MB/rank
 # Input parameters
 dim = 3
-mbperrank = 10
+mbperrank = 100
 ranks = 16
 floatingbytes = 8
 
@@ -35,10 +35,10 @@ fastdg = 1
 sumfact = 1
 vectorize_quad = 1
 vectorize_grads = 1
-instrumentation_level = 4
+instrumentation_level = 2
 opcounter = 1, 0 | expand opcount
 time_opcounter = 0, 1 | expand opcount
 
 [formcompiler.ufl_variants]
 cell = hexahedron
-degree = 2, 3, 4, 5, 6, 7, 8 | expand
+degree = 2, 3, 4, 5, 6, 7, 8, 9, 10 | expand
diff --git a/bin/plot_measurements.py b/bin/plot_measurements.py
index 9de098ddc30d225af351eb5a4d6583359d043d9f..87bda47817f3db5772d6810267dcc9b471c07951 100755
--- a/bin/plot_measurements.py
+++ b/bin/plot_measurements.py
@@ -1,5 +1,6 @@
 #!/usr/bin/env python
 
+import itertools
 import pandas
 import matplotlib
 import sys
@@ -8,24 +9,27 @@ matplotlib.use("PDF")
 
 from matplotlib import pyplot as plt
 
-title = sys.argv[1]
+what = sys.argv[1]
+title = sys.argv[2]
 filename = title.lower().replace(" ", "_") + ".pdf"
 
-flopframe = pandas.read_csv("./poissondg-insn2/floprates.csv", header=None, delimiter=" ", names=("exec", "degree", "what", "GFlops"))
-flopframe = flopframe[flopframe.what == "residual_evaluation"]
+flopframe = pandas.read_csv("./floprates.csv", header=None, delimiter=" ", names=("exec", "degree", "what", "GFlops"))
+flopframe = flopframe[flopframe.what == what]
 
-timeframe = pandas.read_csv("./poissondg-insn2/doftimes.csv", header=None, delimiter=" ", names=("exec", "degree", "what", "DOFs"))
-timeframe = timeframe[timeframe.what == "residual_evaluation"]
+timeframe = pandas.read_csv("./doftimes.csv", header=None, delimiter=" ", names=("exec", "degree", "what", "DOFs"))
+timeframe = timeframe[timeframe.what == what]
 
 fig, ax1 = plt.subplots()
 
 ax2 = ax1.twinx()
-ax1.plot(frame['degree'], flopframe['GFlops'])
-ax2.plot(frame['degree'], timeframe['DOFs'])
+x, y = list(itertools.izip(*sorted(itertools.izip(flopframe['degree'], flopframe['GFlops']))))
+ax1.plot(x, y, "b-", x, y, "bo")
+x, y = list(itertools.izip(*sorted(itertools.izip(timeframe['degree'], timeframe['DOFs']))))
+ax2.plot(x, y, "r-", x, y, "ro")
 
 ax1.set_xlabel("Polynomial degree")
-ax1.set_ylabel("GFlops")
-ax2.set_ylabel("DOFs / s")
+ax1.set_ylabel("GFlops /s", color="b")
+ax2.set_ylabel("MDOFs / s", color="r")
 plt.title(title)
 
 plt.savefig(filename)
diff --git a/python/dune/perftool/pdelab/driver.py b/python/dune/perftool/pdelab/driver.py
index 524ce57967e4c4c42f94d08b30c433ba2aab40ed..e11e659fc0c8a4f2e3045cb13485c492c626bb2b 100644
--- a/python/dune/perftool/pdelab/driver.py
+++ b/python/dune/perftool/pdelab/driver.py
@@ -1471,6 +1471,46 @@ def evaluate_residual_timer():
     return evaluation
 
 
+@preamble
+def apply_jacobian_timer():
+    # Set the matrix_free option to True!
+    from dune.perftool.options import set_option
+    set_option("matrix_free", True)
+
+    formdata = _driver_data['formdata']
+    n_go = name_gridoperator(formdata)
+    v = name_vector(formdata)
+    t_v = type_vector(formdata)
+    setup_timer()
+
+    if get_option('instrumentation_level') >= 2:
+        # Write back times
+        from dune.perftool.generation import post_include
+        post_include("HP_DECLARE_TIMER(apply_jacobian);", filetag="driver")
+        timestream = name_timing_stream()
+        print_times = []
+
+    from dune.perftool.generation import get_global_context_value
+    formdatas = get_global_context_value("formdatas")
+    for formdata in formdatas:
+        lop_name = name_localoperator(formdata)
+        if get_option('instrumentation_level') >= 3:
+            print_times.append("{}.dump_timers({}, argv[0], true);".format(lop_name, timestream))
+
+    if get_option('instrumentation_level') >= 2:
+        evaluation = ["HP_TIMER_START(apply_jacobian);",
+                      "{}.jacobian_apply({}, j);".format(n_go, v),
+                      "HP_TIMER_STOP(apply_jacobian);",
+                      "DUMP_TIMER(apply_jacobian, {}, true);".format(timestream)]
+        evaluation.extend(print_times)
+    else:
+        evaluation = ["{}.jacobian_apply({}, j);".format(n_go, v)]
+
+    evaluation = ["{} j({});".format(t_v, v), "j=0.0;"] + evaluation
+
+    return evaluation
+
+
 @preamble
 def assemble_matrix_timer():
     formdata = _driver_data['formdata']
@@ -1743,6 +1783,7 @@ def generate_driver(formdatas, data):
         # In case of operator conunting we only assemble the matrix and evaluate the residual
         #assemble_matrix_timer()
         evaluate_residual_timer()
+#        apply_jacobian_timer()
     elif is_stationary():
         # We could also use solve if we are not interested in visualization
         vtkoutput()
diff --git a/python/dune/perftool/pdelab/signatures.py b/python/dune/perftool/pdelab/signatures.py
index d2c48409d746278ea35d9abb22d99a9bad6675e8..97560a248a4aa28f3db489859c64dcb48d54c58c 100644
--- a/python/dune/perftool/pdelab/signatures.py
+++ b/python/dune/perftool/pdelab/signatures.py
@@ -274,7 +274,7 @@ def jacobian_apply_skeleton_templates():
     lfsvt = type_testfunctionspace()
     cct = type_coefficientcontainer()
     avt = type_accumulation_variable()
-    return (geot, lfsut, cct, lfsvt, lfsut, cct, lfsvt, avt)
+    return (geot, lfsut, cct, lfsvt, lfsut, cct, lfsvt, avt, avt)
 
 
 def jacobian_apply_skeleton_args():
diff --git a/python/dune/perftool/sumfact/sumfact.py b/python/dune/perftool/sumfact/sumfact.py
index c4890794ba31ae13c3911bda8f7fc9688153e167..45894cc09a5c73e24e22bb78ab5209623ca9eb22 100644
--- a/python/dune/perftool/sumfact/sumfact.py
+++ b/python/dune/perftool/sumfact/sumfact.py
@@ -269,7 +269,7 @@ def generate_accumulation_instruction(visitor, accterm, measure, subdomain_id):
         # variable.
         if get_option('fastdg'):
             ft = get_global_context_value("form_type")
-            if ft == 'residual':
+            if ft == 'residual' or ft == 'jacobian_apply':
                 accum = accum + ".data()"
                 size = basis_functions_per_direction() ** world_dimension()
                 globalarg(accum, dtype=np.float64, shape=(size,), managed=False)