Speeding up the computation in an expression

Question

Hi,

I need to evaluate an integral within an expression class (see below) but it is running extremely slow when using tplquad in scipy.

Do you guys have any suggestions to speed it up?

Would defining the expression in C++ and calling an integrator in C++ speed up the process?

Thanks!

class Integration(Expression):
     def __init__(self, ulimit=0.001, mesh=mesh):
                 self.ulimit = ulimit
                 self.mesh = mesh

    def diff_volume(self, p,t,r):
        return r**2*sin(p)

  def eval(self, value, x):

       upper=project(self.ulimit,FunctionSpace(self.mesh, "CG", 1))
       upper_array = upper.vector().array()

    r1 = 0.
    r2 = 1.
    # limits for theta
    t1 = 0
    t2 = 2*pi
    # limits for phi
    p1 = 0
    p2 = pi

    results = []
    for p in upper_array:
            triple_Res = tplquad(self.diff_volume, r1, r2, lambda r:   t1, lambda r:   t2,
                                       lambda r,t: p1, lambda r,t: p)[0]

            results.append(triple_Res)
    results = np.array(results)


    S = FunctionSpace(mesh, 'CG', 1)
    results_exp = Function(S)
    results_exp.vector()[:] = results;
    value[0] = results_exp(x[0], x[1], x[2]);

Answer 1

There are many ways to speed this up. First, try to reuse FunctionSpace and Function-objects (avoid creation in the eval-method). In general, keep anything that can be computed independently of x out of the eval-function.

In addition, definitely avoid project in the eval-method.

You could convert to a C++-expression, but you have a lot to gain by modifying your python code.

Edit: If you require more specific help, please post a complete working example.

Comments

Thanks a lot, Øyvind !
Below is a complete working example:

from dolfin import *
import numpy as np
from numpy import *
from matplotlib import pylab as plt
from scipy import optimize
from scipy import integrate
from scipy.integrate import tplquad
import os

set_log_active(False)

class Integration(Expression):
    def __init__(self, E=0.001, mesh=mesh):
        self.E = E
    self.mesh = mesh

    def f2(self, E11, E22, E33, b, c):

    p = E11*cos(c)**2 + E22*cos(b)**2*sin(c)**2 + E33*sin(b)**2*sin(c)**2
    invexp = lambda x: 1/(2*pi)**0.5/0.11*exp(-(0.28-x)**2/(2*0.11**2))*1.35*10**6*(p-x)/(1+2*x)

    return integrate.quad(invexp, 0, p, epsabs=1e-2, epsrel=1e-2)[0]

    def Jacobian(self, a,b,c):
        d1 = -((cos(b)*sin(c))/(cos(c)*sin(a) - cos(a)*cos(b)*sin(c)) - (cos(a)*sin(b)*sin(b)*sin(c)*sin(c))/(cos(c)*sin(a) - cos(a)*cos(b)*sin(c))**2)/((sin(b)*sin(b)*sin(c)*sin(c))/(cos(c)*sin(a) - cos(a)*cos(b)*sin(c))**2 + 1)
        d2 = -((cos(c)*sin(b))/(cos(c)*sin(a) - cos(a)*cos(b)*sin(c)) + (sin(b)*sin(c)*(sin(a)*sin(c) + cos(a)*cos(b)*cos(c)))/(cos(c)*sin(a) - cos(a)*cos(b)*sin(c))**2)/((sin(b)*sin(b)*sin(c)*sin(c))/(cos(c)*sin(a) - cos(a)*cos(b)*sin(c))**2 + 1)
        d3 = (sin(a)*sin(b)*sin(c))/(1 - (cos(a)*cos(c) + cos(b)*sin(a)*sin(c))**2)**(1/2)
        d4 = (cos(a)*sin(c) - cos(b)*cos(c)*sin(a))/(1 - (cos(a)*cos(c) + cos(b)*sin(a)*sin(c))**2)**(1/2)
        J=d1*d4-d2*d3
        return J

     # Distribution Function 
    def R(self, a,b,c):

       sigmaphi=0.14
       mphi=0.69
       phi=arccos(sin(a)*cos(b)*sin(c)+cos(a)*cos(c))
       if phi >= 0 and phi <= pi/2:
        return 1/(2*(2*pi)**0.5*sigmaphi)*exp(-(mphi-phi)**2/(2*sigmaphi**2))
       else:
        return 1/(2*(2*pi)**0.5*sigmaphi)*exp(-(pi-mphi-phi)**2/(2*sigmaphi**2))

    def Fn(self, c, b, a, E11, E22, E33):

    return self.f2(E11, E22, E33, b, c)*self.R(a,b,c)*self.Jacobian(a,b,c)*cos(c)**2


    def eval(self, value, x):

    E11=project(self.E[0,0],FunctionSpace(self.mesh, "DG", 0))
    E11 = E11.vector().array()

    E22=project(self.E[1,1],FunctionSpace(self.mesh, "DG", 0))
    E22 = E22.vector().array()

    E33=project(self.E[2,2],FunctionSpace(self.mesh, "DG", 0))
    E33 = E33.vector().array()

        stress11 = []
        print len(E11)
        for k in range(0,len(E11)):
            print k

      #  S11 = lambda c,b,a: self.f2(E11[k], E22[k], E33[k], b, c)*self.R(a,b,c)*self.Jacobian(a,b,c)*cos(c)**2

            print E11[k], E22[k], E33[k]

            res =  tplquad(self.Fn, -pi/24,pi/24,
                                           lambda a : 0, 
                                           lambda a : 2*pi , 
                                           lambda a , b : 0 ,
                                           lambda a , b : pi, args=((E11[k], E22[k], E33[k])), epsabs=1e-2, epsrel=1e-2)[0]#, 

            stress11.append(res)




        stress11 = np.array(stress11)

        S = FunctionSpace(mesh, 'DG', 0)
        results_exp = Function(S)
        results_exp.vector()[:] = stress11;
        value[0] = results_exp(x[0], x[1], x[2]);


def Fmat(u):
    d = u.geometric_dimension()
    I = Identity(d)             # Identity tensor
    F = I + grad(u)             # Deformation gradient
    return F

def Emat(u):
    F = Fmat(u)
    d = u.geometric_dimension()
    I = Identity(d)             # Identity tensor
    return 0.5*(F.T*F - I)


# Define boundary
class Left(SubDomain):
    def inside(self, x, on_boundary):
        return abs(x[0] - 0.0) < DOLFIN_EPS and on_boundary

class Right(SubDomain):
    def inside(self, x, on_boundary):
        return abs(x[0] - 10.0) < DOLFIN_EPS and on_boundary

class Bottom(SubDomain):
    def inside(self, x, on_boundary):
        return abs(x[2] - 0.0) < DOLFIN_EPS and on_boundary

class Top(SubDomain):
    def inside(self, x, on_boundary):
        return abs(x[2] - 1.0) < DOLFIN_EPS and on_boundary

class Lower(SubDomain):
    def inside(self, x, on_boundary):
        return abs(x[1] - 0.0) < DOLFIN_EPS and on_boundary

class Upper(SubDomain):
    def inside(self, x, on_boundary):
        return abs(x[1] - 10.0) < DOLFIN_EPS and on_boundary


os.system('rm nonlinearelasticity*')
parameters["form_compiler"]["quadrature_degree"] = 1

parameters["form_compiler"]["cpp_optimize"] = True
ffc_options = {"optimize": True, \
               "eliminate_zeros": True, \
               "precompute_basis_const": True, \
               "precompute_ip_const": True}


nx = 1  # divisions in r direction
ny = 1  # divisions in theta direction
nz = 1  # divisions in theta direction
mesh = BoxMesh(0, 0, 0, 10, 10, 1, nx, ny, nz)

V = VectorFunctionSpace(mesh, 'CG', 1)


# Initialize sub-domain instances
left = Left()
right = Right()
bottom = Bottom()
top = Top()
lower = Lower()
upper = Upper()

# Initialize mesh function for boundary domains
boundaries = FacetFunction("size_t", mesh)
boundaries.set_all(0)
left.mark(boundaries, 1)
right.mark(boundaries, 2)
bottom.mark(boundaries, 3)
top.mark(boundaries, 4)
lower.mark(boundaries, 5)
upper.mark(boundaries, 6)

ds = Measure("ds")[boundaries]
n = FacetNormal(mesh)

# Set up boundary condition at left end
c = Constant((0.0))
bcleft0 = DirichletBC(V.sub(0), c, boundaries, 1)
bcleft1 = DirichletBC(V.sub(1), c, boundaries, 1)
bcleft2 = DirichletBC(V.sub(2), c, boundaries, 1)
bclower0 = DirichletBC(V.sub(0), c, boundaries, 5)
bclower1 = DirichletBC(V.sub(1), c, boundaries, 5)
bclower2 = DirichletBC(V.sub(2), c, boundaries, 5)
bcbottom2 = DirichletBC(V.sub(2), c, boundaries, 3)
corner = DirichletBC(V, Constant((0,0,0)), boundaries, 7)

uspecified = Expression("u", u = 0.0)
bc_right0 = DirichletBC(V.sub(0), uspecified, boundaries, 2)
bc_upper1 = DirichletBC(V.sub(1), uspecified, boundaries, 6)


# Set up boundary conditions
bcs = [bcleft0, bclower1, bc_right0, bc_upper1, bcbottom2]

## Define variational problem
du = TrialFunction(V)
u = Function(V)
v = TestFunction(V)

kappa = 1e5
b = 0
c = 0
u = Expression(('0.20*x[0]', '0.20*x[1]', '0.20*x[2]'), domain=mesh)
E = Emat(u)
F2 = Integration(E=E,mesh=mesh)

plot(F2, mesh = mesh)
interactive()

---

This is not a working example. Indentation is all messed up.

Besides, please try the suggestions in the answer above.

---

Sorry about it.... I will try your suggestions and let you know if I have any issues...
Thanks !