Implicit model

Table of Contents

• Problem description
• Defining the system of equations
• Solving the system of equations
• Remarks on the implicit function theorem

Problem description

One of the most notable differences between constitutive models and feed-forward neural networks is that updating certain stiff systems with implicit methods is often more computationally efficient compared to explicit algorithms.

A generally nonlinear, recursive, implicit system of equations take the following form

\begin{align*} \mathbf{r}(\tilde{\mathbf{s}}) & = f(\tilde{\mathbf{s}}, \mathbf{g}; \mathbf{p}), \\ \mathbf{s} &= \mathop{\mathrm{root}}\limits_{\tilde{\mathbf{s}}} (\mathbf{r}). \end{align*}

Here \( \mathbf{r} \) represents the residual (for root-fiding) of the system of equations, \( \mathbf{s} \) represents the unknowns to be solved for, and \( \mathbf{g} \) represents the given variables

The Perzyna viscoplasticity model mentioned in a previous tutorial takes this form:

\begin{align} \boldsymbol{\varepsilon}^e & = \boldsymbol{\varepsilon} - \boldsymbol{\varepsilon}^p, \label{1} \\ \boldsymbol{\sigma} & = 3K\operatorname{vol}\boldsymbol{\varepsilon}^e + 2G\operatorname{dev}\boldsymbol{\varepsilon}^e, \label{2} \\ \bar{\sigma} & = J_2(\boldsymbol{\sigma}), \label{3} \\ f^p & = \bar{\sigma} - \sigma_y, \label{4} \\ \boldsymbol{N} & = \pdv{f^p}{\boldsymbol{\sigma}}, \label{5} \\ \dot{\gamma} & = \left( \dfrac{\left< f^p \right>}{\eta} \right)^n, \label{6} \\ \dot{\boldsymbol{\varepsilon}}^p & = \dot{\gamma} \boldsymbol{N}. \label{7} \end{align}

Its residual can be defined using backward-Euler time integration as

\begin{align} \mathbf{r}\left( \boldsymbol{\varepsilon}^p \right) = \boldsymbol{\varepsilon}^p - \boldsymbol{\varepsilon}^p_n - \left( t - t_n \right) \dot{\boldsymbol{\varepsilon}}^p, \label{8} \end{align}

the solution procedure of which is often referred to as the return-mapping algorithm in the solid mechanics community.

This tutorial illustrates the use of ImplicitUpdate in conjunction with a Newton-Raphson solver to perform the implicit update.

Defining the system of equations

After searching among existing models provided by NEML2, we are able to translate each of these equations into a NEML2 model. The input file looks like

[Models]
  [eq1]
    type = SR2LinearCombination
    from = 'strain plastic_strain'
    to = 'elastic_strain'
    weights = '1 -1'
  []
  [eq2]
    type = LinearIsotropicElasticity
    coefficients = '1e5 0.3'
    coefficient_types = 'YOUNGS_MODULUS POISSONS_RATIO'
    strain = 'elastic_strain'
  []
  [eq3]
    type = SR2Invariant
    invariant_type = 'VONMISES'
    tensor = 'stress'
    invariant = 'effective_stress'
  []
  [eq4]
    type = YieldFunction
    yield_stress = 5
  []
  [surface]
    type = ComposedModel
    models = 'eq3 eq4'
  []
  [eq5]
    type = Normality
    model = 'surface'
    function = 'yield_function'
    from = 'stress'
    to = 'flow_direction'
  []
  [eq6]
    type = PerzynaPlasticFlowRate
    reference_stress = 100
    exponent = 2
    yield_function = 'yield_function'
  []
  [eq7]
    type = AssociativePlasticFlow
  []
  [eq8]
    type = SR2BackwardEulerTimeIntegration
    variable = 'plastic_strain'
  []
  [system]
    type = ComposedModel
    models = 'eq1 eq2 surface eq5 eq6 eq7 eq8'
  []
[]

Note the one-to-one correspondance between the models and the equations.

The structure of the system of equations can be summarized using the code below.

• C++

  #include "neml2/neml2.h"
   
  int
  main()
  {
    using namespace neml2;
    set_default_dtype(kFloat64);
    auto system = load_model("input1.i", "system");
    std::cout << *system << std::endl;
  }

  Output:

  Name:       system
  Input:      plastic_strain [SR2]
              plastic_strain~1 [SR2]
              strain [SR2]
              t [Scalar]
              t~1 [Scalar]
  Output:     plastic_strain_residual [SR2]
  Parameters: eq2_E [Scalar][Double][cpu]
              eq2_nu [Scalar][Double][cpu]
              eq4_sy [Scalar][Double][cpu]
              eq6_eta [Scalar][Double][cpu]
              eq6_n [Scalar][Double][cpu]
  Buffers:    eq1_w_0 [Scalar][Double][cpu]
              eq1_w_1 [Scalar][Double][cpu]

• Python

  import neml2
  import torch
   
  torch.set_default_dtype(torch.double)
  system = neml2.load_model("input1.i", "system")
  print(system)

  Output:

  Name:       system
  Input:      plastic_strain [SR2]
              plastic_strain~1 [SR2]
              strain [SR2]
              t [Scalar]
              t~1 [Scalar]
  Output:     plastic_strain_residual [SR2]
  Parameters: eq2_E [Scalar][Double][cpu]
              eq2_nu [Scalar][Double][cpu]
              eq4_sy [Scalar][Double][cpu]
              eq6_eta [Scalar][Double][cpu]
              eq6_n [Scalar][Double][cpu]
  Buffers:    eq1_w_0 [Scalar][Double][cpu]
              eq1_w_1 [Scalar][Double][cpu]

Solving the system of equations

Once the system of equations are properly defined, we can use the ImplicitUpdate to solve the system of equations. The ImplicitUpdate is responsible for the following:

• (Re)declare the solution to the system of equations as output variables.
• Validate the shape of input variables and residual variables to make sure the system is square.
• Invoke a solver to solve the system of equations.
• Apply the implicit function theorem to calculate exact derivatives (up to machine precision).

NEML2 offers two fully vectorized Newton solvers to be used in conjunction with ImplicitUpdate:

• Newton, the (vectorized) vanilla version of the Newton-Raphson algorithm which always takes the "full" step.
• NewtonWithLineSearch, similar to Newton but offers several commonly used (again fully vectorized) line search strategies.

In addition, the assembly routines as well as the application of the implicit function theorem are also implemented in a vectorized fashion.

In the input file, the nonlinear system is defined under the EquationSystems section (e.g., type = NonlinearSystem with a model), and the nonlinear solver references a linear solver (e.g., type = DenseLU).

The additional sections needed in the input file are

[EquationSystems]
  [eq_sys]
    type = NonlinearSystem
    model = 'system'
    unknowns = 'plastic_strain'
  []
[]
 
[Solvers]
  [newton]
    type = Newton
    rel_tol = 1e-08
    abs_tol = 1e-10
    max_its = 50
    verbose = true
    linear_solver = 'lu'
  []
[]
 
 
[Models]
  [model]
    type = ImplicitUpdate
    equation_system = 'eq_sys'
    solver = 'newton'
  []
[]

The ImplicitUpdate model can then be invoked in the same way as regular models.

• C++

  #include "neml2/neml2.h"
  #include "neml2/tensors/Scalar.h"
  #include "neml2/tensors/SR2.h"
   
  int
  main()
  {
    using namespace neml2;
    set_default_dtype(kFloat64);
    auto model = load_model("input2.i", "model");
   
    // Create input variables
    // Unspecified variables are assumed to be zero
    auto strain = SR2::fill(0.01, 0.005, -0.001);
    auto t = Scalar::full(1);
   
    // Solve the implicit model
    auto outputs = model->value({{"strain", strain}, {"t", t}});
   
    // Get the solution
    std::cout << "\nPlastic strain:\n" << outputs["plastic_strain"] << std::endl;
  }

  Output:

  ITERATION   0, |R| = 3.540990e+01, |R0| = 3.540990e+01
  ITERATION   1, |R| = 8.850542e+00, |R0| = 3.540990e+01
  ITERATION   2, |R| = 2.210703e+00, |R0| = 3.540990e+01
  ITERATION   3, |R| = 5.507485e-01, |R0| = 3.540990e+01
  ITERATION   4, |R| = 1.357799e-01, |R0| = 3.540990e+01
  ITERATION   5, |R| = 3.211516e-02, |R0| = 3.540990e+01
  ITERATION   6, |R| = 6.470185e-03, |R0| = 3.540990e+01
  ITERATION   7, |R| = 7.366970e-04, |R0| = 3.540990e+01
  ITERATION   8, |R| = 1.601343e-05, |R0| = 3.540990e+01
  ITERATION   9, |R| = 8.269535e-09, |R0| = 3.540990e+01
   
  Plastic strain:
  0.001 *
   5.2193
   0.3262
  -5.5455
   0.0000
   0.0000
   0.0000
  [ CPUDoubleType{6} ]

• Python

  import neml2
  from neml2.tensors import Scalar, SR2
  import torch
   
  torch.set_default_dtype(torch.double)
  model = neml2.load_model("input2.i", "model")
   
  # Create input variables
  # Unspecified variables are assumed to be zero
  strain = SR2.fill(0.01, 0.005, -0.001)
  t = Scalar.full(1)
   
  # Solve the implicit model
  outputs = model.value({"strain": strain, "t": t})
   
  # Get the solution
  print("\nPlastic strain:")
  print(outputs["plastic_strain"])

  Output:

  ITERATION   0, |R| = 3.540990e+01, |R0| = 3.540990e+01
  ITERATION   1, |R| = 8.850542e+00, |R0| = 3.540990e+01
  ITERATION   2, |R| = 2.210703e+00, |R0| = 3.540990e+01
  ITERATION   3, |R| = 5.507485e-01, |R0| = 3.540990e+01
  ITERATION   4, |R| = 1.357799e-01, |R0| = 3.540990e+01
  ITERATION   5, |R| = 3.211516e-02, |R0| = 3.540990e+01
  ITERATION   6, |R| = 6.470185e-03, |R0| = 3.540990e+01
  ITERATION   7, |R| = 7.366970e-04, |R0| = 3.540990e+01
  ITERATION   8, |R| = 1.601343e-05, |R0| = 3.540990e+01
  ITERATION   9, |R| = 8.269535e-09, |R0| = 3.540990e+01
   
  Plastic strain:
  0.001 *
   5.2193
   0.3262
  -5.5455
   0.0000
   0.0000
   0.0000
  [ CPUDoubleType{6} ]
  <Tensor of shape [][][6]>

Remarks on the implicit function theorem

Unlike other regular models, declaring variables on the correct sub-axes is important because NEML2 relies on the reserved sub-axes (state, forces, etc.)

• To determine whether the variable value and derivatives should be calculated during the assembly of residual and Jacobian. For example, the derivatives with respect to all variables on the forces sub-axis are skipped because they are not required in the assembly of the linear system.
• To efficiently reuse the factorization of the system Jacobian when applying the implicit function theorem.

Note that NEML2 guarantees the correctness of the variable derivatives after one or more implicit updates, up to machine precision. The same guarantee also applies to user-defined custom models.

This is a significant advantage compared to some of the alternative constitutive model libraries, especially in the context of coupling with external PDE solvers. For example, in the context of finite element method, thermodynamic forces (e.g. strain) are calculated at each quadrature point, and the constitutive library (e.g. NEML2) is responsible for updating the thermodynamic state variables (e.g. stress, plastic strain, etc.), which are then used in the residual definition of the discretized PDE. Therefore, the exact derivatives of the state variables with respect to the forces are the key to the assembly of the exact Jacobian of the descretized PDE, which is in turn the fundamental requirement for optimal convergence for many nonlinear solvers.

Previous	Next
Model parameters (revisited)	Transient driver