diff --git a/applications/CMakeLists.txt b/applications/CMakeLists.txt
index 5a75b84e1c0c3b4d4daa19381d7032489d2cbddb..d22155b0c2bc4a1c7e5cdacf2eeba609540e09c5 100644
--- a/applications/CMakeLists.txt
+++ b/applications/CMakeLists.txt
@@ -1,2 +1,3 @@
+add_subdirectory(conv_diff_reac)
 add_subdirectory(poisson_dg)
 add_subdirectory(poisson_dg_tensor)
diff --git a/applications/conv_diff_reac/CMakeLists.txt b/applications/conv_diff_reac/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..5966da0692584e01cf69f3043a00fd4738a0a6a9
--- /dev/null
+++ b/applications/conv_diff_reac/CMakeLists.txt
@@ -0,0 +1,11 @@
+dune_add_formcompiler_system_test(UFLFILE cdr_ernupwind.ufl
+                                  BASENAME app_cdr_ernupwind
+                                  INIFILE cdr_ernupwind.mini
+                                  NO_TESTS
+                                  )
+
+dune_add_formcompiler_system_test(UFLFILE conv_diff_dg_ernupwind.ufl
+                                  BASENAME verify_app_cdr_ernupwind
+                                  INIFILE verify_ernupwind.mini
+                                  NO_TESTS
+                                  )
diff --git a/applications/conv_diff_reac/cdr_condupwind.mini b/applications/conv_diff_reac/cdr_condupwind.mini
new file mode 100644
index 0000000000000000000000000000000000000000..c326e64ca611aec3a8b0badd83110c99ad1dad4b
--- /dev/null
+++ b/applications/conv_diff_reac/cdr_condupwind.mini
@@ -0,0 +1,51 @@
+__name = app_nonlinear_conv_diff_condupwind_{__exec_suffix}
+__exec_suffix = deg{formcompiler.ufl_variants.degree}_{opcount_suffix}_level{formcompiler.instrumentation_level}
+
+opcount_suffix = opcount, nonopcount | expand opcount
+{opcount_suffix} == opcount and {formcompiler.instrumentation_level} != 4 | exclude
+
+# Calculate the size of the grid to equlibritate it to 100 MB/rank
+# Input parameters
+dim = 3
+mbperrank = 100
+ranks = 16
+floatingbytes = 8
+
+# Metaini Calculations
+memperrank = {mbperrank} * 1048576 | eval
+dofsperdir = {formcompiler.ufl_variants.degree} + 1 | eval
+celldofs = {dofsperdir} ** {dim} | eval
+cellsperrank = {memperrank} / ({floatingbytes} * {celldofs}) | eval
+cellsperdir = {cellsperrank} ** (1/{dim}) | eval | toint
+firstdircells = {ranks} * {cellsperdir} | eval
+dimminusone = {dim} - 1 | eval
+ones = 1 | repeat {dimminusone}
+otherdircells = {cellsperdir} | repeat {dimminusone}
+
+# Setup the grid!
+extension = 1.0 | repeat {dim}
+cells = {firstdircells} {otherdircells}
+partitioning = {ranks} {ones}
+
+# Set up the timing identifier
+identifier = nonlinear_convdiff_condupwind_deg{formcompiler.ufl_variants.degree}
+
+[wrapper.vtkcompare]
+name = {__name}
+extension = vtu
+
+[formcompiler]
+fastdg = 1
+sumfact = 1
+vectorize_quad = 1
+vectorize_grads = 1
+instrumentation_level = 2, 3, 4 | expand
+opcounter = 1, 0 | expand opcount
+time_opcounter = 0, 1 | expand opcount
+exact_solution_expression = g
+compare_l2errorsquared = 1e-6
+quadrature_order = {formcompiler.ufl_variants.degree} * 2 | eval
+
+[formcompiler.ufl_variants]
+cell = hexahedron
+degree = 1, 2, 3, 4, 5, 6, 7, 8, 9, 10 | expand
diff --git a/applications/conv_diff_reac/cdr_condupwind.ufl b/applications/conv_diff_reac/cdr_condupwind.ufl
new file mode 100644
index 0000000000000000000000000000000000000000..0071766cb4bff3bb787879b29c555fc6743cd89b
--- /dev/null
+++ b/applications/conv_diff_reac/cdr_condupwind.ufl
@@ -0,0 +1,37 @@
+dim = 3
+x = SpatialCoordinate(cell)
+
+g = x[0]*x[0] + x[1]*x[1] + x[2]*x[2]
+I = Identity(3)
+A = as_matrix([[x[i]*x[j] + I[i,j] for j in range(3)] for i in range(3)])
+f = -6.
+c = 5.
+b = x
+
+V = FiniteElement("DG", cell, degree)
+
+u = TrialFunction(V)
+v = TestFunction(V)
+
+n = FacetNormal(cell)('+')
+
+alpha = 3.0
+h_ext = CellVolume(cell) / FacetArea(cell)
+cse_gamma_ext = (alpha * degree * (degree + dim - 1)) / h_ext
+h_int = Min(CellVolume(cell)('+'), CellVolume(cell)('-')) / FacetArea(cell)
+cse_gamma_int = (alpha * degree * (degree + dim - 1)) / h_int
+upwind = conditional(inner(b, n) > 0, u('+'), u('-'))
+
+theta = -1.0
+
+r = (inner(A*grad(u) - b*u, grad(v)) + (c*u-f)*v)*dx \
+  + (inner(n, A*avg(grad(u))) + inner(b*upwind, n))*jump(v)*dS \
+  + cse_gamma_int*jump(u)*jump(v)*dS \
+  - theta*jump(u)*inner(A*avg(grad(v)), n)*dS \
+  - inner(n, A*grad(u) + b*u)*v*ds \
+  + cse_gamma_ext*u*v*ds \
+  + theta*u*inner(A*grad(v), n)*ds \
+  - theta*g*inner(A*grad(v), n)*ds \
+  - cse_gamma_ext*g*v*ds
+
+forms = [r]
diff --git a/applications/conv_diff_reac/cdr_ernupwind.mini b/applications/conv_diff_reac/cdr_ernupwind.mini
new file mode 100644
index 0000000000000000000000000000000000000000..40638c1159d4747e2cf788c32fd45e11413e02ae
--- /dev/null
+++ b/applications/conv_diff_reac/cdr_ernupwind.mini
@@ -0,0 +1,51 @@
+__name = app_nonlinear_conv_diff_ernupwind_{__exec_suffix}
+__exec_suffix = deg{formcompiler.ufl_variants.degree}_{opcount_suffix}_level{formcompiler.instrumentation_level}
+
+opcount_suffix = opcount, nonopcount | expand opcount
+{opcount_suffix} == opcount and {formcompiler.instrumentation_level} != 4 | exclude
+
+# Calculate the size of the grid to equlibritate it to 100 MB/rank
+# Input parameters
+dim = 3
+mbperrank = 100
+ranks = 16
+floatingbytes = 8
+
+# Metaini Calculations
+memperrank = {mbperrank} * 1048576 | eval
+dofsperdir = {formcompiler.ufl_variants.degree} + 1 | eval
+celldofs = {dofsperdir} ** {dim} | eval
+cellsperrank = {memperrank} / ({floatingbytes} * {celldofs}) | eval
+cellsperdir = {cellsperrank} ** (1/{dim}) | eval | toint
+firstdircells = {ranks} * {cellsperdir} | eval
+dimminusone = {dim} - 1 | eval
+ones = 1 | repeat {dimminusone}
+otherdircells = {cellsperdir} | repeat {dimminusone}
+
+# Setup the grid!
+extension = 1.0 | repeat {dim}
+cells = {firstdircells} {otherdircells}
+partitioning = {ranks} {ones}
+
+# Set up the timing identifier
+identifier = nonlinear_convdiff_ernupwind_deg{formcompiler.ufl_variants.degree}
+
+[wrapper.vtkcompare]
+name = {__name}
+extension = vtu
+
+[formcompiler]
+fastdg = 1
+sumfact = 1
+vectorize_quad = 1
+vectorize_grads = 1
+instrumentation_level = 2, 3, 4 | expand
+opcounter = 1, 0 | expand opcount
+time_opcounter = 0, 1 | expand opcount
+exact_solution_expression = g
+compare_l2errorsquared = 1e-6
+quadrature_order = {formcompiler.ufl_variants.degree} * 2 | eval
+
+[formcompiler.ufl_variants]
+cell = hexahedron
+degree = 1, 2, 3, 4, 5, 6, 7, 8, 9, 10 | expand
diff --git a/applications/conv_diff_reac/cdr_ernupwind.ufl b/applications/conv_diff_reac/cdr_ernupwind.ufl
new file mode 100644
index 0000000000000000000000000000000000000000..6e1763f521cf7fdfb8b22bf6730f0eb732eb811f
--- /dev/null
+++ b/applications/conv_diff_reac/cdr_ernupwind.ufl
@@ -0,0 +1,38 @@
+dim = 3
+x = SpatialCoordinate(cell)
+
+g = x[0]*x[0] + x[1]*x[1] + x[2]*x[2]
+I = Identity(3)
+A = as_matrix([[x[i]*x[j] + I[i,j] for j in range(3)] for i in range(3)])
+f = -6.
+c = 5.
+b = x
+
+V = FiniteElement("DG", cell, degree)
+
+u = TrialFunction(V)
+v = TestFunction(V)
+
+n = FacetNormal(cell)('+')
+
+alpha = 3.0
+h_ext = CellVolume(cell) / FacetArea(cell)
+cse_gamma_ext = (alpha * degree * (degree + dim - 1)) / h_ext
+h_int = Min(CellVolume(cell)('+'), CellVolume(cell)('-')) / FacetArea(cell)
+cse_gamma_int = (alpha * degree * (degree + dim - 1)) / h_int
+
+theta = -1.0
+
+r = (inner(A*grad(u) - b*u, grad(v)) + (c*u-f)*v)*dx \
+  + inner(n, A*avg(grad(u)))*jump(v)*dS \
+  - (avg(u) - 0.5*jump(u))*inner(b, n)*jump(v)*dS \
+  + cse_gamma_int*jump(u)*jump(v)*dS \
+  - theta*jump(u)*inner(A*avg(grad(v)), n)*dS \
+  - inner(n, A*grad(u) + b*g)*v*ds \
+  + cse_gamma_ext*u*v*ds \
+  + theta*u*inner(A*grad(v), n)*ds \
+  - theta*g*inner(A*grad(v), n)*ds \
+  - cse_gamma_ext*g*v*ds
+
+
+forms = [r]
diff --git a/applications/conv_diff_reac/ern_verify.mini b/applications/conv_diff_reac/ern_verify.mini
new file mode 100644
index 0000000000000000000000000000000000000000..531b23ad761ea68ae69b5677c4c1fc0758c2de5c
--- /dev/null
+++ b/applications/conv_diff_reac/ern_verify.mini
@@ -0,0 +1,21 @@
+__name = verify_app_cdr_ernupwind
+
+# Setup the grid!
+extension = 1.0 1.0 1.0
+cells = 8 8 8
+
+[wrapper.vtkcompare]
+name = {__name}
+extension = vtu
+
+[formcompiler]
+fastdg = 1
+sumfact = 1
+vectorize_quad = 1
+vectorize_grads = 1
+exact_solution_expression = g
+compare_l2errorsquared = 1e-6
+
+[formcompiler.ufl_variants]
+cell = hexahedron
+degree = 2