Segmentation fault when assembling forms with the curl operator applied on derivatives on hexahedral elements

Hi all,

I have the following form:

S = \sum_{K \in \mathcal{T}_h} (curl(\mathcal{L}\textbf{u}),curl(\mathcal{L}\textbf{v}))_K

Being \mathcal{T}_h each one of the domain cells and \mathcal{L} \textbf{u}:= \sigma \textbf{u} + (\beta \cdot \nabla)\textbf{u} - \mu \Delta \textbf{u}

When I try to build this stabilization term, I get the error of “Segmentation Fault (core dumped)” after calling the JIT compiler (it is not immediate, it takes some time until the error appears).

In order to make the debugging easier, I have considered the terms of \sigma and \mu equal to 1 and \beta = (1,0,0)^T.

Also, I tried to operate with each term separately and I have observed that the terms where the derivatives appear (advective and viscous terms), the error of Seg. Fault appears too.

The following MWE is replicating the issue (as I mentioned before, I have considered \sigma and \mu equal to 1)

 """
 Mimum Working Example

 Testing the curl operation over derivatives
 """
from dolfin import *

# Mesh generation
mesh = UnitCubeMesh.create(2,2,2,CellType.Type.hexahedron)

# Test and trial functions
V = VectorFunctionSpace(mesh,"Lagrange",3)
u = TrialFunction(V)
v = TestFunction(V)

# Stabilization term generation
dx  =  dx(metadata = {"quadrature_degree":4}) #interior facets
beta = as_vector([Constant(1.0),Constant(0.0),Constant(0.0)])

# Working
print("Trying test 1")
Laplu = u
Laplv = v
S = assemble(inner(curl(Laplu),curl(Laplv))*dx)
print("Test 1 worked")
print("--------------------------------------------------")

# Not Working
print("Trying test 2")
Laplu = grad(u)*beta
Laplv = grad(v)*beta
S = assemble(inner(curl(Laplu),curl(Laplv))*dx)
print("Test 2 worked")
print("--------------------------------------------------")

# Not Working
print("Trying test 3")
Laplu = u + grad(u)*beta
Laplv = v + grad(v)*beta
S = assemble(inner(curl(Laplu),curl(Laplv))*dx)
print("Test 3 worked")
print("--------------------------------------------------")

# Not working
print("Trying test 4")
Laplu = - div(grad(u))
Laplv = - div(grad(v))
S = assemble(inner(curl(Laplu),curl(Laplv))*dx)
print("Test 4 worked")

# Not working
print("Trying test 5")
Laplu = u + grad(u)*beta - div(grad(u))
Laplv = v + grad(v)*beta - div(grad(v))
S = assemble(inner(curl(Laplu),curl(Laplv))*dx)
print("Test 5 worked")

Also I have seen the following post:
fenics-project / DOLFIN / issues / #1072 - Segmentation fault when assembling forms with high-order derivatives on hexahedral elements — Bitbucket
in which the author has a similar problem, but it has no replies.

I really appreciate your time and your help!

If you want to use anything other than simplices I recommend you use dolfin-x. Old dolfin is not in development anymore and it’s highly unlikely anyone will want to fix this.

An immediate short-term workaround is to use the TSFC form representation, as discussed in my answer here. (This can also result in more efficient code than UFLACS.) However, TSFC can fail in other places, like integration over interior facets with quad/hex elements. It is possible to assemble different terms with different form representations, though, e.g.,

print("Trying test 6")

# Need to switch back to UFLACS for interior facets; comment this line for
# an error:
parameters["form_compiler"]["representation"] = 'uflacs'

S = assemble(inner(avg(u),avg(v))*dS)
print("Test 6 worked")
3 Likes

In the case of having the follwing:

S = inner(curl(Laplu),curl(Laplv))*dx
Js = inner(avg(u),avg(v))*dS

How could I assemble these two functions together (S + Js)?

I tried the following but it didn’t work

parameters["form_compiler"]["representation"] = "tsfc"
print("Trying test 7")

Laplu = u + grad(u)*beta - div(grad(u))
Laplv = v + grad(v)*beta - div(grad(v))
S = inner(curl(Laplu),curl(Laplv))*dx

parameters["form_compiler"]["representation"] = "uflacs"
Js = (inner(avg(u),avg(v))*dS)

S_total = assemble(S + Js )

print("Test 7 worked")
print("--------------------------------------------------")