Flagship project · open source
Differentiable FEA — Beyond the Stiffness Matrix
Matrix-free differentiable finite elements: weak-form residuals and Jacobian–vector products are evaluated through automatic differentiation, enabling Newton–Krylov solution without explicitly assembling the global stiffness matrix.
- Problem
- Traditional nonlinear FEA repeatedly forms and solves large linearized systems involving a global stiffness (Jacobian) matrix. Here, large-strain hyperelastic problems are solved without forming it, and with exact derivatives.
- Method
- A PyTorch solver. Automatic differentiation supplies the residual and the Jacobian–vector products; preconditioned conjugate gradients needs only those products, inside an inexact Newton iteration.
- Physics & math
- Weak-form mechanics: equilibrium is a stationary point of the strain energy, so the residual and the tangent are its first and second derivatives. Design gradients follow from the implicit function theorem and one adjoint solve.
- Result
-
- < 10−13 relative agreement with FEniCSx/dolfinx (≤ 9×10−14) on identical discretizations
- 11.6× faster GPU tangent product after profiling (element-wise versus batched matrix products; 40,000 Q4 elements, T4)
- 2.9 s for a complete nonlinear solve of 10,000 Q4 elements on a T4 GPU
Scope. A linear system is still solved at every Newton step; what is avoided is assembling and storing the global matrix. Matrix-free evaluation does not by itself make every simulation faster: on a CPU in 2D an assembled sparse matrix is still faster per product. The advantage is expected for high-order elements, 3D and GPU memory, and has not been measured there yet. v0.1, 2D plane strain.
Completed: released, validated against FEniCSx and analytic solutionsPlanned: 3D elements, multigrid preconditioning, multi-GPU
Matrix-free Newton–Krylov
- Weak form∫ P(F) : ∇v dV − ℓ(v) = 0Equilibrium as a stationary point of the strain energy.
- Residual R(u)R(u) = fint(u) − fextEvaluated element by element at quadrature points.
- Automatic differentiationd/dε R(u + εv) at ε = 0Forward-mode derivative through the constitutive law.
- J(u)v no global matrixJ(u) vThe action of the tangent on a vector, which is all a Krylov method needs.
- Krylov solverJ δu = −R, solved inexactlyPreconditioned conjugate gradients with a loose, adaptive tolerance.
- Newton updateu ← u + α δuLine search, then repeat.
↺ repeat from the residual until ‖R(u)‖ is below tolerance
