diff --git a/python/dune/perftool/options.py b/python/dune/perftool/options.py
index 85efbb6630d537e2d10c3849c10562763a072bd4..85080483c61297415b3b5eb47271baaeb0828a3a 100644
--- a/python/dune/perftool/options.py
+++ b/python/dune/perftool/options.py
@@ -68,6 +68,7 @@ class PerftoolOptionsArray(ImmutableRecord):
     simplify = PerftoolOption(default=False, helpstr="Whether to simplify expressions using sympy")
     precision_bits = PerftoolOption(default=64, helpstr="The number of bits for the floating point type")
     assure_statement_ordering = PerftoolOption(default=False, helpstr="Whether special care should be taken for a good statement ordering in sumfact kernels, runs into a loopy scheduler performance bug, but is necessary for production.")
+    generate_jacobians = PerftoolOption(default=True, helpstr="Whether jacobian_* methods should be generated. This is set to false automatically, when numerical_jacobian is set to true.")
 
     # Arguments that are mainly to be set by logic depending on other options
     max_vector_width = PerftoolOption(default=256, helpstr=None)
@@ -126,6 +127,9 @@ def process_options(opt):
     if opt.sumfact:
         opt = opt.copy(unroll_dimension_loops=True)
 
+    if opt.numerical_jacobian:
+        opt = opt.copy(generate_jacobians=False)
+
     return opt
 
 
diff --git a/python/dune/perftool/pdelab/localoperator.py b/python/dune/perftool/pdelab/localoperator.py
index 23c92da10c0945b6e4196b1ccde44f0cdbdb5c83..76d6b0c162e2b5677905ac6c79cfefc8f9f3e01b 100644
--- a/python/dune/perftool/pdelab/localoperator.py
+++ b/python/dune/perftool/pdelab/localoperator.py
@@ -821,22 +821,23 @@ def generate_localoperator_kernels(formdata, data):
         jacform = preprocess_form(jacform).preprocessed_form
 
         with global_context(form_type="jacobian"):
-            for measure in set(i.integral_type() for i in jacform.integrals()):
-                logger.info("generate_localoperator_kernels: measure {}".format(measure))
-                with global_context(integral_type=measure):
-                    with global_context(kernel=assembler_routine_name()):
-                        kernel = [k for k in get_backend(interface="generate_kernels_per_integral")(jacform.integrals_by_type(measure))]
-                operator_kernels[(measure, 'jacobian')] = kernel
-
-            # Generate dummy functions for those kernels, that vanished in the differentiation process
-            # We *could* solve this problem by using lambda_* terms but we do not really want that, so
-            # we use empty jacobian assembly methods instead
-            alpha_measures = set(i.integral_type() for i in form.integrals())
-            jacobian_measures = set(i.integral_type() for i in jacform.integrals())
-            for it in alpha_measures - jacobian_measures:
-                with global_context(integral_type=it):
-                    from dune.perftool.pdelab.signatures import assembly_routine_signature
-                    operator_kernels[(it, 'jacobian')] = [LoopyKernelMethod(assembly_routine_signature(), kernel=None)]
+            if get_option("generate_jacobians"):
+                for measure in set(i.integral_type() for i in jacform.integrals()):
+                    logger.info("generate_localoperator_kernels: measure {}".format(measure))
+                    with global_context(integral_type=measure):
+                        with global_context(kernel=assembler_routine_name()):
+                            kernel = [k for k in get_backend(interface="generate_kernels_per_integral")(jacform.integrals_by_type(measure))]
+                    operator_kernels[(measure, 'jacobian')] = kernel
+
+                # Generate dummy functions for those kernels, that vanished in the differentiation process
+                # We *could* solve this problem by using lambda_* terms but we do not really want that, so
+                # we use empty jacobian assembly methods instead
+                alpha_measures = set(i.integral_type() for i in form.integrals())
+                jacobian_measures = set(i.integral_type() for i in jacform.integrals())
+                for it in alpha_measures - jacobian_measures:
+                    with global_context(integral_type=it):
+                        from dune.perftool.pdelab.signatures import assembly_routine_signature
+                        operator_kernels[(it, 'jacobian')] = [LoopyKernelMethod(assembly_routine_signature(), kernel=None)]
 
         # Jacobian apply methods for matrix-free computations
         if get_option("matrix_free"):