Adding vanishing boundary terms changes solution

Question

Hello, I've think I've got a general question about the FEM method. When transforming a second order PDE into the weak form we get the on term where the source is integrated over the whole volume and by partial integration and Gauss' law another term corresponding to the Neumann boundary condition.

$$ [...] + \int f dV - \int g dS $$

Both f and g are zero in my problem, so I thought I didn't have to add them to the formulation.
But after adding them, while setting both f and g to zero (Constant('0.0')) I get different results than from before. Is there an explanation for that?

Comments

Hi, I think you should post a minimal example illustrating your problem. It seems unlikely that adding the zero function could have any effect on any equation :/

---

Here's the code I'm using.

from fenics import *
import os
import numpy as np

parameters["form_compiler"]["optimize"]     = True
parameters["form_compiler"]["cpp_optimize"] = True

T = 4.0
Nt = 300
dt = T/Nt
L = 8.0
Nx = 200


mesh = IntervalMesh(Nx, -L, L)
V = FunctionSpace(mesh, 'CG', 1)

u_D = Constant(0.0)

def boundary(x, on_boundary):
    return on_boundary

bc = DirichletBC(V, u_D, boundary)

gamma = 0.281
sigma = gamma/np.sqrt(8*np.log(2.0))
y0 = 2.915;

ic = Expression('1/sqrt(8*pi*pow(sigma, 2.0)) * exp(-pow(x[0]-y0, 2)) / (2*pow(sigma, 2.0)) \
+ 1/sqrt(8*pi*pow(sigma, 2.0)) * exp(-pow(x[0]+y0, 2)) / (2*pow(sigma, 2.0))',
                 y0=y0,sigma=sigma,pi=np.pi,degree=2)

v = TestFunction(V)
u0 = Function(V)
u = Function(V)
f = Constant(0.0)
g = Constant(0.0)

name = "diffusion"
q = 1.25
I = Expression(('(2.0-q)*sinh(x[0])',),q=q,degree=2)
factor = 160.0/1150.0

def operator(u, v):
    return dot(grad(pow(u,2.0-q)),grad(v))*dx 
              + u*dot(I,grad(v))*dx 
              - g*v*ds 
              + f*v*dx


F = inner(u-u0, v)*dx + dt*operator(u, v)

# normalize to one. integrate over the whole domain
u0.interpolate(ic)
S = assemble(u0*dx(mesh))
u0.vector()[:] /= S
u.assign(u0)

solver = NewtonSolver()
solver.parameters["linear_solver"] = "lu"
solver.parameters["convergence_criterion"] = "incremental"
solver.parameters["relative_tolerance"] = 1e-6

solve(F == 0, u, bc)

vtkfile = File('diffusion/diffusion.pvd')
t = 0.0
u.rename('u', 'density')


vtkfile << (u, t)

while t < T:
    solve(F == 0, u, bc)
    vtkfile << (u, t)

    u0.assign(u)
    t += dt

---

Ok, great. Can you specify how the two solutions differ, i.e. how large is the difference?

Thanks!