Why integral of zero function is not zero?

Hi! I’m solving simple parabolic equation and tried to add zero convection to it. However, I see that integral of convective term is not zero. Moreover the value may change from run to run. I’m using Fenics 2019.1.0 under Ubuntu for Windows 10. Could you please clarify what I misunderstand? Thank you!

from dolfin import *

vtkfile_s=File("solution.pvd")

num_steps=2
t=0
T=1.0
dt=T/num_steps

x0=y0=z0=0.0
xe=ye=ze=1.0

mesh = BoxMesh(Point(x0, y0, z0), Point(xe, ye, ze), num_steps, num_steps, num_steps)

V = FunctionSpace(mesh, 'CG', 1)

subd= MeshFunction('size_t',mesh,mesh.topology().dim(),1)
dx=Measure("dx",domain=mesh,subdomain_data=subd)

def boundary_d(x, on_boundary):
  return (near (x[0], x0) or near(x[0], xe)) and on_boundary

u_D = Constant('1.0e0')
bc = DirichletBC(V, u_D, boundary_d)

u_ps=interpolate(u_D,V)
u = TrialFunction(V)
v = TestFunction(V)
f = Constant('1.0e0')

class Velocity(UserExpression):
      def \_\_init\_\_(self, mesh, **kwargs):
            self.mesh = mesh
            super().\_\_init\_\_(**kwargs)
      def eval_cell(self, values, x, ufc_cell):
            values=[0.0,0.0,0.0]
            return   
      def value_shape(self):
            return (3,)
velocity=Velocity(mesh)

F=(u-u_ps)\*v\*dx + dt\*( dot(grad(u), grad(v)) )\*dx - dt\*( f\*v )\*dx + dt\*( dot(grad(u),velocity)\*v )\*dx 
a, L = lhs(F), rhs(F)

u = Function(V)
u.assign(u_ps)
vtkfile_s<<(u,t)

print('Integral is not zero:', assemble(dot(grad(u), velocity)*dx)) 
    
for i in range(num_steps):
&emsp;     t+=dt
&emsp;     solve(a == L, u, bc)
&emsp;     u_ps.assign(u)
&emsp;     vtkfile_s<<(u,t)
&emsp;     print('Integral is not zero:', assemble(dot(grad(u), velocity)*dx))

This is due to the fact that your are not overwriting the input-values in eval cell.
Changing velocity as follows solves your problem:

class Velocity(UserExpression):
    def __init__(self, mesh, **kwargs):
        self.mesh = mesh
        super().__init__(**kwargs)
    def eval_cell(self, values, x, ufc_cell):
        values[:]=[0,0,0]
        return
    def value_shape(self):
        return (3,)

Thank you very much!