Solution nonconvergence

General
1

Hi, I’m currently solving the diffusion problem of a plate, but when I apply a boundary on both sides, I always get a non-convergent result，I hope you can give me some suggestions, here is an example of my code：

from fenics import *
import numpy as np
import matplotlib.pyplot as plt

# 定义已知表达式
kB = Constant(8.6e-5)    
r = Constant(0.3)  
Sext = Constant(0.0)  
fai = Constant(1e24)   
f = 0.002
u_0 = Expression('0', degree=1)

t_total = 10 
num_steps = 500   
dt = t_total / num_steps   

def D(T):    
    return 4.17e-7*exp(-0.39/(kB*T))

def KrW(T):   
    return 3.2e-15*exp(-1.16/(kB*T))

def KrCuCrZr(T):   
    return 2.9e-14*exp(-1.92/(kB*T))

def S(T):   
    return 9e-3*exp(-1.04/(kB*T))

mesh = Mesh('10_10_100s_slab/mesh.xml.gz')

W = FunctionSpace(mesh, 'P', 1)
V = FunctionSpace(mesh, 'P', 1)

class BoundaryX0(SubDomain):    
    def inside(self, x, on_boundary):
        return on_boundary and abs(x[0] - 0.00) < DOLFIN_EPS

class BoundaryX1(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and abs(x[0] - 0.01) < DOLFIN_EPS
    
class BoundaryY0(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and abs(x[1] - 0.00) < DOLFIN_EPS

class BoundaryY1(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and abs(x[1] - 0.01) < DOLFIN_EPS

bx0 = BoundaryX0()    
bx1 = BoundaryX1()
by0 = BoundaryY0()
by1 = BoundaryY1()

boundary_markers = MeshFunction('size_t', mesh, mesh.topology().dim()-1)    # 创建网格函数

#boundary_markers.set_all(9999)  
bx0.mark(boundary_markers, 0)    
bx1.mark(boundary_markers, 1)
by0.mark(boundary_markers, 2)
by1.mark(boundary_markers, 3)

ds = Measure('ds', domain=mesh, subdomain_data=boundary_markers)
        
T = Function(W)
c1 = Function(V)
c_n = project(u_0, V)  
v = TestFunction(V)

K1 = KrW(T)*S(T)**2

F = c1*v*dx + dt*D(T)*dot(grad(c1),grad(v))*dx - c_n*v*dx - dt*Sext*v*dx - dt*((1-r)*fai-KrW(T)*c1**2)*v*ds(3) + dt*(KrCuCrZr(T)*c1**2)*v*ds(2) - dt*(KrW(T)*c1**2-K1*f)*v*ds(0)

timeseries_T = TimeSeries('10_10_100s_slab/temperature_series')

progress = Progress('Time-stepping', num_steps)
t = 0
bcs=[]
J = derivative(F, c1)
problem = NonlinearVariationalProblem(F, c1, bcs, J)
solver = NonlinearVariationalSolver(problem)
prm = solver.parameters
prm['nonlinear_solver'] = 'newton'
prm['newton_solver']['absolute_tolerance'] = 1e-8
prm['newton_solver']['relative_tolerance'] = 1e-5
prm['newton_solver']['maximum_iterations'] = 100
prm['newton_solver']['relaxation_parameter'] = 1.0

plt.figure(figsize=(12,8))
for n in range(num_steps):
   
    t += dt
    timeseries_T.retrieve(T.vector(), t)
    solver.solve()
    plot(c1)
    c_n.assign(c1)
    
    set_log_level(LogLevel.PROGRESS)
    progress += 1
    #set_log_level(LogLevel.ERROR)

Be grateful！

2

Hi, To reproduce the error I need 10_10_100s_slab/mesh.xml.gz and 10_10_100s_slab/temperature_series files.

3

Okay, I’m terribly sorry.

10_10_100s_slab.rar - Google Drive

4

I sent you an access request.

5

Sorry, I have agreed to share, you can try to download

6

Yeah, now I’m able to reproduce the error as -

RuntimeError: 

*** -------------------------------------------------------------------------
*** DOLFIN encountered an error. If you are not able to resolve this issue
*** using the information listed below, you can ask for help at
***
***     fenics-support@googlegroups.com
***
*** Remember to include the error message listed below and, if possible,
*** include a *minimal* running example to reproduce the error.
***
*** -------------------------------------------------------------------------
*** 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.2.0.dev0
*** Git changeset:  ubuntu
*** -------------------------------------------------------------------------

I tried couple of things like increasing maximum iterations and relax the tolerance requirements as done in this post How to increase maximum number of iterations? and I also tried easing the relaxing parameter as suggested here Default absolute tolerance and relative tolerance and also used CustomSolver() given here.
But none of them worked which led me to conclude that there might be a chance of malformed Jacobian therefore I suggest you to kindly read this post Default absolute tolerance and relative tolerance and check Jacobian once.

1 Like
7

Thanks for your answer, I tried to solve the following Jacobian matrix, using the following code:

bcs=[]
J = derivative(F, c1)
J_mat = assemble(J)
J_array = J_mat.array()
print(J_array)
plt.spy(J_array)
plt.show()

Here are the results：
1689656621105
I don’t know how to analyze this result, can you continue to give me some help?

8

I don’t have time to read your code in detail. As @violetus states, rather than looking at numerical values of the matrix, you could consider slowly building up the ‘nonlinearity’ of your problem to see its effects on the Jacobian. Your material coefficients are very difficult to deal with, so a first step would be to set everything to a reasonable constant value. Then, step by step, introduce more complex models until you find the source of the problem. Getting highly nonlinear problems to converge can be quite challenging.

1 Like