### Why does the Solution fail to converge? (PETSc reason DIVERGED_ITS)

667
views
0
13 months ago by
Hello everyone,

I am using Fenics installed via conda-forge on MacOS (Sierra) and I run my code via Terminal.
I am trying to solve a problem (linear elasticity coupled with damage) where I use the following geometry (done with Gmsh; file name is ct_mm.xml for use in the code):

It is supposed to be a CT-Specimen, where the force is applied at the hole in the upper right corner.

Here is the code:
import dolfin as dlfn
import numpy as np
#===============================================================================
def get_Y(disp):
eps = dlfn.as_tensor(1.0/2.0*(disp[i].dx(j)+disp[j].dx(i)),(i,j))
eps_dev = dlfn.as_tensor(eps[i,j]-eps[k,k]/3.*delta[i,j], (i,j))
Y = 1./wf*(heavisideFunction(sigmaH)*K/2.*(eps[k,k])**2.0
+ mu*eps_dev[i,j]*eps_dev[j,i])
F = dlfn.FunctionSpace(mesh, "DG", 0)
Yp = dlfn.project(Y, F)
Yp.rename("Y","damage driving force")
return Yp
#===============================================================================
# auxialiary objects
delta = dlfn.Identity(3)
i,j,k = dlfn.indices(3)
def heavisideFunction(x):
return dlfn.conditional(dlfn.gt(x, 0.0), 1.0, 0.0)
def macaulayFunction(x):
return dlfn.conditional(dlfn.gt(x, 0.0), x, 0.0)
#===============================================================================
# solver parameters
abs_tol = 1e-16
rel_tol = 1e-15
max_iter = 1000
#===============================================================================
# geometry of CT-specimen
ct_H = 24.0# mm
ct_G = 50.0# mm
ct_B = 10.0# mm
clength=ct_G # characteristic length of continuum body
cV=ct_G*ct_B*ct_H # characteristic volume
cA=ct_G*ct_B # characteristic area
mesh=dlfn.Mesh('ct_mm.xml')
#===============================================================================
# material parameters for mineral fiber
nu = 0.23
E = 80.0e3 #MPa
# yield strength
sigmaf = 2000.0 #MPa
mu = E/(2.0*(1.0+nu))
# energy at failure
wf = sigmaf**2/2.0/E
# critical Damage
Dc = .99
# damage evolution parameter
b = 0.02 * wf # dim 1
# initial damage resistance
Yinit_value = 25.0 / wf
# bulk modulus
K = 2.*mu*(1.+nu)/3./(1.-2.*nu)
#===============================================================================
# element surfaces
cells = dlfn.MeshFunction('size_t',mesh, 'ct_mm_physical_region.xml')
facets = dlfn.MeshFunction('size_t',mesh, 'ct_mm_facet_region.xml')
dA = dlfn.Measure('ds', domain=mesh, subdomain_data=facets)
dV = dlfn.Measure('dx', domain=mesh, subdomain_data=cells)
#===============================================================================
# info from ct_mm.msh:
# $PhysicalNames # 5 # 2 1 "test_x" # 2 2 "force" # 2 3 "symmy" # 2 4 "symmz" # 3 5 "volume" #$EndPhysicalNames
#===============================================================================
# traction vector
tr = dlfn.Expression(("0.0","A*time","0.0"), A=40.0, time=0.0, degree=1)# in Mpa
#===============================================================================
# Define elements
P1 = dlfn.FiniteElement("DG", mesh.ufl_cell(), 0)# D
P2 = dlfn.VectorElement("CG", mesh.ufl_cell(), 1)# u
TH = P2 * P1
# mixed (product) space
W = dlfn.FunctionSpace(mesh, TH)
# define spaces for auxiliary fields
ScalarSpace = dlfn.FunctionSpace(mesh, "DG", 0)# D0, Yinit, sigmaH
TensorSpace = dlfn.TensorFunctionSpace(mesh, "DG", 0)# sigma0
# define auxiliary fields
D0 = dlfn.Function(ScalarSpace)
sigma0 = dlfn.Function(TensorSpace)
Yinit = dlfn.project(Yinit_value, ScalarSpace)
Y0 = dlfn.Function(ScalarSpace)
#===============================================================================
# boundary conditions
bc1 = dlfn.DirichletBC(W.sub(0).sub(0), 0.0, facets, 1)
bc2 = dlfn.DirichletBC(W.sub(0).sub(1), 0.0, facets, 3)
bc3 = dlfn.DirichletBC(W.sub(0).sub(2), 0.0, facets, 4)
#===============================================================================
print 'initializing time'
t = 0.0
t_end = 1.0
dt = 0.1
#===============================================================================
while t <= t_end:
print 'time: ', t
tr.time = t
bc = [ bc1, bc2, bc3]
#===========================================================================
# definitions for the variational formulation
w = dlfn.Function(W)
(u,D) = dlfn.split(w)
del_w = dlfn.TestFunction(W)
(del_u, del_D) =  dlfn.split(del_w)
dw = dlfn.TrialFunction(W)
(du,dD) = dlfn.split(dw)
#===========================================================================
# strain tensor
eps = dlfn.as_tensor(1.0/2.0*(u[i].dx(j)+u[j].dx(i)),(i,j))
# strain deviator
eps_dev = dlfn.as_tensor(eps[i,j]-eps[k,k]/3.*delta[i,j], (i,j))
# Cauchy stress tensor
sigmaH = dlfn.project(sigma0[k,k]/3., ScalarSpace)
sigma = dlfn.as_tensor(1./wf*((1.0-
heavisideFunction(sigmaH)*D)*K*eps[k,k]*delta[i,j]+2.0*mu*(1.0-D)*eps_dev[i,j]),(i,j))
# damage driving force
Y = 1./wf*(heavisideFunction(sigmaH)*K/2.*(eps[k,k])**2.0
+ mu*eps_dev[i,j]*eps_dev[j,i])
#===========================================================================
# Form
dw_int = 1./cV*(sigma[j,i]*del_u[i].dx(j)+(D-Dc*(1.-dlfn.exp(-
b*macaulayFunction(Y-Yinit))))*del_D)*dV(5)
dw_ext = 1./(cA*clength*wf)*tr[i]*del_u[i]*dA(2)
Form = dw_int - dw_ext
#===========================================================================
# definition of local tangent (2nd variation)
tangent = dlfn.derivative(Form,w)
#===========================================================================
problem = dlfn.NonlinearVariationalProblem(Form, w, bcs = bc, J = tangent)
solver = dlfn.NonlinearVariationalSolver(problem)
prm = solver.parameters
prm["newton_solver"]["absolute_tolerance"] = abs_tol
prm["newton_solver"]["relative_tolerance"] = rel_tol
prm["newton_solver"]["maximum_iterations"] = 25
prm["newton_solver"]["relaxation_parameter"] = 1.0
prm["newton_solver"]["linear_solver"] = "gmres"
prm["newton_solver"]["preconditioner"] = "ilu"#"hypre_amg"
prm["newton_solver"]["krylov_solver"]["absolute_tolerance"] = abs_tol * 1e-3
prm["newton_solver"]["krylov_solver"]["relative_tolerance"] = rel_tol
prm["newton_solver"]["krylov_solver"]["maximum_iterations"] = max_iter
solver.solve()
#===========================================================================
(u, D) = w.split(True)
#
Y = get_Y(u)
print "Y: ", np.amax(Y.vector().array())
dmgc=Y.vector().array() - Yinit.vector().array()
dmgc[dmgc<0]=0 # set all values in diff which are <0 to zero
differenz=Y.vector().array()-Y0.vector().array()
differenz[differenz<0]=0
if dmgc.any()==True and differenz.any()==True:
print "damage occured"
for p in range(np.size(D.vector().array())):
if D.vector().array()[p] > D0.vector().array()[p]:
D0 = D
print "damage updated"
else:
print "no damage increase, Y < Y0"
Y0 = Y
print "Y0 updated"
sigma0 = dlfn.project(sigma, TensorSpace)
print "stress updated"
t=t+dt
​

I guess the critical point is the applied force (see traction vector in the code).
For very small values (e.g. A=1.0 in the above expression) the code even works, so it might be a scaling problem?
On the other hand the value for Y (see line print "Y: ", np.amax(Y.vector().array())) is far too small and for A=40.0
I get the following Error (as stated in the question title):
time:  0.7
Solving nonlinear variational problem.
Newton iteration 0: r (abs) = 4.988e-04 (tol = 1.000e-16) r (rel) = 1.000e+00 (tol = 1.000e-15)
Traceback (most recent call last):
File "mix_ct_mwe.py", line 150, in <module>
solver.solve()
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
***
***
*** Remember to include the error message listed below and, if possible,
*** include a *minimal* running example to reproduce the error.
***
*** -------------------------------------------------------------------------
*** Error:   Unable to solve linear system using PETSc Krylov solver.
*** Reason:  Solution failed to converge in 1000 iterations (PETSc reason DIVERGED_ITS, residual norm ||r|| = 1.592554e-17).
*** Where:   This error was encountered inside PETScKrylovSolver.cpp.
*** Process: 0
***
*** DOLFIN version: 2017.1.0
*** Git changeset:
*** -------------------------------------------------------------------------

Abort trap: 6​

What can I try to find out where the problem is and how to solve it?

Thanks,
Philipp
Community: FEniCS Project
Hi, If you do not mind, can you please send me the mesh files (ct_mm.xml) and  (ct_mm_physical_region.xml).

Thanks for help
written 11 months ago by hirshikesh

5
13 months ago by
Reason:  Solution failed to converge in 1000 iterations (PETSc reason DIVERGED_ITS, residual norm ||r|| = 1.592554e-17)

It says that the solver failed to converge, but the residual norm is quite small (10-17). May be the tolerance you chose is two small. If you really nead that kind of tolerance, you could always increse the number of max iterations and see if it converges.

Cheers.

Hi Fabien,