Convergence problem with changing initial condition

Question

Hi
I have a code for a time-dependent problem.There is a part in my code that corresponds to the initial condition:

# Previous solution
W_const = Constant(1.0)
C_previous= interpolate(W_const, V_c)

It works perfectly but only for:

W_const = Constant(1.0)

The problem is when I change the the value of the initial condition (W_const) for example:

W_const = Constant(1200.0)

I get an error indicating it does not converge with NewtonSolver. Here is the error I am getting:

*** 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: 1.6.0
*** Git changeset:  unknown

Here is my complete code:

from dolfin import *
import mshr

#Diffusion constant
D = 10 * (10**(-11))

#The gas constant
R = 8.31

#Electrical mobility
mu = 4.11 * (10**(-14.))

#Anion concentration
C_0 = 1200

#Dielectric permittivity
epsilon = 0.025

#Temperature
T = D / (mu*R)

#Faraday number
F = 96485

#length scale
l = 200 * (10**(-6))

#Debye screening length
lambda_d = ((epsilon*R*T)/(2*pow(F,2)*C_0))**(0.5)

eps = lambda_d / l

# Define boundaries
class Left(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[0], 0.0)


class Right(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[0], 1.0)


class Bottom(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[1], 0.0)


class Top(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[1], 1.0)


# Initialize sub-domain instances
left = Left()
top = Top()
right = Right()
bottom = Bottom()

# Domain
domain = mshr.Rectangle(Point(0,0), Point(1,1))

# Mesh
mesh = mshr.generate_mesh(domain, 60)

# Marking boundary domains
boundaries = FacetFunction("size_t", mesh)
boundaries.set_all(0)
left.mark(boundaries, 2)
top.mark(boundaries, 1)
right.mark(boundaries, 4)
bottom.mark(boundaries, 3)

#Defining the function spaces for the concentration
V_c = FunctionSpace(mesh, 'CG', 1)

#Defining the function spaces for the voltage
V_phi = FunctionSpace(mesh, 'CG', 1)

#Defining the mixed function space
W = MixedFunctionSpace((V_c, V_phi))

#Defining the "Trial" functions
u = Function(W)
c, phi = split(u)

#Defining the test  functions
(v_1, v_2) = TestFunctions(W)

# Time variables
dt = 0.05; t = 0; T = 3

#Zero Neumann boundary conditions
g_L=Constant(0.0)
g_R=Constant(0.0)

# Previous solution (equal to C_0)
W_const = Constant(1200.0)
C_previous= interpolate(W_const, V_c)

# Define Dirichlet boundary conditions at top boundary (V=4.0 volt)
Voltage = Constant(4.0)
bc_top = DirichletBC(W.sub(1), Voltage, boundaries, 1)

# Define Dirichlet boundary conditions at bottom boundary (V=0.0 volt)
bc_bottom = DirichletBC(W.sub(1), 0.0, boundaries, 3)

#Collecting the boundary conditions
bcs = [bc_top, bc_bottom]

# Define variational form
a = dot(grad(phi), grad(v_2)) * dx \
    + c*v_1 * dx \
    - (1./(2 * (pow(eps,2)))) * c * v_2 * dx \
    + dt * eps * dot(grad(c),grad(v_1)) * dx \
    + dt * eps * c * dot(grad(phi),grad(v_1)) * dx

L = C_previous * v_1 * dx - (1./(2 * (pow(eps,2)))) * v_2 * dx + eps*g_L*v_1*ds(2) + eps*g_L*v_1*ds(4)

# For nonlinear solver
F = a - L
f = File("conscen3.pvd")

while t <= T:
    solve(F == 0, u, bcs)

    (c, phi) = u.split(True)
    C_previous.assign(c)
    #print phi(0.8,0.2)
    t += dt
    plot(c,key="c", interactive=False)
    f<<c

Could you please help me figure out this issue?
Thanks in advance!

Answer 1

I solved this convergence problem by changing the dimensions of the domain. When we increase the initial condition from 1 to 1200 we have to change the dimensions of the domain for it to converge. I did this:

domain = mshr.Rectangle(Point(0,0), Point(10**(-3),10**(-3)))

Now it works right.

Answer 2

Hi there,

completely ignoring the physics of your problem for the moment, have you tried choosing a different solver and/or increasing the number of iterations steps?

See e.g. here for how this can be done.

Does this help?

Comments

Thank you for your response. I had never tried changing the solver. I gave it a try according to the link you provided. In the first step I just tried to change the tolerance and number of time steps. My code would fail after 50 iterations and by increasing the number of iterations to 100 it still fails. Maybe I am doing somwthing wrong here.
This is what I did:

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

I think I have to change other parameteres like: solver type (a nonlinear solver except newton)

Could you please suggest me something helpul for my code as I am not that familiar with changing the solver parameteres.
Thanks again!

---

Hi Leo,

can you, during each timestep, monitor the progress? E.g. via the set_log_level statement:

#CRITICAL  = 50, // errors that may lead to data corruption and suchlike
#ERROR     = 40, // things that go boom
#WARNING   = 30, // things that may go boom later
#INFO      = 20, // information of general interest
#PROGRESS  = 16, // what's happening (broadly)
#TRACE     = 13, // what's happening (in detail)
#DBG       = 10  // sundry

set_log_level(INFO)

Also, what happens if you do 1000 iterations of the Newton solver?

If none of that helps, then your variational form is likely not correct.

---

The variational formulation looks incorrect. Have you checked it? E.g: why is there no \Delta t on the boundary integral terms?

---

I changed the the number of iterations to 1000 but it failed to converge again. This is really weired because it works fine when I put:

# Previous solution
W_const = Constant(1.0)

This is the message I get:

  Newton iteration 1: r (abs) = 3.337e-09 (tol = 1.000e-08) r (rel) = 1.689e-03 (tol = 1.000e-08)

  Newton solver finished in 1 iterations and 1 linear solver iterations.

It fails to converge when I increase it to :

# Previous solution
W_const = Constant(1000.0)

This is the message I get at the last iteration:

 Newton iteration 50: r (abs) = 4.096e+08 (tol = 1.000e-08) r (rel) = 3.481e+06 (tol = 1.000e-08)

Do you know of any idea to fix this convergence issue? Is it good to change the solver from Newton solver to another one? Because increasing the number of time steps did not work.

---

Hi nate
Do you mean theses terms?

eps*g_L*v_1*ds(2) + eps*g_L*v_1*ds(4)

Actually I should not have included these 2 terms in the variational problem because they are equal to zero. It does not matter if delta-t is multiplied in them or not.