assemble step takes a long time relative to solver

Question

I've been doing some time tests and found that the vast majority of the time involved in running my code goes into assembling the matrix from the bilinear form. For example, the code below which solves a 3D elasticity problem on a 30x30x30 unit cube mesh takes almost two minutes to assemble the matrix and only a few seconds to solve the linear system. Is this to be expected? Is there anyway to decrease the assembly time? I've tried including

parameters['form_compiler']['optimize'] = True

but that doesn't increase the speed. Here is my code:

from dolfin import *
import numpy, copy, time

parameters['form_compiler']['optimize'] = True
parameters['form_compiler']['cpp_optimize'] = True

start_time = time.time()


# define periodic boundary
class PeriodicBoundary(SubDomain):

def inside(self, x, on_boundary):
    # return True if on left or bottom boundary AND NOT on one of the two slave edges
    return bool((near(x[0], 0.) or near(x[1], 0.) or near(x[2], 0.)) and
                (not ((near(x[1], 1.) and near(x[2], 0.)) or
                      (near(x[1], 0.) and near(x[2], 1.)) or
                      (near(x[0], 0.) and near(x[2], 1.)) or
                      (near(x[0], 1.) and near(x[2], 0.)) or
                      (near(x[0], 0.) and near(x[1], 1.)) or
                      (near(x[0], 1.) and near(x[1], 0.)))) and on_boundary)

def map(self, x, y):
    if near(x[0], 1.) and near(x[1], 1.) and near(x[2], 1.):
        y[0] = x[0] - 1.
        y[1] = x[1] - 1.
        y[2] = x[2] - 1.
    elif near(x[1], 1.) and near(x[2], 1.):
        y[0] = x[0]
        y[1] = x[1] - 1.
        y[2] = x[2] - 1.
    elif near(x[0], 1.) and near(x[1], 1.):
        y[0] = x[0] - 1.
        y[1] = x[1] - 1.
        y[2] = x[2]
    elif near(x[0], 1.) and near(x[2], 1.):
        y[0] = x[0] - 1.
        y[1] = x[1]
        y[2] = x[2] - 1.
    elif near(x[0], 1.):
        y[0] = x[0] - 1.
        y[1] = x[1]
        y[2] = x[2]
    elif near(x[1], 1.):
        y[0] = x[0]
        y[1] = x[1] - 1.
        y[2] = x[2]
    elif near(x[2], 1):
        y[0] = x[0]
        y[1] = x[1]
        y[2] = x[2] - 1.
    else:
        y[0] = -1000.
        y[1] = -1000.
        y[2] = -1000.

# Create mesh
nn = 30
mesh = UnitCubeMesh(nn, nn, nn)
V = VectorFunctionSpace(mesh, "Lagrange", 1, constrained_domain=PeriodicBoundary())
R = VectorFunctionSpace(mesh, "R", 0, constrained_domain=PeriodicBoundary())
P = MixedFunctionSpace([V, R])

class Fiber(SubDomain):
    def inside(self, x, on_boundary):
        return True if ((x[0]-0.5)**2+(x[1]-0.5)**2 <= 0.25**2) else False
fib = Fiber()

domains = CellFunction("size_t", mesh)
domains.set_all(0)
fib.mark(domains, 1)
dy = Measure("dx")[domains]

# define tensors
bbEm = numpy.zeros(81)
bbEm = bbEm.reshape(3,3,3,3)
bbEf = numpy.zeros(81)
bbEf = bbEf.reshape(3,3,3,3)
delta = [[1,0,0],[0,1,0],[0,0,1]]

# Elastic parameters
Em, Ef, nu = 10.5, 150000.0, 0.3
mum, lmbdam = Em/(2*(1 + nu)), Em*nu/((1 + nu)*(1 - 2*nu))
muf, lmbdaf = Ef/(2*(1 + nu)), Ef*nu/((1 + nu)*(1 - 2*nu))
# define bbEm
for ii in range(3):
    for jj in range(3):
        for kk in range(3):
            for ll in range(3):
                bbEm[ii,jj,kk,ll] = mum*(delta[ii][kk]*delta[jj][ll]+delta[ii][ll]*delta[kk][jj])+lmbdam*delta[ii][jj]*delta[kk][ll]
# define bbEf
for ii in range(3):
    for jj in range(3):
        for kk in range(3):
            for ll in range(3):
                bbEf[ii,jj,kk,ll] = muf*(delta[ii][kk]*delta[jj][ll]+delta[ii][ll]*delta[kk][jj])+lmbdaf*delta[ii][jj]*delta[kk][ll]
bbEmC = Constant(bbEm)
bbEfC = Constant(bbEf)

# define basis for V
e0 = Constant([1,0,0])
e1 = Constant([0,1,0])
e2 = Constant([0,0,1])
e = [e0, e1, e2]

# Define functions
(u, c) = TrialFunctions(P)
(v, d) = TestFunctions(P)

# Set up linear problem
a = bbEmC[i,j,k,l]*grad(u[k])[l]*grad(v[i])[j]*dy(0) \
    + bbEfC[i,j,k,l]*grad(u[k])[l]*grad(v[i])[j]*dy(1) \
    + (inner(c, v) + inner(u, d))*dy


print "time before assemble", time.time() - start_time
A = assemble(a)
print "time after assemble", time.time() - start_time
z = Function(P)
b = None

# solve the problem
ii, jj = 0, 0
L = - bbEmC[i,j,k,l]*outer(e[ii],e[jj])[k,l]*grad(v[i])[j]*dy(0) \
    - bbEfC[i,j,k,l]*outer(e[ii],e[jj])[k,l]*grad(v[i])[j]*dy(1)
b = assemble(L, tensor=b)
solve(A, z.vector(), b, 'cg', 'default')

print "end time", time.time() - start_time

Answer 1

Call list_timings() at the end to get a summary of where the time is spent. Looks like the culprit is sparsity pattern building:

end time 18.3907630444
Summary of timings                |  Average time  Total time  Reps
-------------------------------------------------------------------
Build sparsity                    |        7.2892      14.578     2

Comments

It's because of the mixed space with Real. I believe this is a known issue.

from dolfin import *

nn = 30
mesh = UnitCubeMesh(nn, nn, nn)

V = VectorFunctionSpace(mesh, "Lagrange", 1)
R = VectorFunctionSpace(mesh, "R", 0)
P = MixedFunctionSpace([V, R])

p = TrialFunction(P)
q = TestFunction(P)
a = inner(p,q)*dx
A = assemble(a)

list_timings()

---

Thanks for the answer. Do you know a way to deal with pure periodic boundary conditions without using the mixed function space?

---

I've never used the periodic boundary conditions myself. But it's the use of Real that kills the sparsity pattern builder, so maybe try to see if you can avoid that?

---

I just found that the method in

http://fenicsproject.org/documentation/dolfin/dev/python/demo/documented/singular-poisson/python/documentation.html

solved my problem. It is much faster.