Applying DirichletBC to DG0 in 1D burgers equation not working

Question

I want to solve a simple Burgers like equation in a 1D domain.

$\frac{\partial u}{\partial t} + 0.001\frac{\partial(u(1-u/0.2)}{\partial x} = 0$

I am using DG0 and implicit Euler. I want to impose a Dirichlet BC of 0.1 at the left side. This seems to only work if I switch to DG1, which it is unstable for my problem. So there are several questions about this:

-How can I implement DirichletBC for DG 0?
-Is a simple upwind scheme appropiate for this nonlinear equation?
-How can I stabilize DG1?

I tried to make the code as simple as possible:

from dolfin import *
import numpy as np

# Parameters
t_end = 100
dt = 0.01
nel = 200

mesh = UnitIntervalMesh(nel)
tdim = mesh.topology().dim()


V = FunctionSpace(mesh, "DG", 0)
initial_condition = Expression("0.1*x[0]+0.1", degree=1, domain=mesh)

class LeftBoundary(SubDomain):
    def inside(self, x, on_boundary):
        return x[0]<DOLFIN_EPS

bc_left = DirichletBC(V, Constant(0.1), LeftBoundary(), method="pointwise")

# Build the measure
measure_boundary = FacetFunction("size_t", mesh)
LeftBoundary().mark(measure_boundary, 1)
ds = ds(subdomain_data=measure_boundary)


n = FacetNormal(mesh)

v = TestFunction(V)
u = TrialFunction(V)

u0 = interpolate(initial_condition,V )


def a(u,v) :
    # Flux
    b = u*(1-u/0.2)*Constant((0.001,))
    # Upwind
    bn = (dot(b,n) + abs(dot(b, n)))/2.0
    # Bilinear form
    a_int = dot(grad(v), - b*u)*dx
    a_vel = dot(jump(v),( bn('+') - bn('-')  ))*dS + dot(v, bn*u)*ds(1)

    a = a_int + a_vel
    return a



# Define variational forms
a0=a(u0,v)
a1=a(u,v)

theta = Constant(1.0)
A = (1/dt)*inner(u, v)*dx - (1/dt)*inner(u0,v)*dx + theta*a1 + (1-theta)*a0

F = A
u_ = Function(V)
F = action(F,u_)    
J = derivative(F, u_, u)

# Create files for storing results
file_output = File("results/u.pvd")


u_.assign(u0)
u_.rename("u", "u")

# Time-stepping
t = 0.0

file_output << u_

while t < t_end:

    print "t =", t, "end t=", t_end

    # Compute
    solve(F==0, u_, bc_left, J=J)
    plot(u_)
    # Save to file
    file_output << u_

    # Move to next time step
    u0.assign(u_)
    t += dt