Compute gradient of scalar field on BoundaryMesh

[BenjaminRueth]

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

[kamensky]

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
f = -div(grad(u_exact))
def R(u,v):
    return (inner(grad(u),grad(v)) - f*v)*dx

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)
consistencyTerm = inner(grad(u_h),n)*v*ds_marked(COMPLEMENT_FLUX_BDRY)
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):
    err = inner(grad(u_exact),n) - 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:
#
#  https://link.springer.com/article/10.1007/BF01385871
#
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:
print("Error in direct flux: "+str(fluxErr(inner(grad(u_h),n))))

# Slow for large meshes; also first-order:
#print("Error in projection to linears: "
#      +str(fluxErr(inner(project(grad(u_h),
#                                 VectorFunctionSpace(mesh,"CG",1)),n))))

[BenjaminRueth]

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.

@BenjaminRueth

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.

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.

[dg0609]

Is there an updated version of this for dolfinx?