How to integrate normal flux on boundary

7 weeks ago by

The following implements a natural convection model applied to a heat-driven cavity benchmark in FEniCS 2017.2.0, and solves until an approximately steady state:

import fenics

""" Solution space """
mesh = fenics.UnitSquareMesh(20, 20)

P1, P2 = fenics.FiniteElement('P', mesh.ufl_cell(), 1), fenics.VectorElement('P', mesh.ufl_cell(), 2)

mixed_element = fenics.MixedElement([P1, P2, P1])

W = fenics.FunctionSpace(mesh, mixed_element)

psi_p, psi_u, psi_theta = fenics.TestFunctions(W)

w = fenics.Function(W)

p, u, theta = fenics.split(w)

""" Benchmark parameters """
mu, Pr, Ra = fenics.Constant(1.), fenics.Constant(0.71), fenics.Constant(1.e6)

ghat = fenics.Constant((0., -1.))

theta_h, theta_c  = fenics.Constant(0.5), fenics.Constant(-0.5)

""" Initial values """
w_n = fenics.interpolate(
        ("0.", "0.", "0.", "theta_h + x[0]*(theta_c - theta_h)"), 
        theta_h = theta_h, theta_c = theta_c, element = mixed_element),

p_n, u_n, theta_n = fenics.split(w_n)

""" Finite difference time discretization: Backward Euler"""
Delta_t = fenics.Constant(0.001)

u_t, theta_t = (u - u_n)/Delta_t, (theta - theta_n)/Delta_t

""" Variational form of the governing equations """
inner, dot, grad, div, sym = \
    fenics.inner,, fenics.grad, fenics.div, fenics.sym
mass = -psi_p*div(u)

f_B = Ra/Pr*theta*ghat

momentum = dot(psi_u, u_t + dot(grad(u), u) + f_B) - div(psi_u)*p \
    + 2.*mu*inner(sym(grad(psi_u)), sym(grad(u)))

energy = psi_theta*theta_t + dot(grad(psi_theta), 1./Pr*grad(theta) - theta*u)
F = (mass + momentum + energy)*fenics.dx

""" Penalty method stabilization """
gamma = fenics.Constant(1.e-7)

F += -psi_p*gamma*p*fenics.dx

""" Boundary conditions """
class HotWall(fenics.SubDomain):

    def inside(self, x, on_boundary):

        return on_boundary and fenics.near(x[0], 0.)
class ColdWall(fenics.SubDomain):

    def inside(self, x, on_boundary):

        return on_boundary and fenics.near(x[0], 1.)

class Walls(fenics.SubDomain):

    def inside(self, x, on_boundary):

        return on_boundary

W_u, W_theta = W.sub(1), W.sub(2)

boundary_conditions = [
    fenics.DirichletBC(W_u, (0., 0.), Walls()),
    fenics.DirichletBC(W_theta, theta_h, HotWall()),
    fenics.DirichletBC(W_theta, theta_c, ColdWall())]
""" Solve until approximate steady state """

fenics.solve(F == 0., w, boundary_conditions)

for timestep in range(4):
    fenics.solve(F == 0., w, boundary_conditions)

For reference, the temperature solution look like this:


Now I wish to compute something like the heat flux on the cold wall, the difficult parts of which are captured in the following:

           $f=\int_{\Gamma_c}^{ }\partial_n\theta ds=\int_{\Gamma_c}\left(\mathbf{n}\cdot\nabla\right)\theta ds=\int_{\Gamma_c}\nabla\theta\cdot\mathbf{n}ds$ƒ =Γcnθds=Γc(n·)θds=Γcθ·nds   

How does one compute this integral?

With minor modification from the FEniCS book, I tried

boundary_parts = fenics.MeshFunction("size_t", mesh, mesh.topology().dim() - 1)

cold_wall_marker = 2

ColdWall().mark(boundary_parts, cold_wall_marker)

nhat = fenics.FacetNormal(mesh)

integrate, ds, dot, grad = fenics.assemble, fenics.ds,, fenics.grad

cold_wall_heat_flux = integrate(dot(grad(theta), nhat)*ds(cold_wall_marker))


which runs but prints "0.0", while I expect a non-zero value.


Community: FEniCS Project
I might be completely wrong, but if theta is constant along the wall and you are calculating an integral along the wall, maybe FEniCS is calculating the gradient using the shape functions along that edge, and the values along that edge. Because these are constant, it will be zero. This is just my assumption, hopefully a FEniCS author can confirm or dismiss my point.
written 7 weeks ago by Miguel  

2 Answers

7 weeks ago by
Probably a bug in the user code along the lines of following

from dolfin import *
mesh = UnitSquareMesh(3, 3)
ff = MeshFunction('size_t', mesh, mesh.topology().dim()-1)
print(assemble(1*ds(subdomain_id=0, domain=mesh)))
print(assemble(1*ds(subdomain_id=0, domain=mesh, subdomain_data=ff)))​

which prints

Thanks for the demonstration. Between your answer and looking at help(fenics.ds) again, I came up with the following:

mesh_function = fenics.MeshFunction("size_t", mesh, mesh.topology().dim() - 1)

cold_wall_id = 2

hot_wall_id = 3

ColdWall().mark(mesh_function, cold_wall_id)

HotWall().mark(mesh_function, hot_wall_id)

nhat = fenics.FacetNormal(mesh)

integrate, dot, grad = fenics.assemble,, fenics.grad

ds = fenics.ds(
    domain = mesh, subdomain_data = mesh_function, subdomain_id = cold_wall_id)

cold_wall_heat_flux = integrate(dot(grad(theta), nhat)*ds)


ds = fenics.ds(
    domain = mesh, subdomain_data = mesh_function, subdomain_id = hot_wall_id)
hot_wall_heat_flux = integrate(dot(grad(theta), nhat)*ds)


which prints


which has the qualities I expected, e.g. approximately opposite flux for the hot and cold walls.
written 7 weeks ago by Alexander G. Zimmerman  
7 weeks ago by
I had a similar problem (with electric potential/current).
Besides the definition of the measure ds which Jan Blechta poitned out, there is an other problem.
The problem is (in my understanding) if you use grad on the solution you differentiate the the basis functions (going from linear to constant and therefore discontinous). This gives you some interpolation errors. You will this if you solve a poisson equation with two different Dirichlet BCs and a zero source term, and than integrate the in and ouflow. They will differ (for my problem about 20%).

I solved this by imposing the Dirichlet BCs through Lagrange Multipliers, the Lagrange Multipliers than have the meaning of the flow through that boundary. This worked for me and had far less performance impact than using higher degree elements (which also gives you better results for the flow interpolation).

Thanks for the warning. In the foreseeable future I will only be using this for regression testing, so thankfully I am okay with the interpolation errors in this case.
written 7 weeks ago by Alexander G. Zimmerman  
Please login to add an answer/comment or follow this question.

Similar posts:
Search »