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,