# Compute gradient of scalar field on BoundaryMesh

Dear community,

I am currently trying to correctly compute the (normal) heat flux of a given temperature field on the boundary of a square domain. I am trying to extend the FEniCS Heat Equation Tutorial with the computation of the heat flux as a postprocessing. The temperature field is already known. The tutorial case I am discussion has the analytical solution u = 1 + x^2 + \alpha y^2 + \beta t. This gives the heat flux q_x = 2x and q_y = 2\alpha x. For my experiments I picked \alpha = 3.0 and \beta = 1.3.

How can I accurately compute the flux? I would like to obtain the exact solution for the flux on the whole boundary. I hope for such a high accuracy since with the numerical tools that are used in the tutorial, I can also compute the exact solution for the temperature and I would like to achieve the same for the flux.

In the following there is a quite lengthy discussion on what I have tried so far. Unfortunately, none of the approaches fits perfectly:

The thread Project gradient on BoundaryMesh explains how to compute the gradient for a Poisson equation. I also got this working for the Heat Equation, by manipulating print("Source term: "+str(assemble(f*dx(domain=mesh)))) to also include the time derivative print("Source term: "+str(assemble(f * dx(domain=mesh) - dudt * dx()))). The balance of source term and flux is fitting.

However, when I look closer, I see that the heat flux only fits in an integral sense. At the corners of my square domain I observe a Gibbâ€™s Phenomenon-like jump (see pictures section below). This is from my perspective understandable, since at the corners a normal flux cannot be determined. I would rather obtain a flux in each coordinate direction.

## Computing Derivatives

This tutorial explains how to compute the gradient in an alternative way. However, the balances from the example above disagree:

Flux: 8.000000000000012
Source term: -8.006720724414592


if I plot the solution (see pictures section below), I can see that the flux is not accurate enough. The flux across the right boundary should be equal to 2 according to the analytical solution. Here we get ~2.08.

## Do it manually

I used an approach described in Toselli, Andrea, and Olof Widlund. Domain decomposition methods-algorithms and theory . Vol. 34. Springer Science & Business Media, 2006., p.3, where Greenâ€™s identity is used to compute the normal flux. I think the code explains it better than I can:

F = u * v / dt * dx + dot(grad(u), grad(v)) * dx - (u_n / dt + f) * v * dx


is the weak form of the problem. I compute u_np1 via

solve(lhs(F) == rhs(F), u_np1, bcs)


Next I define

F_known_u = u_np1 * v / dt * dx + dot(grad(u_np1), grad(v)) * dx - (u_n / dt + f) * v * dx


which is identical to the weak form, but u = TrialFunction(V) being replaced with the known solution u_np1. Finally, I run

fluxes_vector = assemble(F_known_u)  # assemble weak form -> evaluate integral
v = TestFunction(V)
fluxes = Function(V)  # create function for flux
area = assemble(v * ds).get_local()
for i in range(area.shape[0]):
if area[i] != 0:  # put weight from assemble on function
fluxes.vector()[i] = fluxes_vector[i] / area[i]  # scale by surface area
else:
assert(abs(fluxes_vector[i]) < 10**-10)  # for non surface parts, we expect zero flux
fluxes.vector()[i] = fluxes_vector[i]


The resulting flux is so far the most accurate one (see pictures section below). At the corners there is still a problem with the singularity, since my flux is not a VectorFunction. The main problem with this approach is, that from my perspective the approach is very low level and not straightforward. This leads to trouble, if I want to run the code on multiple processes (mpirun -n np python3 heat.py; this might also originate from my limited understanding of the parallel computing capabilities of FEniCS). This was the main motivation for starting to look for an alternative approach.

## Pictures

Short explanation: I am plotting the normal heat flux over the arc length along the boundary. I traverse the square domain [0,1] x [0,1] starting at (0,0) in clockwise fashion.

Left: Project gradient on BoundaryMesh; Center: Computing Derivatives; Right: Do it manually

## edit: Mixed Finite Elements

I can also solve the heat equation using mixed finite elements. However, this take me quite far away from the basic tutorial for the heat equation that I mentioned above and that it not really fit together with the â€śpostprocessingâ€ť idea that I have in mind.

from fenics import Function, SubDomain, RectangleMesh, FunctionSpace, Point, Expression, Constant, DirichletBC, \
TrialFunction, TestFunction, File, solve, plot, lhs, rhs, grad, inner, dot, dx, ds, assemble, interpolate, \
project, near, VectorFunctionSpace, BoundaryMesh, Measure, FacetNormal, FiniteElement, VectorElement, TrialFunctions, TestFunctions

from matplotlib import pyplot as plt

x_left, x_right = 0, 1
y_bottom, y_top = 0, 1
p0 = Point(x_left, y_bottom)
p1 = Point(x_right, y_top)
nx, ny = 20, 20
alpha = 3.0  # parameter alpha
beta = 1.0  # parameter beta
dt = 1

u_D = Expression('1 + x[0]*x[0] + alpha*x[1]*x[1] + beta*t', degree=2, alpha=alpha, beta=beta, t=0)
f = Constant(beta - 2 - 2 * alpha)
mesh = RectangleMesh(p0, p1, nx, ny, diagonal='crossed')
P = FiniteElement("Lagrange", mesh.ufl_cell(), 2)  # use quadratic temperature triangle
Q = VectorElement("Lagrange", mesh.ufl_cell(), 1)  # use linear flux
V = FunctionSpace(mesh, P * Q)

u_n = interpolate(u_D, V.sub(0).collapse())
u_D.t = dt

(v, w) = TestFunctions(V)
(u, q) = TrialFunctions(V)
mixed_np1 = Function(V)
u_np1, q_np1 = mixed_np1.split()

def pde(u, q, v, w):
dudt = (u - u_n) / dt
F = 0
F += dot(q, grad(v)) * dx
F += dudt * v * dx
F += -f * v * dx
F += dot(grad(u), w) * dx
F += -dot(q, w) * dx
return F

bcStr = "near(x[0]," + str(x_left) + ")" \
" || " \
"near(x[0]," + str(x_right) + ")" \
" || " \
"near(x[1]," + str(y_bottom) + ")" \
" || " \
"near(x[1]," + str(y_top) + ")"

bcStr = "("+bcStr+")"
bc = DirichletBC(V.sub(0), u_D, bcStr)
F = pde(u, q, v, w)
solve(lhs(F) == rhs(F), mixed_np1, bc)

dudt = (u_np1 - u_n) / dt

normal = FacetNormal(mesh)

print("Flux: "+str(assemble(dot(q_np1, normal) * ds)))
print("Source term: "+str(assemble((f - dudt) * dx(domain=mesh))))

import numpy as np
from matplotlib import pyplot as plt
plt.figure(1)
plt.plot(np.linspace(0, 1), [q_np1(x_left, y)[0] for y in np.linspace(0, 1)])
plt.plot(np.linspace(1, 2), [q_np1(x, y_top)[1] for x in np.linspace(0, 1)])
plt.plot(np.linspace(2, 3), [q_np1(x_right, y)[0] for y in np.linspace(1, 0)])
plt.plot(np.linspace(3, 4), [q_np1(x, y_bottom)[1] for x in np.linspace(1, 0)])
plt.show()


1 Like

For domains with corners, I would recommend extracting the flux from each smooth section of the boundary separately. When extracting flux from only a portion of the boundary, there is an extra â€śtrickâ€ť that needs to be added on top of what I wrote up in the earlier post that you linked to.

Take a look at the following example, which I prepared some time ago for a research collaborator of mine; it allows for 3/2 rate of convergence in L^2 on the boundary, which is, I think, about as good as you can hope to get with linear elements:

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

# The full boundary, on which we apply a Dirichlet BC:
class Bdry(SubDomain):
def inside(self,x,on_boundary):
return on_boundary

# Interior of the domain:
class Interior(SubDomain):
def inside(self,x,on_boundary):
return (not (near(x[0],0) or near(x[0],1)
or near(x[1],0) or near(x[1],1)))

# The boundary on the right side of the domain, on which we want to
# extract the flux:
class BdryOfInterest(SubDomain):
def inside(self,x,on_boundary):
return (x[0] > 1.0 - DOLFIN_EPS) and on_boundary

# Everything BUT the right side of the domain from which we want to
# extract the boundary flux; this must include all nodes not used to
# approximate flux on the boundary of interest, including those in the
# interior of the domain (so on_boundary should not be used in the
# return value).
class AntiBdry(SubDomain):
def inside(self,x,on_boundary):
return x[0] < 1.0 - DOLFIN_EPS

# Mark parts of the boundary to which integrals will need to be restricted.
# (Parts not explicitly marked are flagged with zero.)
FLUX_BDRY = 1
COMPLEMENT_FLUX_BDRY = 2
boundaryMarkers = MeshFunction("size_t",mesh,mesh.topology().dim()-1,
COMPLEMENT_FLUX_BDRY)
BdryOfInterest().mark(boundaryMarkers,FLUX_BDRY)
ds_marked = ds(subdomain_data=boundaryMarkers)

# Desired exact solution:
u_exact = Expression("x[0]*x[0] + x[1]*x[1]", degree=2, domain=mesh)

# Boundary conditions to apply to solution, u_h:
BCs = [DirichletBC(V,u_exact,Bdry()),]

# Boundary conditions to apply when solving for flux on whole boundary:
interiorBCs = [DirichletBC(V,Constant(0.0),Interior()),]

# Boundary conditions to apply to the flux solution when we are only
# interested in flux on the right side of the domain:
antiBCs = [DirichletBC(V,Constant(0.0),AntiBdry()),]

# Residual of the Poisson equation
def R(u,v):

u = TrialFunction(V)
v = TestFunction(V)
F = R(u,v)
a = lhs(F)
L = rhs(F)

# Solve and put solution in u_h
u_h = Function(V)
solve(a==L,u_h,BCs)

# Solve in the trace space for what Neumann data, i.e.,
#
#  \mathbf{q}\cdot\mathbf{n} = \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.
qn = TrialFunction(V)

################################################################################

# The trick:  Since we want to use the corner nodes to approximate the
# flux on our boundary of interest, test functions will end up being
# nonzero on an $O(h)$ part of the complement of the boundary of interest.
# Thus we need to integrate a consistency term on that part of the boundary.

n = FacetNormal(mesh)
FBdry = qn*v*ds_marked(FLUX_BDRY) - R(u_h,v) + consistencyTerm

FBdry_inconsistent = qn*v*ds_marked(FLUX_BDRY) - R(u_h,v)

# Applying the consistent flux extraction on the full boundary
# and then restricting the result is also sub-optimal; this flux extraction
# technique doesn't appear to play nicely with corners.
FBdry_full = qn*v*ds - R(u_h,v)

################################################################################

# Get $\mathbf{q}\cdot\mathbf{n}$ on the boundary of interest with and
# without the consistency term:
def solveFor_qn_h(FBdry, BCs):
aBdry = lhs(FBdry)
LBdry = rhs(FBdry)
ABdry = assemble(aBdry,keep_diagonal=True)
bBdry = assemble(LBdry)
[BC.apply(ABdry,bBdry) for BC in BCs]
qn_h = Function(V)
solve(ABdry,qn_h.vector(),bBdry)
return qn_h

qn_h = solveFor_qn_h(FBdry, antiBCs)
qn_h_inconsistent = solveFor_qn_h(FBdry_inconsistent, antiBCs)
qn_h_full = solveFor_qn_h(FBdry_full, interiorBCs)

# Compare fluxes with the exact solution:
import math
def fluxErr(qn):
return math.sqrt(assemble(err*err*ds_marked(FLUX_BDRY)))

# Converges at 3/2 order in $L^2$, as expected for smooth solutions, according
# to this paper:
#
#
print("Error in consistent flux on boundary of interest: "
+str(fluxErr(qn_h)))

# Converges at only 1/2 order:
print("Error in restricted consistent flux from full boundary: "
+str(fluxErr(qn_h_full)))

# Converges at only 1/2 order:
print("Error in inconsistent flux on boundary of interst: "
+str(fluxErr(qn_h_inconsistent)))

# Converges at first order:

# Slow for large meshes; also first-order:
#print("Error in projection to linears: "

Thanks for the detailed answer and the code. This approach increases the accuracy a lot. I created the same kind of picture that I already used above. The dotted line shows the flux, if the full domain is considered (qn_h_full), while the solid line shows the flux for a single edge (qn_h).
However, the solution for my problem that I used in the end is actually a different one: You mentioned that a 3/2 rate of convergence is as good as I can hope for linear elements. Therefore, I switched to quadratic elements for temperature (V = FunctionSpace(mesh, 'P', 2)) and used this approach.
With quadratic elements I am able to get the exact solution and since I am using a VectorFunctionSpace the corners do not cause any problems.