Computing derivative of residual vector w.r.t material parameter for finite element model updating

dolfinx
,
1

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

  1. Given the assembled matrix A above, what is the recommended way in FEniCSx to obtain the “physical” sensitivity vector dR/dm, for a constant material parameter represented as a quadrature Function?

  2. Is there any supported way (with FEMExternalOperator) to obtain dR/dm directly 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