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

**URL:** <https://fenicsproject.discourse.group/t/phase-field-the-newtons-method-does-not-converge-issue/12148>\
**Category:** mathematics\
**Created:** [September 8, 2023, 9:06am UTC](https://fenicsproject.discourse.group/t/phase-field-the-newtons-method-does-not-converge-issue/12148 "2023-09-08T09:06:50Z")\
**Posts on this page:** 1\
**Page:** 1

<div class="post-metadata">

**Author:** ![CHANG\_KUAN-HUI](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/chang_kuan-hui/32/5487_2.png) [@CHANG\_KUAN-HUI](https://fenicsproject.discourse.group/u/CHANG_KUAN-HUI)\
**Post date:** [September 8, 2023, 9:06am UTC](https://fenicsproject.discourse.group/t/phase-field-the-newtons-method-does-not-converge-issue/12148/1 "2023-09-08T09:06:50Z")

</div>

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:

 ![image](https://global.discourse-cdn.com/free1/uploads/fenicsproject1/original/2X/e/ed3817da190a71fe9c1c36fd4edfdbeef58bfdea.png)  
**φ** represents phase field.  
U represents concentration field.

Here is my code:

```auto
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

mesh = UnitSquareMesh.create(96, 96, CellType.Type.quadrilateral)
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)

dphidx = grad(phi)[0]
dphidy = grad(phi)[1]
thetaa = atan(dphidy/dphidx)

AA = 1 + epsilon_4 * cos(fold * thetaa)
dAAdphi = diff(AA, phi)
f = abs(dot(1,phi))**2 * AA * dAAdphi
dfdx = grad(f)[0]
dfdy = grad(f)[1]
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
L1 = uu1*q2*dx - AA**2*dot(grad(phi), grad(q2))*dx + dot(grad(phi), grad(phi))*AA*dAAdphi*q2*dx
L2 = 0.5*((1+k)-(1-k)*phi)*((mu-mu0)/dt)*v1*dx - (1+(1-k)*mu)*0.5*((phi - phi0)/dt)*v1*dx - dot(grad(uu2_mid), grad(v1))*dx
L3 = uu2*v2*dx - D*(1-phi)*0.5*dot(grad(mu), grad(v2))*dx - ((pow(2,0.5))/4)*(1+(1-k)**mu)*((phi - phi0)/dt)*dot(grad(phi), grad(v2))/sqrt(dot(grad(phi), grad(phi)))*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