PETSc problem when solving nonlinear variational problem

Question

I'm solving a mixed nonlinear variational problem and the iteration stops after the very first one, giving the following output:

Solving nonlinear variational problem.
Newton iteration 0: r (abs) = 7.503e+01 (tol = 1.000e-10) r (rel) = 1.000e+00 (tol = 1.000e-09)
Traceback (most recent call last):
File "lmain.py", line 28, in
solver.solve()
RuntimeError:
*** Error: Unable to successfully call PETSc function 'MatSetValuesLocal'.
*** Reason: PETSc error code is: 63 (Argument out of range).
*** Where: This error was encountered inside /build/dolfin-4SStI2/dolfin-2016.2.0/dolfin/la/PETScMatrix.cpp.

Below is the minimal working code snippet.

from dolfin import *
import numpy as np

mesh = UnitSquareMesh(32,32)
U = VectorElement("Lagrange",mesh.ufl_cell(), 1, dim=2) # displacement
M = FiniteElement("Lagrange",mesh.ufl_cell(), 1)       # mu
W = VectorElement("Lagrange",mesh.ufl_cell(), 1, dim=2) # Lagrangian multipliers
S = FunctionSpace(mesh, MixedElement([U,M,W]))

s = Function(S)
s = interpolate(Constant((1.0,1.0,1.0,1.0,1.0)), S)
u, mu, lamd = s.split()
L = inner(u,u)*dx + inner(grad(mu),grad(mu))*dx + \
inner(mu*(grad(lamd)+grad(lamd).T), grad(u)+grad(u).T)*dx + inner(div(lamd),div(u))*dx

ds = TestFunction(S)
du = TrialFunction(S)
F = derivative(L, s, ds)
JF = derivative(F, s, du)

bc0 = DirichletBC(S.sub(0), Constant([0.0,0.0]), "on_boundary")
bc1 = DirichletBC(S.sub(1), Constant(0.0), "on_boundary")
bc2 = DirichletBC(S.sub(2), Constant([0.0,0.0]), "on_boundary")
bc=[bc0, bc1, bc2];

problem = NonlinearVariationalProblem(F,s,bc,JF)
solver = NonlinearVariationalSolver(problem)
solver.solve()
u_s, mu_s, lamd_s = s.split()

Now the only related message I'm able to locate is this one, yet it's still not clear to me how to apply that solution to my problem. Currently I'm not even sure how to approach or comprehend this error message. Any help is greatly appreciated!

Answer 1

The error is due to the line

u, mu, lamd = s.split()

The unfortunate naming of the two functions s.split() and split(s) suggests they are the same, but this not the case. Here, it is split(s) that should be used, i.e.

u, mu, lamd = split(s)

Generally, you want to use split(s) to get subfunctions for use in variational forms, and use s.split() to get subfunctions for plotting or other postprocessing.

Comments

Thanks so much for your suggestion! It works now. This could take me forever to find out...