# Flux integral is different from flux in the Form!

**URL:** <https://fenicsproject.discourse.group/t/flux-integral-is-different-from-flux-in-the-form/8131>\
**Category:** variational formulation\
**Created:** [April 24, 2022, 10:14am UTC](https://fenicsproject.discourse.group/t/flux-integral-is-different-from-flux-in-the-form/8131 "2022-04-24T10:14:23Z")\
**Posts on this page:** 1\
**Page:** 1

<div class="post-metadata">

**Author:** ![ilia145](https://avatars.discourse-cdn.com/v4/letter/i/ee7513/32.png) [@ilia145](https://fenicsproject.discourse.group/u/ilia145)\
**Post date:** [April 24, 2022, 10:14am UTC](https://fenicsproject.discourse.group/t/flux-integral-is-different-from-flux-in-the-form/8131/1 "2022-04-24T10:14:23Z")

</div>

Hi  
I have a simple question in a time-dependent diffusion problem.  
Why flux integration on the interior domain has different result from that defined in the From equation as a constant flux? This difference is such that in through time steps, computed flux reaches a constant (~0.17 in this setting) but different from Neumann boundary condition (0.2) that specified in the right hand side Form. Also, computed flux trend is affected by diffusion constant.

The complete code is

```auto
from fenics import *
from mshr import *
import matplotlib.pyplot as plt

circle1= Circle(Point(30, 30), 30)
circle2= Circle(Point(30, 30), 15)
Domain= circle1 - circle2
mesh= generate_mesh(Domain, 15)

# Neumann b.c. will defined in inner boundary
B= CompiledSubDomain('pow(x[0]-30, 2) + pow(x[1]-30, 2) < pow(15+0.5, 2)')
B_Marker= MeshFunction('size_t', mesh, 1)
B.mark(B_Marker, 1)
ds= Measure('ds', domain= mesh, subdomain_data= B_Marker)

V= FunctionSpace(mesh, 'P', 1)
u0= interpolate(Constant(1), V)
u= TrialFunction(V)
v= TestFunction(V)
D= Constant(1)
g= Constant(2)
dt= Constant(0.1)
a= u*v*dx + dt*D*dot(grad(u), grad(v))*dx
L= u0*v*dx + dt*g*v*ds(1)
n= FacetNormal(mesh)
flux= dt*D*dot(grad(u0), n)*ds(1)
DS= assemble(1*ds(1))

A= assemble(a)
u= Function(V)
flux_array= []
for i in range(200):
    b= assemble(L)
    solve(A, u.vector(), b)
    u0.assign(u)
    flux_array.append(assemble(flux)/DS)
plt.plot(flux_array)

```