Using SNESVI solver without bounds should throw an error

The code below does not raise DOLFIN error but rather cryptic PETSc error.

from dolfin import *

# Load mesh and subdomains
mesh = Mesh("stokes-taylor-hood/dolfin_fine.xml.gz")
sub_domains = MeshFunction("size_t", mesh, "stokes-taylor-hood/dolfin_fine_subdomains.xml.gz")

# Define function spaces
V = VectorFunctionSpace(mesh, "CG", 2)
Q = FunctionSpace(mesh, "CG", 1)
S = VectorFunctionSpace(mesh, "DG", 1, dim=2)
W = MixedFunctionSpace([V, Q, S])
print W.dim()

# No-slip boundary condition for velocity
noslip = Constant((0, 0))
bc0 = DirichletBC(W.sub(0), noslip, sub_domains, 0)

# Inflow boundary condition for velocity
inflow = Expression(("-sin(x[1]*pi)", "0.0"))
bc1 = DirichletBC(W.sub(0), inflow, sub_domains, 1)

# Boundary condition for pressure at outflow
zero = Constant(0)
bc2 = DirichletBC(W.sub(1), zero, sub_domains, 2)

# Collect boundary conditions
bcs = [bc0, bc1, bc2]

# Define variational problem
w = Function(W)
(u, p, s) = split(w)
(v, q, t) = TestFunctions(W)
s = as_tensor(((s[0], s[1]), (s[1], -s[0])))
t = as_tensor(((t[0], t[1]), (t[1], -t[0])))
d = sym(grad(u))

# Linear Stokes for initial guess
F = (inner(s, grad(v)) - div(v)*p + q*div(u) + inner(s-d, t))*dx
solve(F == 0, w, bcs)

# Stokes with non-linear constitutive law
#g = s - Constant(1.0)*inner(d, d)*d
g = s - Constant(1.0)*(1.0+inner(d, d))*d
F = (inner(s, grad(v)) - div(v)*p + q*div(u) + inner(g, t))*dx

# Define the solver parameters
snes_solver_parameters = {"nonlinear_solver": "snes",
                          "linear_solver"   : "lu",
                          "snes_solver"     : { "method": "vinewtonssls",
                                                "line_search": "basic",
                                                "maximum_iterations": 20,
                                                "report": True,
                                                "error_on_nonconvergence": True,
                                              }}

# Set up the non-linear problem
J = derivative(F, w)
problem = NonlinearVariationalProblem(F, w, bcs, J=J)

# Set up the non-linear solver
solver  = NonlinearVariationalSolver(problem)
solver.parameters.update(snes_solver_parameters)
info(solver.parameters, True)
solver.solve()

Alternatively can be reproduced using "snes_solver": {"method": "vinewtonssls"} in demo_contact-vi-snes.py and calling solver.solve() (without bounds).

Comments (5)