Application of Laplace operator very slow

1

Hi all,

for post-processing, I’m repeatedly applying the Laplace operator to a function like

u = Function(V)
out = Function(V)

lap = div(grad(u))

for i in range(1e4):
  # some computations to get u ...
  out.vector()[:] = project(lap, V).vector()[:]

Running this for large meshes in parallel is incredbly slow (much slower than the computation of u itself). Is there a smarter way to do this?

Thank you very much!

2

Hi, a call to project includes assembly of the mass matrix of the space, the right hand side vector and finally a solve of the resulting linear system. For repeated projections this is inefficient. Note that the matrix stays needs to be assembled only once. Similarly, you can reuse the solver.

from dolfin import *

mesh = UnitSquareMesh(128, 128)

V = FunctionSpace(mesh, 'CG', 2)
u = TrialFunction(V)
v = TestFunction(V)

f = Function(V)
# Projection lhs and rhs
lhs = inner(u, v)*dx
rhs = inner(div(grad(f)), v)*dx

A, b = map(assemble, (lhs, rhs))
solver = PETScKrylovSolver('cg', 'amg')
# Build solver once
solver.set_operator(as_backend_type(A))
solver.parameters['relative_tolerance'] = 1E-12
solver.parameters['monitor_convergence'] = True

laplace = Function(V)
x = laplace.vector()
for i in range(10):
    # Update
    f.assign(interpolate(Constant(i+1), V))
    # Only b need to be recomputed
    assemble(rhs, tensor=b)
    niters = solver.solve(x, b)

    print(i, x.norm('l2'), 'with', niters)
3

Thank you very much! Great idea!