diff --git a/test/sumfact/hyperbolic/CMakeLists.txt b/test/sumfact/hyperbolic/CMakeLists.txt
index 8399811dfe665e2b7d096c4d2ab6800a07bce657..25388b32023e5f20a5439212b4a1c2a92a525fad 100644
--- a/test/sumfact/hyperbolic/CMakeLists.txt
+++ b/test/sumfact/hyperbolic/CMakeLists.txt
@@ -2,3 +2,8 @@ dune_add_formcompiler_system_test(UFLFILE lineartransport.ufl
                                   BASENAME sumfact_lineartransport
                                   INIFILE lineartransport.mini
                                   )
+
+dune_add_formcompiler_system_test(UFLFILE linearacoustics.ufl
+                                  BASENAME sumfact_linearacoustics
+                                  INIFILE linearacoustics.mini
+                                  )
diff --git a/test/sumfact/hyperbolic/linearacoustics.mini b/test/sumfact/hyperbolic/linearacoustics.mini
new file mode 100644
index 0000000000000000000000000000000000000000..9ac464c09493274efcad320ac2e4c969b59447f0
--- /dev/null
+++ b/test/sumfact/hyperbolic/linearacoustics.mini
@@ -0,0 +1,25 @@
+__name = sumfact_linearacoustics_{__exec_suffix}
+__exec_suffix = {diff_suffix}
+
+diff_suffix = numdiff, symdiff | expand diff
+
+extension = 1.0 1.0
+cells = 20 20
+elementType = hexahedral
+
+[instat]
+T = 1.5
+dt = 0.0025
+
+[wrapper.vtkcompare]
+name = {__name}
+extension = vtu
+
+[formcompiler]
+numerical_jacobian = 1, 0 | expand diff
+explicit_time_stepping = 0
+sumfact = 1
+# This tests that all mass is transported out of the domain.
+# While this is not the best of tests, it is something easily checked for.
+#exact_solution_expression = g
+#compare_l2errorsquared = 1e-10
diff --git a/test/sumfact/hyperbolic/linearacoustics.ufl b/test/sumfact/hyperbolic/linearacoustics.ufl
new file mode 100644
index 0000000000000000000000000000000000000000..09386752df35e0593693feefe96fbcb29b775d2d
--- /dev/null
+++ b/test/sumfact/hyperbolic/linearacoustics.ufl
@@ -0,0 +1,24 @@
+cell = quadrilateral
+x = SpatialCoordinate(cell)
+
+#V = FiniteElement("DG", cell, 1)
+
+MV = VectorElement("DG", cell=cell, degree=1, dim=3)
+
+n = FacetNormal(cell)('+')
+
+u = TrialFunction(MV)
+v = TestFunction(MV)
+
+rho, q0, q1 = split(u)
+rhot, q0t, q1t = split(v)
+
+mass = inner(u, v)*dx
+
+flux = as_matrix([[q0, q1], [rho, 0.], [0., rho]])
+
+r = -1.*inner(flux, grad(v))*dx \
+  + inner(dot(avg(flux), n) + 0.5*jump(u), jump(v))*dS \
+  + inner(u, v)*ds
+
+forms = [mass, r]