Model composition

Table of Contents

• Problem description
• Writing the input file
• Running the models: the hard way
• Running the models: the easy way

Problem description

We have been working with the linear, isotropic elasticity model in the previous tutorials. We started with that example because it is arguably the simplest possible material model in the context of solid mechanics. It is simple not just because of the simplicity in the description of the material behavior, but also due to the fact that its mathematical formulation only involves one linear equation.

Much more complicated, nonlinear models can be created using NEML2.

Using a Perzyna-type viscoplasticity model as an example, it can be formulated as

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

The above formulation makes a series of constitutive choices:

• The strain is small and can be additively decomposed into elastic and plastic strains.
• The elastic response is linear and isotropic.
• The plastic flow is isochoric.
• There is no isotropic hardening associated with plasticity.
• There is no kinematic hardening associated with plasticity.
• There is no back stress associated with plasticity.
• The plastic flow associative.
• The plastic rate-sensitivity follows a power-law relation.

Any change in one of the constitutive choices will result in a new material model. Suppose there are a total of $N$ constitutive choices, each having $k$ variants, the total number of possible material models would be $k^N$.

In other words, the number of possible material models grows exponentially with the number of constitutive choices, and implementing all combinations is practically infeasible.

To address such challenge, NEML2 introduces a model composition mechanism which allows multiple models to be "stitched" together in a flexible, modular manner.

This tutorial demonstrates model composition using a much simplified model (without loss of generality). The model can be written as

$$\begin{array}{}\text{(1)}& \overline{a}& =I_1(\mathit{a}),\text{(2)}& \overline{b}& =J_2(\mathit{b}),\text{(3)}& \dot{b}& =\overline{b}\mathit{a}+\overline{a}\mathit{b},\end{array}$$

where $\mathit{a}$ and $\mathit{b}$ are symmetric second order tensors.

Writing the input file

Let us first search for available models describing this set of equations:

• $\text{(1)}\mathrm{\&}\text{(2)}$ correspond to SR2Invariant;
• $\text{(3)}$ corresponds to LinearCombination.

The input file then looks like

[Models]
  [eq1]
    type = SR2Invariant
    tensor = 'forces/a'
    invariant = 'state/a_bar'
    invariant_type = I1
  []
  [eq2]
    type = SR2Invariant
    tensor = 'state/b'
    invariant = 'state/b_bar'
    invariant_type = VONMISES
  []
  [eq3]
    type = SR2LinearCombination
    from_var = 'forces/a state/b'
    to_var = 'state/b_rate'
    coefficients = '1 1'
    coefficient_as_parameter = true
  []
[]

Running the models: the hard way

Now that all three models are defined in the input file, we can load and evaluate them in sequence, with a bit of effort:

C++Python

• #include "neml2/models/Model.h"
  #include "neml2/tensors/SR2.h"
   
  int
  main()
  {
    using namespace neml2;
    set_default_dtype(kFloat64);
    load_input("input.i");
    auto & eq1 = get_model("eq1");
    auto & eq2 = get_model("eq2");
    auto & eq3 = get_model("eq3");
   
    // Create the input variables
    auto a_name = VariableName("forces", "a");
    auto b_name = VariableName("state", "b");
    auto a = SR2::fill(0.1, 0.05, -0.03, 0.02, 0.06, 0.03);
    auto b = SR2::fill(100, 20, 10, 5, -30, -20);
   
    // Evaluate the first model to get a_bar
    auto a_bar_name = VariableName("state", "a_bar");
    auto a_bar = eq1.value({{a_name, a}})[a_bar_name];
   
    // Evaluate the second model to get b_bar
    auto b_bar_name = VariableName("state", "b_bar");
    auto b_bar = eq2.value({{b_name, b}})[b_bar_name];
   
    // Evaluate the third model to get b_rate
    eq3.set_parameter("c_0", b_bar);
    eq3.set_parameter("c_1", a_bar);
    auto b_rate_name = VariableName("state", "b_rate");
    auto b_rate = eq3.value({{a_name, a}, {b_name, b}})[b_rate_name];
   
    std::cout << "b_rate: \n" << b_rate << std::endl;
  }

  Output:

  b_rate: 
   22.6184
    7.7092
   -1.9855
    3.8519
    3.9188
    1.1109
  [ CPUDoubleType{6} ]

• import neml2
  from neml2.tensors import SR2
  import torch
   
  torch.set_default_dtype(torch.double)
  neml2.load_input("input.i")
  eq1 = neml2.get_model("eq1")
  eq2 = neml2.get_model("eq2")
  eq3 = neml2.get_model("eq3")
   
  # Create the input variables
  a = SR2.fill(0.1, 0.05, -0.03, 0.02, 0.06, 0.03)
  b = SR2.fill(100, 20, 10, 5, -30, -20)
   
  # Evaluate the first model to get a_bar
  a_bar = eq1.value({"forces/a": a})["state/a_bar"]
   
  # Evaluate the second model to get b_bar
  b_bar = eq2.value({"state/b": b})["state/b_bar"]
   
  # Evaluate the third model to get b_rate
  eq3.c_0 = b_bar
  eq3.c_1 = a_bar
  b_rate = eq3.value({"forces/a": a, "state/b": b})["state/b_rate"]
   
  print("b_rate:")
  print(b_rate)

  Output:

  b_rate:
   22.6184
    7.7092
   -1.9855
    3.8519
    3.9188
    1.1109
  [ CPUDoubleType{6} ]
  <Tensor of shape [][6]>

Running the models: the easy way

We were able to successfully calculate $\dot{\mathit{b}}$ by

1. calculating $\overline{a}$ by evaluating $\text{(1)}$,
2. calculating $\overline{b}$ by evaluating $\text{(2)}$,
3. setting the two coefficients of $\text{(3)}$ to be $\overline{b}$ and $\overline{a}$ respectively,
4. calculating $\dot{\mathit{b}}$ by evaluating $\text{(3)}$.

However, that is not ideal because we had to

• Manually evaluate the equations and figure out the evaluation order, and
• Manually set the parameters in $\text{(3)}$ as outputs from $\text{(1)}\mathrm{\&}\text{(2)}$.

This manual method is not scalable when the number of equations, variables, and parameters increase.

Using NEML2's model composition capability can address these issues without sacrificing modularity. ComposedModel allows us to compose a new model from the three existing models:

[Models]
  [eq1]
    type = SR2Invariant
    tensor = 'forces/a'
    invariant = 'state/a_bar'
    invariant_type = I1
  []
  [eq2]
    type = SR2Invariant
    tensor = 'state/b'
    invariant = 'state/b_bar'
    invariant_type = VONMISES
  []
  [eq3]
    type = SR2LinearCombination
    from_var = 'forces/a state/b'
    to_var = 'state/b_rate'
    coefficients = 'eq2 eq1'
    coefficient_as_parameter = true
  []
  [eq]
    type = ComposedModel
    models = 'eq1 eq2 eq3'
  []
[]

Note

The names of the other two models are used to specify the coefficients of $\text{(3)}$, i.e. ‘coefficients = 'eq2 eq1’`. This syntax is different from what was covered in the previous tutorial on model parameters and will be explained in more details in the next tutorial.

Let us first inspect the composed model and compare it against the three sub-models:

C++Python

• #include "neml2/models/Model.h"
   
  int
  main()
  {
    using namespace neml2;
   
    load_input("input_composed.i");
    auto & eq1 = get_model("eq1");
    auto & eq2 = get_model("eq2");
    auto & eq3 = get_model("eq3");
    auto & eq = get_model("eq");
   
    std::cout << "eq1:\n" << eq1 << std::endl << std::endl;
    std::cout << "eq2:\n" << eq2 << std::endl << std::endl;
    std::cout << "eq3:\n" << eq3 << std::endl << std::endl;
    std::cout << "eq:\n" << eq << std::endl << std::endl;
  }

  Output:

  eq1:
  Name:       eq1
  Input:      forces/a [SR2]
  Output:     state/a_bar [Scalar]
   
   
  eq2:
  Name:       eq2
  Input:      state/b [SR2]
  Output:     state/b_bar [Scalar]
   
   
  eq3:
  Name:       eq3
  Input:      forces/a [SR2]
              state/a_bar [Scalar]
              state/b [SR2]
              state/b_bar [Scalar]
  Output:     state/b_rate [SR2]
   
   
  eq:
  Name:       eq
  Input:      forces/a [SR2]
              state/b [SR2]
  Output:     state/b_rate [SR2]

• import neml2
  from neml2.tensors import SR2
   
  neml2.load_input("input_composed.i")
  eq1 = neml2.get_model("eq1")
  eq2 = neml2.get_model("eq2")
  eq3 = neml2.get_model("eq3")
  eq = neml2.get_model("eq")
   
  print("eq1:")
  print(eq1, "\n")
  print("eq2:")
  print(eq2, "\n")
  print("eq3:")
  print(eq3, "\n")
  print("eq:")
  print(eq, "\n")

  Output:

  eq1:
  Name:       eq1
  Input:      forces/a [SR2]
  Output:     state/a_bar [Scalar]
   
   
  eq2:
  Name:       eq2
  Input:      state/b [SR2]
  Output:     state/b_bar [Scalar]
   
   
  eq3:
  Name:       eq3
  Input:      forces/a [SR2]
              state/a_bar [Scalar]
              state/b [SR2]
              state/b_bar [Scalar]
  Output:     state/b_rate [SR2]
   
   
  eq:
  Name:       eq
  Input:      forces/a [SR2]
              state/b [SR2]
  Output:     state/b_rate [SR2]

Note that the composed model "eq" automatically:

• Identified the input variables $\mathit{a}$ and $\mathit{b}$,
• Identified the output variable $\dot{\mathit{b}}$,
• Registered the parameters of $\text{(3)}$ as input variables, and
• Sorted out the evaluation order.

The composed model can be evaluated in the same way as regular models:

C++Python

• #include "neml2/models/Model.h"
  #include "neml2/tensors/SR2.h"
   
  int
  main()
  {
    using namespace neml2;
    set_default_dtype(kFloat64);
    auto & eq = load_model("input_composed.i", "eq");
   
    // Create the input variables
    auto a_name = VariableName("forces", "a");
    auto b_name = VariableName("state", "b");
    auto a = SR2::fill(0.1, 0.05, -0.03, 0.02, 0.06, 0.03);
    auto b = SR2::fill(100, 20, 10, 5, -30, -20);
   
    // Evaluate the composed model
    auto b_rate_name = VariableName("state", "b_rate");
    auto b_rate = eq.value({{a_name, a}, {b_name, b}})[b_rate_name];
   
    std::cout << "b_rate: \n" << b_rate << std::endl;
  }

  Output:

  b_rate: 
   22.6184
    7.7092
   -1.9855
    3.8519
    3.9188
    1.1109
  [ CPUDoubleType{6} ]

• import neml2
  from neml2.tensors import SR2
  import torch
   
  torch.set_default_dtype(torch.double)
  eq = neml2.load_model("input_composed.i", "eq")
   
  # Create the input variables
  a = SR2.fill(0.1, 0.05, -0.03, 0.02, 0.06, 0.03)
  b = SR2.fill(100, 20, 10, 5, -30, -20)
   
  # Evaluate the composed model
  b_rate = eq.value({"forces/a": a, "state/b": b})["state/b_rate"]
   
  print("b_rate:")
  print(b_rate)

  Output:

  b_rate:
   22.6184
    7.7092
   -1.9855
    3.8519
    3.9188
    1.1109
  [ CPUDoubleType{6} ]
  <Tensor of shape [][6]>

Previous	Next
Vectorization	Model parameters (revisited)