Newton Solver did not converge phase field model

1

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

dw/dt = (1/chi(xi,w))*[ div * [ (D(xi,w)/R*T)*(grad(w) + Fa*grad(phi)) ]
- (dxi/dt)*ft(w) ]

div*(Le1(xi)*grad(phi)) = Fa*Csm*dxi/dt

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



# Constants (no usadas, referencia)
Fa = 96485.3365  # C/mol
R = 8.314472     # J/(molxK)
T = 293.0        # K
R = Constant(R)
T = Constant(T)
Fa = Constant(Fa)

#Dimensional Parameters normalizados
#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

F1 = (xi-xi_n)/k*v_1*dx + L*kappa*dot(grad(xi),grad(v_1))*dx + L*dg(xi)*v_1*dx + Ls*(exp(phi*AA/Constant(2.0))-Constant(14.89)*cl(w)*(Constant(1.0)-h(xi))*exp(-phi*AA/Constant(2.0)))*dh(xi)*v_1*dx

F2 = chi(xi,w)*(w-w_n)/k*v_2*dx + D(xi, w)*dot(grad(w),grad(v_2))*dx + D(xi, w)*AA*dot(grad(phi),grad(v_2))*dx + (h(xi)-h(xi_n))/k*ft(w)*v_2*dx

F3 = Le1(xi)*dot(grad(phi),grad(v_3))*dx + ft2*(xi-xi_n)/k*v_3*dx
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

2

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?)

3

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