<figure><img src="image_1.png" /></figure>

HAL
open science

# Variational Solvers for Irreversible Evolutionary Systems

Andrés a León Baldelli, Pierluigi Cesana

► To cite this version:

Andrés a León Baldelli, Pierluigi Cesana. Variational Solvers for Irreversible Evolutionary Systems. Journal of Open Source Software, 2024, 10.21105/joss.06897. hal-04795478

HAL Id: hal-04795478

https://hal.sorbonne-universite.fr/hal-04795478v1

Submitted on 9 Dec 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Variational Solvers for Irreversible Evolutionary Systems

Andrés A León Baldelli and Pierluigi Cesana

1 Sorbonne Université, CNRS, Institut Jean Le Rond d'Alembert, F-75005 Paris, France 2 Institute of Mathematics for Industry, Kyushu University, 744 Motooka, Nishi-ku, Fukuoka, 819-0395, Japan. ¶
Corresponding author

DOI: 10.21105/joss.06897

## Software

- Review

- Repository

- Archive

Editor: Mark A. Jensen

## Reviewers:

- @jorgensd

- @jhale

Submitted: 12 April 2024

Published: 05 December 2024

## License

Authors of papers retain copyright and release the work under a Creative Commons Attribution 4.0 International License (CC BY 4.0).

## Summary

We study irreversible evolutionary processes with a general energetic notion of stability. With this contribution, we release three nonlinear variational solvers as modular components (based on FEniCSx/dolfinx) that address three mathematical optimisation problems. They are general enough to apply, in principle, to evolutionary systems with instabilities, jumps, and emergence of patterns. Systems with these qualities are commonplace in diverse arenas spanning from quantum to continuum mechanics, economy, social sciences, and ecology. Our motivation proceeds from fracture mechanics, with the ultimate goal of deploying a transparent numerical platform for scientific validation and prediction of large scale natural fracture phenomena. Our solvers are used to compute one solution to a problem encoded in a system of two inequalities: one (pointwise almost-everywhere) constraint of irreversibility and one global energy statement.

## Statement of need

Quasi-static evolution problems arising in fracture are strongly nonlinear (Marigo, 2023), (Bourdin et al., 2008). They can admit multiple solutions, or none (León Baldelli & Maurini, 2021). This demands both a functional theoretical framework and practical computational tools for real case scenarios. Due to the lack of uniqueness of solutions, it is fundamental to leverage the full variational structure of the problem and investigate solutions up to second order, to detect nucleation of stable modes and transitions of unstable states. The stability of a multiscale system along its nontrivial evolutionary paths in phase space is a key property that is difficult to check: numerically, for real case scenarios with several length scales involved, and analytically, in the infinite-dimensional setting. Despite the concept of unilateral stability is classical in the variational theory of irreversible systems (Mielke & Roubíček, 2015) and the mechanics of fracture (Francfort & Marigo, 1998) (see also Nguyen (2000)), few studies have explored second-order criteria for crack nucleation and evolution. Although sporadic, these studies are significant, including (Pham et al., 2011), (Pham & Marigo, 2013), (Sicsic et al., 2014), (León Baldelli & Maurini, 2021), and (Zolesi & Maurini, 2024). The current literature in computational fracture mechanics predominantly focuses on unilateral first-order criteria, systematically neglecting the exploration of higher-order information for critical points. To the best of our knowledge, no general numerical tools are available to address second-order criteria in evolutionary nonlinear irreversible systems and fracture mechanics.

To fill this gap, our nonlinear solvers offer a flexible toolkit for advanced stability analysis of systems which evolve with constraints.
solver = {Hybrid,Bifurcation,Stability}Solver(
    E,                 # An energy (dolfinx.fem.form)
    state,             # A dictionary of fields describing the system
    bcs,               # A list of boundary conditions
    [bounds],          # A list of bounds (lower and upper) for the order parameter
    parameters        # A dictionary of numerical parameters
)

where [bounds]=[lower, upper] are required for the HybridSolver. Calling solver-solving(<args>)
triggers the solution of the corresponding variational problem. Here, <args> depend on the
solver (see the documentation for details).

HybridSolver solves a (first order) constrained nonlinear variational inequality, implementing
a two-phase hybrid strategy which is *ad hoc* for energy models typical of applications in
damage and fracture mechanics. The first phase (iterative alternate minimisation) is based on
a de-facto industry standard, conceived to exploit the (partial, directional) convexity of the
underlying mechanical models (Bourdin et al., 2000). Once an approximate-solution enters
the attraction set around a critical point, the solver switches to perform a fully nonlinear step
solving a block-matrix problem via Newton's method. This guarantees a precise estimation of
the convergence of the first-order nonlinear problem based on the norm of the (constrained)
residual.

BifurcationSolver is a variational eigenvalue solver which uses SLEPc.EPS to explore the lower part of the spectrum of the Hessian of the energy, automatically computed performing two directional derivatives. Constraints are accounted for by projecting the full Hessian onto the subspace of inactive constraints (Jorge Nocedal, 1999). The relevance of this approach is typical of systems with threshold laws. Thus, the solve method returns a boolean value indicating whether the restricted Hessian is positive definite. Internally, the solver stores the lower part of the operators' spectrum as an array.

StabilitySolver solves a constrained variational eigenvalue inequality in a convex cone, to
check whether the (restricted) nonlinear Hessian operator is positive therein. Starting from an
initial guess z₀*, it iteratively computes (eigenvalue, eigenvector) pairs (λₖ, zₖ) converging to a
limit (λ*, z*) (as k → ∞), by implementing a simple projection and scaling algorithm (Moreau,
1962), (Pinto da Costa & Seeger, 2010). The positivity of λ* (the smallest eigenvalue) allows
to conclude on the stability of the current state (or lack thereof), hence effectively solving
P(0). Notice that, if the current state is unstable (λ* < 0), the minimal eigenmode indicates
the direction of energy decrease.

We dedicate a separate contribution to illustrate how the three solvers are algorithmically
combined to solve problem P(0) in the case of fracture. Figure 1 illustrates the numerical
convergence properties of the StabilitySolver in a 1d verification test.

In a supplementary document, we perform a thorough verification of the code through
parametric benchmark for investigating the stability of a 1D mechanical system, providing
analytical expressions used for comparison with numerical solutions, as well as all parameters
(numerical and physical) employed in the calculations. Accuracy and reliability of the solvers
is shown by the close agreement between numerical and analytic solutions in a benchmark
minimisation of (a constrained) Rayleigh ratio, a key problem for applications in structural
mechanics and stability analysis.
<figure><img src="image_3.png" /><figcaption>**Figure 1:** Rate of convergence for StabilitySolver in 1d (cf. benchmark problem in the Addendum). Targets are the eigenvalue $\lim_k \lambda_k$ := $\lambda^*$ (pink) and the associated eigen-vector $x^*$ (error curve in blue). Note that the residual vector (green) for the cone problem need not be zero at a solution.</figcaption></figure>

## Acknowledgements

ALB acknowledges the students of MEC647 (Complex Crack Propagation in Brittle Materials) of the Modélisation Multiphysique Multiéchelle des Matériaux et des Structures master program at ENSTA Paris Tech/École Polytechnique for their contributions, motivation, and feedback; Yves Capdeboscq, Jean-Jacques Marigo, Sebastien Neukirch, and Luc Nguyen, for constructive discussions and one key insight that was crucial for this project. The work of PC was supported by the JSPS Innovative Area grant JP21H00102 and JSPS Grant-in-Aid for Scientific Research (C) JP24K06797. PC holds an honorary appointment at La Trobe University and is a member of GNAMPA.

## References

Balay, S., Abhyankar, S., Adams, M. F., Benson, S., Brown, J., Brune, P., Buschelman, K., Constantinescu, E., Dalcin, L., Dener, A., Eijkhout, V., Faibussowitsch, J., Gropp, W. D., Hapla, V., Isaac, T., Jolivet, P., Karpeev, D., Kaushik, D., Knepley, M. G., ... Zhang, J. (2023). *PETSc/TAO users manual* (ANL-21/39 - Revision 3.20). Argonne National Laboratory. https://doi.org/10.2172/2205494

Baratta, I. A., Dean, J. P., Dokken, J. S., Habera, M., Hale, J. S., Richardson, C. N., Rognes, M. E., Scroggs, M. W., Sime, N., & Wells, G. N. (2023). *DOLFINx: The next generation FEniCS problem solving environment*. Zenodo. https://doi.org/10.5281/zenodo.10447666

Bažant, Z. P. (1988). *Stable States and Paths of Structures with Plasticity or Damage*. *Journal of Engineering Mechanics, 114*. https://doi.org/10.1061/(ASCE)0733-9399(1988)114:12(2003)
Bourdin, B., Francfort, G. A., & Marigo, J.-J. (2000). Numerical experiments in revisited brittle fracture. *Journal of the Mechanics and Physics of Solids*, **48**(4), 797-826. https://doi.org/10.1016/S0022-5096(99)00028-9

Bourdin, B., Francfort, G. A., & Marigo, J.-J. (2008). The variational approach to fracture. In *Journal of Elasticity* (Vol. 91, pp. 5-148). Springer. https://doi.org/10.1007/978-1-4020-6395-4

Dalcin, L. D., Paz, R. R., Kler, P. A., & Cosimo, A. (2011). Parallel distributed computing using Python. *Advances in Water Resources*, **34**(9), 1124-1139. https://doi.org/10.1016/j.advwatres.2011.04.013

Francfort, G. A., & Marigo, J.-J. (1998). Revisiting brittle fracture as an energy minimization problem. *Journal of the Mechanics and Physics of Solids*, **46**(8), 1319-1342. https://doi.org/10.1016/S0022-5096(98)00034-9

Habera, M., & Zilian, A. (2024). *Dolfiny: Python wrappers for DOLFINx*. https://github.com/michalhabera/dolfiny

Hernandez, V., Roman, J. E., & Vidal, V. (2005). SLEPC: A scalable and flexible toolkit for the solution of eigenvalue problems. *ACM Trans. Math. Software*, **31**(3), 351-362. https://doi.org/10.1145/1089014.1089019

Jorge Nocedal, S. W. (1999). *Numerical optimization*. Springer. https://doi.org/10.1007/978-0-387-40065-5

León Baldelli, A. A., & Maurini, C. (2021). Numerical bifurcation and stability analysis of variational gradient-damage models for phase-field fracture. *Journal of the Mechanics and Physics of Solids*, **152**, 104424. https://doi.org/10.1016/j.jmps.2021.104424

Marigo, J.-J. (2023). La mécanique de l'endommagement au secours de la mécanique de la rupture : L'évolution de cette idée en un demi-siècle. *Comptes Rendus Mécanique*. https://doi.org/10.5802/crmeca.156

Mielke, A., & Roubíček, A. (2015). *Rate-independent systems*. Springer New York, NY. https://doi.org/10.1007/978-1-4939-2706-7

Moreau, J. J. (1962). Décomposition orthogonale d'un espace hilbertien selon deux cônes mutuellement polaires. *Comptes Rendus Hebdomadaires Des séances de l'Académie Des Sciences*, **255**, 238-240.

Nguyen, Q. S. (1994). Bifurcation and stability in dissipative media (plasticity, friction, fracture). *Applied Mechanics Reviews*, **47**(1), 1-31. https://doi.org/10.1115/1.3111068

Nguyen, Q. S. (2000). Standard dissipative systems and stability analysis. In G. Maugin, R. Drouot, & F. Sidoroff (Eds.), *Continuum thermodynamics: The art and science of modeling matter's behavior* (pp. 343-354). Springer. https://doi.org/10.1007/0-306-46946-4

Petryk, H., & Thermann, K. (1985). Second-order bifurcation in elastic-plastic solids. *Journal of the Mechanics and Physics of Solids*, **33**(6), 577-593. https://doi.org/10.1016/0022-5096(85)90004-3

Pham, K., & Marigo, J.-J. (2013). Stability of Homogeneous States with Gradient Damage Models: Size Effects and Shape Effects in the Three-Dimensional Setting. *Journal of Elasticity*, **110**(1), 63-93. https://doi.org/10.1007/s10659-012-9382-5

Pham, K., Marigo, J.-J., & Maurini, C. (2011). The issues of the uniqueness and the stability of the homogeneous response in uniaxial tests with gradient damage models. *Journal of the Mechanics and Physics of Solids*, **59**(6), 1163-1190. https://doi.org/10.1016/j.jmps.2011.03.010

Pinto da Costa, A., & Seeger, A. (2010). Cone-constrained eigenvalue problems: theory and algorithms. *Computational Optimization and Applications*, **45**(1), 25-57. https://doi.org/10.1007/s10659-010-0387-9
//doi.org/10.1007/s10589-008-9167-8

Sicsic, P., Marigo, J.-J., & Maurini, C. (2014). Initiation of a periodic array of cracks in the thermal shock problem: A gradient damage modeling. *Journal of the Mechanics and Physics of Solids*, **63**, 256-284. https://doi.org/10.1016/j.jmps.2013.09.003

Zolesi, C., & Maurini, C. (2024). Stability and crack nucleation in variational phase-field models of fracture: Effects of length-scales and stress multi-axiality. *Journal of the Mechanics and Physics of Solids*, **192**, 105802.