diff --git a/tests/regression/test_upwind_flux.py b/tests/regression/test_upwind_flux.py
new file mode 100644
index 0000000000..9cce699345
--- /dev/null
+++ b/tests/regression/test_upwind_flux.py
@@ -0,0 +1,98 @@
+"""
+This code does the following.
+
+First, obtains an upwind DG0 approximation to div(u*D).
+Then, tries to find a BDM1 flux F such that div F is
+equal to this upwind approximation.
+
+we have
+
+\int_e phi D_1 dx = -\int_e grad phi . u D dx
+                    + \int_{\partial e} phi u.n \tilde{D} ds
+
+where \tilde{D} is the value of D on the upwind face. For
+DG0, grad phi = 0.
+
+Then, if we define F such that
+
+\int_f phi F.n ds = \int_f phi u.n \tilde{D} ds
+
+then
+
+\int_e phi div(F) ds = \int_{\partial e} phi u.n \tilde{D} ds
+as required.
+"""
+from firedrake import *
+import pytest
+
+
+def run_test():
+    mesh = UnitIcosahedralSphereMesh(refinement_level=2)
+    global_normal = Expression(("x[0]/sqrt(x[0]*x[0]+x[1]*x[1]+x[2]*x[2])",
+                                "x[1]/sqrt(x[0]*x[0]+x[1]*x[1]+x[2]*x[2])",
+                                "x[2]/sqrt(x[0]*x[0]+x[1]*x[1]+x[2]*x[2])"))
+    mesh.init_cell_orientations(global_normal)
+
+    # Define function spaces and basis functions
+    V_dg = FunctionSpace(mesh, "DG", 0)
+    M = FunctionSpace(mesh, "RT", 1)
+
+    # advecting velocity
+    u0 = Expression(('-x[1]', 'x[0]', '0'))
+    u = Function(M).project(u0)
+
+    # Mesh-related functions
+    n = FacetNormal(mesh)
+
+    # ( dot(v, n) + |dot(v, n)| )/2.0
+    un = 0.5*(dot(u, n) + abs(dot(u, n)))
+
+    # D advection equation
+    phi = TestFunction(V_dg)
+    D = TrialFunction(V_dg)
+    a_mass = phi*D*dx
+    a_int = dot(grad(phi), -u*D)*dx
+    a_flux = (dot(jump(phi), un('+')*D('+') - un('-')*D('-')))*dS
+
+    arhs = (a_int + a_flux)
+
+    D1 = Function(V_dg)
+
+    D0 = Expression("exp(-pow(x[2],2) - pow(x[1],2))")
+    D = Function(V_dg).interpolate(D0)
+
+    D1problem = LinearVariationalProblem(a_mass, action(arhs, D), D1)
+    D1solver = LinearVariationalSolver(D1problem)
+    D1solver.solve()
+
+    # Surface Flux equation
+    V1 = FunctionSpace(mesh, "RT", 1)
+    w = TestFunction(V1)
+    Ft = TrialFunction(V1)
+    Fs = Function(V1)
+
+    aFs = (inner(w('+'), n('+'))*inner(Ft('+'), n('+')) +
+           inner(w('-'), n('-'))*inner(Ft('-'), n('-')))*dS
+    LFs = 2.0*(inner(w('+'), n('+'))*un('+')*D('+')
+               + inner(w('-'), n('-'))*un('-')*D('-'))*dS
+
+    Fsproblem = LinearVariationalProblem(aFs, LFs, Fs)
+    Fssolver = LinearVariationalSolver(Fsproblem,
+                                       solver_parameters={'ksp_type': 'preonly'})
+    Fssolver.solve()
+
+    divFs = Function(V_dg)
+
+    solve(a_mass == phi*div(Fs)*dx, divFs)
+
+    assert errornorm(divFs, D1, degree_rise=0) < 1e-12
+
+
+def test_upwind_flux():
+    run_test()
+
+
+@pytest.mark.xfail(reason="No support for interior facet integrals in MPI")
+@pytest.mark.parallel
+def test_upwind_flux_parallel():
+    run_test()