diff --git a/python/dune/perftool/pdelab/driver/instationary.py b/python/dune/perftool/pdelab/driver/instationary.py
index 1d4b9e87be8c0e5f55a55556c2a2abb3ec7c62d4..173c1b6e7fa237ff4b422942408b9cdb02ed12a0 100644
--- a/python/dune/perftool/pdelab/driver/instationary.py
+++ b/python/dune/perftool/pdelab/driver/instationary.py
@@ -36,21 +36,11 @@ from dune.perftool.options import get_option
 
 
 def solve_instationary():
-    # Test if form is linear in ansatzfunction
-    linear = is_linear()
-
-    # Test wether we want to do matrix free operator evaluation
-    matrix_free = get_option('matrix_free')
-
     # Create time loop
-    if linear and matrix_free:
-        assert False
-    elif linear and not matrix_free:
+    if get_option('matrix_free'):
+        raise NotImplementedError("Instationary matrix free not implemented!")
+    else:
         time_loop()
-    if not linear and matrix_free:
-        assert False
-    elif not linear and not matrix_free:
-        assert False
 
     print_residual()
     print_matrix()
diff --git a/test/hyperbolic/CMakeLists.txt b/test/hyperbolic/CMakeLists.txt
index 53148a9ed8881e9252bd687222c91795ceab14fd..2b00da61d3925c4d7e514136fdf0c7af4c169286 100644
--- a/test/hyperbolic/CMakeLists.txt
+++ b/test/hyperbolic/CMakeLists.txt
@@ -7,3 +7,8 @@ dune_add_formcompiler_system_test(UFLFILE linearacoustics.ufl
                                   BASENAME linearacoustics
                                   INIFILE linearacoustics.mini
                                   )
+
+dune_add_formcompiler_system_test(UFLFILE shallowwater.ufl
+                                  BASENAME shallowwater
+                                  INIFILE shallowwater.mini
+                                  )
diff --git a/test/hyperbolic/shallowwater.mini b/test/hyperbolic/shallowwater.mini
new file mode 100644
index 0000000000000000000000000000000000000000..374b2e89a41010a46817fee8ded947555e35354d
--- /dev/null
+++ b/test/hyperbolic/shallowwater.mini
@@ -0,0 +1,22 @@
+__name = shallowwater_{__exec_suffix}
+__exec_suffix = {diff_suffix}
+
+diff_suffix = numdiff, symdiff | expand diff
+
+extension = 10.0 10.0
+upperright = 10.0 10.0
+cells = 20 20
+elementType = simplical
+elements = 20 20
+
+[instat]
+T = 0.5
+dt = 0.005
+
+[wrapper.vtkcompare]
+name = {__name}
+extension = vtu
+
+[formcompiler]
+numerical_jacobian = 1, 0 | expand diff
+explicit_time_stepping = 1
diff --git a/test/hyperbolic/shallowwater.ufl b/test/hyperbolic/shallowwater.ufl
new file mode 100644
index 0000000000000000000000000000000000000000..9fca496e9fab29d58e91e79f15fb6fed3b4265af
--- /dev/null
+++ b/test/hyperbolic/shallowwater.ufl
@@ -0,0 +1,39 @@
+cell = quadrilateral
+x = SpatialCoordinate(cell)
+
+f = 0.2 * exp(-(pow((x[0] - 5), 2) / 2 + pow(x[1] - 5, 2) / 2)) + 1.
+
+V = FiniteElement("DG", cell, 1)
+MV = MixedElement(V, V, V)
+
+g = 10.0
+
+n = FacetNormal(cell)('+')
+
+u = TrialFunction(MV)
+v = TestFunction(MV)
+
+mass = inner(u, v)*dx
+
+h, q0, q1 = split(u)
+
+flux = as_matrix([[q0,                  q1],
+                  [h*q0*q0 + 0.5*g*h*h, h*q0*q1],
+                  [h*q0*q1,             h*q1*q1 + 0.5*g*h*h]])
+
+bflux = as_matrix([[-q0,                  -q1],
+                  [h*q0*q0 + 0.5*g*h*h, h*q0*q1],
+                  [h*q0*q1,             h*q1*q1 + 0.5*g*h*h]])
+
+
+# Define numerical fluxes to choose from
+llf_flux = dot(avg(flux), n) - 0.5*jump(u)
+boundary_flux = 0.5*dot(flux + bflux, n) + as_vector([0., q0, q1])
+numerical_flux = llf_flux
+
+r = -1. * inner(flux, grad(v))*dx \
+  - inner(numerical_flux, jump(v))*dS \
+  + inner(boundary_flux, v)*ds
+
+forms = [mass, r]
+dirichlet_expression = f, 0.0, 0.0
diff --git a/test/sumfact/hyperbolic/CMakeLists.txt b/test/sumfact/hyperbolic/CMakeLists.txt
index 25388b32023e5f20a5439212b4a1c2a92a525fad..bca871950a0a01a07f20c8663ddadc7627dff6ae 100644
--- a/test/sumfact/hyperbolic/CMakeLists.txt
+++ b/test/sumfact/hyperbolic/CMakeLists.txt
@@ -7,3 +7,8 @@ dune_add_formcompiler_system_test(UFLFILE linearacoustics.ufl
                                   BASENAME sumfact_linearacoustics
                                   INIFILE linearacoustics.mini
                                   )
+
+dune_add_formcompiler_system_test(UFLFILE shallowwater.ufl
+                                  BASENAME sumfact_shallowwater
+                                  INIFILE shallowwater.mini
+                                  )
diff --git a/test/sumfact/hyperbolic/shallowwater.mini b/test/sumfact/hyperbolic/shallowwater.mini
new file mode 100644
index 0000000000000000000000000000000000000000..8346ebf34a742fcd8afa4c19dec08701c4f81fbe
--- /dev/null
+++ b/test/sumfact/hyperbolic/shallowwater.mini
@@ -0,0 +1,17 @@
+__name = sumfact_shallowwater
+
+extension = 10.0 10.0
+cells = 20 20
+elementType = hexahedral
+
+[instat]
+T = 0.5
+dt = 0.005
+
+[wrapper.vtkcompare]
+name = {__name}
+extension = vtu
+
+[formcompiler]
+explicit_time_stepping = 1
+sumfact = 1
diff --git a/test/sumfact/hyperbolic/shallowwater.ufl b/test/sumfact/hyperbolic/shallowwater.ufl
new file mode 100644
index 0000000000000000000000000000000000000000..9fca496e9fab29d58e91e79f15fb6fed3b4265af
--- /dev/null
+++ b/test/sumfact/hyperbolic/shallowwater.ufl
@@ -0,0 +1,39 @@
+cell = quadrilateral
+x = SpatialCoordinate(cell)
+
+f = 0.2 * exp(-(pow((x[0] - 5), 2) / 2 + pow(x[1] - 5, 2) / 2)) + 1.
+
+V = FiniteElement("DG", cell, 1)
+MV = MixedElement(V, V, V)
+
+g = 10.0
+
+n = FacetNormal(cell)('+')
+
+u = TrialFunction(MV)
+v = TestFunction(MV)
+
+mass = inner(u, v)*dx
+
+h, q0, q1 = split(u)
+
+flux = as_matrix([[q0,                  q1],
+                  [h*q0*q0 + 0.5*g*h*h, h*q0*q1],
+                  [h*q0*q1,             h*q1*q1 + 0.5*g*h*h]])
+
+bflux = as_matrix([[-q0,                  -q1],
+                  [h*q0*q0 + 0.5*g*h*h, h*q0*q1],
+                  [h*q0*q1,             h*q1*q1 + 0.5*g*h*h]])
+
+
+# Define numerical fluxes to choose from
+llf_flux = dot(avg(flux), n) - 0.5*jump(u)
+boundary_flux = 0.5*dot(flux + bflux, n) + as_vector([0., q0, q1])
+numerical_flux = llf_flux
+
+r = -1. * inner(flux, grad(v))*dx \
+  - inner(numerical_flux, jump(v))*dS \
+  + inner(boundary_flux, v)*ds
+
+forms = [mass, r]
+dirichlet_expression = f, 0.0, 0.0