Projection of mixed elements over functionspace : jit error

I/O
1

Hi, I have been trying to project a vector of 20 numbers over a functionspace. it works till 2nd time steps but at 3rd time steps the following error is happening.
*Error: Unable to perform just-in-time compilation of form.
***** Reason: ffc.jit failed with message:

I have been trying to run following code:

#**********************************************************************************************************
#**********************************************************************************************************

import sys, pdb

sys.setrecursionlimit(100000)
from dolfin import *
import os as os
#import numpy as np
from math import *
import vtk

from mpi4py import MPI as pyMPI

from matplotlib import pylab as plt


#************************************************************************************

deg = 4
flags = ["-O3", "-ffast-math", "-march=native"]
parameters["form_compiler"]["representation"]="uflacs"
parameters["form_compiler"]["quadrature_degree"]=deg

print(parameters.linear_algebra_backend)



le = 1.0  #length
wd = 1.0  #width
ht = 0.1  #height


#Mesh create
mesh = BoxMesh( Point(0.0, 0.0, 0.0),Point( le, wd, ht),4,4,2)
print "Plotting a BoxMesh"
File ("test/slabmesh"+".pvd") << mesh

print mesh.num_cells()
print mesh.num_vertices()

comm = mesh.mpi_comm()



class Top(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[2], ht) and on_boundary

class Bottom(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[2], 0) and on_boundary


class Left(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[0], 0) and on_boundary

class Right(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[0], le) and on_boundary

class Front(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[1], wd) and on_boundary

class Back(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[1], 0) and on_boundary



#Faceboundaries 
facetboundaries = MeshFunction("size_t", mesh, mesh.topology().dim() - 1)
facetboundaries.set_all(0)
left = Left()
left.mark(facetboundaries, 1) # Mark left id as 1
right = Right()
right.mark(facetboundaries, 2) # Mark right id as 2
back = Back()
back.mark(facetboundaries, 4) # Mark back id as 4





#..........................................................................................
dx = dolfin.dx(mesh)
ds = dolfin.ds(mesh) 


#VectorFunctionSpace for fiber vector
VFS = VectorFunctionSpace(mesh, 'CG', 2)
FS = FunctionSpace(mesh, 'CG', 1)
TFS = TensorFunctionSpace(mesh, "CG", 1) #, shape=(3,3))

Telem = TensorElement('CG', mesh.ufl_cell(), 2, quad_scheme="default")
Telem._quad_scheme = 'default'
Velem = VectorElement('CG', mesh.ufl_cell(), 1, quad_scheme="default")
Velem._quad_scheme = 'default'
Qelem = FiniteElement('DG', mesh.ufl_cell(), 1, quad_scheme="default")
Qelem._quad_scheme = 'default'
Relem = FiniteElement("Real", mesh.ufl_cell(), 0, quad_scheme="default")
Relem._quad_scheme = 'default'
VRelem = MixedElement([Relem, Relem, Relem, Relem, Relem, Relem])

W = FunctionSpace(mesh, MixedElement([Velem,VRelem])) 


w = Function(W)
dw = TrialFunction(W)
wtest = TestFunction(W)

du,dc = TrialFunctions(W)
(u, c) = split(w)
(v, cq) = TestFunctions(W)


tTot = 100
T1st = tTot /2
value = 0.0

## Boundary condition declaration 
ux = Expression(("m0"), m0 = value, degree = 0)
cl = Expression(("0.0", "0.0", "0.0"), degree = 0)
#bcleft = DirichletBC(W, cl, facetboundaries, 1)
bcback = DirichletBC(W.sub(0).sub(1), Expression(("0.0"), degree = 0), facetboundaries, 4)
bcleft = DirichletBC(W.sub(0).sub(0),  Expression(("0.0"), degree = 0), facetboundaries, 1)
bcright = DirichletBC(W.sub(0).sub(0), ux, facetboundaries, 2)

bcs = [bcleft, bcback, bcright] 

#########################################################################
#COntinuum formulation

N = FacetNormal(mesh )

d = u.ufl_domain().geometric_dimension()
I = Identity(d)
F = I + grad(u)
#print F

C = F.T*F
J = det(F)
n0 = J*inv(F.T)*N 
E = (F.T*F - I)*0.5

Ufile = File("test/Utensor.pvd")
Enh = 1.0E6 #Pa
nu = 0.19

sigmaMax = 2.0E5


k = Enh/(3*(1-2*nu))
G = Enh/(2*(1+nu))


sigmaNH = k*ln(J)*I/J + G*(F*F.T - I*J**2/3)/J


	
#################################################
#FIber formulation
n = 20
k0 = 1.5e-06
k1 = 7.0e-07
kd = 1.0e-03

Qelem_array = MixedElement([Qelem]*n)
Velem_array = MixedElement([Velem]*n)
#Telem_array = MixedElement([Telem]*n)
WQ = FunctionSpace(mesh, Qelem_array)
wq = FunctionSpace(mesh, Qelem)
WV = FunctionSpace(mesh, Velem_array)
#WT = FunctionSpace(mesh, Telem_array)
Wphim = FunctionSpace(mesh, Qelem)


Q = Function(WQ)
Q.rename("fibervol","fibervol")
V = Function(WV)

Qm = Function(wq)

Qm.vector().array()[:] = 0.05
Q.vector().array()[:] = 0.0

thetain = -pi/2 
thetafin = pi/2 -pi/n
dtheta = pi/n 
K = 0.0
E0 = 0.12


THETA = []
for i in range(0,n):
	THETA.append(thetain)

	thetain = thetain + dtheta

knu = 5.0E01
sigmax = 2.0E05
phitotal = 0.05
zer = Constant((0.0))
phi = as_vector([(zer*k) for k in range (0,n)])

phivt = phi[0]
for k in range (0,n-1):
	phivt = phivt + phi[k+1]
phim =(phitotal -(phivt/n))




Vcell = as_vector([[cos(THETA[k]), sin(THETA[k]), 0.0] for k in range(0,n)])

Vnew = as_vector([(F*Vcell[k]) for k in range(0,n)])

lamda = as_vector([(dolfin.sqrt(inner(C*Vcell[k],Vcell[k]))) for k in range(0,n)])


strain = as_vector([((lamda[k]**2 - 1)*0.5) for k in range (0,n)])


strainOld = strain
straindot = (strain - strainOld)/1.0


fact = as_vector([(dolfin.exp((strain[k]/0.12)**2)) for k in range(0,n)])


fpas = as_vector([(conditional(ge(0.0, strain[k]), (strain[k]/0.1)**2, 0.0)) for k in range (0,n)])


fe = fact + fpas



fedot = as_vector([((knu*straindot[k] +2.0)/dolfin.sqrt((knu*straindot[k] + 2.0)**2 +1.0))/ (1+ 2/dolfin.sqrt(5.0)) for k in range (0,n)])


sigmatheta = as_vector([(sigmax*fe[k]*fedot[k]) for k in range(0,n)])


MCell =as_tensor([(outer(Vnew[k], Vnew[k])) for k in range(0,n)])


sigmacell = [(phi[k]*sigmatheta[k]*MCell[k]) for k in range(0,n)]



SigmaX = sigmacell[0]
for k in range (0,n-1):
	SigmaX = SigmaX + sigmacell[k+1]  

SigmaCell = as_tensor(SigmaX)


cauchy = sigmaNH +  SigmaCell/n


piola = J*cauchy*inv(F.T)


F1 = - inner(grad(v),piola)*dx

X = SpatialCoordinate(mesh)
Wrigid = inner(as_vector([c[0], c[1], c[2]]), u) + \
	 inner(as_vector([0.0, 0.0, c[3]]), cross(X, u)) + \
	 inner(as_vector([c[4], 0.0, 0.0]), cross(X, u)) + \
	 inner(as_vector([0.0, c[5], 0.0]), cross(X, u))

F2 = derivative(Wrigid, w, wtest)*dx

Ftotal = F1 + F2
Jac1 = derivative(F1, w, dw)
Jac2 = derivative(F2, w, dw)


Jac = Jac1 +Jac2



problem = NonlinearVariationalProblem(Ftotal, w, bcs, Jac)
solver = NonlinearVariationalSolver(problem)

prm = solver.parameters
prm['newton_solver']['absolute_tolerance'] = 1E-7
prm['newton_solver']['relative_tolerance'] = 1E-6
prm['newton_solver']['maximum_iterations'] = 40
prm['newton_solver']['relaxation_parameter'] = 1.0
prm['newton_solver']['linear_solver']='umfpack'




fileu = File("test/disp.pvd")
fileu << w.sub(0)
track = 0.0
tsteps = 0.0
totalT = 100
dt = 1.0




while(track<=totalT):

	solver.solve() 


	
	phiold = phi
	phi = as_vector([(phiold[k] + (k0 + k1*sigmax*fact[k]*fedot[k])*phim - kd*phiold[k]) for k in range(0,n)])


	phivt = phi[0]
	for k in range (0,n-1):
		phivt = phivt + phi[k+1]
	phim = (phitotal -(phivt/n))

	PHIv =  project(phi, WQ) ##The error happens here

	fileu << w.sub(0)

	tsteps = tsteps + dt
	track = track+dt

	if tsteps < totalT:
		ux.m0 = ux.m0+0.01

#***************************************************************************************************************
#***************************************************************************************************************

Can anybody shade some light on this issue? I need to extract the phi value over a node or a specific point in a mesh.

2

Could you reduce this problem to a minimal working example which reproduces your error? In its current form it’s not immediately compatible with the latest version of DOLFIN.

1 Like
3

Hi Nate, thank you for prompt reply.

I tried to shorten the code. The error happens at following line under while loop:

PHIv =  project(phi, WQ) #Jit Error happens here

I tried to shorten the code in following way. will it work?

import sys
sys.setrecursionlimit(100000)
from dolfin import *


#************************************************************************************

deg = 4
flags = ["-O3", "-ffast-math", "-march=native"]
parameters["form_compiler"]["representation"]="uflacs"
parameters["form_compiler"]["quadrature_degree"]=deg

le = 1.0  #length
wd = 1.0  #width
ht = 0.1  #height

#Mesh create
mesh = BoxMesh( Point(0.0, 0.0, 0.0),Point( le, wd, ht),4,4,2)
print "Plotting a BoxMesh"
File ("test/slabmesh"+".pvd") << mesh

print mesh.num_cells()
print mesh.num_vertices()

comm = mesh.mpi_comm()

#..........................................................................................
dx = dolfin.dx(mesh)
n = 20


Velem = VectorElement('CG', mesh.ufl_cell(), 1, quad_scheme="default")
Velem._quad_scheme = 'default'
Qelem = FiniteElement('DG', mesh.ufl_cell(), 1, quad_scheme="default")
Qelem._quad_scheme = 'default'
Relem = FiniteElement("Real", mesh.ufl_cell(), 0, quad_scheme="default")
Relem._quad_scheme = 'default'
VRelem = MixedElement([Relem, Relem, Relem, Relem, Relem, Relem])


Qelem_array = MixedElement([Qelem]*n)
WQ = FunctionSpace(mesh, Qelem_array)
W = FunctionSpace(mesh, MixedElement([Velem,VRelem])) 


w = Function(W)
dw = TrialFunction(W)
wtest = TestFunction(W)

du,dc = TrialFunctions(W)
(u, c) = split(w)
(v, cq) = TestFunctions(W)


tTot = 100
T1st = tTot /2
value = 0.0

## Boundary condition declaration 

class Left(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[0], 0) and on_boundary

class Right(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[0], le) and on_boundary


class Back(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[1], 0) and on_boundary



#Faceboundaries 
facetboundaries = MeshFunction("size_t", mesh, mesh.topology().dim() - 1)
facetboundaries.set_all(0)
left = Left()
left.mark(facetboundaries, 1) # Mark left id as 1
right = Right()
right.mark(facetboundaries, 2) # Mark right id as 2
back = Back()
back.mark(facetboundaries, 4) # Mark back id as 4
ux = Expression(("m0"), m0 = value, degree = 0)
bcback = DirichletBC(W.sub(0).sub(1), Expression(("0.0"), degree = 0), facetboundaries, 4)
bcleft = DirichletBC(W.sub(0).sub(0),  Expression(("0.0"), degree = 0), facetboundaries, 1)
bcright = DirichletBC(W.sub(0).sub(0), ux, facetboundaries, 2)

bcs = [bcleft, bcback, bcright] 

#########################################################################
#Continuum formulation

N = FacetNormal(mesh )

d = u.ufl_domain().geometric_dimension()
I = Identity(d)
F = I + grad(u)
C = F.T*F
J = det(F)
n0 = J*inv(F.T)*N 
E = (F.T*F - I)*0.5

Enh = 1.0E6 #Pa
nu = 0.19

k = Enh/(3*(1-2*nu))
G = Enh/(2*(1+nu))

sigmaNH = k*ln(J)*I/J + G*(F*F.T - I*J**2/3)/J

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

piola = J*sigmaNH*inv(F.T)
F1 = - inner(grad(v),piola)*dx

X = SpatialCoordinate(mesh)
Wrigid = inner(as_vector([c[0], c[1], c[2]]), u) + \
	 inner(as_vector([0.0, 0.0, c[3]]), cross(X, u)) + \
	 inner(as_vector([c[4], 0.0, 0.0]), cross(X, u)) + \
	 inner(as_vector([0.0, c[5], 0.0]), cross(X, u))

F2 = derivative(Wrigid, w, wtest)*dx

Ftotal = F1 + F2
Jac1 = derivative(F1, w, dw)
Jac2 = derivative(F2, w, dw)

Jac = Jac1 +Jac2

problem = NonlinearVariationalProblem(Ftotal, w, bcs, Jac)
solver = NonlinearVariationalSolver(problem)

thetain = -pi/2 
thetafin = pi/2 -pi/n
dtheta = pi/n 
K = 0.0
E0 = 0.12

THETA = []
for i in range(0,n):
	THETA.append(thetain)
	thetain = thetain + dtheta
phitotal = Constant((0.05))
zer = Constant((0.0))
phi = as_vector([(zer*k) for k in range (0,n)])

phivt = phi[0]
for k in range (0,n-1):
	phivt = phivt + phi[k+1]
phim =(phitotal -(phivt/n))

Vcell = as_vector([[cos(THETA[k]), sin(THETA[k]), 0.0] for k in range(0,n)])
lamda = as_vector([(dolfin.sqrt(inner(C*Vcell[k],Vcell[k]))) for k in range(0,n)])
strain = as_vector([((lamda[k]**2 - 1)*0.5) for k in range (0,n)])

track = 0.0
tsteps = 0.0
totalT = 100
dt = 1.0

while(tsteps<=totalT):

	solver.solve() 

	
	phiold = phi

	## Updating Phi
	phi = as_vector([(phiold[k] + (7.0E-7*strain[k])*phim - 1.0E-3*phiold[k]) for k in range(0,n)])
	PHIv =  project(phi, WQ)     #Jit Error happens here


	phivt = phi[0]
	for k in range (0,n-1):
		phivt = phivt + phi[k+1]
	phim = (phitotal -(phivt/n))

	tsteps = tsteps + dt
	ux.m0 = ux.m0+0.002
4

You are using Python2, indicating that you are using a version of FEniCS prior to the 2018.1.0 release. This is quite outdated, and I would strongly encourage you to update to a newer version with python3.

5

Another thing that seems strange to me is that you do not update w in your while loop. However, I must admit it is very hard for me to tell what you want to do, as your example is far from minimal.
To me, it seems like your projection loop is fully decoupled from the problem you are solving.

1 Like
6

Hi, thanks for replying.

The error was updating phi under while loop. Since it was declared under while loop, the form was getting bigger at each time step and taking longer time to project. I solved that issue.

import sys, pdb
sys.setrecursionlimit(100000)
from dolfin import *


#************************************************************************************

deg = 4
flags = ["-O3", "-ffast-math", "-march=native"]
parameters["form_compiler"]["representation"]="uflacs"
parameters["form_compiler"]["quadrature_degree"]=deg

le = 1.0  #length
wd = 1.0  #width
ht = 0.1  #height

#Mesh create
mesh = BoxMesh( Point(0.0, 0.0, 0.0),Point( le, wd, ht),4,4,2)


print mesh.num_cells()
print mesh.num_vertices()

comm = mesh.mpi_comm()

#..........................................................................................
dx = dolfin.dx(mesh)
n = 20


Velem = VectorElement('CG', mesh.ufl_cell(), 1, quad_scheme="default")
Velem._quad_scheme = 'default'
Qelem = FiniteElement('DG', mesh.ufl_cell(), 1, quad_scheme="default")
Qelem._quad_scheme = 'default'
Relem = FiniteElement("Real", mesh.ufl_cell(), 0, quad_scheme="default")
Relem._quad_scheme = 'default'
VRelem = MixedElement([Relem, Relem, Relem, Relem, Relem, Relem])


Qelem_array = MixedElement([Qelem]*n)
WQ = FunctionSpace(mesh, Qelem_array)
W = FunctionSpace(mesh, MixedElement([Velem,VRelem])) 


w = Function(W)
dw = TrialFunction(W)
wtest = TestFunction(W)

du,dc = TrialFunctions(W)
(u, c) = split(w)
(v, cq) = TestFunctions(W)


tTot = 100
T1st = tTot /2
value = 0.0

## Boundary condition declaration 

class Left(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[0], 0) and on_boundary

class Right(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[0], le) and on_boundary


class Back(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[1], 0) and on_boundary



#Faceboundaries 
facetboundaries = MeshFunction("size_t", mesh, mesh.topology().dim() - 1)
facetboundaries.set_all(0)
left = Left()
left.mark(facetboundaries, 1) # Mark left id as 1
right = Right()
right.mark(facetboundaries, 2) # Mark right id as 2
back = Back()
back.mark(facetboundaries, 4) # Mark back id as 4
ux = Expression(("m0"), m0 = value, degree = 0)
bcback = DirichletBC(W.sub(0).sub(1), Expression(("0.0"), degree = 0), facetboundaries, 4)
bcleft = DirichletBC(W.sub(0).sub(0),  Expression(("0.0"), degree = 0), facetboundaries, 1)
bcright = DirichletBC(W.sub(0).sub(0), ux, facetboundaries, 2)

bcs = [bcleft, bcback, bcright] 

#########################################################################
#Continuum formulation

N = FacetNormal(mesh )

d = u.ufl_domain().geometric_dimension()
I = Identity(d)
F = I + grad(u)
C = F.T*F
J = det(F)
n0 = J*inv(F.T)*N 
E = (F.T*F - I)*0.5

Enh = 1.0E6 #Pa
nu = 0.19

k = Enh/(3*(1-2*nu))
G = Enh/(2*(1+nu))

sigmaNH = k*ln(J)*I/J + G*(F*F.T - I*J**2/3)/J

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

piola = J*sigmaNH*inv(F.T)
F1 = - inner(grad(v),piola)*dx

X = SpatialCoordinate(mesh)
Wrigid = inner(as_vector([c[0], c[1], c[2]]), u) + \
	 inner(as_vector([0.0, 0.0, c[3]]), cross(X, u)) + \
	 inner(as_vector([c[4], 0.0, 0.0]), cross(X, u)) + \
	 inner(as_vector([0.0, c[5], 0.0]), cross(X, u))

F2 = derivative(Wrigid, w, wtest)*dx

Ftotal = F1 + F2
Jac1 = derivative(F1, w, dw)
Jac2 = derivative(F2, w, dw)

Jac = Jac1 +Jac2

problem = NonlinearVariationalProblem(Ftotal, w, bcs, Jac)
solver = NonlinearVariationalSolver(problem)

thetain = -pi/2 
thetafin = pi/2 -pi/n
dtheta = pi/n 
K = 0.0
E0 = 0.12

THETA = []
for i in range(0,n):
	THETA.append(thetain)
	thetain = thetain + dtheta
phiold = Function(WQ)
phinew = Function(WQ)
phiold.vector()[:] = 0.0
phinew.vector()[:] = 0.0

Vcell = as_vector([[cos(THETA[k]), sin(THETA[k]), 0.0] for k in range(0,n)])
lamda = as_vector([(dolfin.sqrt(inner(C*Vcell[k],Vcell[k]))) for k in range(0,n)])
strain = as_vector([((lamda[k]**2 - 1)*0.5) for k in range (0,n)])
phi = as_vector([(phiold[k] + (7.0E-7*strain[k]) - 1.0E-3*phiold[k]) for k in range(0,n)])

track = 0.0
tsteps = 0.0
totalT = 100
dt = 1.0

while(tsteps<totalT):

	solver.solve() 

	
	phinew.vector()[:]= project(phi,WQ).vector().array()[:]
	phiold.vector()[:] = phinew.vector().array()[:]

	print tsteps
	print phinew.vector().array()
	tsteps = tsteps + dt
	ux.m0 = ux.m0+0.002