Fixing the convergence issue by changing the "absolute_tolerance"

Question

I have a code for solving a time-dependent and non-linear problem. The thing is my code does not converge for some reason. Here is my complete code:

from dolfin import *
import mshr
import numpy as np
from numpy import cos, pi

#Diffusion constant
D = 1. * (10**(-9))

#The gas constant
R = 8.31

#Anion concentration
C_0 = 1200.

#Dielectric permittivity
epsilon = 0.025

#Temperature
Temp = 293.0

#Electrical mobility
mu = D / (R*Temp)

#Faraday number
F_N = 96500.

#Charge number
Charge_num=Constant(1.)


# 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], 0.03)


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


class Top(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[1], (200 * (10**(-6))))

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

# Domain
class Mesh:
    def __init__(self, lx,ly, Nx,Ny):
        self.lx = lx
        self.ly = ly
        self.Nx = Nx
        self.Ny = Ny

    def mesh(self):
        m = RectangleMesh(Point(0,0), Point(1, 1),self.Nx, self.Ny)
        x = m.coordinates()
         #Refine on top and bottom sides
        x[:,:] = (x[:,:] - 0.5) * 2.
        x[:,:] = 0.5 * (cos(pi * (x[:,:] - 1.) / 2.) + 1.)
        #Scale
        x[:,0] = x[:,0]*self.lx
        x[:,1] = x[:,1]*self.ly
        return m

mesh0=Mesh(0.03,200.*(10 ** (-6)), 1000,15)
mesh1 = mesh0.mesh()
# Marking boundary domains
boundaries = FacetFunction("size_t", mesh1)
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(mesh1, 'CG', 1)

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

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

# Defining the "Trial" functions

z = Function(W)
dz=TrialFunction(W)

c, phi = split(z)

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

# Time variables
dt = 0.005
t = 0
T = 0.2

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

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

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

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

# Variational form For nonlinear solver
F = dot(grad(phi), grad(v_2)) * dx \
    + c * v_1 * dx \
    - (F_N / epsilon) * c * v_2 * dx \
    + dt * D * dot(grad(c), grad(v_1)) * dx \
    + dt * Charge_num * mu * F_N * c * dot(grad(phi), grad(v_1)) * dx \
    - C_previous * v_1 * dx + (F_N / epsilon) * C_0 * v_2 * dx

while t <= T:
    #solve(F == 0, u, bcs)
    J = derivative(F, z, dz) #Define the Jacobian
    problem = NonlinearVariationalProblem(F, z, bcs,J)
    solver = NonlinearVariationalSolver(problem)    
    solver.solve()
    (c, phi) = z.split(True)
    C_previous.assign(c)
    t += dt
    plot(c,key="c", interactive=False)

interactive()

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

I fixed this issue by changing the "absolute_tolerance". I just increased this value from 1.000e-10 (default value) to 5e-7. I just added this line:

solver.parameters["newton_solver"]["absolute_tolerance"]=5e-7

When I run this code I am not getting the same error and it gives me some results. This time I got this message:

Solving nonlinear variational problem.
  Newton iteration 0: r (abs) = 4.925e-07 (tol = 5.000e-07) r (rel) = 1.000e+00 (tol = 1.000e-09)
  Newton solver finished in 0 iterations and 0 linear solver iterations.

I just want to know if these results are reliable or not. Is it a right approach to fix the convergence problem by increasing the absolute_tolerance? If not, do you have any idea how to resolve this issue in my code?
Thanks!

Answer 1

There are two issues here:

You seem to be getting a solution with really small values for the units used. You can try using the relative tolerance instead which is a tolerance normalized by the magnitude of your solution. Try:

solver.parameters["newton_solver"]["relative_tolerance"] = 1.0e-3

This might work but also to have a solution less dependent on machine precision errors (when dealing with large and small constants for example), you might consider changing your units.

For example, instead of using all SI units you can use microns for length and change all constants accordingly. After doing this you probably won't need to change the Newton solver settings and your solution will be more accurate.

To check you did this correctly, all the terms in your differential equation should end up being multiplied by the same constant after doing the change of units. For example, if going to meters from microns all the terms should be scaled by 1.0e6 to some power. And remember some units like watts have a lengh scale as well. Get all units in terms of mass length and time units.

Comments

Thanks for your response. Unfortunately using:

solver.parameters["newton_solver"]["relative_tolerance"] = 1.0e-3

was not helpful. Anyway, I have not still gotten the answer of my question. Does it make sense to solve the convergence issue by increasing the absolute tolerance of the newton solver?
In addition, when I am getting this message:

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

can I trust the given results?