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:

## Project gradient on BoundaryMesh

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