Update:
With the help of this post, the problem has been partially resolved.
However, there are still some discrepancies:
The residual values do not match exactly, but the final calculated force is correct. So it is not a concern anymore.
When using iterative convergence based on the dx-norm, the solution converges successfully. However, when using residual-based convergence with the rhs-norm, the solution fails to converge.
Could you please help identify any potential issues in the implementation? Below is the updated code snippet for reference:
import os
import dolfinx.fem.petsc
import dolfinx.nls.petsc
from dolfinx import fem, io, log, mesh
import basix
import ufl
import numpy as np
from petsc4py.PETSc import ScalarType
from petsc4py import PETSc
from mpi4py import MPI
from myfunctionstrial import mtrl_parameters, fracture_parameters, degradation_fn, \
bracket_operator_plus, bracket_operator_minus, log_info, log_setup
meshfilename = "senp.msh"
if meshfilename.endswith(".msh"):
# dmesh -> DOLFINxmesh, ct -> cell_tags, ft -> facet_tags
dmesh, ct, ft = io.gmshio.read_from_msh(meshfilename, MPI.COMM_WORLD, 0, gdim=2)
else:
raise ValueError("Only .msh or .xdmf file formats for meshes are supported!")
# function elements P->mechanics field (vector), S->fracture field (scalar), R-> history field (scalars)
P = basix.ufl.element("CG", dmesh.basix_cell(), 1, shape=(dmesh.geometry.dim,))
S = basix.ufl.element("CG", dmesh.basix_cell(), 1)
R = basix.ufl.element("DG", dmesh.basix_cell(), 0)
# mixed element
mix = basix.ufl.mixed_element([P,S])
# create a function space consisting of a mixed element at each node
V = fem.functionspace(dmesh,mix)
# create a separate function space for history variable
W = fem.functionspace(dmesh,R)
# test functions du->mechanics field, dd->fracture field
du , dd = ufl.TestFunctions(V)
# mixed solution vector creation, sol-> n+1 th load increment step, solold-> n th time step (previous)
sol = fem.Function(V)
solold = fem.Function(V)
# history variable vector creation
H = fem.Function(W)
Hold = fem.Function(W) # history field at the previous time step
# individual solution vectors at n+1 step, u-> mechanics field, d-> fracture field
u, d = ufl.split(sol)
# individual solution vectors at n th step, uold-> mechanics field, dold-> fracture field
uold, dold = ufl.split(solold)
lmda, mu, K = mtrl_parameters()
G_cr, l, total_u, delta_u, k =fracture_parameters()
# function for returning facets in the top edge
def top(x):
return np.isclose(x[1],1)
# function for returning facets in the bottom edge
def bot(x):
return np.isclose(x[1],0)
# function for returning facets in the crack face
def crack(x):
#crack is at (0,0.5) to (0.5,0.5)
return np.logical_and(abs(x[1]-0.5) <= 0.001, x[0] <= 0.5)
# subspace of mechanics field, accessing all components of u
V0,_ = V.sub(0).collapse()
# subspace of mechanics field, accessing u along x-direction
V01, _ = V.sub(0).sub(0).collapse()
# subspace of mechanics field, accessing u along y-direction
V02, _ = V.sub(0).sub(1).collapse()
# subspace of fracture field, accessing d
V1,_ = V.sub(1).collapse()
# obtain corresponding facets as per the requirement
bot_facets = mesh.locate_entities_boundary(dmesh, 1, bot)
top_facets = mesh.locate_entities_boundary(dmesh, 1, top)
crack_facets = mesh.locate_entities_boundary(dmesh, 1, crack)
# obtain dofs on the corresponding facets using the subspaces
# on bottom face all components of u are chosen
bot_dofs = fem.locate_dofs_topological((V.sub(0),V0), 1, bot_facets)
# on top face only x-component of u is chosen
top_dofs = fem.locate_dofs_topological((V.sub(0).sub(0),V01), 1, top_facets)
# on top face only y-component of u is chosen
top_dofs_y = fem.locate_dofs_topological((V.sub(0).sub(1),V02), 1, top_facets)
# the phase field parameter d dofs is chosen on the crack facets
crack_dofs = fem.locate_dofs_topological((V.sub(1),V1), 1, crack_facets)
# expression for defining values to the corresponding boundary conditions
# boundary values values needed for bottom face
u_bot = fem.Function(V0)
u_botval = fem.Constant(dmesh,ScalarType([0.0,0.0]))
u_botexpr = fem.Expression(u_botval,V0.element.interpolation_points())
u_bot.interpolate(u_botexpr)
# boundary values values needed for top face
u_top = fem.Function(V01)
u_topval = fem.Constant(dmesh,ScalarType(0.001))
u_topexpr = fem.Expression(u_topval,V01.element.interpolation_points())
u_top.interpolate(u_topexpr)
# boundary values values needed for top face
u_topy = fem.Function(V02)
u_topyval = fem.Constant(dmesh,ScalarType(0))
u_topyexpr = fem.Expression(u_topyval,V02.element.interpolation_points())
u_topy.interpolate(u_topyexpr)
# expression for defining the value of d at crack facets
d_crack = fem.Function(V1)
d_crack.x.array[:] = 1
# defining dirichlet boundary condition on the bottom face, all components of u=0 (fixed)
bc_ubot = fem.dirichletbc(u_bot,bot_dofs,V.sub(0))
# defining dirichlet boundary condition on the top face, for applying load increments along x
bc_utop = fem.dirichletbc(u_top, top_dofs, V.sub(0).sub(0))
# defining dirichlet boundary condition on the top face, for applying uy=0
bc_utopy = fem.dirichletbc(u_topy, top_dofs_y, V.sub(0).sub(1))
# defining dirichlet boundary condition on the pre-crack face to be d=1
bc_dcrack = fem.dirichletbc(d_crack, crack_dofs, V.sub(1))
# defining tractions - (Neumann Boundary Conditions) - zero tractions on the boundary
t_bar = fem.Constant(dmesh,np.array([0.0,0.0],dtype=ScalarType))
# defining body forces - zero body forces
b = fem.Constant(dmesh,np.array([0.0,0.0],dtype=ScalarType))
# facet normal (will be required for computation of forces)
n = ufl.FacetNormal(dmesh)
ds_top = ufl.Measure(integral_type='ds', domain=dmesh, subdomain_data=ft, subdomain_id=2)
###############################################################################
# defining constitutive and kinematical relations
###############################################################################
def To2D(A):
return ufl.as_tensor([[ A[0,0], A[0,1]],
[ A[1,0], A[1,1]] ])
def To3D(A):
return ufl.as_tensor([[ A[0,0], A[0,1], 0],
[ A[1,0], A[1,1], 0],
[ 0, 0, 0]])
# infinitesimal strain tensor
def epsilon(x_u):
return ufl.sym(ufl.grad(x_u))
# stress tensor
def sigma(x_u, x_d):
eps = epsilon(x_u)
eps=To3D(eps)
dev_epsilon = eps - ( (1/3)*ufl.tr(eps)*To3D(ufl.Identity(len(x_u))) )
stress3d= (degradation_fn(x_d)+ k) * ( (lmda+mu)*bracket_operator_plus(ufl.tr(epsilon(x_u)))*To3D(ufl.Identity(len(x_u))) + 2*mu*dev_epsilon ) +\
(lmda+mu)* bracket_operator_minus(ufl.tr(epsilon(x_u))) *To3D(ufl.Identity(len(x_u))) # Deviator split
return To2D(stress3d)
#return (degradation_fn(x_d)+ k)*(2.0*mu*epsilon(x_u)+lmda*ufl.tr(epsilon(x_u))*ufl.Identity(len(x_u))) # alternative - for hybrid implementation (ref. thesis)
###############################################################################
# defining fracture related expressions
###############################################################################
# tensile part of elastic strain energy density
def psi_e(x_u):
eps = epsilon(x_u)
eps=To3D(eps)
dev_epsilon = eps - ( (1/3)*ufl.tr(eps)*To3D(ufl.Identity(len(x_u))) )
return 0.5*(lmda+mu)*(bracket_operator_plus(ufl.tr(epsilon(x_u))))**2 + mu*ufl.inner(dev_epsilon,dev_epsilon)
# History field variable
def H(x_uold,x_u,x_Hold):
return ufl.conditional(ufl.operators.gt(psi_e(x_u),x_Hold),psi_e(x_u),x_Hold)
# mechanics field variational form
#mech = ufl.inner(sigma(u, dold),epsilon(du))*ufl.dx - ufl.dot(b,du)*ufl.dx - ufl.dot(t_bar,du)*ds_top
mech = ufl.inner(sigma(u, d),epsilon(du))*ufl.dx - ufl.dot(b,du)*ufl.dx - ufl.dot(t_bar,du)*ds_top
# fracture field variational form
frac = G_cr*l*ufl.dot(ufl.grad(d),ufl.grad(dd)) *ufl.dx + ((G_cr/l) + 2*H(uold,u,Hold))*ufl.inner(d,dd) * ufl.dx - 2*dd*H(uold,u,Hold) *ufl.dx
# combined variational form
F = mech + frac
# uncomment below lines to know about the residual and convergence status at each solution step
# visible only in console (will not be logged in log file)
#if (MPI.COMM_WORLD.rank ==0):
# log.set_log_level(log.LogLevel.INFO)
###############################################################################
# general initializations
###############################################################################
residual = dolfinx.fem.form(F)
dsol = ufl.TrialFunction(V)
J = ufl.derivative(F, sol, dsol)
jacobian = dolfinx.fem.form(J)
# Next, we define the matrix `A`, right hand side vector `L` and the correction function `du`
delsol = fem.Function(V)
delu, deld = ufl.split(delsol)
delu = delsol.sub(0)
deld = delsol.sub(1)
A = dolfinx.fem.petsc.create_matrix(jacobian)
L = dolfinx.fem.petsc.create_vector(residual)
solver = PETSc.KSP().create(dmesh.comm)
solver.setType(PETSc.KSP.Type.PREONLY)
solver.getPC().setType(PETSc.PC.Type.LU)
opts = PETSc.Options()
prefix = f"solver_{id(solver)}"
solver.setOptionsPrefix(prefix)
option_prefix = solver.getOptionsPrefix()
opts[f"{option_prefix}ksp_type"] = "preonly"
opts[f"{option_prefix}pc_type"] = "lu"
opts[f"{option_prefix}pc_factor_mat_solver_type"] = "mumps"
solver.setFromOptions()
solver.setOperators(A)
load = 0
# initializing u, d, uold, dold, aold, vold for problem solving and results saving
u = sol.sub(0)
d = sol.sub(1)
uold = solold.sub(0)
dold = solold.sub(1)
while load <= 1:
# incrementing displacement
u_topval.value = load*total_u
u_top.interpolate(u_topexpr)
#solver:
i = 0
max_iterations=10
converged = False
while i < max_iterations:
# Assemble Jacobian and residual
with L.localForm() as loc_L:
loc_L.set(0)
A.zeroEntries()
dolfinx.fem.petsc.assemble_matrix(A, jacobian, bcs=[bc_ubot, bc_utop, bc_dcrack, bc_utopy])
A.assemble()
dolfinx.fem.petsc.assemble_vector(L, residual)
L.ghostUpdate(addv=PETSc.InsertMode.ADD, mode=PETSc.ScatterMode.REVERSE)
L.scale(-1)
# Compute b - J(u_D-u_(i-1))
dolfinx.fem.petsc.apply_lifting(L, [jacobian], [[bc_ubot, bc_utop, bc_dcrack, bc_utopy]], x0=[sol.x.petsc_vec],scale=1)
# Set du|_bc = u_{i-1}-u_D
dolfinx.fem.petsc.set_bc(L, [bc_ubot, bc_utop, bc_dcrack, bc_utopy], sol.x.petsc_vec, 1.0)
L.ghostUpdate(addv=PETSc.InsertMode.INSERT_VALUES, mode=PETSc.ScatterMode.FORWARD)
# Solve linear problem
solver.solve(L, delsol.x.petsc_vec)
delsol.x.scatter_forward()
# Update u_{i+1} = u_i + delta u_i
sol.x.array[:] += delsol.x.array
i += 1
dolfinx.fem.petsc.assemble_vector(L, residual)
residual_norm=L.norm(1)
correction_norm = delsol.x.petsc_vec.norm(0)
MPI.COMM_WORLD.Barrier()
if(MPI.COMM_WORLD.rank ==0):
print(f"Iteration {i-1}: residual norm {residual_norm}")
print(f"Iteration {i-1}: Correction norm {correction_norm}")
if correction_norm < 1e-6:
converged=True
break
#
if not converged:
if MPI.COMM_WORLD.rank == 0:
print("Nonlinear solver did not converge.")
break
# bounding d value to be between 0 and 1 # not necessary given the problem formulation ?????
interimp = sol.sub(1).collapse()
interimp.x.array [interimp.x.array < 0] = 0
interimp.x.array [interimp.x.array > 1] = 1
sol.sub(1).interpolate(interimp)
sol.x.scatter_forward()
# updating u, d for next iteration
solold.x.array[:] = sol.x.array
solold.x.scatter_forward()
# incrementing load step
load += delta_u
log_info("Simulation completed.",my_logger)