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{array}{rl}\mathbf{r}(\mathbf{s})& =f(\mathbf{s},\mathbf{f},{\mathbf{s}}_n,{\mathbf{f}}_n;\mathbf{p}),\\ \mathbf{s}& =\underset{\mathbf{s}}{\mathrm{root}}(\mathbf{r}).\end{array}$$

Here $\mathbf{r}$ represents the residual (for root-fiding) of the system of equations, and $\mathbf{s}$, $\mathbf{f}$, ${\mathbf{s}}_n$, ${\mathbf{f}}_n$ are defined by four reserved sub-axes in NEML2:

• The state sub-axis hosts input variables in set $\mathbf{s}$. The variables on this sub-axis are the primary unknowns to be solved for.
• The forces sub-axis hosts prescribed input variables in set $\mathbf{f}$. These variables are prescribed and, by definition, do not change while the system is being solved.
• The old_state and old_forces sub-axes respectively correspond to ${\mathbf{s}}_n$ and ${\mathbf{f}}_n$. These variables correspond to the previous solution to the system to facilitate the recursive definition of internal variables in history-dependent models. The equivalent plastic strain in plasticity models is a well-known example.

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

$$\begin{array}{}\text{(1)}& {\mathit{\epsilon }}^e& =\mathit{\epsilon }-{\mathit{\epsilon }}^p,\text{(2)}& \mathit{\sigma }& =3K\mathrm{vol}{\mathit{\epsilon }}^e+2G\mathrm{dev}{\mathit{\epsilon }}^e,\text{(3)}& \overline{\sigma }& =J_2(\mathit{\sigma }),\text{(4)}& f^p& =\overline{\sigma }-{\sigma }_y,\text{(5)}& \mathit{N}& =\frac{\partial f^p}{\partial \mathit{\sigma }},\text{(6)}& \dot{\gamma }& ={\left({\displaystyle \frac{⟨f^p⟩}{\eta }}\right)}^n,\text{(7)}& {\dot{\mathit{\epsilon }}}^p& =\dot{\gamma }\mathit{N}.\end{array}$$

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

$$\begin{array}{}\text{(8)}& \mathbf{r}\left({\mathit{\epsilon }}^p\right)={\mathit{\epsilon }}^p-{\mathit{\epsilon }}_n^p-(t-t_n){\dot{\mathit{\epsilon }}}^p,\end{array}$$

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_var = 'forces/E state/Ep'
    to_var = 'state/Ee'
    coefficients = '1 -1'
  []
  [eq2]
    type = LinearIsotropicElasticity
    coefficients = '1e5 0.3'
    coefficient_types = 'YOUNGS_MODULUS POISSONS_RATIO'
    strain = 'state/Ee'
    stress = 'state/S'
  []
  [eq3]
    type = SR2Invariant
    invariant_type = 'VONMISES'
    tensor = 'state/S'
    invariant = 'state/s'
  []
  [eq4]
    type = YieldFunction
    yield_stress = 5
    yield_function = 'state/fp'
    effective_stress = 'state/s'
  []
  [surface]
    type = ComposedModel
    models = 'vonmises yield'
  []
  [eq5]
    type = Normality
    model = 'surface'
    function = 'state/fp'
    from = 'state/S'
    to = 'state/N'
  []
  [eq6]
    type = PerzynaPlasticFlowRate
    reference_stress = 100
    exponent = 2
    yield_function = 'state/fp'
    flow_rate = 'state/gamma_rate'
  []
  [eq7]
    type = AssociativePlasticFlow
    flow_rate = 'state/gamma_rate'
    flow_direction = 'state/N'
    plastic_strain_rate = 'state/Ep_rate'
  []
  [eq8]
    type = SR2BackwardEulerTimeIntegration
    variable = 'state/Ep'
  []
  [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++Python

• #include "neml2/models/Model.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:      forces/E [SR2]
              forces/t [Scalar]
              old_forces/t [Scalar]
              old_state/Ep [SR2]
              state/Ep [SR2]
  Output:     residual/Ep [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_c_0 [Scalar][Double][cpu]
              eq1_c_1 [Scalar][Double][cpu]

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

  Output:

  Name:       system
  Input:      forces/E [SR2]
              forces/t [Scalar]
              old_forces/t [Scalar]
              old_state/Ep [SR2]
              state/Ep [SR2]
  Output:     residual/Ep [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_c_0 [Scalar][Double][cpu]
              eq1_c_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.
• Assemble the residual vector and Jacobian matrix of the underlying linear system.
• Invoke a solver to solve the system of equations.
• Apply the implicit function theorem to calculate exact derivatives (up to machine precision).

NEML2 offers three 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.
• NewtonWithTrustRegion, the trust-region variant of the Newton-Raphson algorithm, where a scalar-valued, fully vectorized quadratic sub-problem is solved to identify the search direction in each update.

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

The resulting input file looks like

[Models]
  # ...
  # Other models defining the system of equations
  # ...
  [system]
    type = ComposedModel
    models = 'eq1 eq2 surface eq5 eq6 eq7 eq8'
  []
  [model]
    type = ImplicitUpdate
    implicit_model = 'system'
    solver = 'newton'
  []
[]
 
[Solvers]
  [newton]
    type = Newton
    rel_tol = 1e-08
    abs_tol = 1e-10
    max_its = 50
    verbose = true
  []
[]

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

C++Python

• #include "neml2/models/Model.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 E = SR2::fill(0.01, 0.005, -0.001);
    auto t = Scalar::full(1);
   
    // Solve the implicit model
    auto outputs = model.value({{VariableName("forces", "E"), E},
                                {VariableName("forces", "t"), t}});
   
    // Get the solution
    std::cout << "\nPlastic strain:\n" << outputs[VariableName("state", "Ep")] << 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} ]

• 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
  E = SR2.fill(0.01, 0.005, -0.001)
  t = Scalar.full(1)
   
  # Solve the implicit model
  outputs = model.value({"forces/E": E, "forces/t": t})
   
  # Get the solution
  print("\nPlastic strain:")
  print(outputs["state/Ep"])

  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.

As long as models are defined using the correct sub-axis definitions and satisfy some mild continuity requirements, 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