Time-dependent Neumann problem does not conserve mass
I try to solve a time-dependent problem (using a backward Euler scheme) with Neumann-Neumann boundary conditions
a = (c1*testfunction*dx - c0*testfunction*dx) # standard for backward Euler time stepping L = dt*(inner(flux_c, grad(testfunction))*dx # standard part of continuity equation - j1*testfunction*ds_(id1) # interface 1 + j2*testfunction*ds_(id2)) # interface 2 F = a-L
c1 is the new solution, c0 the old one
dt the time step
flux_c a user-defined function that calculates the flux in the volume
id_1 and id_2 mark the facets for the two interfaces
j1 and j2 are the imposed flux (expressed in something per square meter)
The imposed flux is balanced, the two terms below have the exact same value:
integral_j1 = assemble(j1*ds_(id1)) integral_j2 = assemble(j2*ds_(id2))
So, in this case, since exit flux = entry flux, i should keep the term below constant :
integral_c = assemble(c1*dx)
Although, it's not the case, i get a significant deviation of the mass within the system (few percents of the mean).
What could induce such deviation?