### Time-dependent Neumann problem does not conserve mass

166
views
0
5 months ago by
Hello everyone,

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?
Any idea?

Best,
Community: FEniCS Project
1
I agree, this is very puzzling; as long as your test space contains testfunction identically equal to 1 (implying grad(testfunction) = 0), then, regardless of what's going on in flux_c, the integral of (c1 - c0) should exactly balance the boundary term, which you already verified was zero.  The only ways I can see this logic breaking down are

1.  There is a strongly-enforced Dirichlet BC on some third part of the boundary, so testfunction must go to zero somewhere, and testfunction=1 is not a valid choice.

2.  The algebraic system of equations is not being solved completely (e.g., due to an iterative solver with nonzero tolerance).
written 5 months ago by David Kamensky