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
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)