JuliaCon 2022 (Times are UTC)

The standard implementation of the Finite Element Method for solid mechanics is based on discretizing the domain into elements and solving the weak form of the point-wise Cauchy's equilibrium equations. This process involves the evaluation of the components of complex tensorial quantities, such as the stress and the stiffness tensor, that are required to calculate the residual force vector and the tangent stiffness matrix. On the other hand the residual force vector and the tangent stiffness matrix coincides, both formally and numerically, with the gradient and the hessian of the free energy of the system, therefore they can be evaluated directly by taking the automatic derivatives of this quantity. The advantage, in this case, is that the free energy is a scalar quantity, which is significantly simpler to evaluate.

In particular, forward mode AD seems particularly suited for the solution of solid mechanics FE problems. In fact, even if FE models can have a very large number of degrees of freedom (DoFs), the free energy of the system in a given configuration is obtained as the sum over the elements of the mesh, and only the degrees of freedom of a single element are involved in the calculation of its contribution to the global residual force vector and tangent stiffness matrix. Therefore we only deal with a limited number of independent variables at time when evaluating individual elements contributions. In this situation a forward mode automatic differentiation scheme implemented through hyper-dual numbers is very efficient for the calculation of higher order derivatives of, however complicate scalar expressions.

The definition of a hyper-dual number system in Julia is particularly straightforward, in fact a structure capable of storing the gradient and the hessian of a quantity, alongside its value, can be simply defined as

struct D2{T,N,M} <:Number 
  v::T
  g::NTuple{N,T}
  h::NTuple{M,T}
end

with
- v the value of the variable,
- g the components of the gradient of v,
- h the Hessian,

and where the type parameters are
- T the type of the values,
- N the number of independent variables controlling the gradient,
- M = N(N+1)/2 is the number of independent elements in the Hessian, we remark that since the hessian is symmetric we only store and operate on half of it.

Thanks to the multiple dispatch feature of Julia, it is sufficient to implement the needed mathematical operator and function over this newly numerical type and the same code that evaluates an expression will also evalute its derivatives, without any change to the original source code. In addition, the usage of macros for implementation of the operations on the gradient and hessian tuples makes the resulting code particualarly efficient.

In this talk we will present AD4SM.jl, a package that implements a second order hyper-dual number system for evaluating both the gradient and hessian of the free energy of a deformable body, allowing the solution of non linear equilibrium problems. The implementation of the dual number system was inspired by the ForwardDiff.jl package, but here the hessian is explicitly introduced, since it is essential for the evaluation of the tangent stiffness matrix, and its symmetry exploited to maximize efficiency. A number of examples are also presented that illustrate how complex non linear problems, with non trivial constraints and boundary conditions, both in two and three dimensions can be numerically stated and solved [1].

References

[1] Vigliotti, A., Auricchio, F., Automatic Differentiation for Solid Mechanics, Archives of Computational Methods in Engineering, Volume 28, Issue 3, Pages 875 - 895 May 2021