Hello,
I am beginner with FEniCSx and am currently working with version 0.8.0.
I am implementing a phase-field brittle fracture problem using a monolithic approach and attempting to integrate a damped Newton method. As part of this process, firstly I am trying to develop a custom Newton solver that completely replicates the results of dolfinx.nls.petsc.NewtonSolver.
Using Dokken’s tutorial as a guide, I have created a custom Newton solver. In my MWE, the commented sections use dolfinx.nls.petsc.NewtonSolver, while the remaining code is for the custom Newton solver.
However, I am encountering the following issues:
- The residuals from the custom solver do not match those from
dolfinx.nls.petsc.NewtonSolver. Additionally, the custom solver requires more iterations and is significantly slower. - The custom solver does not support parallel execution. When I attempt to run it in parallel, the code hangs and does not proceed.
Here is the code snippet:
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
MPI.COMM_WORLD.Barrier()
start_time = MPI.Wtime()
###############################################################################
# loading a 2D mesh - Single Edge Notched shear (SENP) specimen discretized with quadrilateral elements
###############################################################################
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!")
###############################################################################
# creating mixed element function space, test functions and solution vectors
###############################################################################
# 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)
#botom line - id - 1
#top line - id - 2
#crack upper line - id - 31
#crack bot line - id -1- 32
#surface - id - 4
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)
###############################################################################
# defining the variational form (monolithic)
###############################################################################
# mechanics field variational form
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
'''
# problem creation
problem = dolfinx.fem.petsc.NonlinearProblem(F,sol,[bc_ubot, bc_utop, bc_dcrack, bc_utopy])
###############################################################################
# Initialize solver and set parameters
###############################################################################
# solver creation - nonlinear solver - newton solver
solver = dolfinx.nls.petsc.NewtonSolver(MPI.COMM_WORLD,problem)
# newton solver parameters
# convergence criterion - incremental tolerance (or residual tolerance)
solver.convergence_criterion = "incremental"
# set relative tolerance
solver.rtol = 1e-6
# set absolute tolerance
solver.atol = 1e-6
# set max iterations
solver.max_it = 30
# krylov subspace solver parameters
ksp = solver.krylov_solver
opts = PETSc.Options()
option_prefix = ksp.getOptionsPrefix()
# using direct solver with LU preconditioner
opts[f"{option_prefix}ksp_type"] = "preonly"
opts[f"{option_prefix}pc_type"] = "lu"
opts[f"{option_prefix}pc_factor_mat_solver_type"] = "mumps"
ksp.setFromOptions()
# 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)
log_info("PETSc problem created and solver parameters set.",my_logger)
'''
###############################################################################
# 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.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)
'''# exception handling for situations when newton solver does not converge
try:
r = solver.solve(sol) # r = [convergence iterations, converge status(true or false)]
assert r[1]
except:
if(MPI.COMM_WORLD.rank ==0):
my_logger.error("Newton Solver did not converge!.")
print("Newton Solver did not converge!.")
break'''
#solver:
i = 0
max_iterations=30
L2_error = []
L_norm = []
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])
# 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
# Compute norm of update
correction_norm = np.linalg.norm(L.array)
L_norm.append(correction_norm)
print(f"Iteration {i}: Correction norm {correction_norm}")
if correction_norm < 1e-6:
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
#
# end stopwatch
MPI.COMM_WORLD.Barrier()
end_time = MPI.Wtime()
print(f"Elapsed time: {end_time-start_time} s.")
I would be grateful if anyone can point out where i am making mistake or the way whereby I can solve the problems.
(If the problem can be resolved this way, I would prefer not to use the SNES solver, as my damped Newton method is specifically customized for my problem. However, if using SNES provides a significantly more convenient solution, I am open to exploring that approach—please guide me on how to implement it.)