Setting Dirichlet boundary conditions in Cahn-Hilliard example

General
1

I’m pretty much following the Cahn-Hilliard example, but trying to set Dirichlet boundary conditions for the concentration, initially only on one side. I would have naively thought that zero-flux for the chemical potential would work, so only setting the concentration condition:

# Class representing the intial conditions
class InitialConditions(UserExpression):
    def __init__(self, **kwargs):
        random.seed(2 + MPI.rank(MPI.comm_world))
        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,)

mesh = UnitSquareMesh.create(96, 96, CellType.Type.quadrilateral)

# Class for interfacing with the Newton solver
class CahnHilliardEquation(NonlinearProblem):
    def __init__(self, a, L, bcs):
        NonlinearProblem.__init__(self)
        self.L = L
        self.a = a
        self.bcs = bcs  # Store boundary conditions
    def F(self, b, x):
        assemble(self.L, tensor=b)
        for bc in self.bcs:  # Apply boundary conditions to the residual
            bc.apply(b)
    def J(self, A, x):
        assemble(self.a, tensor=A)
        for bc in self.bcs:  # Apply boundary conditions to the Jacobian
            bc.apply(A)

# Model parameters
kappa  = 0.25*10e-04  # surface parameter
dt     = 5e-05  # time step
X = 0.99*7/2 # Chi Flory Huggins parameter
# time stepping family, e.g.: 
# theta=1 -> backward Euler, theta=0.5 -> Crank-Nicolson
theta  = 0.5      
# Form compiler optionsmai
parameters["form_compiler"]["optimize"]     = True
parameters["form_compiler"]["cpp_optimize"] = True

P1 = FiniteElement("Lagrange", mesh.ufl_cell(), 1)
ME = FunctionSpace(mesh, P1*P1)

# Define boundary condition
def boundary(x, on_boundary):
    return on_boundary

def left_boundary(x, on_boundary):
    return on_boundary and near(x[0], 0)

def right_boundary(x, on_boundary):
    return on_boundary and near(x[0], 1)


# Define the boundary condition for each component if needed
# bc_c = DirichletBC(ME.sub(0), Constant(0), boundary)  # For 'c'
# bc_mu = DirichletBC(ME.sub(1), Constant(0), boundary)  # For 'mu'
# bc_mu_left = DirichletBC(ME.sub(1), Constant(1), left_boundary)  # For 'mu'
# bc_mu_right = DirichletBC(ME.sub(1), Constant(1), right_boundary)  # For 'mu'
# bc_c_left = DirichletBC(ME.sub(0), Constant(1), left_boundary)  # Concentration 0.9 on the left edge
bc_c_right = DirichletBC(ME.sub(0), Constant(0.0005), right_boundary)


# You can apply both conditions if both fields need the same Dirichlet condition
# bcs = [bc_c_left, bc_mu_left, bc_c_right, bc_mu_right]
bcs = [bc_c_right]

# 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
tis = time.time()
# Split mixed functions
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)
print(time.time()-tis)
# Compute the chemical potential df/dc
c = variable(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*c*(1-c)*dot(grad(mu_mid), grad(q))*dx
L1 = mu*v*dx - (ln(c/(1-c))+X*(1-2*c))*v*dx - kappa*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, bcs)
solver = NewtonSolver()
# solver.parameters["linear_solver"] = "lu"
solver.parameters["linear_solver"] = 'gmres'
solver.parameters["preconditioner"] = 'ilu'
# solver.parameters["convergence_criterion"] = "incremental"
solver.parameters["relative_tolerance"] = 1e-6

# Output file
# file = File("output.pvd", "compressed")
file_c = XDMFFile('2comp.xdmf')

# Step in time
t = 0.0
T = 200*dt
ti = time.time()
while (t < T):
#     file << (u.split()[0], t)
    file_c.write(u.split()[0], t)
    t += dt
    u0.vector()[:] = u.vector()
    solver.solve(problem, u.vector())
    print(time.time() - ti)
file_c.close()

However, the solver crashes in this case after a few iterations, unless I specify zero as the constant. Do I have some problem in the setup or is this a maths issue?
Thanks!