Can't solve Cahn-Hillard equation with gmres hypre_euclid

Errors
1

Hi everyone,
Not sure if this is a problem, but I’m wondering why I can’t use gmres with the hypre_euclid preconditioner to solve the problem defined in the Cahn-Hillard example.

Indeed, if I run the following (which is just a slight modification of the original script that is compliant with the current fenics version and modified to use gmres and hypre_eucluid):

import random
from fenics import *


# Class representing the intial conditions
class InitialConditions(UserExpression):
    def __init__(self):
        super(InitialConditions, self).__init__()
        random.seed(2 + MPI.comm_world.Get_rank())

    def eval(self, values, x):
        values[0] = 0.63 + 0.02*(0.5 - random.random())
        values[1] = 0.0

    def value_shape(self):
        return (2,)


# Model parameters
lmbda  = 1.0e-02  # surface parameter
dt     = 5.0e-06  # time step
theta  = 0.5      # time stepping family, e.g. theta=1 -> backward Euler, theta=0.5 -> Crank-Nicolson

# Create mesh and define function spaces
mesh = UnitSquareMesh(96, 96)
CG1 = FiniteElement("CG", triangle, 1)
mixed_element = MixedElement([CG1, CG1])
ME = FunctionSpace(mesh, mixed_element)

# Define trial and test functions
du    = TrialFunction(ME)
q, v  = TestFunctions(ME)

# Define functions
u   = Function(ME)  # current solution
u0  = Function(ME)  # solution from previous converged step

# Split mixed functions
dc, dmu = split(du)
c,  mu  = split(u)
c0, mu0 = split(u0)

# Create intial conditions and interpolate
u_init = InitialConditions()
u.interpolate(u_init)
u0.interpolate(u_init)

# Compute the chemical potential df/dc
c = variable(c)
f    = 100*c**2*(1-c)**2
dfdc = diff(f, c)

# mu_(n+theta)
mu_mid = (1.0-theta)*mu0 + theta*mu

# Weak statement of the equations
L0 = c*q*dx - c0*q*dx + dt*dot(grad(mu_mid), grad(q))*dx
L1 = mu*v*dx - dfdc*v*dx - lmbda*dot(grad(c), grad(v))*dx
L = L0 + L1

# Compute directional derivative about u in the direction of du (Jacobian)
a = derivative(L, u, du)

# Create nonlinear problem and Newton solver
problem = NonlinearVariationalProblem(L, u, J=a)
solver = NonlinearVariationalSolver(problem)
solver.parameters['newton_solver']['linear_solver'] = "gmres"
solver.parameters['newton_solver']['preconditioner'] = "hypre_euclid"

# Step in time
t = 0.0
T = 50*dt
while (t < T):
    t += dt
    u0.vector()[:] = u.vector()
    solver.solve()

I get the following error:

RuntimeError: 

*** -------------------------------------------------------------------------
*** DOLFIN encountered an error. If you are not able to resolve this issue
*** using the information listed below, you can ask for help at
***
***     fenics-support@googlegroups.com
***
*** Remember to include the error message listed below and, if possible,
*** include a *minimal* running example to reproduce the error.
***
*** -------------------------------------------------------------------------
*** Error:   Unable to solve linear system using PETSc Krylov solver.
*** Reason:  Solution failed to converge in 0 iterations (PETSc reason DIVERGED_NANORINF, residual norm ||r|| = 0.000000e+00).
*** Where:   This error was encountered inside PETScKrylovSolver.cpp.
*** Process: 0
*** 
*** DOLFIN version: 2019.2.0.dev0
*** Git changeset:  unknown
*** -------------------------------------------------------------------------

Thank you in advance.