Newton Solver issue with an non zero initial guess

lmersel · February 17, 2022, 5:42pm

Hello everyone,
I posted previously on ‘variational formulation’ section an issue about a phase field problem on single edge shear test (see the link ) . But I don’t know if my issue comes from my formulation or the Newton solver I use.

I added few update in the previous code where I forced the initial guess value with the following expression :
u_ = interpolate( Expression( ('0.2', '0.2'), degree = 0 ), V) (u_ the unknow function), avoiding an zero initial value but the Newton solver don’t seem to solve the problem. So I wonder, if finally the problem wouldn’t come from the nonlinear stress formulation I use or I just don’t initialize correctly the problem.

The residual gives an nan value even if I put an non zero inital guess.

Newton iteration 0: r (abs) = -nan (tol = 1.000e-08) r (rel) = -nan (tol = 1.000e-08)

May anyone help me with this?

Thanks for your help !

Lamia

dparsons · February 17, 2022, 10:44pm

By default the solver uses a zero initial guess. To use nonzero initial guess, you need to specify that’s what you want,

prm['newton_solver']['krylov_solver']['nonzero_initial_guess'] = True

lmersel · February 18, 2022, 3:45pm

Thank you very much for your answer ! Indeed I had already tried this command but the issue remains. I think the problem comes from the expression of the stress tensor based on a spectral decomposition.

with

If I initialize with a constant value of u, the initial strain tensor will become zero and leads to an NaN value for some componants of the eigenvectors because it tries to compute initially zero by zero.

def eig_val_2D(E):
    delta =sqrt( E[0,0]**2 +E[1,1]**2 - 2 * E[0,0]* E[1,1] + 4 * E[0,1]+E[1,0])
    E00 = (E[0,0]+E[1,1])/2. + 0.5 * delta
    E01 = 0.0
    E10 = 0.0
    E11 = (E[0,0]+E[1,1])/2. - 0.5 * delta
    return E00, E01, E11, E10

def eig_vec_2D(E):
    delta = sqrt(E[0, 0]**2 + E[1, 1]**2 - 2*E[0, 0]*E[1, 1] + 4*E[0, 1]*E[1, 0] )
    v0x = 1.0 # associated to 0.5*tr(E)+0.5 * delta
    v1y = 1.0 # associated to 0.5*tr(E)-0.5 * delta
    try:
        w = (E[1, 1] - E[0, 0] + delta) / (E[0, 1] )
        v0y = 0.5 * w
        v1x = -0.5 * w
    except ufl.log.UFLException as err:
        v0y, v1x= 0., 0.
    return v0x, v0y,v1x,v1y

An open source code from Tianju Xue seems also to meet the same trouble with the same stress model and the use of the Newton solver. They propose to compute the initial value with a linear model to start with a non zero strain tensor then switch with the model based an the spectral decomposition.
I will try to do this.

Thanks for the help !

Lmersel

nate · February 18, 2022, 3:49pm

Have you considered the following tips and ticks?

You have the line:

w = (E[1, 1] - E[0, 0] + delta) / (E[0, 1] )

If E is some kind of rate of strain tensor then your initial guess needs to be \in C^1 in order to not yield a singularity.

lmersel · February 18, 2022, 4:30pm

Thanks for the link ! It will be really helpful.
My issue is illustrated with this point :

Does your Newton solver fail on the first step producing `NaN` s?

You likely have an initial guess leading to singularities in the Jacobian/residual. E.g: for a solution variable u

Zero initial guess with coefficients of the type 1/u or √u
Piecewise constant initial guess when invoking ∇u−1 or √εII(u)

Indeed I never try an initial displacement being ∈C1. I will try this.

lmersel · February 25, 2022, 4:56pm

Hello !

After some long research, I finally didn’t manage to solve my issue.
I decided to simplify my problem and remove the damage equation. I have just keeped the linear elastic equation and I do a pure traction test applied on the top edge of a square, discretized with 2*2 elements.

But my problem keeps diverging with a Newton solver. At the first time step, I use a linear elastic stress formulation (hooke law in plane strain) to avoid trouble. Then for the following time step, I switch to the spectral decompositon of the stress and I finally got divergence issue.
Fortunately, I managed to highlight the issue. Insteed of having the value 0., I got 10e-18 order value. And when I divided by this, I don’t have an error but the results are huge (with 10e+25 order that is logical). I wanted an error to force the code with a try/except but currently the value belongs to [-1e-18; 1e18].

Just for the information, I also tested a guess value belongs to C1 class but the code diverged too.

@bleyerj Bleyer had already treated this kind of formulation with the package mgis.fenics but I would like to manage this formulation in 2D with fenics first before used mgis package for more complexe cases (3D).
I joined are the results I got when I run the code. I joined the code too to show how I compute the eigen values and vectors.

I will be grateful for any advice !

Thanks a lot for the help.

Lamia

Results :

Display U_ unknow function
__________________________________________________________

 ux.vector().get_local()    :  [ 0.00000000e+00  0.00000000e+00 -1.20492473e-17 -2.25378701e-17
  0.00000000e+00  6.59121014e-18 -1.78618593e-18 -1.40755132e-18
  0.00000000e+00  5.58999480e-18 -2.81510265e-18  0.00000000e+00
  0.00000000e+00]

 uy.vector().get_local()    :  [0.22222222 0.11111111 0.16666667 0.22222222 0.         0.05555556
 0.11111111 0.16666667 0.22222222 0.         0.05555556 0.11111111
 0.        ]

Strain component in cartesian basis
__________________________________________________________

eps00    : [ 8.23315702e-17  1.99907200e-18  3.48658126e-17 -2.79660615e-17
 -4.26654172e-17 -3.50600942e-18 -1.23538316e-18 -3.48658126e-17
 -2.63994473e-17  2.54952951e-17  3.50600942e-18  4.71694314e-19
 -8.32517308e-18]
 eps11    : [0.22222222 0.22222222 0.22222222 0.22222222 0.22222222 0.22222222
 0.22222222 0.22222222 0.22222222 0.22222222 0.22222222 0.22222222
 0.22222222]
 eps01    : [ 3.75936277e-18  2.55270551e-17  2.26150643e-17 -2.40409858e-17
  6.97678205e-19  3.27099423e-17 -1.21789723e-17 -3.94820142e-17
  7.99832289e-18 -7.85679934e-18 -3.00930695e-17 -7.87934707e-18
  7.76037125e-18]

Discriminant (check issue with sqrt()
__________________________________________________________

delta :  [0.04938272 0.04938272 0.04938272 0.04938272 0.04938272 0.04938272
 0.04938272 0.04938272 0.04938272 0.04938272 0.04938272 0.04938272
 0.04938272]
lambda0  :  [0.33333333 0.33333333 0.33333333 0.33333333 0.33333333 0.33333333
 0.33333333 0.33333333 0.33333333 0.33333333 0.33333333 0.33333333
 0.33333333]

Eigen value lambda0 et lamdba1
__________________________________________________________

lambda0 =   [0.22222222 0.22222222 0.22222222 0.22222222 0.22222222 0.22222222
 0.22222222 0.22222222 0.22222222 0.22222222 0.22222222 0.22222222
 0.22222222]
lambda1 =  [ 7.24453566e-17 -1.49764460e-16 -9.16925266e-17  9.88792381e-17
 -5.04721926e-17 -1.92802123e-16  6.76542156e-17  1.85367594e-16
 -7.59148035e-17  7.11236625e-17  1.75454889e-16  5.60893924e-17
 -6.40195569e-17]

Eigen vector components v0 = [v0x, v0y], v1= [v1x, v1y]
__________________________________________________________

v0x =   [1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1.]
v0y =  [ 1.16585344e+17 -1.42508396e+16  6.18573817e+16  5.78184618e+16
  3.43871879e+16  2.76925921e+16 -4.31071142e+15 -2.26450483e+16
 -1.58441840e+16  2.44384602e+15 -1.51142633e+16  5.77919406e+15
 -1.43041773e+16]
v1x =   [-1.16585344e+17  1.42508396e+16 -6.18573817e+16 -5.78184618e+16
 -3.43871879e+16 -2.76925921e+16  4.31071142e+15  2.26450483e+16
  1.58441840e+16 -2.44384602e+15  1.51142633e+16 -5.77919406e+15
  1.43041773e+16]
v1y =  [1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1.]

Implementation :

"""
Usefull function
"""
from fenics import *
import numpy as np
import matplotlib.pyplot as plt
import sys
import time

# Define strain-rate and stress tensor
def epsilon(u):
    return sym(nabla_grad(u))


def eig_val_2D(m):
    #DOLFIN_EPS = 1E-14
    delta = sqrt( m[0,0]**2 + m[1,1]**2 - 2 * m[0,0]* m[1,1] + 4 * m[0,1]+m[1,0]) 
    m00   = (m[0,0]+m[1,1])/2. + 0.5 * delta
    m01   = 0.0
    m10   = 0.0
    m11   = (m[0,0]+m[1,1])/2. - 0.5 * delta
    return m00, m01, m11, m10

import ufl
def eig_vec_2D(m):
    delta = m[0, 0]**2 + m[1, 1]**2 - 2*m[0, 0]*m[1, 1] + 4*m[0, 1]*m[1, 0] 
    v0x = 1.0 # associated to 0.5*tr(m)+0.5 * delta
    v1y = 1.0 # associated to 0.5*tr(m)-0.5 * delta

    #w = (m[1, 1] - m[0, 0] + sqrt(delta)) / (m[0, 1]) #+DOLFIN_EPS )
    #v0y =  0.5 * w
    #v1x = -0.5 * w
    v0y = conditional(ufl.eq(m[0,1], 0.), 0., 0.5 * (m[1, 1] - m[0, 0] + sqrt(delta)) / (m[0, 1]) )
    v1x = conditional(ufl.eq(m[0,1], 0. ), 0., -0.5 * (m[1, 1] - m[0, 0] + sqrt(delta)) / (m[0, 1]) )
    return v0x, v0y,v1x,v1y

def strain_p(m):
    #m = epsilon(u)
    m00, m01, m11, m10 = eig_val_2D(m) 
    w00 = conditional(gt(m00, 0.), m00, 0.)
    w11 = conditional(gt(m11, 0.), m11, 0.)
    D_pos_eigval = as_tensor([[w00, 0.],[0., w11]])
    v1x, v1y,v2x,  v2y = eig_vec_2D(m)
    P_pos_eigval = as_tensor([[v1x, v2x],[v1y, v2y]])
    return dot(P_pos_eigval, dot(D_pos_eigval,P_pos_eigval.T ))


def sigma_miehe(u, G , lmbda, ndim):
    eps             = epsilon(u)
    eps_plus        = strain_p(eps)
    eps_moins       = eps - eps_plus
    sig_plus  = lmbda * tr(eps_plus) * Identity(ndim) + 2*G * eps_plus
    sig_moins = lmbda * tr(eps_moins) * Identity(ndim) + 2*G * eps_moins
    return  sig_plus + sig_moins

def sigma_bourdin(u, G, lmbda,ndim):
    return  2 * G * epsilon(u) + lmbda*tr(epsilon(u))*Identity(ndim)

# Bulk energy ----------------------------------
def psi_plus(u, lmbda, G ):

  eps = epsilon(u)
  eps1, m01, eps2, m10= eig_val_2D(eps)
  tr_eps_plus = ppos(tr(eps))
  eps_plus_doubledot_eps_plus = ppos(eps1)**2 + ppos(eps2)**2
  energy = (1./2.)*lmbda*(tr_eps_plus**2) + G * eps_plus_doubledot_eps_plus
  return energy


# Historical variable -----------------------------------
def H(u_old, u_new,  H_old, lmbda, G):
    psi_i_new = psi_plus(u_new, lmbda, G) 
    psi_i_old = psi_plus(u_old, lmbda, G) 
    psi_new = ppos(psi_i_new) 
    psi_old = ppos(psi_i_old )
    return conditional(lt(psi_old, psi_new), psi_new, H_old)


 print("======================================================================")

ERROR = 40
set_log_level(ERROR) # log level

PROGRESS = 16
#set_log_level(PROGRESS) # log level

parameters.parse()   # read paramaters from command line
# set some dolfin specific parameters
parameters["form_compiler"]["optimize"]     = True
parameters["form_compiler"]["cpp_optimize"] = True
parameters["form_compiler"]["representation"] = "uflacs" # quadrature is deprecated


# ------------------
# folder name : save file Create VTK file for saving solution
# ------------------

savedir        = "results_CubeLinElas_Miehe_nx10/"
vtkfileU       = File(savedir+'displacement.pvd')
vtkfileSigma   = File(savedir+'stress.pvd')



"""
=========================================================
                        Problem
=========================================================
"""

E, nu     = Constant(210000), Constant(0.3) #[MPa]
G         = Constant(E / (2.0*(1.0 + nu))) # [MPa]
lmbda     = Constant (E * nu / ((1.0 + nu)* (1.0 - 2.0*nu)))
K0, Gc    = Constant(lmbda+2./3.*G), Constant(2.7)
#cw, lw    = Constant(2), Constant(0.01) # [mm]
#psi_cr    = 0#Constant(Gc/(2*lw))



T            = 1.           # final time
num_steps    = 10    # number of time steps
#dt           = T / num_steps # time step size
time_step    = np.linspace(0, T, num_steps)

n = 2 # ele :  20201 cellsize h = 0.01 mm l/h = 1
mesh = UnitSquareMesh(n, n,'crossed')
Nnode = len(mesh.coordinates())
coord = mesh.coordinates()
ndim = mesh.topology().dim()

V         = VectorFunctionSpace(mesh,'Lagrange', 1)
V_alpha   = FunctionSpace(mesh,'Lagrange', 1)
V_H       = FunctionSpace(mesh,'DG', 0)
Tens      = TensorFunctionSpace(mesh, 'Lagrange', 1)


"""
=========================================================
                        Boundary condition
=========================================================
"
Disp       = Expression( 't', t = 0.0, degree = 1 )
boundary_bottom  = 'near(x[1], 0)'
boundary_up      = 'near(x[1], 1)'
boundary_left    = 'near(x[0], 0)'
boundary_right   = 'near(x[0], 1)'

bottom     = CompiledSubDomain(boundary_bottom)
up         = CompiledSubDomain(boundary_up)
left       = CompiledSubDomain(boundary_left)
right      = CompiledSubDomain(boundary_right)
boundaries = MeshFunction("size_t", mesh,1)
boundaries.set_all(0)
bottom.mark(boundaries, 1) # mark left as 1
up.mark(boundaries, 2) # mark right as 2
left.mark(boundaries, 3) # mark right as 3
right.mark(boundaries, 4) # mark right as 4
ds = Measure("ds",subdomain_data=boundaries) # left: ds(1), right: ds(2), crack : ds(4), ds(0) all the edges

bcu_bottom = DirichletBC( V.sub(1) , Constant((0.)), boundary_bottom)
bcu_right  = DirichletBC( V.sub(0) , Constant(0.), boundary_right)
bcu_left   = DirichletBC( V.sub(0) , Constant(0.), boundary_left)
bcu_up     = DirichletBC( V.sub(1) , Disp, boundary_up )
bcu        = [bcu_bottom, bcu_up, bcu_right, bcu_left]

du,v,u_        = TrialFunction(V), TestFunction(V), Function(V)
sig_funct      = Function(Tens, name="Stress")


"""
=========================================================
                        Define Variational problem
=========================================================
"""

#u_guess = interpolate(Expression( ('x[1]', 'x[0]*x[0]'),  degree = 2 ),V)
#u_.assign(u_guess)

it = 0
DOLFIN_EPS = 1E-14
for t in time_step[1:]:

    it +=1
    Disp.t = t

    print("\n")
    print("======================================================================")
    print('\n -> Time Step : %d, Time [s] : %.3g, Ux_given [mm] :  %.3g mm ' %(it, t, Disp.t))
    print("======================================================================")
    
    if it < 4:
        sigma = sigma_bourdin(u_  , G, lmbda,ndim)
    else:
        sigma = sigma_miehe(u_ , G , lmbda, ndim)

    #sigma = sigma_miehe(u_ , G , lmbda, ndim)
    f = inner( sigma, epsilon(v) ) * dx + dot( Constant((0., 0.)), v ) * ds
    J = derivative(f,  u_, du)
    
    problem = NonlinearVariationalProblem(f, u_, bcu, J)
    solver  = NonlinearVariationalSolver(problem)
    newton_prm = solver.parameters['newton_solver']
    newton_prm['relative_tolerance'] = 1e-8
    newton_prm['absolute_tolerance'] = 1e-8
    newton_prm['maximum_iterations'] = 100
    newton_prm['krylov_solver']['nonzero_initial_guess'] = True
    
    #---------------- RESOLUTION
    iter_u = solver.solve()

    #---------------- Dispay - Check errors
    ux = u_.sub(0, deepcopy=True)
    uy = u_.sub(1, deepcopy=True)
    print("\n")
    print(" Display U_ unknow function")
    print("__________________________________________________________")
    print("\n ux.vector().get_local()    : ", ux.vector().get_local()[:])
    print("\n uy.vector().get_local()    : ", uy.vector().get_local()[:])

    
    eps = epsilon(u_)
    print("\nStrain component in cartesian basis")
    print("__________________________________________________________")
    print(' \neps00    :',project(epsilon(u_)[0, 0], V_alpha).vector().get_local() )
    print(' eps11    :',project(epsilon(u_)[1, 1], V_alpha).vector().get_local() )
    print(' eps01    :',project(epsilon(u_)[0, 1], V_alpha).vector().get_local() )

    print("\nDiscriminant (check issue with sqrt() ")
    print("__________________________________________________________")
    delta = (eps[0,0]**2 + eps[1,1]**2 - 2 * eps[0,0]* eps[1,1] + 4 * eps[0,1]+eps[1,0])
    print('\ndelta : ', project(delta, V_alpha).vector().get_local())
    #print('lmbPart 1 : ', project((eps[0,0]+eps[1,1])/2. , V_alpha).vector().get_local())
    #print('lmbPart 2 : ', project( sqrt(delta), V_alpha).vector().get_local())
    print('lambda0  : ', project((eps[0,0]+eps[1,1])/2. + sqrt(delta) , V_alpha).vector().get_local())

    print("\nEigen value lambda0 et lamdba1")
    print("__________________________________________________________")
    lambda00, m01, lambda11, m10 = eig_val_2D(eps)
    print('\nlambda0 =  ', project(lambda00, V_alpha).vector().get_local())
    print('lambda1 = ' , project(lambda11, V_alpha).vector().get_local())

    print("\nEigen vector components v0 = [v0x, v0y], v1= [v1x, v1y]")
    print("__________________________________________________________")
    v0x, v0y, v1x, v1y = eig_vec_2D(eps)
    print('\nv0x =  ', project(v0x, V_alpha).vector().get_local())
    print('v0y = ' , project(v0y, V_alpha).vector().get_local())
    print('v1x =  ', project(v1x, V_alpha).vector().get_local())
    print('v1y = ' , project(v1y, V_alpha).vector().get_local())

    print(" \nPositive strain compoenents in cartesian basis")
    print(" eps_plus = P*macaulay(D)*P^T = [ [eps00_p,eps01_p ], [eps11_p],eps10_p ]")
    print("__________________________________________________________")
    eps_p = strain_p(eps)
    print(' \neps00_p    :',project(eps_p[0, 0], V_alpha).vector().get_local() )
    print(' eps11_p    :',project(eps_p[1, 1], V_alpha).vector().get_local() )
    #print(' eps01_p    :',project(eps_p[0, 1], V_alpha).vector().get_local() )

    print(" \nStress components in cartesian basis")
    print(" \nDisplay stress component project in cartesian basis  ")
    print(" sig = sig_plus + sig_minus =  [ [s00, s01 ], [s11], s10 ]")
    print("__________________________________________________________")
    print(' \nsig00    :',project(sigma[0, 0], V_alpha).vector().get_local() )
    print(' sig11    :',project(sigma[1, 1], V_alpha).vector().get_local() )
    #print(' sig01    :',project(sigma[0, 1], V_alpha).vector().get_local() )
    

    #---------------- PostProcess
    sig_funct.assign(project(sigma, Tens))
    elastic_energy         = 0.5*inner(sigma, epsilon(u_))*dx
    elastic_energy_density = project(0.5*inner(sigma, epsilon(u_)), V_alpha)
    elastic_energy_value   = assemble(elastic_energy)

    Fn = assemble(sigma[0,0]*ds(2))
    Ft = assemble(sigma[1,1]*ds(2))
    ux = u_.sub(0, deepcopy=True)
    uy = u_.sub(1, deepcopy=True)

    vtkfileU      << u_
    vtkfileSigma  << sig_funct

    print( "\nEel_density_max: %.8g" % (elastic_energy_density.vector().max()))
    print("\n(Fn, Ft): (%.8g,%.8g)"%(Fn,Ft))
    print("\nElastic energies : (%g)"%(elastic_energy_value)) 
    print("-------------------------------------------------------------")

Topic		Replies	Views
Getting residual "nan" from the start variational formulation	5	621	October 4, 2023
Spectral decomposition - Nan value issues variational formulation	0	731	February 16, 2022
Newtonsolver does not converge - Hyperelasticity	5	1026	May 2, 2020
Newton Solver did not converge phase field model	2	835	October 28, 2020
Newton solver, nonlinear problem convergence	4	974	November 28, 2020

Newton Solver issue with an non zero initial guess

Does your Newton solver fail on the first step producing NaN s?

Related topics

Does your Newton solver fail on the first step producing `NaN` s?