Hello there,
I am trying to replicate an analysis in a paper. I need to load a hyperelastic body until a certain instability. To this end, I would like to check the stiffness matrix during NR iterations, to check when it loses positive definiteness.
As a little test to make sure I was getting the right quantity out, I compare the Jacobian matrix I am able to extract to the stiffness matrix of the linear problem. As the applied bc are “very small” (and one iterations only), the two should be very close (and indeed, the solutions are very close). But for the stiffness matrices there is a significant difference for some elements.
Here a MWE.
import numpy as np
import ufl
from petsc4py import PETSc
from mpi4py import MPI
import dolfinx
from dolfinx import fem, mesh, plot, io
from petsc4py.PETSc import ScalarType
from dolfinx.fem import locate_dofs_geometrical, locate_dofs_topological
L = 1.0
H = 1.0
domain = mesh.create_unit_square(MPI.COMM_WORLD, 5,5, mesh.CellType.quadrilateral)
V = fem.VectorFunctionSpace(domain, ("CG", 1))
def clamped_bottom(x):
return np.isclose(x[1], 0)
u_zero = np.array((0,)*domain.geometry.dim, dtype=ScalarType)
bc_bottom = fem.dirichletbc(u_zero, locate_dofs_geometrical(V, clamped_bottom), V)
def top(x):
return np.isclose(x[1], 1.0)
V_x = V.sub(0).collapse()[0]
V_y = V.sub(1).collapse()[0]
blocked_dofs_top_y = dolfinx.fem.locate_dofs_geometrical((V.sub(1), V_y), top)
non_zero_u_y = dolfinx.fem.Function(V_y)
with non_zero_u_y.vector.localForm() as bc_local:
bc_local.set(0.0)
bc_top = dolfinx.fem.dirichletbc(non_zero_u_y, blocked_dofs_top_y, V.sub(1))
boundary_top = dolfinx.mesh.locate_entities_boundary(domain, domain.topology.dim-1, top)
boundary_dofs_x = locate_dofs_topological(V.sub(0), domain.topology.dim-1, boundary_top)
bc_top_x = fem.dirichletbc(ScalarType(0.), boundary_dofs_x, V.sub(0))
bcs = [bc_bottom, bc_top_x, bc_top]
v = ufl.TestFunction(V)
u = fem.Function(V)
# Spatial dimension
d = len(u)
# Identity tensor
I = ufl.variable(ufl.Identity(d))
# Deformation gradient
F = ufl.variable(I + ufl.grad(u))
# Right Cauchy-Green tensor
C = ufl.variable(F.T * F)
# Invariants of deformation tensors
Ic = ufl.variable(ufl.tr(C))
J = ufl.variable(ufl.det(F))
# Elasticity parameters
E = PETSc.ScalarType(1.0e4)
nu = PETSc.ScalarType(0.3)
mu = fem.Constant(domain, E/(2*(1 + nu)))
lmbda = fem.Constant(domain, E*nu/((1 + nu)*(1 - 2*nu)))
lmbda = fem.Constant(domain,PETSc.ScalarType(8000))
# Stored strain energy density (compressible neo-Hookean model)
psi = (mu / 2) * (Ic - 3) - mu * ufl.ln(J) + (lmbda / 2) * (ufl.ln(J))**2
# Stress
# Hyper-elasticity
P = ufl.diff(psi, F)
metadata = {"quadrature_degree": 4}
#ds = ufl.Measure('ds', domain=domain, subdomain_data=facet_tag, metadata=metadata)
dx = ufl.Measure("dx", domain=domain, metadata=metadata)
F = ufl.inner(ufl.grad(v), P)*dx
problem = fem.petsc.NonlinearProblem(F, u, bcs)
from dolfinx import nls
solver = nls.petsc.NewtonSolver(domain.comm, problem)
# Set Newton solver options
solver.atol = 1e-8
solver.rtol = 1e-8
solver.convergence_criterion = "incremental"
from dolfinx import log
log.set_log_level(log.LogLevel.INFO)
d_max = 0.005
iter = 2
for d_i, d in enumerate(np.linspace(0, d_max, iter)):
with non_zero_u_y.vector.localForm() as bc_local:
bc_local.set(d)
# bc_top_z = fem.dirichletbc(ScalarType(d), boundary_dofs_z, V.sub(2))
# problem = fem.petsc.NonlinearProblem(F, u, [bc_bottom, bc_top_x, bc_top_y, bc_top_z])
# solver = nls.petsc.NewtonSolver(domain.comm, problem)
bc_top = dolfinx.fem.dirichletbc(non_zero_u_y, blocked_dofs_top_y, V.sub(1))
num_its, converged = solver.solve(u)
assert(converged)
u.x.scatter_forward()
print(f"Time step {d_i}, Number of iterations {num_its}, Displacement {d }")
########## LINEAR PROBLEM
u_lin = ufl.TrialFunction(V)
v_lin = ufl.TestFunction(V)
def P_lin(u_lin):
# return lmbda* ufl.nabla_div(u_lin) * ufl.Identity(len(u_lin)) + 2*ufl.sym(ufl.grad(u))
return 2.0 * mu * ufl.sym(ufl.grad(u_lin)) + lmbda * ufl.tr(ufl.sym(ufl.grad(u_lin))) * I
def epsilon(u_lin):
return ufl.sym(ufl.grad(u_lin)) # Equivalent to 0.5*(ufl.nabla_grad(u) + ufl.nabla_grad(u).T)
F_lin = ufl.inner(P_lin(u_lin),epsilon(v_lin))*ufl.dx
f_zero = fem.Constant(domain, ScalarType((0, 0)))
L_zero = ufl.dot(f_zero, v_lin) * ufl.dx
problem_lin = fem.petsc.LinearProblem(F_lin, L_zero, bcs, petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
u_lin = problem_lin.solve()
So one can check that indeed u.x.array and u_lin.x.array are very close.
Now the stiffness matrices.
For the nonlinear problem
A = solver._A
A.getValuesCSR()
A.convert('dense')
A_mat = A.getDenseArray()
an for the linear case
A_lin = fem.petsc.assemble_matrix(fem.form(F_lin), bcs = bcs)
A_lin.assemble()
A_lin.getValuesCSR()
A_lin.convert('dense')
A_lin_mat = A_lin.getDenseArray()
There are elements that are quite different, for example A_lin_mat[-3,-4] and A_mat[-3,-4] .
I checked the rhs of the linear systems are also very close.
I am not grasping why the stiffness matrices have so different elements. The linear systems (internally constructed in both cases) yield the same results, is it some numbering subtlety I am missing, or what? Tried to check on smaller meshes but the problem does not appear.