Should periodic boundary conditions with DG function space and ghosted mesh work?

Question

Hi!

I have heard rumors that the new support for ghosted meshes and DG in parallel should be enough to enable periodic BCs for DG, but the toy example below of linear advection in 1D seems to show that it is not so. Are there some more knobs to turn except parameters["ghost_mode"] to make periodic BCs work, or is this an expected failure?

Results for function spaces CG and DG are included below the code. This was tested with FEniCS 1.5.

FSPACE = 'DG'
vel = 0.1
dt = 0.01
tmax = 10
u0code = '(x[0] > 0.2 && x[0] < 0.4) ? pow(cos(16*(x[0]-0.3)), 2) : 0.0'

# ---------------------------------------
from dolfin import *
import numpy
from matplotlib import pyplot

import logging
logging.getLogger('FFC').setLevel(logging.WARNING)
logging.getLogger('UFL').setLevel(logging.WARNING)

# Make mesh ghosted to permit evaluation of DG terms
# ("shared_facet", "shared_vertex" or "none")
# http://fenicsproject.org/pipermail/fenics-support/2014-July/000778.html
parameters["ghost_mode"] = "shared_vertex"

class PeriodicDomain(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[0], 0)

    def map(self, x, y):
        y[0] = x[0] - 1.0

constrained_domain = PeriodicDomain()

mesh = UnitIntervalMesh(1000)
V = FunctionSpace(mesh, FSPACE, 1, constrained_domain=constrained_domain)
u = TrialFunction(V)
v = TestFunction(V)

uf = Function(V)  # New value
upf = Function(V) # Previous value

u0 = Expression(u0code)
project(u0, V, function=upf)

uc0 = Constant(vel)
u_conv = Constant([vel])
if FSPACE == 'CG':
    eq = (Constant(1/dt)*(u - upf)*v + uc0*u.dx(0)*v)*dx
elif FSPACE == 'DG':
    n = FacetNormal(mesh)
    flux_nU = u*(dot(u_conv, n) + abs(dot(u_conv, n)))/2
    flux = flux_nU('+') - flux_nU('-')
    eq = Constant(1/dt)*(u - upf)*v*dx - u*uc0*v.dx(0)*dx + flux*jump(v)*dS

a, L = lhs(eq), rhs(eq)

pyplot.figure(figsize=(10,6))
for t in numpy.arange(numpy.ceil(tmax/dt))*dt:
    solve(a==L, uf)
    upf.assign(uf)

    if near(t, round(t)):
        pyplot.plot(mesh.coordinates()[:,0], upf.compute_vertex_values(), label='t=%g' % t)

pyplot.ylim(-0.1, 1.1)
pyplot.legend(loc='best')
pyplot.show()

Result with FSPACE='CG':

Result with FSPACE='DG':