Solving Cahn-Hilliard equation with a multi-grid method

Linear Algebra
1

Solving Cahn-Hilliard equation with a multi-grid method.

Hi, I am a beginner at FEniCS.
I have tried solving the Chan-Hilliard equation with a multi-grid solver.
But It is not stable, I must set up a very smaller time-step than the standard Krylov-solver.
I think my implementation is wrong. But I have no idea. Please help me, thanks.
Following is my source code.

import random
from dolfin import *
from mpi4py import MPI
comm = MPI.COMM_WORLD

N = 60
mesh_1 = UnitSquareMesh(N, N)
mesh_2 = refine(mesh_1)
mesh_3 = refine(mesh_2)
meshes = [mesh_1, mesh_2, mesh_3]
P1 = FiniteElement("Lagrange", mesh_1.ufl_cell(), 1)
P2 = FiniteElement("Lagrange", mesh_2.ufl_cell(), 1)
P3 = FiniteElement("Lagrange", mesh_3.ufl_cell(), 1)
MEs = [FunctionSpace(mesh_1, P1*P1), FunctionSpace(mesh_2, P2*P2), FunctionSpace(mesh_3, P3*P3)]
print(mesh_3.num_vertices())

# Create collection of PETSc DM objects
dm_collection = PETScDMCollection(MEs)

class InitialConditions(UserExpression):
    def __init__(self, **kwargs):
        random.seed(2)
        super().__init__(**kwargs)
    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,)

class CahnHilliardEquation(NonlinearProblem):
    def __init__(self, a, L):
        NonlinearProblem.__init__(self)
        self.L = L
        self.a = a
    def F(self, b, x):
        assemble(self.L, tensor=b)
    def J(self, A, x):
        assemble(self.a, tensor=A)

class KrylovSolver(NewtonSolver):
    def __init__(self):
        NewtonSolver.__init__(self, comm, PETScKrylovSolver(), PETScFactory.instance())
    def solver_setup(self, A, P, problem, iteration):
        self.linear_solver().set_operator(A)

        # Set PETSc solver type
        PETScOptions.set("ksp_type", "gmres")
        PETScOptions.set("pc_type", "mg")

        # Set PETSc MG type and levels
        PETScOptions.set("pc_mg_levels", 3)
        PETScOptions.set("pc_mg_galerkin")

        # Set smoother
        PETScOptions.set("mg_levels_ksp_type", "chebyshev")
        PETScOptions.set("mg_levels_pc_type", "jacobi")

        # Set tolerance and monitor residual
        PETScOptions.set("ksp_monitor_true_residual")
        PETScOptions.set("ksp_atol", 1.0e-12)
        PETScOptions.set("ksp_rtol", 1.0e-12)

        self.linear_solver().set_from_options()
        self.linear_solver().set_dm(dm_collection.get_dm(-1))
        self.linear_solver().set_dm_active(False)


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

parameters["linear_algebra_backend"] = "PETSc"

du    = TrialFunction(MEs[-1])
q, v  = TestFunctions(MEs[-1])

u   = Function(MEs[-1])  # current solution
u0  = Function(MEs[-1])  # solution from previous converged step

dc, dmu = split(du)
c,  mu  = split(u)
c0, mu0 = split(u0)

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

c = variable(c)
f    = 100*c**2*(1-c)**2
dfdc = diff(f, c)

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 = CahnHilliardEquation(a, L)
solver = KrylovSolver()

# Output file
file = File("output.pvd", "compressed")

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