Comparison between MOOSE and pyzag runs of NEML2 viscoplastic model

[fusmatrs]

Hello all,

I am trying to use NEML2 as the constitutive model in a project. As part of this, I am trying to compare a single element MOOSE model run using a NEML2 material model. For simplicity I am using the viscoplasticity_isoharden.i MOOSE input (modified to also output strains & stresses) and the corresponding NEML2 model from the MOOSE tests.

After MOOSE has solved, I am taking the strains (converted to Mandel notation) and using those to drive the same material model in python using pyzag.
I was hoping that there would be agreement between the stresses given by MOOSE and those from pyzag.
For elastic models this seems to work perfectly, the pyzag and MOOSE outputs agree. However, when I have a viscoplastic model there is a discrepancy.

I assumed that the MOOSE and pyzag outputs should always agree as they are both derived from the same NEML2 model. Unlike MOOSE, the pyzag run of the model gives stresses in the Y and Z directions, and lowers the stress in the X direction (see attached image). The trace of the stress tensors from MOOSE and pyzag are identical.

Assuming these two approaches should agree, is there something I am not understanding?

Thanks

[fusmatrs]

Slight correction - I noticed I was doing something odd with my times, which when fixed gives a better result:

However, if we make the loading a little more complex, but applying a vertical rather than horizontal displacement (to generate some shear), the problem becomes worse:

[hugary1995]

Sorry for the late reply -- I forgot to set up notification for this discussion forum. I just did.

My first guess would be that the MOOSE simulation might be using large deformation kinematics. But I just checked the MOOSE input file, and it seems like it's using small strain. So that's not it.

The fact that the uniaxial loading case yields a pretty good agreement, whereas the multiaxial loading does not, leads me to think that the discrepancy might be stemming from how you define the strain in MOOSE. I have two suggestions:

1. Compare the strain in MOOSE and the strain you passed into the pyzag model.
2. Try turning off volumetric_locking_correction in the MOOSE model. We use the F-bar approach in MOOSE, and so the stabilized strain tensor might not exactly match what you expect it to be.

[fusmatrs]

No problem, thanks for getting back to me.

1. If I increase the number of elements in the model, there's good agreement except on the faces with shear. Is there a rotation applied to the strains before MOOSE sends them to NEML2? From global to local coordinates?
2. The volumetric_locking_correction doesn't seem to make a difference.

[hugary1995]

There's no special treatment to rotation in small strain. Could you elaborate on the boundary conditions you applied? I can try reproducing it. Alternatively, just post the [BCs] block here.

You shouldn't need to increase the number of elements to see agreement. In fact, they should only match if and only if the MOOSE boundary conditions match the pyzag input strain, and if you have one single element with one quadrature point. When you use the default quadrature order, i.e., 8 quadrature points for HEX8, and don't apply Dirichlet BC on all nodes, there will be slight differences as you are technically solving for a different equation in MOOSE, i.e., div(stress) = 0, which is not applicable in pyzag.

[fusmatrs]

Ok, I may have misunderstood how both MOOSE and pyzag are working then. I had assumed that given the complete strain tensor history at a node, pyzag would be able to reproduce the same stress tensor history as MOOSE.

The BCs are here:

[BCs]
[xfix]
type = DirichletBC
variable = disp_x
boundary = left
value = 0
[]
[yfix]
type = DirichletBC
variable = disp_y
boundary = bottom
value = 0
[]
[zfix]
type = DirichletBC
variable = disp_z
boundary = back
value = 0
[]
[ydisp]
type = FunctionDirichletBC
variable = disp_y
boundary = right
function = t
preset = false
[]
[]

[hugary1995]

There's nothing wrong with your understanding. You just need to be careful with the boundary conditions since they need to fully constrain (with Dirichlet-type BC) all nodes of the HEX8 element. That's really the only way of making sure balance of linear momentum is not at play.

See my input file below:

Modified input file

[GlobalParams]
  displacements = 'disp_x disp_y disp_z'
[]

[Mesh]
  [gmg]
    type = GeneratedMeshGenerator
    dim = 3
    nx = 1
    ny = 1
    nz = 1
  []
[]

[Physics]
  [SolidMechanics]
    [QuasiStatic]
      [all]
        strain = SMALL
        new_system = true
        add_variables = true
        formulation = TOTAL
        volumetric_locking_correction = false
      []
    []
  []
[]

[NEML2]
  input = 'isoharden_neml2.i'
  [all]
    model = 'model'
    verbose = true
    device = 'cpu'

    moose_input_types = 'MATERIAL     MATERIAL     POSTPROCESSOR POSTPROCESSOR MATERIAL     MATERIAL'
    moose_inputs = '     neml2_strain neml2_strain time          time          neml2_stress equivalent_plastic_strain'
    neml2_inputs = '     forces/E     old_forces/E forces/t      old_forces/t  old_state/S  old_state/internal/ep'

    moose_output_types = 'MATERIAL     MATERIAL'
    moose_outputs = '     neml2_stress equivalent_plastic_strain'
    neml2_outputs = '     state/S      state/internal/ep'

    moose_derivative_types = 'MATERIAL'
    moose_derivatives = 'neml2_jacobian'
    neml2_derivatives = 'state/S forces/E'
  []
[]

[Postprocessors]
  [time]
    type = TimePostprocessor
    outputs = none
    execute_on = 'INITIAL TIMESTEP_BEGIN'
  []
[]

[Materials]
  [convert_strain]
    type = RankTwoTensorToSymmetricRankTwoTensor
    from = 'mechanical_strain'
    to = 'neml2_strain'
  []
  [stress]
    type = ComputeLagrangianObjectiveCustomSymmetricStress
    custom_small_stress = 'neml2_stress'
    custom_small_jacobian = 'neml2_jacobian'
  []
[]

[BCs]
  [xfix]
    type = DirichletBC
    variable = disp_x
    boundary = 'left right'
    value = 0
  []
  [yfix]
    type = DirichletBC
    variable = disp_y
    boundary = 'left'
    value = 0
  []
  [zfix]
    type = DirichletBC
    variable = disp_z
    boundary = 'left right'
    value = 0
  []
  [ydisp]
    type = FunctionDirichletBC
    variable = disp_y
    boundary = right
    function = t
    preset = false
  []
[]

[Executioner]
  type = Transient
  solve_type = NEWTON
  petsc_options_iname = '-pc_type'
  petsc_options_value = 'lu'
  automatic_scaling = true
  dt = 1e-3
  dtmin = 1e-3
  num_steps = 5
  residual_and_jacobian_together = true
[]

[AuxVariables]
  [strain_xx]
    order = CONSTANT
    family = MONOMIAL
    [AuxKernel]
      type = MaterialRankTwoTensorAux
      property = 'mechanical_strain'
      i = 0
      j = 0
      execute_on = 'INITIAL TIMESTEP_END'
    []
  []
  [strain_yy]
    order = CONSTANT
    family = MONOMIAL
    [AuxKernel]
      type = MaterialRankTwoTensorAux
      property = 'mechanical_strain'
      i = 1
      j = 1
      execute_on = 'INITIAL TIMESTEP_END'
    []
  []
  [strain_zz]
    order = CONSTANT
    family = MONOMIAL
    [AuxKernel]
      type = MaterialRankTwoTensorAux
      property = 'mechanical_strain'
      i = 2
      j = 2
      execute_on = 'INITIAL TIMESTEP_END'
    []
  []
  [strain_yz]
    order = CONSTANT
    family = MONOMIAL
    [AuxKernel]
      type = MaterialRankTwoTensorAux
      property = 'mechanical_strain'
      i = 1
      j = 2
      execute_on = 'INITIAL TIMESTEP_END'
    []
  []
  [strain_xz]
    order = CONSTANT
    family = MONOMIAL
    [AuxKernel]
      type = MaterialRankTwoTensorAux
      property = 'mechanical_strain'
      i = 0
      j = 2
      execute_on = 'INITIAL TIMESTEP_END'
    []
  []
  [strain_xy]
    order = CONSTANT
    family = MONOMIAL
    [AuxKernel]
      type = MaterialRankTwoTensorAux
      property = 'mechanical_strain'
      i = 0
      j = 1
      execute_on = 'INITIAL TIMESTEP_END'
    []
  []
  [stress_xx]
    order = CONSTANT
    family = MONOMIAL
    [AuxKernel]
      type = MaterialRankTwoTensorAux
      property = 'cauchy_stress'
      i = 0
      j = 0
      execute_on = 'INITIAL TIMESTEP_END'
    []
  []
  [stress_yy]
    order = CONSTANT
    family = MONOMIAL
    [AuxKernel]
      type = MaterialRankTwoTensorAux
      property = 'cauchy_stress'
      i = 1
      j = 1
      execute_on = 'INITIAL TIMESTEP_END'
    []
  []
  [stress_zz]
    order = CONSTANT
    family = MONOMIAL
    [AuxKernel]
      type = MaterialRankTwoTensorAux
      property = 'cauchy_stress'
      i = 2
      j = 2
      execute_on = 'INITIAL TIMESTEP_END'
    []
  []
  [stress_yz]
    order = CONSTANT
    family = MONOMIAL
    [AuxKernel]
      type = MaterialRankTwoTensorAux
      property = 'cauchy_stress'
      i = 1
      j = 2
      execute_on = 'INITIAL TIMESTEP_END'
    []
  []
  [stress_xz]
    order = CONSTANT
    family = MONOMIAL
    [AuxKernel]
      type = MaterialRankTwoTensorAux
      property = 'cauchy_stress'
      i = 0
      j = 2
      execute_on = 'INITIAL TIMESTEP_END'
    []
  []
  [stress_xy]
    order = CONSTANT
    family = MONOMIAL
    [AuxKernel]
      type = MaterialRankTwoTensorAux
      property = 'cauchy_stress'
      i = 0
      j = 1
      execute_on = 'INITIAL TIMESTEP_END'
    []
  []
[]

[Postprocessors]
  [strain_xx]
    type = ElementAverageValue
    variable = 'strain_xx'
    execute_on = 'INITIAL TIMESTEP_END'
  []
  [strain_yy]
    type = ElementAverageValue
    variable = 'strain_yy'
    execute_on = 'INITIAL TIMESTEP_END'
  []
  [strain_zz]
    type = ElementAverageValue
    variable = 'strain_zz'
    execute_on = 'INITIAL TIMESTEP_END'
  []
  [strain_yz]
    type = ElementAverageValue
    variable = 'strain_yz'
    execute_on = 'INITIAL TIMESTEP_END'
  []
  [strain_xz]
    type = ElementAverageValue
    variable = 'strain_xz'
    execute_on = 'INITIAL TIMESTEP_END'
  []
  [strain_xy]
    type = ElementAverageValue
    variable = 'strain_xy'
    execute_on = 'INITIAL TIMESTEP_END'
  []
  [stress_xx]
    type = ElementAverageValue
    variable = 'stress_xx'
    execute_on = 'INITIAL TIMESTEP_END'
  []
  [stress_yy]
    type = ElementAverageValue
    variable = 'stress_yy'
    execute_on = 'INITIAL TIMESTEP_END'
  []
  [stress_zz]
    type = ElementAverageValue
    variable = 'stress_zz'
    execute_on = 'INITIAL TIMESTEP_END'
  []
  [stress_yz]
    type = ElementAverageValue
    variable = 'stress_yz'
    execute_on = 'INITIAL TIMESTEP_END'
  []
  [stress_xz]
    type = ElementAverageValue
    variable = 'stress_xz'
    execute_on = 'INITIAL TIMESTEP_END'
  []
  [stress_xy]
    type = ElementAverageValue
    variable = 'stress_xy'
    execute_on = 'INITIAL TIMESTEP_END'
  []
[]

[Outputs]
  csv = true
[]

I saved the results into a csv:

Results

time,strain_xx,strain_xy,strain_xz,strain_yy,strain_yz,strain_zz,stress_xx,stress_xy,stress_xz,stress_yy,stress_yz,stress_zz
0,0,0,0,0,0,0,0,0,0,0,0,0
0.001,0,0.0005,0,0,2.1175823681358e-20,0,0,30.317859510235,0,0,1.2840112947698e-15,0
0.002,0,0.001,0,0,4.2351647362715e-20,0,0,47.436202229863,0,0,2.0090013090656e-15,0
0.003,0,0.0015,0,0,-6.098637220231e-20,0,0,55.907575170376,0,0,-3.8661704691976e-15,0
0.004,0,0.002,0,0,8.470329472543e-20,0,0,59.949680136601,0,0,4.6633291115342e-15,0
0.005,0,0.0025,0,0,1.0842021724855e-19,0,0,61.927815102818,0,0,4.08255860225e-15,0

The MOOSE strain tensor is then passed to pyzag to calculate the stress. See below for an example python script:

Example pyzag script

import pandas as pd
import torch
from pyzag import nonlinear, chunktime
import neml2
from neml2.tensors import Scalar, SR2, Tensor

torch.set_default_dtype(torch.float64)

# moose results
data = pd.read_csv("isoharden_out.csv")

# model
neml2_model = neml2.load_model("isoharden_neml2.i", "implicit_rate")
pyzag_model = neml2.pyzag.NEML2PyzagModel(neml2_model)

# times
times = torch.tensor(data["time"]).unsqueeze(-1)

# strain
exx = Scalar(torch.tensor(data["strain_xx"]))
eyy = Scalar(torch.tensor(data["strain_yy"]))
ezz = Scalar(torch.tensor(data["strain_zz"]))
eyz = Scalar(torch.tensor(data["strain_yz"]))
exz = Scalar(torch.tensor(data["strain_xz"]))
exy = Scalar(torch.tensor(data["strain_xy"]))
e = SR2.fill(exx, eyy, ezz, eyz, exz, exy).torch()

# pyzag solver
solver = nonlinear.RecursiveNonlinearEquationSolver(
    pyzag_model,
    step_generator=nonlinear.StepGenerator(block_size=1),
    predictor=nonlinear.PreviousStepsPredictor(),
    nonlinear_solver=chunktime.ChunkNewtonRaphson(rtol=1.0e-8, atol=1.0e-10),
)

# solve
forces = pyzag_model.forces_asm.assemble_by_variable({"forces/E": e, "forces/t": times})
state0 = torch.zeros(pyzag_model.nstate)
with torch.no_grad():
    results = nonlinear.solve(solver, state0, 6, forces.torch())


# stress calculated by moose
sxx = Scalar(torch.tensor(data["stress_xx"]))
syy = Scalar(torch.tensor(data["stress_yy"]))
szz = Scalar(torch.tensor(data["stress_zz"]))
syz = Scalar(torch.tensor(data["stress_yz"]))
sxz = Scalar(torch.tensor(data["stress_xz"]))
sxy = Scalar(torch.tensor(data["stress_xy"]))
s_moose = SR2.fill(sxx, syy, szz, syz, sxz, sxy)
print(s_moose)

# stress calculated by pyzag
states = pyzag_model.state_asm.split_by_variable(Tensor(results, 1))
s_pyzag = states["state/S"]
print(s_pyzag)

The MOOSE and pyzag results are identical.

[fusmatrs]

Thanks, that's really helpful.
Does this mean that for more complex geometry models I shouldn't expect the stress calculated at nodes / elements away from the BCs to agree between the two?
My hope was to use NEML2 constitutive models to predict stresses from a measured strain history over a surface. The MOOSE part was to generate some ground truth to assess the approach.