Dear Community,
I am working on a Phase Field problem and to compute the stress tensor I am using a spectral decomposition of the strain. I simplyfied the problem in working in 2D so that I can compute by hand the eigenvalue and eigenvector. But I get an NaN value issue at the first timestep of my resolution concerning the mechanical formulation. I tried to avoid this with ‘isnan’ method or a " try / except" but my attempt failed. Do you any idea of how to solve this ?
I joined my code below.
Thanks for your help !
Lm
# Single notched shear test with spectral decomposition of the strain
from fenics import *
import numpy as np
import ufl
from math import isnan
def w(alpha):
return alpha**2
def w_prime(alpha):
return 2*alpha
def gk( alpha, k= 1e-6):
return (1-alpha)**2 + k
def gk_prime(alpha, k= 1e-6):
return 2*(alpha-1)
ppos = lambda x: (abs(x)+x)/2.
pneg = lambda x: x-ppos(x)
# Define strain-rate and stress tensor
def epsilon(u):
return sym(nabla_grad(u))
def eig_val_2D(m):
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
def eig_vec_2D(m):
delta = sqrt(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
#ifNaN = isnan(float("(m[1, 1] - m[0, 0] + delta) / m[0, 1]"))
#ValueError: could not convert string to float: '(m[1, 1] - m[0, 0] + delta) / m[0, 1]'
#if ifNan == True:
# v0y, v1x= 0., 0.
DOLFIN_EPS = 1E-14
w = (m[1, 1] - m[0, 0] + delta )/(m[0, 1]+DOLFIN_EPS)
v0y = 0.5 * w
v1x = - 0.5 * w
"""
try:
v0y = 0.5 * (m[1, 1] - m[0, 0] + delta) / m[0, 1]
v1x = - 0.5 * (m[1, 1] - m[0, 0] + delta) / m[0, 1]
except ufl.log.UFLException as err:
v0y, v1x= 0., 0.
"""
return v0x, v0y,v1x,v1y
def strain_p(m):
m00, m01, m11, m10 = eig_val_order_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_order_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 strain_n(m):
m00, m01, m11, m10 = eig_val_order_2D(m)
w00 = conditional(lt(m00, 0.), m00, 0.)
w11 = conditional(lt(m11, 0.), m11, 0.)
D_neg_eigval = as_tensor([[w00, 0.],[0., w11]])
v1x, v1y,v2x, v2y = eig_vec_order_2D(m)
P_neg_eigval = as_tensor([[v1x, v2x],[v1y, v2y]])
return dot(P_neg_eigval, dot(D_neg_eigval,P_neg_eigval.T ))
# Stress tensor ----------------------------------
def sigma_miehe(u, g, G , lmbda, ndim):
eps = epsilon(u)
eps_plus = strain_p(eps)
eps_moins = strain_n(eps)
sig_plus = lmbda * tr(eps_plus) * Identity(ndim) + 2*G * eps_plus
sig_moins = lmbda * tr(eps_moins) * Identity(ndim) + 2*G * eps_moins
#print(project(g * sig_plus + sig_moins, V_alpha).vector().get_local())
return g * sig_plus + sig_moins
# Positive part of the energy density ------------------
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)
import sys
import time
tic = time.time()
"""
=========================================================
Physical parameters
=========================================================
"""
E, nu = 210000, 0.3 #[MPa]
G = E / (2.0*(1.0 + nu)) # [MPa]
lmbda = E * nu / ((1.0 + nu)* (1.0 - 2.0*nu))
K0, Gc = lmbda+2./3.*G, 2.7
cw, lw = 2, 0.01 # [mm]
"""
=========================================================
Time discretization
=========================================================
"""
T = 0.01 # final time
num_steps = 300 # number of time steps
dt = T / num_steps # time step size
time_step = np.linspace(0, T, num_steps)
"""
=========================================================
Create mesh
=========================================================
"""
nx, ny = 10, 10
#nx, ny = 150, 150 # 160 400 cell h = 0.005 mm
mesh = UnitSquareMesh(nx, ny,'crossed')
Nnode = len(mesh.coordinates())
coord = mesh.coordinates()
ndim = mesh.topology().dim()
"""
=========================================================
Space discretization
=========================================================
"""
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_shear = Expression( ('t', '0.'), t = 0.0, degree = 1 ) # Constant((0.,0.))
# boolean boundary function
boundary_bottom = 'near(x[1], 0)'
boundary_up = 'near(x[1], 1)'
boundary_left = 'near(x[0], 0)'
boundary_right = 'near(x[0], 1)'
class Precrack(SubDomain):
def inside(self, x, on_boundary):
return (between(x[0], (0., 0.5)) and near(x[1], 0.5) ) # and : &&
crack = Precrack()
#Up, bottom, left, right and crack domain
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
crack.mark(boundaries, 5) # mark right as 5
ds = Measure("ds",subdomain_data=boundaries) # left: ds(1), right: ds(2), crack : ds(4), ds(0) all the edges
bcu_bottom = DirichletBC( V , Constant((0.,0.)), boundary_bottom)
bcu_right = DirichletBC( V.sub(1) , Constant(0.), boundary_right)
bcu_left = DirichletBC( V.sub(1) , Constant(0.), boundary_left)
bcu_up = DirichletBC( V , Disp_shear, boundary_up )
bcu = [bcu_bottom, bcu_up, bcu_right, bcu_left]
bcp_bottom = DirichletBC(V_alpha, Constant(0.), boundary_bottom)
bcp_up = DirichletBC(V_alpha, Constant(0.), boundary_up)
bcp_crack = DirichletBC(V_alpha, Constant(1.), boundaries,5)
bc_alpha = [bcp_bottom, bcp_up, bcp_crack]
"""
=========================================================
Define Trial and Test Function
Define functions for solutions at current times steps
=========================================================
"""
u,v,u_ = TrialFunction(V), TestFunction(V), Function(V)
alpha,beta,alpha_ = TrialFunction(V_alpha), TestFunction(V_alpha), Function(V_alpha)
u_old, alpha_old = Function(V), Function(V_alpha)
H_old = interpolate(Expression("0.0",degree=0), V_alpha)
g, g_prime = gk(alpha_,1.0e-6), gk_prime(alpha_,1.0e-6)
sig_funct = Function(Tens, name="Stress")
eps_funct = Function(Tens, name="Strain")
"""
=========================================================
Define Variational problem
=========================================================
"""
#---------------- Mechanical equation
tn = Constant((0., 0.)) # external stress
f_u = inner( sigma_miehe(u_, alpha_, lmbda, G, ndim ) , epsilon(v) ) * dx + dot( tn, v ) * ds
J_u = derivative(f_u, u_, u) # jacobi
## Check value : Nan issue
#eps = epsilon(u_)
#e00, e01, e11, e10 = eig_val_2D(eps)
#print(" e00 : " ,project(e00, V_alpha).vector().get_local())
#v0x, v0y,v1x,v1y = eig_vec_2D(eps)
#print(' v0y : ' ,project(v0y, V_alpha).vector().get_local() )
#print("\n u_.vector().get_local() : ", u_.vector().get_local())
#--------------- Phase Field equation
f_alpha = H(u_old, u_, H_old, lmbda, G)*inner(beta, g_prime) * dx \
+ (Gc/lw* 1/cw) * (beta*w_prime(alpha_) + (2.*lw**2)*inner(grad(beta), grad(alpha_))) * dx
J_alpha = derivative(f_alpha, alpha_, alpha) # jacobian
# Constraints for the phase field
alpha_min = interpolate(Constant(DOLFIN_EPS), V_alpha) # lower bound
alpha_max = interpolate(Constant(1.0), V_alpha) # upper bound
# Solver parametrization
problem_u = NonlinearVariationalProblem(f_u,u_,bcu , J_u)
problem_alpha = NonlinearVariationalProblem(f_alpha,alpha_ ,bc_alpha, J_alpha)
problem_alpha.set_bounds(alpha_min, alpha_max)
solver_u = NonlinearVariationalSolver(problem_u)
solver_alpha = NonlinearVariationalSolver(problem_alpha)
newton_prm = solver_u.parameters['newton_solver']
newton_prm['relative_tolerance'] = 1e-8
newton_prm['absolute_tolerance'] = 1e-8
newton_prm['maximum_iterations'] = 2000
newton_prm['error_on_nonconvergence'] = False
snes_prm = {"nonlinear_solver": "snes",
"snes_solver" : { "method": "vinewtonssls",
"maximum_iterations": 1e-8,
"relative_tolerance": 1e-8,
"absolute_tolerance": 2000,
"report": True,
"error_on_nonconvergence": False,
}}
solver_alpha.parameters.update(snes_prm)
#---------------- RESOLUTION
t= 0.
it = 0
tol = 0.001 # 1e-6
maxiter = 300
stop_tt = 3
import ufl
for i in time_step[:stop_tt]:
Disp_shear.t = t
loading[it]= np.array([it,Disp_shear.t])
print("\n -> TimeLoop : %d "%(it))
print('Time : %.3g s, Ux_given : %.3g mm ' %(t, Disp_shear.t))
k = 1
res = 1.
print("\n ----- > Start Staggered Loop")
while (res > tol) and (k < maxiter):
#----Step 1 : Displacement resolution
iter_u = solver_u.solve()
print("\n u_.vector().get_local() : ", u_.vector().get_local())
## ----> Nan issue
#eps = epsilon(u_)
#e00, e01, e11, e10 = eig_val_2D(eps)
#v1x, v1y,v2x,v2y = eig_vec_2D(eps)
#print(" eps[0, 0] : " ,project( e00, V_alpha).vector().get_local())
## ---------
#----Step 2 : Damage resolution
iter_alpha = solver_alpha.solve()
#----Step 3 : Residual
err_u = errornorm(u_, u_old, norm_type = 'l2', mesh = None)
err_alpha = errornorm(alpha_, alpha_old, norm_type = 'l2', mesh = None)
res = max(err_u, err_alpha)
print( "\n Iteration: %d, Residual: %.8g" % (k, res) )
#---- Step 4 : Update previous solution
H_old.assign(project(H(u_old, u_, H_old, lmbda, G), V_alpha))
u_old.assign(u_)
alpha_old.assign(alpha_)
#g = gk(alpha_,1.0e-6)
k +=1
iter = [k]
t += dt
it += 1
toc = time.time() - tic
print('Elapsed CPU time: ', toc, '[sec]')