Automatic assembly of the Jacobian for problems with mixed function spaces and constrained_domain != None fails

[Matthias Liertzer]

The Jacobian of a form is incorrectly assembled for a problem involving mixed function_spaces and periodic boundary conditions. As a test case I could reproduce the problem for a simple 1D Helmholtz problem with periodic boundary conditions:

#!/usr/bin/env python
from  __future__ import print_function
from dolfin import *

mesh = UnitIntervalMesh(50)

# Sub domain for Periodic boundary condition
class PeriodicBoundary(SubDomain):
    def inside(self, x, on_boundary):
        return bool(near(x[0], 0) and on_boundary)
    def map(self, x, y):
        y[0] = x[0] - 1.0

pbc = PeriodicBoundary()

V = FunctionSpace(mesh, "CG", 1, constrained_domain=pbc)
R = FunctionSpace(mesh, "R", 0, constrained_domain=pbc)
W = V * R

u_mixed = interpolate(Expression(("1.5*sin(k*x[0])","k"), k=4*pi), W)
uw, kw = u_mixed.split()
uv, kv = TestFunctions(W)

F = (inner(nabla_grad(uw), nabla_grad(uv)) - uw * uv * kw**2) * dx + (uw**2 - 1.0) * kv * dx

J = derivative(F, (uw, kw))
solve(F == 0, u_mixed, J=J) # fails with PETSc error 62 (wrong arguments)

The main issue lies in the defective assembly from the form J to a corresponding matrix Jmat, which yields a non-square, and therefore non-invertible matrix:

Jmat = PETScMatrix()
assemble(J, tensor=Jmat)
print("Size of matrix J: ({}, {})".format(Jmat.size(0), Jmat.size(1)))
# yields: Size of matrix J: (51, 52)

For now, the problem can be circumvented by explicitly passing the correct mixed function space to the assemble routine.

class MyNonlinearProblem(NonlinearProblem):
    def __init__(self, F, J, W):
        NonlinearProblem.__init__(self)
        self.formF = F
        self.formJ = J
        self.W = W

    def F(self, vecF, x):
        assemble(self.formF, tensor=vecF)

    def J(self, matJ, x):
        assemble(self.formJ, tensor=matJ, function_spaces=self.W)

problem = MyNonlinearProblem(F, J, W)
solver = NewtonSolver()
solver.solve(problem, u_mixed.vector())

[Martin Sandve Alnæs]

Try specifying the variation to derivative() explicitly:

duk = TrialFunction(W)
J = derivative(F, (uw, kw), du)

The constrained_space property isn't known to the symbolic layer.

[Martin Sandve Alnæs]

What do you expect this to do?

u_mixed = interpolate(Expression(("1.5*sin(k*x[0])","k"), k=4*pi), W)
uw, kw = u_mixed.split()

Now u_mixed and uw,kw are three separate functions. You should do:

u_mixed = interpolate(Expression(("1.5*sin(k*x[0])","k"), k=4*pi), W)
uw, kw = split(u_mixed)

and then later

du_mixed = TrialFunction(W)
J = derivative(F, u_mixed, du_mixed)

[Matthias Liertzer]

Thanks for the fast response!

In the end it all boils down to the unfortunate fact that

u_mixed.split() != split(u_mixed),

which is not really obvious from first sight to someone starting using FEniCS. Note, that the code above actually worked flawlessly except for the case when constrained domains were used.