Hi all,
I am doing finite-element model updating for hyperelastic material models where I need analytic sensitivities of the residual w.r.t. constant material parameters:
\frac{\partial \mathbf{R}}{\partial m} \in \mathbb{R}^{N_u}
I do not use the adjoint method, so I need explicit parameter derivatives (implicit-function-theorem based optimization) and afterward solve a linear system with the Jacobian matrix to obtain the derivative of the displacement vector with repect to the material parameter:
\mathbf{J} \frac{\partial \mathbf{u}}{\partial m} = - \frac{\partial \mathbf{R}}{\partial m} , \quad \mathbf{K} = \frac{\partial \mathbf{R}}{\partial \mathbf{u}} \in \mathbb{R}^{N_u \times N_u}
Because my material model uses spline-based energies, I rely on dolfinx_external_operator. This introduces the constraint that material parameters (although constant) must be quadrature Functions, since external operators evaluate operands at quadrature points.
As a consequence, the only working pattern is:
Qe = basix.ufl.quadrature_element(
domain.topology.cell_name(),
degree=quadrature_degree,
value_shape=()
)
Q = fem.functionspace(domain, Qe)
m = fem.Function(Q) # Constant material parameter represented as a field
m.x.array[:] = 2.0
# Define external operator Psi and residual R in terms of m
# ...
dR_dm = ufl.derivative(R, m)
dR_dm_exp = ufl.algorithms.expand_derivatives(dR_dm)
R_repl, ops = replace_external_operators(dR_dm_exp)
A = assemble_matrix(fem.form(R_repl)) # <-- Matrix, not vector
Here, the assembly for dR/dm therefore gives a matrix. However, for a constant scalar parameter the mathematical object I actually want is a vector.
Questions
-
Given the assembled matrix
Aabove, what is the recommended way in FEniCSx to obtain the “physical” sensitivity vectordR/dm, for a constant material parameter represented as a quadratureFunction? -
Is there any supported way (with
FEMExternalOperator) to obtaindR/dmdirectly as a vector, without assembling the full Jacobian in parameter space? In my calibration problems I typically have 20–30 material parameters; storing 20–30 large matrices for all parameter sensitivities is quite expensive.
Minimal Code Example
The code below defines a variational problem for hyperelasticity with a simple linear energy density function expressed as external operator.
from mpi4py import MPI
import numpy as np
from petsc4py import PETSc
from dolfinx import *
import basix
import ufl
from petsc4py import PETSc
import dolfinx.fem as femx
from dolfinx.fem.petsc import NonlinearProblem
from dolfinx.nls.petsc import NewtonSolver
from dolfinx_external_operator import (
FEMExternalOperator,
replace_external_operators,
evaluate_operands,
evaluate_external_operators,
)
default_scalar_type = PETSc.ScalarType
def make_linear_iso_energy_external():
"""
External operator representing dPsi/dI1 for a very simple
isochoric energy Psi(I1_bar; slope) = slope * I1_bar
"""
def psi_external(derivatives):
# derivatives is (k_I1, k_slope)
if derivatives == (0, 0):
return lambda I1_vals, slope_vals: slope_vals.reshape(-1)
elif derivatives == (1, 0):
return lambda I1_vals, slope_vals: np.zeros_like(I1_vals).reshape(-1)
elif derivatives == (0, 1):
return lambda I1_vals, slope_vals: np.ones_like(slope_vals).reshape(-1)
else:
raise NotImplementedError(f"No implementation for derivatives={derivatives}")
return psi_external
if __name__ == "__main__":
comm = MPI.COMM_WORLD
# Mesh and spaces
domain = mesh.create_unit_cube(comm, 1, 1, 1)
V = fem.functionspace(domain, ("CG", 1, (domain.geometry.dim,))) # vector P1
u = fem.Function(V)
v = ufl.TestFunction(V)
# Measure (quadrature of degree 3)
quadrature_degree = 2
dx = ufl.Measure(
"dx",
domain=domain,
metadata={"quadrature_scheme": "default", "quadrature_degree": quadrature_degree},
)
# Kinematics
d = domain.geometry.dim
I = ufl.Identity(d)
F = ufl.variable(I + ufl.grad(u))
C = F.T * F
J = ufl.det(F)
I1 = ufl.tr(C)
# Quadrature space for external operator and material parameter (slope)
Qe = basix.ufl.quadrature_element(
domain.topology.cell_name(),
degree=quadrature_degree,
value_shape=(),
)
Q = fem.functionspace(domain, Qe)
# Slope parameter
slope = fem.Function(Q)
slope.x.array[:] = 5.0 # constant in space, but stored in quadrature field
# External operator for dPsi/dI1_bar depending on two parameters
dPsi_dI1 = FEMExternalOperator(I1, slope, function_space=Q)
dPsi_dI1.external_function = make_linear_iso_energy_external()
# Simple Piola stress: isochoric part via external operator + volumetric penalty
P_iso = ufl.diff(I1, F) * dPsi_dI1
P_vol = ufl.diff(J, F) * (10.0 * (J - 1.0)) # volumetric energy: 5*(J - 1)^2
P = P_iso + P_vol
# Weak form
R = ufl.inner(P, ufl.grad(v)) * dx
# solve the forward problem (not shown here)
# ...
# Residual w.r.t slope parameter
dR_ds = ufl.derivative(R, slope)
dR_ds_exp = ufl.algorithms.expand_derivatives(dR_ds) # Needed to expand derivatives of external operators
dR_ds_replaced, dR_ds_ops = replace_external_operators(dR_ds_exp)
dR_ds_form = femx.form(dR_ds_replaced)
s_ops = evaluate_operands(dR_ds_ops)
evaluate_external_operators(dR_ds_ops, s_ops)
dR_dS_mat = femx.petsc.assemble_matrix(dR_ds_form)
dR_dS_mat.assemble() # <-- Matrix