Assembler and SystemAssembler give different results

Code to reproduce:

from dolfin import *

# Create mesh and define function space
mesh = UnitSquare(24, 24)
V = FunctionSpace(mesh, "DG", 1)

# Define test and trial functions
v = TestFunction(V)
u = TrialFunction(V)

# Define normal component, mesh size and right-hand side
n = V.cell().n
h = CellSize(mesh)
h_avg = (h('+') + h('-'))/2
f = Function(V, "500.0*exp(-(pow(x[0] - 0.5, 2) + pow(x[1] - 0.5, 2)) / 0.02)")

# Define parameters
alpha = 4.0
gamma = 8.0

# Define bilinear form
a = dot(grad(v), grad(u))*dx \
   - dot(avg(grad(v)), jump(u, n))*dS \
   - dot(jump(v, n), avg(grad(u)))*dS \
   + alpha/h_avg*dot(jump(v, n), jump(u, n))*dS \
   - dot(grad(v), u*n)*ds \
   - dot(v*n, grad(u))*ds \
   + gamma/h*v*u*ds

# Define linear form
L = v*f*dx

# Compute solution
problem = VariationalProblem(a, L)
u = problem.solve()

_problem = VariationalProblem(a, L)
_problem.parameters["symmetric"] = True
_u = _problem.solve()

print "Assembler result: ", u.vector().norm("l2")
print "SystemAssembler result: ", _u.vector().norm("l2")

Looks like problem could be with RHS vector.