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!