# Phase Field - The Newton's method does not converge issue

I wrote a program to solve a problem similar to the Cahn-Hilliard equation,
but I encountered an issue where the results do not converge.

Here is the formula:

φ represents phase field.
U represents concentration field.

Here is my code:

``````import math
import random
from dolfin import *

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.2
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)

# Model parameters

dt     = 2.0e-04  # time step
theta  = 0.5      # time stepping family, e.g. theta=1 -> backward Euler, theta=0.5 -> Crank-Nicolson

epsilon_4 = 0.05
fold = 4
alpha = 0
beta = 0
M_PI = math.pi
Cinf = 0.25
k = 0.1
Cl = 0.74
a1 = 5.0 * math.sqrt(2.0) / 8.0
a2 = 47.0 / 75.0
epsilon = 100
Vp = 5.0e-6
d0 = 1.3e-8
D = 1.0e-9
mCinf = 2
m = mCinf / Cinf
G = 100.0e2
lT = m*(1 - k)*Cl / G
theta_inf = 1.0

Dtelda = a1 * a2 * epsilon
Vptelda = a1 * a2 * Vp * d0 / D * pow(epsilon, 2.0)
lTtelda = lT / d0 / epsilon
lamda = a1 * epsilon
tau = a1 * a2 * pow(epsilon, 3.0) * pow(d0, 2.0)

T0 = -mCinf / k
Tinf = T0 + G * Vp * theta_inf * lTtelda / Vptelda

# Step in time
t = 0.0
T = 50*dt

parameters["form_compiler"]["optimize"]     = True
parameters["form_compiler"]["cpp_optimize"] = True

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

du     = TrialFunction(ME)
q1,q2  = TestFunctions(ME)
v1,v2  = TestFunctions(ME)

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

# Split mixed functions
dmu,  duu2 = split(du)
dphi, duu1 = split(du)

phi, uu1 = split(u)
mu, uu2  = split(u)

phi0, uu10 = split(u0)
mu0, uu20 = split(u0)

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

phi = variable(phi)

thetaa = atan(dphidy/dphidx)

AA      = 1 + epsilon_4 * cos(fold * thetaa)
f       = abs(dot(1,phi))**2 * AA * dAAdphi
temp    = (grad(phi)[1] - Vp * t) / lT

uu1_mid = (1.0 - theta) * uu10 + theta * uu1
uu2_mid = (1.0 - theta) * uu20 + theta * uu2

L0 = (1-(1-k)*temp)*(AA**2)*((phi - phi0)/dt)*q1*dx - phi*q1*dx + phi**3*q1*dx + lamda*(1-phi**2)*(mu + temp)*q1*dx - dot(grad(uu1_mid), grad(q1))*dx
L = L0 + L1 + L2 + L3

a = derivative(L, u, du)

problem = CahnHilliardEquation(a, L)
solver = NewtonSolver()
solver.parameters["linear_solver"] = "lu"
solver.parameters["convergence_criterion"] = "incremental"
solver.parameters["relative_tolerance"] = 1e-6

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

while (t < T):
t += dt
u0.vector()[:] = u.vector()
solver.solve(problem, u.vector())
file << (u.split()[0], t)
``````

thank you all