# Newton Solver did not converge phase field model

Dear all,

I´m trying to solve the following nonlineal problem:

dxi/dt = -L*(g'(xi) -kappa * Laplacian(xi))
- L*h'(xi)*(exp(0.5*Fa*phi/R*T)
- (c(w)/c0)*(1-h(xi))*exp(-0.5*Fa*phi/R*T))

- (dxi/dt)*ft(w) ]

Where xi, w and phi are the variables to find.

The functions are defined in the code:

``````from fenics import *
from dolfin import MPI
import matplotlib.pyplot as plt
import numpy as np
from datetime import datetime

# Domain dimensions
lox= 200
loy= 200
#co = 0.068             # mol/m^3

Fa = 96485.3365  # C/mol
R = 8.314472     # J/(molxK)
T = 293.0        # K
R = Constant(R)
T = Constant(T)
Fa = Constant(Fa)

#Define expressions used in variational form
L = Constant(6.25)
kappa = Constant(0.3) # =R*T
Ls = Constant(0.001)
alpha = Constant(0.5)
AA = Constant(38.69)   # nF/RT
W = Constant(2.4)
es = Constant(-13.8)
el = Constant(2.631)  # mu0 = -R*T*ln(c0),, el = ln(cl)
A = Constant(1.0) # =R*T/R*T
c0 = Constant(1.0/14.89)
dv = Constant(5.5)   # Csm/Clm
M0 = Constant(317.9)
S1 = Constant(1000000)
S2 = Constant(1.19)
ft2 = Constant(0.0074)

#Mesh
# Create mesh and define function space
nx, ny = 2000, 2000
mesh = RectangleMesh(Point(0, 0), Point(lox, loy), nx, ny)
P1 = FiniteElement('P', triangle, 2)
V = FunctionSpace(mesh, MixedElement([P1,P1,P1]))

#Define trial and  test functions
v_1, v_2, v_3 = TestFunctions(V)

# Define functions for solutions at previous and at current time
u = Function(V)                # At current time
u_n = Function(V)

#Split system function to access the components
xi, w, phi = split(u)
xi_n, w_n, phi_n = split(u_n)

#Parameter for initial conditions
phie= -0.45

ruido=Expression('-0.04*sin(0.1*x[1])', degree= 2)

#Create initial conditions
u_init=Expression(('0.5*(1.0-1.0*tanh((x[0]-20.0-ruido)*2))','0.0','-0.225*(1.0-tanh((x[0]-20.0-ruido)*2))'), degree= 2, ruido=ruido)

u_n = project(u_init,V)

# Define boundary condition

def boundary0(x, on_boundary):
return on_boundary and near(x[0], 0)
def boundaryL(x, on_boundary):
return on_boundary and near(x[0], lox)

bc_xi1 = DirichletBC(V.sub(0), Constant(1.0), boundary0)

bc_xi2 = DirichletBC(V.sub(0), Constant(0.0), boundaryL)

bc_c2 = DirichletBC(V.sub(1), Constant(0.0), boundaryL)

bc_phi1 = DirichletBC(V.sub(2), Constant(phie), boundary0)

bc_phi2 = DirichletBC(V.sub(2), Constant(0.0), boundaryL)

bcs = [bc_xi1, bc_xi2, bc_c2, bc_phi1, bc_phi2 ] # Condiciones dirichlet

# SwitchingFunctionMaterial
def h(_x):
return _x
def dh(_x):
return 1

# BarrierFunctionMaterial
def g(_x):
return W*_x**2.0*(Constant(1.0) - _x)**2
def dg(_x):
return W*Constant(2.0)*(_x * (Constant(1.0) - _x) ** 2 - (Constant(1.0) - _x) * _x ** 2)

# Concentrarion
def cl(_x):
return exp((_x-el)/A)/(Constant(1.0)+exp((_x-el)/A))
def dcldw(_x):
return exp((_x+el)/A)/(A*(exp(_x/A)+exp(el/A))**2)
def cs(_x):
return exp((_x-es)/A)/(Constant(1.0)+exp((_x-es)/A))
def dcsdw(_x):
return exp((_x+es)/A)/(A*(exp(_x/A)+exp(es/A))**2)

def chi(_xi,_w):
return dcldw(_w)*(Constant(1.0)-h(_xi))+dcsdw(_w)*h(_xi)*dv

# Mobility defined by D*c/(R*T), whereR*T is normalized by the chemical potential
# M0*(1-h) is the effective diffusion coefficient; cl*(1-h) is the ion concentration
def D(_xi,_w):
return M0*(Constant(1.0)-h(_xi))*cl(_w)*(Constant(1.0)-h(_xi))

def ft(_w):
return cs(_w)*dv-cl(_w)

# Effective conductivity, DerivativeParsedMaterial
def Le1(_xi):
return S1*h(_xi)+S2*(Constant(1.0)-h(_xi))

#Numerical Parameters
Tf= 108.0
dt= 0.02
num_steps=int(Tf/dt)
k = Constant(dt)

#Define variational (dimensional) problem

F=F1+F2+F3

#Jacobian
J = derivative(F,u)

t=0.0
for n in range(num_steps):

t+=dt
#print('mpi rank', mpi_rank, 'time', t)

solve(F == 0, u, bcs, J=J, solver_parameters={"newton_solver":{"absolute_tolerance":1.0e-6, "relative_tolerance":1.0e-6, "maximum_iterations":100}})

u_n.assign(u)

``````

the message error is the following:

File “/usr/lib/python3/dist-packages/dolfin/fem/solving.py”, line 266, in _solve_varproblem
solver.solve()

RuntimeError:

*** Error: Unable to solve nonlinear system with NewtonSolver.
*** Reason: Newton solver did not converge because maximum number of iterations reached.
*** Where: This error was encountered inside NewtonSolver.cpp.
*** Process: 0

*** DOLFIN version: 2019.1.0
*** Git changeset: unknown

I have tried to increase the mesh but it doesn´t work
Anyone who can help me?
Thanks

Please format your post with markdown syntax, such that the code is properly formatted, using ``` for encapsulation. You can format Math with latex syntax using \$ encapsulation

Additionally, please reduce the code example, by removing any variable that is not needed to reproduce the error. (Follow the guidelines: Read before posting: How do I get my question answered?)

Thank you for your adivce, It’s the first time that I use it.
I’ve removed the lines that are not necessary in the code.
I hope it’s clearer now