Fast projection in case of explicit relations

Question

I have a CG2 function, the "natural" gradient and Hessian of which is a DG1 vector and DG0 tensor, respectively. I.e. one should not have to solve any variational problem to calculate these quantities.

When I calculate the Hessian and project it onto a DG0 field, it however takes a long time (considering it is an explicit operation). I have constructed an example below; it is the last line that I want to speed up.

from dolfin import *
mesh = RectangleMesh(0,0,1,1,40,40) 
u = interpolate(Expression('sin(x[1])*cos(x[0])'),FunctionSpace(mesh,'CG',2))
H = project(grad(grad(u)),TensorFunctionSpace(mesh,'DG',0))

Answer 1

Hi,

There are quite many ways to speed this up. If you are going to compute this projection many times, then much work can be preassembled. Some suggestions are:

from dolfin import *
mesh = RectangleMesh(0,0,1,1,40,40) 
V = FunctionSpace(mesh,'CG',2)
u = interpolate(Expression('sin(x[1])*cos(x[0])'), V)
S = TensorFunctionSpace(mesh,'DG',0)

t0 = Timer("Original")
H = project(grad(grad(u)), S)
print norm(H)
H.vector().zero()

t0 = Timer("Compress diagonal coefficient matrix  for faster solve")
A = assemble(inner(TrialFunction(S), TestFunction(S))*dx)
b = assemble(inner(grad(grad(u)), TestFunction(S))*dx)
A.compress()
solve(A, H.vector(), b, "cg", "default")
print norm(H)
H.vector().zero()

t0 = Timer("Extract diagonal from coefficient matrix - solve directly")
A = assemble(inner(TrialFunction(S), TestFunction(S))*dx)
b = assemble(inner(grad(grad(u)), TestFunction(S))*dx)
ones = Function(S)
ones.vector()[:] = 1
A_diag = A * ones.vector()
A_diag.set_local(1.0/A_diag.array())
H = Function(S)
H.vector()[:] = b * A_diag
t0.stop()
print norm(H)
H.vector().zero()

###
M = assemble(inner(grad(grad(TrialFunction(V))), TestFunction(S))*dx)
M.compress()
A = assemble(inner(TrialFunction(S), TestFunction(S))*dx)
ones = Function(S)
ones.vector()[:] = 1
A_diag = A * ones.vector()
A_diag.set_local(1.0/A_diag.array())

t0 = Timer("Extract diagonal + precompute all required")
b = M*u.vector()
H.vector()[:] = b * A_diag
t0.stop()
print norm(H)

t0 = Timer("Local Solver")
ls = LocalSolver()
ls.solve(H.vector(), Form(inner(TrialFunction(S), TestFunction(S))*dx), Form(inner(grad(grad(u)), TestFunction(S))*dx))
print norm(H)

list_timings()

Comments

I tried the "direct solve"-option, but it seems like the problem is

S = TensorFunctionSpace(mesh,'DG',0)

This basically takes all of the time, and the code is used for (anisotropic) mesh adaptation, so it has to be recomputed. I think it is related to the "Init dofmap", and I guess there is no way to get around this?

If I run your code with

RectangleMesh(0,0,1,1,160,160) 

and time the definition of V and S, there is a factor of ~50 difference. There is only a factor of 2 difference in the number of degrees of freedom.

---

Ok. I think you need to reformulate your question then, or perhaps ask another one. I don't really know how to speed up the creation of a TensorFunctionSpace.