Limited-memory BFGS phasefield solver for crack simulations

L-BFGS-phasefield-solver/CMakeLists.txt

@@ -0,0 +1,41 @@
+##
+#  CMake script for the linear thermoelastic monolithic approach:
+##
+
+# Set the name of the project and target:
+SET(TARGET "main")
+
+# Declare all source files the target consists of. Here, this is only
+# the one step-X.cc file, but as you expand your project you may wish
+# to add other source files as well. If your project becomes much larger,
+# you may want to either replace the following statement by something like
+#  FILE(GLOB_RECURSE TARGET_SRC  "source/*.cc")
+#  FILE(GLOB_RECURSE TARGET_INC  "include/*.h")
+#  SET(TARGET_SRC ${TARGET_SRC}  ${TARGET_INC})
+# or switch altogether to the large project CMakeLists.txt file discussed
+# in the "CMake in user projects" page accessible from the "User info"
+# page of the documentation.
+SET(TARGET_SRC
+  ${TARGET}.cc
+  SpectrumDecomposition.h
+  SpectrumDecomposition.cc
+  )
+
+# Usually, you will not need to modify anything beyond this point...
+
+CMAKE_MINIMUM_REQUIRED(VERSION 3.3.0)
+
+FIND_PACKAGE(deal.II 9.4.0
+  HINTS ${deal.II_DIR} ${DEAL_II_DIR} ../ ../../ $ENV{DEAL_II_DIR}
+  )
+IF(NOT ${deal.II_FOUND})
+  MESSAGE(FATAL_ERROR "\n"
+    "*** Could not locate a (sufficiently recent) version of deal.II. ***\n\n"
+    "You may want to either pass a flag -DDEAL_II_DIR=/path/to/deal.II to cmake\n"
+    "or set an environment variable \"DEAL_II_DIR\" that contains this path."
+    )
+ENDIF()
+
+DEAL_II_INITIALIZE_CACHED_VARIABLES()
+PROJECT(${TARGET})
+DEAL_II_INVOKE_AUTOPILOT()

L-BFGS-phasefield-solver/README.md

@@ -0,0 +1,130 @@
+## Phasefield L-BFGS monolithic solver
+A monolithic solver based on the limited-memory BFGS (L-BFGS) method for phase-field fracture simulations.
+
+### Purpose
+This repository provides the source code and the input files for the numerical examples used in the paper titled “A novel phase-field monolithic scheme for brittle crack propagation based on the limited-memory BFGS method with adaptive mesh refinement”. The L-BFGS monolithic solver has the following features:
+1. It uses the limited-memory BFGS (L-BFGS) method to overcome the non-convexity of the total energy functional of the phase-field fracture formulation.
+2. It uses the history variable (the maximum of the positive strain energy in history) approach to enforce the phase-field irreversibility.
+3. It adopts an adaptive mesh refinement technique to reduce computational cost.
+4. It works for both 2D and 3D phase-field fracture simulations.
+
+### Content
+The repository contains the following content:
+1. the source code of the L-BFGS method for the phase-field monolithic solver.
+2. the input files for several 2D and 3D phase-field fracture simulations included in the aforementioned manuscript.
+
+### Latest update:
+1. (Sept. 10th, 2025) Add a compiler macro so that the code works for the older version of deal.ii due to the interface change of the function `interpolate()` in the `SolutionTransfer` class.
+2. (Sept. 1st, 2025) Add a new gradient-based line search method. Comparing with the previously implemented Strong Wolfe line search, the new line search method could reduce the wall-clock time by 30 to 50 percent.
+
+### Representative results
+A series of widely adopted phase-field crack benchmark problems are included in this code. Here are some examples:
+1. Simple tension test (pre-refined mesh):
+<p align="center">
+<img src="./doc/Tension_test_prerefine_mesh.png" width="250"><img src="./doc/Tension_test_step_15.png" width="250">
+</p>
+2. Simple shear test (adpatively refined mesh):
+<p align="center">
+<img src="./doc/shear_test_adaptive_refine_step_14_mesh.png" width="250"><img src="./doc/shear_test_adaptive_refine_step_14.png" width="250">
+</p>
+3. 3D torsion test (adpatively refined mesh):
+<p align="center">
+<img src="./doc/torsion_mesh_3.png" width="450"><img src="./doc/torsion_fracture_surface_view_3.png" width="450">
+</p>
+
+### Phase-field model adopted in this work
+The phase-field fracture approach aims to minimize the following energy functional:
+<p align="center">
+<img src="./doc/eq1.png" width="450">
+</p>
+where the approximated crack surface is defined as
+<p align="center">
+<img src="./doc/eq2.png" width="360">
+</p>
+the strain energy density function is based on the isotropic linear elasticity and the additive decomposition
+<p align="center">
+<img src="./doc/eq3.png" width="280">
+</p>
+and the phase-field degradation function adopts the following qudratic form
+<p align="center">
+<img src="./doc/eq4.png" width="150">
+</p>
+
+Using the divergence theorem and the technique of integration by parts, the corresponding Euler-Lagrange equations are written as
+<p align="center">
+<img src="./doc/eq5.png" width="200">
+</p>
+In the actual implementation, however, the governing equations of the cracked solid system are modified as
+<p align="center">
+<img src="./doc/eq6.png" width="250">
+</p>
+where a viscosity regularization term is introduced to stabilize the numerical treatment. Even though this regularization term might not be necessary, we still include it in the formulation for flexibility. This term can always be turned off by setting the viscosity coefficient as zero. Among various approaches to enforce the phase-field irreversibility, the approach based on the history variable of the positive strain energy is adopted due to its relative simplicity,
+<p align="center">
+<img src="./doc/eq8.png" width="200">
+</p>
+
+### Main idea of the limited-memory BFGS (L-BFGS) solver
+The classical BFGS method (a type of quasi-Newton method) involves the following approximate Hessian matrix update during each iteration:
+<p align="center">
+<img src="./doc/eq7.png" width="220">
+</p>
+
+**The problem of this update in the context of finite element simulations is that the second term and the third term both generate a fully dense matrix of n by n needs to be stored, where n represents the number of degrees of freedom.** The above Hessian matrix update is too restrictive even for a mid-size finite element problem due to the required memory for the storage of the fully dense matrix. This limitation motivated this work to introduce the limited-memory feature for the phase-field crack simulations. The limited-memory BFGS method implemented in this work follows the algorithm represented in Chapter 7.2 (page 176) of the following textbook by *Nocedal J, Wright SJ. Numerical optimization (2nd edition). Springer New York, NY, 2006*.
+
+### How to compile
+The L-BFGS finite element procedure is implemented in [deal.II](https://www.dealii.org/) (originally with version 9.4.0 and also works for 9.5.1, it is also tested with the develop branch as Sept. 10th, 2025). In order to use the code (**main.cc**) provided here, deal.II should be configured with MPI and at least with the interfaces to BLAS, LAPACK, Threading Building Blocks (TBB), and UMFPACK. For optional interfaces to other software packages, see https://www.dealii.org/developer/readme.html.
+
+Once the deal.II library is compiled, for instance, to "~/dealii-dev/bin/", follow the steps listed below:
+1. cmake -DDEAL_II_DIR=~/dealii-dev/bin/  .
+2. make debug or make release
+3. make
+
+### How to run
+1. Go into one of the examples folders.
+2. For instance, to run a 2D test case: go into examples/simple_tension_test/
+3. Run via ./../../main 2
+4. For instance, to run a 3D test case: go into examples/3D_torsion/
+5. Run via ./../../main 3
+
+### How to expand this code
+If you want to use the current code to solve new 2D or 3D phase-field crack problems, you need to do the following:
+1. Add a new mesh under the function `void make_grid()`.
+2. Add the boundary conditions for your new mesh in the function `void make_constraints(const unsigned int it_nr)`.
+3. Modify the text file `timeDataFile` for the load/time step sizes and the text file `materialDataFile` for the material properties.
+4. Modify the input file `parameters.prm` accordingly.
+
+If you want to use a new phase-field degradation function (the current code uses the standard quadratic degradation function), you can modify the following functions `double degradation_function(const double d)`, `double degradation_function_derivative(const double d)`, and `double degradation_function_2nd_order_derivative(const double d)`.
+
+If you want to modify the phase-field model completely but still use the L-BFGS monolithic solver, you need to modify the calculations of the initial BFGS matrix
+```
+void assemble_system_B0_one_cell(
+      const typename DoFHandler<dim>::active_cell_iterator &cell,
+      ScratchData_ASM &                                     scratch,
+      PerTaskData_ASM &                                     data) const;
+```
+and the residuals
+```
+void assemble_system_rhs_BFGS_one_cell(
+      const typename DoFHandler<dim>::active_cell_iterator &cell,
+      ScratchData_ASM_RHS_BFGS &                           scratch,
+      PerTaskData_ASM_RHS_BFGS &                           data) const;
+```
+
+### How to cite this work:
+Jin T, Li Z, Chen K. A novel phase-field monolithic scheme for brittle crack propagation based on the limited-memory BFGS method with adaptive mesh refinement. Int J Numer Methods Eng. 2024;e7572. doi: 10.1002/nme.7572
+```
+@Article{2024:jin.li.ea:novel,
+  author  = {Jin, Tao and Li, Zhao and Chen, Kuiying},
+  title   = {A novel phase-field monolithic scheme for brittle crack propagation based on the limited-memory BFGS method with adaptive mesh refinement},
+  journal = {International Journal for Numerical Methods in Engineering},
+  year    = 2024,
+  volume  = 125,
+  number  = 22,
+  pages   = {e7572},
+  month   = nov,
+  issn    = {1097-0207},
+  url     = {https://onlinelibrary.wiley.com/doi/10.1002/nme.7572},
+  doi     = {10.1002/nme.7572},
+  publisher = {Wiley}
+}
+```

L-BFGS-phasefield-solver/SpectrumDecomposition.cc

@@ -0,0 +1,39 @@
+#include "SpectrumDecomposition.h" // use double quote
+#include <deal.II/base/tensor.h>
+#include <deal.II/base/symmetric_tensor.h>
+#include <deal.II/base/exceptions.h>
+
+
+#include <iostream>
+
+namespace usr_spectrum_decomposition
+{
+  void print()
+  {
+    Assert(3 >= 2, dealii::StandardExceptions::ExcNotImplemented());
+    dealii::Tensor<1, 3> Br_tilde;
+    std::cout << Br_tilde << std::endl;
+    std::cout << "Hello world!" << std::endl;
+  }
+
+  double positive_ramp_function(const double x)
+  {
+    return std::fmax(x, 0.0);
+  }
+
+  double negative_ramp_function(const double x)
+  {
+    return std::fmin(x, 0.0);
+  }
+
+  double heaviside_function(const double x)
+  {
+    if (std::fabs(x) < 1.0e-16)
+      return 0.5;
+
+    if (x > 0)
+      return 1.0;
+    else
+      return 0.0;
+  }
+}  

L-BFGS-phasefield-solver/SpectrumDecomposition.h

@@ -0,0 +1,147 @@
+#ifndef usrcodes_spectrum_decomposition_h
+#define usrcodes_spectrum_decomposition_h
+#include <deal.II/base/tensor.h>
+#include <deal.II/base/symmetric_tensor.h>
+#include <deal.II/lac/vector.h>
+#include <deal.II/lac/full_matrix.h>
+#include <deal.II/lac/lapack_full_matrix.h>
+#include <deal.II/base/patterns.h>
+#include <fstream>
+#include <iostream>
+
+
+namespace usr_spectrum_decomposition
+{
+  using namespace dealii;
+  void print();
+
+  double positive_ramp_function(const double x);
+
+  double negative_ramp_function(const double x);
+
+  double heaviside_function(const double x);
+
+  // templated function has to be defined in the header file
+  // perform a spectrum decomposition of a symmetric tensor
+  // input: a symmetric tensor (SymmetricTensor<2, matrix_dimension>)
+  // output: eigenvalues  (Vector<double>)
+  //         eigenvectors (std::vector<Tensor<1, dim>>)
+  template <int dim>
+  void spectrum_decomposition(SymmetricTensor<2, dim> const & symmetric_tensor,
+			      Vector<double> & myEigenvalues,
+			      std::vector<Tensor<1, dim>> & myEigenvectors)
+  {
+
+    const std::array< std::pair< double, Tensor< 1, dim > >, dim >
+      myEigenSystem = eigenvectors(symmetric_tensor);
+
+    for (int i = 0; i < dim; i++)
+      {
+        myEigenvalues[i] = myEigenSystem[i].first;
+        myEigenvectors[i] = myEigenSystem[i].second;
+      }
+  }
+
+  template <int dim>
+  SymmetricTensor<2, dim> positive_tensor(Vector<double> const & eigenvalues,
+					  std::vector<Tensor<1, dim>> const & eigenvectors)
+  {
+    SymmetricTensor<2, dim> positive_part_tensor;
+    positive_part_tensor = 0;
+    for (int i = 0; i < dim; i++)
+      positive_part_tensor += positive_ramp_function(eigenvalues[i])
+                            * symmetrize(outer_product(eigenvectors[i],
+                                                       eigenvectors[i]));
+    return positive_part_tensor;
+  }
+
+  template <int dim>
+  SymmetricTensor<2, dim> negative_tensor(Vector<double> const & eigenvalues,
+					  std::vector<Tensor<1, dim>> const & eigenvectors)
+  {
+    SymmetricTensor<2, dim> negative_part_tensor;
+    negative_part_tensor = 0;
+    for (int i = 0; i < dim; i++)
+      negative_part_tensor += negative_ramp_function(eigenvalues[i])
+                            * symmetrize(outer_product(eigenvectors[i],
+                                                       eigenvectors[i]));
+    return negative_part_tensor;
+  }
+
+  template <int dim>
+  void positive_negative_projectors(Vector<double> const & eigenvalues,
+                                    std::vector<Tensor<1, dim>> const & eigenvectors,
+			            SymmetricTensor<4, dim> & positive_projector,
+				    SymmetricTensor<4, dim> & negative_projector)
+  {
+    Assert(dim <= 3,
+	   ExcMessage("Project tensors only work for dim <= 3."));
+
+    std::array<SymmetricTensor<2, dim>, dim> M;
+    for (int a = 0; a < dim; a++)
+      M[a] = symmetrize(outer_product(eigenvectors[a], eigenvectors[a]));
+
+    std::array<SymmetricTensor<4, dim>, dim> Q;
+    for (int a = 0; a < dim; a++)
+      Q[a] = outer_product(M[a], M[a]);
+
+    std::array<std::array<SymmetricTensor<4, dim>, dim>, dim> G;
+    for (int a = 0; a < dim; a++)
+      for (int b = 0; b < dim; b++)
+	for (int i = 0; i < dim; i++)
+	  for (int j = 0; j < dim; j++)
+	    for (int k = 0; k < dim; k++)
+              for (int l = 0; l < dim; l++)
+        	G[a][b][i][j][k][l] = M[a][i][k] * M[b][j][l]
+				    + M[a][i][l] * M[b][j][k];
+
+    positive_projector = 0;
+    for (int a = 0; a < dim; a++)
+      {
+	double lambda_a = eigenvalues[a];
+	positive_projector += heaviside_function(lambda_a)
+			    * Q[a];
+	for (int b = 0; b < dim; b++)
+	  {
+	    if (b != a)
+	      {
+		double lambda_b = eigenvalues[b];
+		double v_ab = 0.0;
+		if (std::fabs(lambda_a - lambda_b) > 1.0e-12)
+		  v_ab = (positive_ramp_function(lambda_a) - positive_ramp_function(lambda_b))
+		       / (lambda_a - lambda_b);
+		else
+		  v_ab = 0.5 * (  heaviside_function(lambda_a)
+		                + heaviside_function(lambda_b) );
+		positive_projector += 0.5 * v_ab * 0.5 * (G[a][b] + G[b][a]);
+	      }
+	  }
+      }
+
+    negative_projector = 0;
+    for (int a = 0; a < dim; a++)
+      {
+	double lambda_a = eigenvalues[a];
+	negative_projector += heaviside_function(-lambda_a)
+			    * Q[a];
+	for (int b = 0; b < dim; b++)
+	  {
+	    if (b != a)
+	      {
+		double lambda_b = eigenvalues[b];
+		double v_ab = 0.0;
+		if (std::fabs(lambda_a - lambda_b) > 1.0e-12)
+		  v_ab = (negative_ramp_function(lambda_a) - negative_ramp_function(lambda_b))
+		       / (lambda_a - lambda_b);
+		else
+		  v_ab = 0.5 * (  heaviside_function(-lambda_a)
+		                + heaviside_function(-lambda_b) );
+		negative_projector += 0.5 * v_ab * 0.5 * (G[a][b] + G[b][a]);
+	      }
+	  }
+      }
+
+  }
+
+}
+#endif

L-BFGS-phasefield-solver/Utilities.h

@@ -0,0 +1,33 @@
+#ifndef usrcodes_utilities_h
+#define usrcodes_utilities_h
+#include <deal.II/grid/tria.h>
+#include <deal.II/dofs/dof_handler.h>
+
+#include <fstream>
+#include <iostream>
+
+
+namespace usr_utilities
+{
+  using namespace dealii;
+
+  template <int dim>
+  std::vector<types::global_dof_index> get_vertex_dofs(
+    const typename Triangulation<dim>::active_vertex_iterator &vertex,
+    const DoFHandler<dim> &dof_handler)
+  {
+    DoFAccessor<0, dim, dim, false> vertex_dofs(
+        &(dof_handler.get_triangulation()),
+        vertex->level(),
+        vertex->index(),
+        &dof_handler);
+    const unsigned int n_dofs = dof_handler.get_fe().dofs_per_vertex;
+    std::vector<types::global_dof_index> dofs(n_dofs);
+    for (unsigned int i = 0; i < n_dofs; ++i)
+    {
+      dofs[i] = vertex_dofs.vertex_dof_index(0, i);
+    }
+    return dofs;
+  }
+}
+#endif

L-BFGS-phasefield-solver/doc/author

@@ -0,0 +1,2 @@
+Tao Jin <[email protected]>
+