Newton solver did not converge because maximum number of iterations reached

1

Hello,
I’m trying to build a multi component Cahn Hilliard system and I can’t seem to get the the newton solver to converge. I’m not sure what the problem is. Any help is welcome.

Sub domain for Periodic boundary condition

class PeriodicBoundary(SubDomain):

def inside(self, x, on_boundary):
    # return True if on left or bottom boundary AND NOT
    # on one of the two corners (0, 1) and (1, 0)
    return bool((near(x[0], 0) or near(x[1], 0)) and
            (not ((near(x[0], 0) and near(x[1], 1)) or
                    (near(x[0], 1) and near(x[1], 0)))) and on_boundary)

def map(self, x, y):
    if near(x[0], 1) and near(x[1], 1):
        y[0] = x[0] - 1.
        y[1] = x[1] - 1.
    elif near(x[0], 1):
        y[0] = x[0] - 1.
        y[1] = x[1]
    else:   # near(x[1], 1)
        y[0] = x[0]
        y[1] = x[1] - 1.

Create mesh and define function space with periodic boundary conditions

mesh = UnitSquareMesh.create(32,32,CellType.Type.quadrilateral)
A = VectorFunctionSpace(mesh, ‘Lagrange’, 2, dim = 6)

Define functions

u_new = Function(A) # solution for the next step
u_old = Function(A) # solution from previous step
u_new.rename(“fields”,“”)

Split mixed functions

c1_new, c2_new, c3_new, mu1_new = split(u_new)
c1_old, c2_old, c3_old, mu1_old = split(u_old)

Define test functions

tf = TestFunction(A)
q_1, q_2, q_3, v_1 = split(tf)

Class representing the intial conditions

class InitialConditions(UserExpression):
def init(self, **kwargs):
random.seed(6 + MPI.rank(MPI.comm_world))
super().init(kwargs)
def eval(self, values, x):
values[0] = 0.6 + 0.01(0.5 - random.random()) # concentration
values[1] = 0.3 + 0.01(0.5 - random.random()) # concentration
values[2] = 1 - values[0] - values[1]
values[3] = 0.0 # chemical potential
def value_shape(self):
return (6,)

Create intial conditions and interpolate

u_init = InitialConditions(degree=1)
u_new.interpolate(u_init)
u_old.interpolate(u_init)

#ij, jk, ik
chi = np.array([2.5,2.5,5.5])
c1_new = variable(c1_new)
c2_new = variable(c2_new)
c3_new = variable(c3_new)

f_1 = c1_new* ln(c1_new) + chi[0] * c1_newc2_new
f_2 = c2_new ln(c2_new) + chi[1] * c2_newc3_new
f_3 = c3_new ln(c3_new) + chi[2] * c3_new*c1_new

dfdc1 = diff(f_1, c1_new)
dfdc2 = diff(f_2, c2_new)
dfdc3 = diff(f_2, c3_new)

lmbda = 1.0e-02 # surface parameter
dt = 5.0 # time step
theta = 0.5

F_1 = inner(c1_new, q_1) * dx - inner(c1_old, q_1) * dx + inner(grad(c1_old), grad(q_1))*dx
F_2 = inner(c2_new, q_2) * dx - inner(c2_old, q_2) * dx + inner(grad(c2_old), grad(q_2))*dx
F_3 = inner(c3_new, q_3) * dx - inner(c3_old, q_3) * dx + inner(grad(c3_old), grad(q_3))dx
Fm_1 = inner(mu1_new, v_1)dx - inner(dfdc1,v_1)dx - chi[0].5(lmbda**2)(inner(grad(c1_new), grad(v_1) ) )dx
Fm_2 = inner(mu1_new, v_1)dx - inner(dfdc1,v_1)dx - chi[1].5(lmbda**2)(inner(grad(c2_new), grad(v_1) ) )dx
Fm_3 = inner(mu1_new, v_1)dx - inner(dfdc1,v_1)dx - chi[2].5(lmbda**2)(inner(grad(c3_new), grad(v_1) ) )*dx

Conc1 = F_1 + Fm_2
Conc2 = F_2 + Fm_3
Conc3 = F_3 + Fm_1
ConcT = Conc1 + Conc2 + Conc3

Output files for concentration and chemical potential

fileC = File(“data/concentration.pvd”, “compressed”)
fileM = File(“data/chem_potential.pvd”, “compressed”)

Step in time

t = 0
T = 50*dt
while (t < T):
t += dt
u_old.vector()[:] = u_new.vector()
solve(ConcT == 0, u_new)
fileC << (u_new.split()[0], t)
fileM << (u_new.split()[1], t)