Non Linear mixed variational problem- error PETSc Krylov solver

1

Hello,
We are using dolfin for simulate mixedVariationalProblem (1D-2D).
We attach a minimum example showing our problem. Indeed, the resolution of our problem by a linear solver works correctly. However, for the continuation of our work we would like to use the non-linear solver but it does not work with the following message “Error: Unable to solve linear system using PETSc Krylov solver”. I wonder if dolfin manages to calculate the jocabian well or if it is not possible to solve this type of mixed problem with Krylov solver.


from dolfin import *


# Create meshes
n = 10
xmax= 10.0
ymax=5.0
nx = int(n*xmax)
ny = int(n*ymax)
p0 = Point([0.0,0.0])
p1 = Point([xmax,ymax])
mesh = RectangleMesh(p0, p1, nx, ny, 'right')


class Bottom(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[1], 0.0)

class Left(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[0], 0.0)

class Top(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[1], ymax)

class Right(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[0], xmax)


def newton_solver_parameters():
    return{"nonlinear_solver": "newton",
           "newton_solver": {"linear_solver": "gmres"}}

def linear_solver_parameters():
    return {"linear_solver": "gmres"}



bottom = Bottom()
left = Left()
top = Top()
right = Right()

domain_mesh = MeshFunction("size_t", mesh, mesh.topology().dim(), 0)
bottom.mark(domain_mesh,1)
boundaries = MeshFunction("size_t", mesh, mesh.topology().dim() - 1, 0)
boundaries.set_all(0)



bottom.mark(boundaries,1)
left.mark(boundaries,1)
top.mark(boundaries,1)
right.mark(boundaries,1)


submesh2 = MeshView.create(boundaries, 1)

boundaries.set_all(0)
bottom.mark(boundaries,1)
left.mark(boundaries,2)
top.mark(boundaries,3)
right.mark(boundaries,4)

submesh2_Measures = MeshFunction("size_t", submesh2, submesh2.topology().dim(), 0)
map1dto2d=submesh2.topology().mapping()[mesh.id()].cell_map()

for i in range(len(submesh2_Measures.array())):
    submesh2_Measures.array()[i]=boundaries.array()[map1dto2d[i]]


# Create function spaces
W1 = FunctionSpace(mesh, "Lagrange", 1)
W2 = FunctionSpace(submesh2, "Lagrange", 1)
# Define the mixed function space W = W1 x W2
W = MixedFunctionSpace( W1, W2)


# Define mixed-domains variational problem
(v1,v2) = TestFunctions(W)
(u1,u2) = TrialFunctions(W)
#u1=u[0]
#u2=u[1]


u_old = Function(W)
u1_old = u_old.sub(0)
u2_old = u_old.sub(1)

u1_old = interpolate(Constant(1.0),W.sub_space(0))
u2_old = interpolate(Constant(0.0),W.sub_space(1))

dx1 = Measure("dx", domain=W.sub_space(0).mesh())
ds1 = Measure("ds",domain=W.sub_space(0).mesh(), subdomain_data=boundaries)
dx2 = Measure("dx", domain=W.sub_space(1).mesh(), subdomain_data=submesh2_Measures)

dt=0.01
T=50
D2D=10
D1D=10

F1 = u1*v1*dx1  + dt*D2D*inner(grad(u1),grad(v1))*dx1 + dt*u1*v1*ds1 - dt*u2*v1*dx2
F2 = u2*v2*dx2  + dt*D1D*inner(grad(u2),grad(v2))*dx2 + dt*(u2-u1)*v2*dx2
F = F1 + F2
L = u1_old*v1*dx1+u2_old*v2*dx2

sol = Function (W)
print("Linear:")
solve(F == L, sol,solver_parameters=linear_solver_parameters())
print("Non Linear")

u_nl=Function(W)
u_nl1 = u_nl.sub(0)
u_nl2 = u_nl.sub(1)

FNL1 = u_nl1*v1*dx1  + dt*D2D*inner(grad(u_nl1),grad(v1))*dx1 + dt*u_nl1*v1*ds1 - dt*u_nl2*v1*dx2
FNL2 = u_nl2*v2*dx2  + dt*D1D*inner(grad(u_nl2),grad(v2))*dx2 + dt*(u_nl2-u_nl1)*v2*dx2
FNL = FNL1 + FNL2 - L

solve(FNL == 0, sol, solver_parameters=newton_solver_parameters())