Dear all,
I am currently working on an incompressible hyperelastic problem on a complicated mesh, which at the moment I cannot get to converge. As part of the debugging process I tried the following code on a simple cube mesh, with Dirichlet boundary conditions set to 0 at the bottom plane x = 0:
import dolfinx
import numpy as np
import ufl
from mpi4py import MPI
from petsc4py import PETSc
from dolfinx.fem.petsc import NonlinearProblem
from dolfinx.nls.petsc import NewtonSolver
from basix.ufl import quadrature_element, element, mixed_element
body = dolfinx.mesh.create_box(MPI.COMM_WORLD,[[0.0,0.0,0.0], [1, 1, 1]], [2, 2, 2], dolfinx.mesh.CellType.hexahedron)
d = body.geometry.dim
def bttm(x):
tol = 1e-10
collapsedV, collapsedMap = V.sub(0).collapse()
zmin = collapsedV.tabulate_dof_coordinates()[:, 0].min() #Z coordinate of the bottom plane
print(zmin)
return np.isclose(x[0],zmin, atol = tol)
U_el = element("CG", body.basix_cell(), 1, shape=(d,))
P_el = element("CG", body.basix_cell(), 1)
V_el = mixed_element([U_el,P_el])
V = dolfinx.fem.functionspace(body, V_el) #Define vector element for displacement field u
up = dolfinx.fem.Function(V)
u, p = ufl.split(up)
v, dp = ufl.TestFunctions(V)
# Deformation gradient
I = ufl.variable(ufl.Identity(d))
F = ufl.variable(I + ufl.grad(u))
gram_A = ufl.variable(F.T*F)
I1 = ufl.variable(ufl.tr(gram_A))
J = ufl.variable(ufl.det(F))
w = (I1 - 3.0)/2.0
metadata = {"quadrature_degree": 4}
dx = ufl.Measure("dx", domain=body, metadata=metadata)
dE = (ufl.inner(ufl.diff(w,F)+p*ufl.diff(J,F),ufl.grad(v))+dp*(J-1.0))*dx
V1, _ = V.sub(0).collapse()
u_D = dolfinx.fem.Function(V1)
u_D.x.array[:] = 0.0
fdim = body.topology.dim - 1
boundary_facets_bttm = dolfinx.mesh.locate_entities_boundary(body, fdim, bttm)
boundary_dofs_b = dolfinx.fem.locate_dofs_topological((V.sub(0), V1), fdim, boundary_facets_bttm)
bc1 = dolfinx.fem.dirichletbc(u_D, boundary_dofs_b, V.sub(0))
dolfinx.log.set_log_level(dolfinx.log.LogLevel.INFO)
problem = NonlinearProblem(dE, up, [bc1])
solver = NewtonSolver(body.comm, problem)
solver.atol = 1e-8
solver.rtol = 1e-8
solver.relaxation_parameter = 1.0
solver.convergence_criterion = "residual"
solver.solve(up)
u, p = up.split()
So in principle nothing should happen (since there are no forces and the only prescribed displacement is 0). However, the following occurs:
[2025-03-04 10:44:14.550] [info] Column ghost size increased from 0 to 0
[2025-03-04 10:44:14.551] [info] Newton iteration 0: r (abs) = 0.7806247497997992 (tol = 1e-08), r (rel) = inf (tol = 1e-08)
[2025-03-04 10:44:14.551] [info] PETSc Krylov solver starting to solve system.
[2025-03-04 10:44:14.552] [info] Newton iteration 1: r (abs) = -nan (tol = 1e-08), r (rel) = -nan (tol = 1e-08)
[2025-03-04 10:44:14.552] [info] PETSc Krylov solver starting to solve system.
[2025-03-04 10:44:14.553] [info] Newton iteration 2: r (abs) = -nan (tol = 1e-08), r (rel) = -nan (tol = 1e-08)
[2025-03-04 10:44:14.553] [info] PETSc Krylov solver starting to solve system.
[2025-03-04 10:44:14.554] [info] Newton iteration 3: r (abs) = -nan (tol = 1e-08), r (rel) = -nan (tol = 1e-08)
...
[2025-03-04 10:44:14.587] [info] PETSc Krylov solver starting to solve system.
[2025-03-04 10:44:14.587] [info] Newton iteration 49: r (abs) = -nan (tol = 1e-08), r (rel) = -nan (tol = 1e-08)
[2025-03-04 10:44:14.588] [info] PETSc Krylov solver starting to solve system.
[2025-03-04 10:44:14.588] [info] Newton iteration 50: r (abs) = -nan (tol = 1e-08), r (rel) = -nan (tol = 1e-08)
Traceback (most recent call last):
File "/home/fenics/shared/Underconstrained_Networks/FEM_Geometries/Hyperelastic_square.py", line 61, in <module>
solver.solve(up)
File "/usr/local/dolfinx-real/lib/python3.12/dist-packages/dolfinx/nls/petsc.py", line 51, in solve
n, converged = super().solve(u.x.petsc_vec)
Why is the 0th iteration resulting in a non zero value and why is this diverging so quickly? I also tried changing the convergence criterion to incremental but that didnt help.
Thanks for any help that you can provide!