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