Flux in and out does not equal.

5 months ago by
Hello There,
I am using Fenics to solve a simple Poisson's equation to compute steady temperature disstribution. Where my domain is like:

The red boundary is labelled as 1 and top green boundary is labelled as 2. They are both Dirichlet boundaries(fix temperature). Others are just Neumann boundaries where the flux is assumed as 0.

I solved the PDE and compute the heat flux through red boundary and green boundary.
Well, there are like 3% error between them no matter how I refined my mesh.
Theoretically they should equal to each other.
But if I compute the flux at the Neumman boundary(which should be 0 but not) and add it to the flux out, then the flux out and flux in equal to each other.
It seems that my Neumman boundary was not strictly applied.
P.S. the way I compute my flux on the boundary:
grad_T = project(grad(T), Q)      //where T is the solved temperature, Q is the VecterFunctionSpace
n = FacetNormal(mesh)
flux_in = assemble( dot(grad_T, n) * ds(1) )
flux_out = assemble( dot(grad_T, n) * ds(2) )
error = abs(flux_in - flux_out) / flux_out

Thank you!
Community: FEniCS Project
By the way, even up to 500,000 vertices the error is still big.
written 5 months ago by SICHENG SUN  

1 Answer

5 months ago by
It's hard to say exactly what might be the problem without going through the full code in detail. However, it's worth noting that the approach of directly evaluating a boundary flux (e.g., \(\nabla T\cdot\mathbf{n}\)) at a Dirichlet boundary is not really best practice for flux extraction from finite element solutions. A more accurate way to extract fluxes is discussed here:


and also in Exercise 8 of Section 2.12 of Tom Hughes's FEM book.  (The same basic reasoning applies to tensor-valued fluxes, like Cauchy stress.)  Here is a short FEniCS program to demonstrate the general idea: 

from dolfin import *
mesh = UnitSquareMesh(100,100)
n = FacetNormal(mesh)
V = FunctionSpace(mesh,"Lagrange",1)

# The Dirichlet boundary
class Bdry(SubDomain):
    def inside(self,x,on_boundary):
        return on_boundary
bcs = [DirichletBC(V,Constant(0.0),Bdry()),]
u = TrialFunction(V)
v = TestFunction(V)
f = Constant(1000.0)

# Residual of the Poisson equation
def R(u,v):
    return (inner(grad(u),grad(v)) - f*v)*dx
F = R(u,v)
a = lhs(F)
L = rhs(F)

# Solve and put solution in Function u
u = Function(V)

# Everything BUT the Dirichlet boundary
class AntiBdry(SubDomain):
    def inside(self,x,on_boundary):
        return x[0] > DOLFIN_EPS and x[0] < 1.0-DOLFIN_EPS \
            and x[1] > DOLFIN_EPS and x[1] < 1.0-DOLFIN_EPS

# Solve in the trace space for what Neumann data, i.e.,
#  h = \nabla u\cdot\mathbf{n}
# would have produced the solution from the Dirichlet problem; this is the
# notion of flux that satisfies the underlying conservation law.
h = TrialFunction(V)
antiBCs = [DirichletBC(V,Constant(0.0),AntiBdry()),]
FBdry = h*v*ds - R(u,v)
aBdry = lhs(FBdry)
LBdry = rhs(FBdry)
ABdry = assemble(aBdry,keep_diagonal=True)
bBdry = assemble(LBdry)
[bc.apply(ABdry,bBdry) for bc in antiBCs]
h = Function(V)

# compare integral of flux and integral of source term (compatibility).  I need
# a lot of elements for the direct flux to converge to the correct value that
# balances the volume source term.
print("\\int_\\Gamma h = "+str(assemble(h*ds)))
print("\\int_\\Gamma \\nabla u\cdot\mathbf{n} = "
print("\\int_\\Omega f = "+str(assemble((f+Constant(0.0)*u)*dx)))
Thanks so much for your wonderful response!
written 5 months ago by SICHENG SUN  
Please login to add an answer/comment or follow this question.

Similar posts:
Search »