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):
  t+=dt
  solve(a == L, u, bc)
  u_ps.assign(u)
  vtkfile_s<<(u,t)
  print('Integral is not zero:', assemble(dot(grad(u), velocity)*dx))