PETSc error code is: 73, Object is in wrong state for simple hyperelastic problem

General
1

Hello Everyone,

I have adopted the hyperelastic Tutorial from Dokken to incorporate another free energy.
My goal is to take the derivative of this free energy \Psi to C:
\frac{d \Psi}{d \boldsymbol{C}} with \boldsymbol{C} = \boldsymbol{F}^T \boldsymbol{F} and \boldsymbol{F} being the deformation gradient.

However, when I do this and then try to solve the system I get the error message:

File “/usr/local/dolfinx-real/lib/python3.10/dist-packages/dolfinx/nls/petsc.py”, line 47, in solve
n, converged = super().solve(u.x.petsc_vec)
RuntimeError: Failed to successfully call PETSc function ‘KSPSolve’. PETSc error code is: 73, Object is in wrong state
2024-06-10 15:42:15.528 ( 0.380s) [main ] loguru.cpp:526 INFO| atexit

I use dolfinx v:0.8

However, when I take the derivative with the deformation gradient, this error does not appear, even when both variables are tensors with the same shape.

I appreciate any advice on what I’m missing, as it seems to be a pretty simple problem but I can not figure it out.

Here is a MME:

from dolfinx import log, default_scalar_type
from dolfinx.fem.petsc import NonlinearProblem
from dolfinx.nls.petsc import NewtonSolver
import numpy as np
import ufl

from mpi4py import MPI
from dolfinx import fem, mesh


def free_energy_mech(_F_def):
    """
    free Energy of the mechanical part
    """
    detF  = ufl.det(_F_def)          # J = det(F)
    _C = ufl.variable(_F_def.T * _F_def)            # Right Cauchy-Green tensor
    Ic = detF**(-2/3)*ufl.tr(_C)    # Invariants of deformation tensors

    # Elasticity parameters
    c1 = default_scalar_type(0.189439) # MPa
    c2 = default_scalar_type(-5.197e-6) # MPa
    c3 = default_scalar_type(0.0024059) # MPa
    D1 = default_scalar_type((0.189439*10)**(-1)) # MPa for nu = 0.452 
    # Stored strain energy density (nearly incompressible Yeoh Model)
    _Psi = c1*(Ic - 3) + c2*(Ic - 3)**2 + c3*(Ic - 3)**3  + 1/D1*(detF-1)**2
    return _Psi


L = 20.0
domain = mesh.create_box(MPI.COMM_WORLD, [[0.0, 0.0, 0.0], [L, 1, 1]], [20, 5, 5], mesh.CellType.hexahedron)
V = fem.functionspace(domain, ("Lagrange", 2, (domain.geometry.dim, )))

def left(x):
    return np.isclose(x[0], 0)


def right(x):
    return np.isclose(x[0], L)


fdim = domain.topology.dim - 1
left_facets = mesh.locate_entities_boundary(domain, fdim, left)
right_facets = mesh.locate_entities_boundary(domain, fdim, right)

# Concatenate and sort the arrays based on facet indices. Left facets marked with 1, right facets with two
marked_facets = np.hstack([left_facets, right_facets])
marked_values = np.hstack([np.full_like(left_facets, 1), np.full_like(right_facets, 2)])
sorted_facets = np.argsort(marked_facets)
facet_tag = mesh.meshtags(domain, fdim, marked_facets[sorted_facets], marked_values[sorted_facets])

u_bc = np.array((0,) * domain.geometry.dim, dtype=default_scalar_type)

left_dofs = fem.locate_dofs_topological(V, facet_tag.dim, facet_tag.find(1))
bcs = [fem.dirichletbc(u_bc, left_dofs, V)]

B = fem.Constant(domain, default_scalar_type((0, 0, 0)))
T = fem.Constant(domain, default_scalar_type((0, 0, 0)))

v = ufl.TestFunction(V)
u = fem.Function(V)

# Spatial dimension
d = len(u)

# Identity tensor
I = ufl.variable(ufl.Identity(d))

# Deformation gradient
F_def = ufl.variable(I + ufl.grad(u))

# Right Cauchy-Green tensor
C = ufl.variable(F_def.T * F_def)

# Invariants of deformation tensors
Ic = ufl.variable(ufl.tr(C))
J = ufl.variable(ufl.det(F_def))

psi = free_energy_mech(F_def)
dPsidC = ufl.diff(psi, C)
# dPsidC = ufl.diff(psi, F_def)

P2nd = 2*dPsidC - 1/3*(ufl.inner(dPsidC,C)*ufl.inv(C)) #2nd Pioloa-Kirchoff
P1st = P2nd* F_def.T

metadata = {"quadrature_degree": 4}
ds = ufl.Measure('ds', domain=domain, subdomain_data=facet_tag, metadata=metadata)
dx = ufl.Measure("dx", domain=domain, metadata=metadata)
# Define form F (we want to find u such that F(u) = 0)
F = ufl.inner(ufl.grad(v), P1st) * dx - ufl.inner(v, B) * dx - ufl.inner(v, T) * ds(2)

problem = NonlinearProblem(F, u, bcs)

solver = NewtonSolver(domain.comm, problem)

# Set Newton solver options
solver.atol = 1e-8
solver.rtol = 1e-8
solver.convergence_criterion = "incremental"

log.set_log_level(log.LogLevel.INFO)
tval0 = -1.5
for n in range(1, 2):
    T.value[2] = n * tval0
    num_its, converged = solver.solve(u)
    assert (converged)
    u.x.scatter_forward()
    print(f"Time step {n}, Number of iterations {num_its}, Load {T.value}")
2

First of all, you are differentiating with respect to the wrong C.
A rewrite like:

from dolfinx import log, default_scalar_type
from dolfinx.fem.petsc import NonlinearProblem
from dolfinx.nls.petsc import NewtonSolver
import numpy as np
import ufl

from mpi4py import MPI
from dolfinx import fem, mesh


def free_energy_mech(_F_def):
    """
    free Energy of the mechanical part
    """
    detF = ufl.det(_F_def)  # J = det(F)
    _C = ufl.variable(_F_def.T * _F_def)  # Right Cauchy-Green tensor
    Ic = detF ** (-2 / 3) * ufl.tr(_C)  # Invariants of deformation tensors

    # Elasticity parameters
    c1 = default_scalar_type(0.189439)  # MPa
    c2 = default_scalar_type(-5.197e-6)  # MPa
    c3 = default_scalar_type(0.0024059)  # MPa
    D1 = default_scalar_type((0.189439 * 10) ** (-1))  # MPa for nu = 0.452
    # Stored strain energy density (nearly incompressible Yeoh Model)
    _Psi = (
        c1 * (Ic - 3)
        + c2 * (Ic - 3) ** 2
        + c3 * (Ic - 3) ** 3
        + 1 / D1 * (detF - 1) ** 2
    )
    return _Psi, _C


L = 20.0
domain = mesh.create_box(
    MPI.COMM_WORLD, [[0.0, 0.0, 0.0], [L, 1, 1]], [20, 5, 5], mesh.CellType.hexahedron
)
V = fem.functionspace(domain, ("Lagrange", 2, (domain.geometry.dim,)))


def left(x):
    return np.isclose(x[0], 0)


def right(x):
    return np.isclose(x[0], L)


fdim = domain.topology.dim - 1
left_facets = mesh.locate_entities_boundary(domain, fdim, left)
right_facets = mesh.locate_entities_boundary(domain, fdim, right)

# Concatenate and sort the arrays based on facet indices. Left facets marked with 1, right facets with two
marked_facets = np.hstack([left_facets, right_facets])
marked_values = np.hstack([np.full_like(left_facets, 1), np.full_like(right_facets, 2)])
sorted_facets = np.argsort(marked_facets)
facet_tag = mesh.meshtags(
    domain, fdim, marked_facets[sorted_facets], marked_values[sorted_facets]
)

u_bc = np.array((0,) * domain.geometry.dim, dtype=default_scalar_type)

left_dofs = fem.locate_dofs_topological(V, facet_tag.dim, facet_tag.find(1))
bcs = [fem.dirichletbc(u_bc, left_dofs, V)]

B = fem.Constant(domain, default_scalar_type((0, 0, 0)))
T = fem.Constant(domain, default_scalar_type((0, 0, 0)))

v = ufl.TestFunction(V)
u = fem.Function(V)

# Spatial dimension
d = len(u)

# Identity tensor
I = ufl.Identity(d)

# Deformation gradient
F_def = ufl.variable(I + ufl.grad(u))

# Right Cauchy-Green tensor

# Invariants of deformation tensors
J = ufl.det(F_def)

psi, C = free_energy_mech(F_def)
dPsidC = ufl.diff(psi, C)
# dPsidC = ufl.diff(psi, F_def)

P2nd = 2 * dPsidC - 1 / 3 * (ufl.inner(dPsidC, C) * ufl.inv(C))  # 2nd Pioloa-Kirchoff
P1st = P2nd * F_def.T

metadata = {"quadrature_degree": 10}
ds = ufl.Measure("ds", domain=domain, subdomain_data=facet_tag, metadata=metadata)
dx = ufl.Measure("dx", domain=domain, metadata=metadata)
# Define form F (we want to find u such that F(u) = 0)
F = ufl.inner(ufl.grad(v), P1st) * dx - ufl.inner(v, B) * dx - ufl.inner(v, T) * ds(2)

gets you further.

However, as I am not a solid mechanician, I cannot comment on the validity of the form, or why it does not converge.

For me the simple case doesn’t even converge, which concerns me.

2 Likes
3

Dear Dokken,

thank you for your very fast reply, now it works!

You are right, I forgot the brackets in the equation. With

P2nd = 2 * (dPsidC - 1 / 3 * (ufl.inner(dPsidC, C) * ufl.inv(C)))

it converges, for the initial step.
Furthermore, the free energy was split into a volumetric and isochoric part here. However, I have only posted the isochoric part here, to keep it more simple.