Deadlock with constant function space

[James Maddison]

Mass matrix assembly deadlocks when using constant function spaces, as illustrated by the following code:

#!/usr/bin/env python
from dolfin import *
set_log_level(DEBUG)
mesh = Rectangle(0.0, 0.0, 1.0, 1.0, 16, 16)
space = FunctionSpace(mesh, "R", 0)
a = assemble(inner(TestFunction(space), TrialFunction(space)) * dx)
L = assemble(inner(TestFunction(space), Constant(1.0)) * dx)
F = Function(space)
sv = LUSolver("default")
sv.solve(a, F.vector(), L)
print assemble(F * dx)

The final log output when run on two processes is:

$ mpirun -np 2 ./test
...
Process 0: Initializing PETSc (ignoring command-line arguments).
Process 1: Initializing PETSc (ignoring command-line arguments).
Process 0: Matrix of size 1 x 1 has 1 (100%) nonzero entries.
Process 0: Elapsed time: 0.000280142 (Build sparsity)
Process 0: Elapsed time: 0.00020504 (Init tensor)
Process 0: Elapsed time: 5.96046e-06 (Delete sparsity)
Process 0: Elapsed time: 9.799e-05 (Assemble cells)
Process 0: Elapsed time: 5.96046e-06 (Apply (matrix))

[Garth Wells (garth-wells)]

The problem is that you're assembling a 1x1 matrix in parallel. It's dubious whether this should be considered a 'bug', but in any case it's unlikely to be fixed for some time since it's not a sensible usage case.

[James Maddison (jamesmadd)]

I'm using an assembly of this kind of part of a larger model. For the purposes of defining the whole model as a variational problem (specifically for use with dolfin-adjoint) one is required to define a very simple solve in this manner (in my case arising from boundary condition lifting). This can obviously be achieved by assembling the relevant quantities and then computing the constant directly, but this breaks the abstraction.

[Patrick Farrell (pefarrell)]

This problem arises because we want to solve a PDE where we don't know the boundary condition, but have an integral constraint on the solution (like the neumann-poisson.py demo, but a bit more complicated).

In that demo, a MixedFunctionSpace is created to solve for the real-valued Lagrange multiplier and the solution together. We can do this, and have done it, but it breaks the symmetry of the problem, which has a massive measured performance penalty.

Instead, the approach we're taking is to solve the PDE with homogeneous boundary conditions, then compute the Lagrange multiplier c (which lies in R0), and then reconstruct the actual solution. We have a variational equation for the Lagrange multiplier, and so pass that equation over to solve(), despite the fact that it's 1x1. However, as James has pointed out, it deadlocks dolfin in parallel, in lots of places.

By the way, if you ever intend to do on-the-fly repartitioning of the domain (for adaptivity, for example), it's quite likely this will happen in some corner case: in my experience of Zoltan, it sometimes decides that it's worth leaving a processor with no degrees of freedom. So this problem will eventually arise in another context.

If we were to go through and fix all the deadlocks for this problem, would that be a candidate for merging? Or is there a parallel-safe way to do this that preserves symmetry?

[Garth Wells (garth-wells)]

On Thursday, 12 July 2012, Patrick Farrell <email address hidden>
wrote:

> This problem arises because we want to solve a PDE where we don't know
> the boundary condition, but have an integral constraint on the solution
> (like the neumann-poisson.py demo, but a bit more complicated).
>
> In that demo, a MixedFunctionSpace is created to solve for the real-
> valued Lagrange multiplier and the solution together. We can do this,
> and have done it, but it breaks the symmetry of the problem, which has a
> massive measured performance penalty.
>
> Instead, the approach we're taking is to solve the PDE with homogeneous
> boundary conditions, then compute the Lagrange multiplier c (which lies
> in R0), and then reconstruct the actual solution. We have a variational
> equation for the Lagrange multiplier, and so pass that equation over to
> solve(), despite the fact that it's 1x1. However, as James has pointed
> out, it deadlocks dolfin in parallel, in lots of places.
>
> By the way, if you ever intend to do on-the-fly repartitioning of the
> domain (for adaptivity, for example), it's quite likely this will happen
> in some corner case: in my experience of Zoltan, it sometimes decides
> that it's worth leaving a processor with no degrees of freedom. So this
> problem will eventually arise in another context.
>
> If we were to go through and fix all the deadlocks for this problem,
> would that be a candidate for merging?

Yes.

I expect that it will always be somewhat fragile because not all back ends
support having no data on a process, e.g. SCOTCH does support this but
ParMETIS does not.

Garth

> Or is there a parallel-safe way
> to do this that preserves symmetry?
>
> --
> You received this bug notification because you are a member of DOLFIN
> Core Team, which is subscribed to DOLFIN.
> https://bugs.launchpad.net/bugs/1020889
>
> Title:
> Deadlock with constant function space
>
> To manage notifications about this bug go to:
> https://bugs.launchpad.net/dolfin/+bug/1020889/+subscriptions
>

--
Garth N. Wells
Department of Engineering, University of Cambridge
http://www.eng.cam.ac.uk/~gnw20