Computational Time of Integrals

1

I’m trying to perform an area integration of my fenics-solution multiplied with an additional given function f using “assemble”. Therefore I compared three different implementation possibilities:

  • List item f as Expression

  • List item f as UserExpression

  • List item f as UserExpression with external function call

Here a minimal working example:

from dolfin import *
from mshr import *
import time
##########################


N_fem=30
l=20.0
w=10.0


time0 = time.time()
rectangle=Rectangle(Point(-l/2,-w/2),Point(l/2,w/2))
mesh = generate_mesh(rectangle , N_fem)
time1 = time.time()

#solve PDE


V = FunctionSpace(mesh,"CG", 1)


u0 = Expression('x[0]',degree=1)
def LeftRight(x, on_boundary):
	tol = 1e-14
	return on_boundary and (near(x[0], -l/2) or near(x[0], l/2))


bc = DirichletBC(V, u0, LeftRight)
u = TrialFunction(V)
v = TestFunction(V)
f = Constant(0)

a = inner(grad(u),grad(v))*dx
L = -f*v*dx

# solve variational problem
u = Function(V)
solve(a == L, u, bc)

time2 = time.time()


#########################
#Integration


class MyFunctionExpression(UserExpression):

	def eval(self, values, x):
		values[0] = x[0]+x[1]

	def value_shape(self):
		return(1,)

def fx(pos):
	return pos[0]+pos[1]

class MyFunctionExpression_ext(UserExpression):

	def eval(self, values, x):
		values[0] = fx(x)

	def value_shape(self):
		return(1,)



# three possible implementations of the
f_1=Expression('x[0] + x[1]',degree=1)
f_2=MyFunctionExpression()
f_3=MyFunctionExpression_ext()


time3 = time.time()
B=assemble(u*dx(mesh))
time4 = time.time()
B=assemble(u*f_1*dx(mesh))
time5 = time.time()
B=assemble(u*f_2[0]*dx(mesh))
time6 = time.time()
B=assemble(u*f_3[0]*dx(mesh))
time7 = time.time()

print('mesh:', time1-time0)
print('FEM:', time2-time1)
print('integration0:', time4-time3)
print('integration1:',time5-time4)
print('integration2:',time6-time5)
print('integration3:',time7-time6)

As a result I get on my computer:

mesh: 0.00998234748840332
FEM: 0.03625941276550293
integration0: 0.00445103645324707
integration1: 0.0037686824798583984
integration2: 0.06827163696289062
integration3: 0.07056236267089844

My questions:
I was wondering that the integration via expression (integration1) is a factor of ~20 (!!) times faster than integration2 and integration3. I also tried to use a cython function f in integration3 which made the code by a factor of ~3 times faster, but that’s nevertheless still far away from integration1. Is it possible to accelerate the variants 2&3 somehow else? In my specific application I need an if-else-statement within the function and am therefore not able to write it as Expression. Is there a possibility to avoid the huge loss in performance in such a case?

Thanks,
Florian

2

Are you sure that the if/else statement you need cannot be implemented either as a UFL Conditional, or using the C++ ternary operator in an Expression? (In case you don’t know, the spatial coordinates can be accessed within UFL via x = SpatialCoordinate(mesh).)

3

Chapter 10 of the FEniCS book notes that with variants of UserExpression the calls to eval needed for integration (done in DOLFIN, i.e. C++) will be done in Python, i.e. involve callbacks which are the source of inefficiency. This is in contrast to Expression which gets compiled to C++ and hence there are no callbacks. Note that Expression constructor is not the only way to get C++ compiled expressions, as shown in @kamensky example here.

4

Hi guys,
thanks a lot for your hints - the reason for the time differences became much clearer now to me! The final solution in my use case was btw. a pure one lined C+±definition of the expression (with the help of some ?-operators), which doesn’t look very elegant within the source code but the performance is nevertheless great! In the future I’ll also keep in mind the multiline C++ expressions in such a case, which for me seem much more convenient than the python-UserExpressions.