Assemble_scalar of hyperelasticity PDE yields NaN

Hello everybody,

my goal is to calculate a single, scalar residual of a PDE. Therefore, create a dolfinx.fem.form from my weak form, provide explicit functions as a “trial” and “test” function, and finally call dolfinx.fem.assemble_scalar on the form to yield a single value for my residual.

This worked fine, as long as I used a Poisson equation (see this entry from the forum, which I recommend to read first). Now, I want to switch to a hyperelastic displacement PDE and for simplicity I use the basic example from dolfinx tutorial.

If you only use zeros as input for the “trial” and “test” functions, the residual is 0 (as one would expect). If I, however, try to calculate the single residual using a random input, I receive NaN, and I do not understand why or what I can change to make it work as intended. Any help is most welcome.

Here is a MWE:

import numpy as np
import dolfinx as dfx
import ufl
from mpi4py import MPI
from petsc4py.PETSc import ScalarType


# Create mesh and define function space
mesh = dfx.mesh.create_unit_square(MPI.COMM_WORLD, 8, 8)
V = dfx.fem.VectorFunctionSpace(mesh, ('CG', 2))
Vc = dfx.fem.FunctionSpace(mesh, ('DG', 0))

# load for the right hand side --> zero for simplicity
f_load = dfx.fem.Constant(mesh, ScalarType((0, 0)))

# "Trial" function (with initial random numbers)
u = dfx.fem.Function(V)
u.x.array[:] = np.random.rand(len(u.x.array))

# "Test" function (with initial random numbers)
v = dfx.fem.Function(V)
v.x.array[:] = np.random.rand(len(v.x.array))

# define material fields for Lamée constants mu and lambda
mu = dfx.fem.Constant(mesh, 0.385e4)
lmbda = dfx.fem.Constant(mesh, 0.577e4)

# Spatial dimension
d = len(u)

# Identity tensor
I = 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))

# formula for hyper-elastic potential (psi)
psi = (mu / 2) * (Ic - 3) - mu * ufl.ln(J) + (lmbda / 2) * (ufl.ln(J)) ** 2

# left hand side of the PDE
a = ufl.inner(ufl.grad(v), ufl.diff(psi, F)) * ufl.dx

# right hand side of the PDE
L = ufl.dot(f_load, v) * ufl.dx

# form for my residual
res_ufl = a - L
res_form = dfx.fem.form(res_ufl)

# calculate my residual
residual = dfx.fem.assemble_scalar(res_form)

# This gives nan for some reason
print(residual)

Thanks a lot,
Vincent

You have a negative J. When you take ln(J) (in the expression for psi), you get NaN.

You get a negative Jacobian because you are prescribing a highly non-physical displacement field. In fact, it looks like this:

import numpy as np
import dolfinx as dfx
import ufl
from mpi4py import MPI
from petsc4py.PETSc import ScalarType


# Create mesh and define function space
n = 8
mesh = dfx.mesh.create_unit_square(MPI.COMM_WORLD, n, n)
V = dfx.fem.VectorFunctionSpace(mesh, ('CG', 2))
Vc = dfx.fem.FunctionSpace(mesh, ('DG', 2))

# load for the right hand side --> zero for simplicity
f_load = dfx.fem.Constant(mesh, ScalarType((0, 0)))

# "Trial" function (with initial random numbers)
np.random.seed(196800207)
u = dfx.fem.Function(V)
u.x.array[:] = np.random.rand(len(u.x.array))

# "Test" function (with initial random numbers)
v = dfx.fem.Function(V)
v.x.array[:] = np.random.rand(len(v.x.array))

# define material fields for Lamée constants mu and lambda
mu = dfx.fem.Constant(mesh, 0.385e4)
lmbda = dfx.fem.Constant(mesh, 0.577e4)

# Spatial dimension
d = len(u)

# Identity tensor
I = 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))

# formula for hyper-elastic potential (psi)
psi = (mu / 2) * (Ic - 3) - mu * ufl.ln(J) + (lmbda / 2) * (ufl.ln(J)) ** 2

# left hand side of the PDE
a = ufl.inner(ufl.grad(v), ufl.diff(psi, F)) * ufl.dx

# right hand side of the PDE
L = ufl.dot(f_load, v) * ufl.dx

# form for my residual
res_ufl = a - L
res_form = dfx.fem.form(res_ufl)

# calculate my residual
residual = dfx.fem.assemble_scalar(res_form)

# This gives nan for some reason
print(residual)

jac = dfx.fem.Function(Vc)
jac.interpolate(dfx.fem.Expression(J, Vc.element.interpolation_points()))

u.name = "u"
jac.name = "J"
with dfx.io.XDMFFile(MPI.COMM_WORLD, "u.xdmf", "w") as xdmf:
    xdmf.write_mesh(mesh)
    xdmf.write_function(u)
with dfx.io.XDMFFile(MPI.COMM_WORLD, "J.xdmf", "w") as xdmf:
    xdmf.write_mesh(mesh)
    xdmf.write_function(jac)

2 Likes