overhead of AdaptiveLinearVariationalSolver seems way to high

Question

The code below is a small modification of the FEniCS demo probgram demo_auto-adaptive-poisson.py, with the major difference being that it solves the Dirichlet problem for Laplace's equation on the unit ball in 3D. Here are the timings from its output:

Time spent for adaptive solve (in seconds):

  Level  |  solve_primal  estimate_error  compute_indicators    mark_mesh  adapt_mesh    update
  ---------------------------------------------------------------------------------------------
  0      |    0.00477398        0.240659           0.0669882  0.000335102   0.0395709  0.272424
  1      |     0.0317322         1.88749            0.451885   0.00264225    0.241248   1.55086
  2      |      0.220263         11.1167             2.48592     0.020308     1.34969   8.18865
  3      |       1.35856         58.3366               12.65     0.177529     6.81677   37.6127

I find the overhead for the adaptively astonishingly large. On the finest mesh, with 75,449 DOFs, the solve takes 1.4 seconds, but the adaptivity takes nearly two minutes, so the cost of adaptivity is about 80 times the cost of the solve. The error estimation alone on the final mesh takes 58 seconds.

Am I doing something wrong? Is there a way to substantially reduce the time for the adaptivity? Otherwise, I think one would almost always be better off just refining the mesh non-adaptively and solving on the finer mesh.

from dolfin import *
from mshr import *
# Dirichlet problems for Laplace equation on unit ball
Omega = Sphere(Point(0., 0., 0.), 1.)
mesh = generate_mesh(Omega, 10)
V = FunctionSpace(mesh, 'Lagrange', 1)  
u = TrialFunction(V)
v = TestFunction(V)
a = (inner(grad(u), grad(v)) + u*v)*dx
L = Constant(0.)*v*dx
g = Expression('pow(x[0], 2)', element=V.ufl_element())    
bc = DirichletBC(V, g, DomainBoundary())
u = Function(V)
problem = LinearVariationalProblem(a, L, u, bc)
# define a goal functional and solve adaptively
M = u*dx
solver = AdaptiveLinearVariationalSolver(problem, M)
solver.parameters["max_iterations"] = 4
solver.parameters["linear_variational_solver"]["linear_solver"] = "cg"
solver.parameters["linear_variational_solver"]["preconditioner"] = "amg"
solver.parameters["linear_variational_solver"]["symmetric"] = True
solver.parameters["error_control"]["dual_variational_solver"]["linear_solver"] = "cg"
solver.parameters["error_control"]["dual_variational_solver"]["preconditioner"] = "amg"
solver.parameters["error_control"]["dual_variational_solver"]["symmetric"] = True
tol= 1.E-5
solver.solve(tol)
# output
solver.summary()
print("computed integral: {:7.5f}  true: {:7.5f}".format(assemble(u*dx), 1.22692))

Answer 1

I don't think you are doing anything wrong. This is simply an
illustration of the fact that some of the adaptivity algorithms do not
scale very well, and this is particularly noticeable in 3D.

I can say a little more about what takes time with the error
estimation. The error estimate is based on computing an dual/adjoint
solution of the same order as the original solution and then
extrapolating that solution into a higher order space and then some
assembly. This extrapolation algorithm is the main culprit in your
example, its run-time really dominates everything else. Sample timings:

Level | A | B | C >= A + B

0 | 0.00595564 | 0.272692 | 0.279996
1 | 0.0383623 | 2.17221 | 2.21893
2 | 0.259767 | 12.2055 | 12.5085
3 | 1.98911 | 62.2913 | 64.4823

where

A = Time to compute dual solution of same order (s)
B = Time to extrapolate dual solution
C = Time to estimate error (involves A, then B, then a minor things)

Please feel free to file an issue regarding scalability of adaptivity
algorithms in general.

But in any case, if you have rather smooth solution and a rather
global goal functional (which seems to be the case in this example), I
would not recommend an adaptive approach.

Some further comments:

• I would use

  g = Expression('pow(x[0], 2)', degree=2)

  rather than

  g = Expression('pow(x[0], 2)', element=V.ufl_element())

• You probably want to evaluate the functional at the final node in
  the adaptive hierarchy. You can access the final solution via
  u.leaf_node(), and so

  print("computed integral: {:7.5f} true: {:7.5f}".format(assemble(u.leaf_node()*dx), 1.22692))