Crystal Plasticity preconditioner and executioner

[p-antonioni]

Hi!

I am performing some Crystal Plasticity simulations and I would have to run a simulation with around some milions of dofs.
Even for my test case, I saw that full Newton solver is faster than PJFNK (with Newton I have a quadratic convergence, with PJFNK superlinear). Maybe because I have not set it properly.
Do you have some suggestions or experience with some solver/preconditioner to have good performances on such simulation? I mean for Newton and PJFNK

Thank you,
Paolo

The set up I used for PJFNK:

[Preconditioning]
    [smp]
    type = SMP
    full = true
    []
[]

[Executioner]
  type = Transient
  solve_type = 'PJFNK'
  petsc_options_iname = '-pc_type -pc_asm_overlap -sub_pc_type -ksp_type -ksp_gmres_restart'
  petsc_options_value = ' asm      2              lu            gmres     200'
  nl_abs_tol = 1e-4
 ...
[]

NEWTON:

[Preconditioning]
  [smp]
    type = SMP
    full = true
    petsc_options_iname = '-pc_type -pc_factor_mat_solver_package '
    petsc_options_value = ' lu       mumps'
  []
[]

[Executioner]
  type = Transient
  petsc_options_iname = '-pc_type -pc_asm_overlap -sub_pc_type -ksp_type -ksp_gmres_restart'
  petsc_options_value = ' asm      2              lu            gmres     200'
  solve_type = 'NEWTON'
  use_pre_SMO_residual = false

[GiudGiud]

@bwspenc

[sapitts]

You are correct, the current crystal plasticity classes do not use AD.

You might consider running the HCP specific class, CrystalPlasticityHCPDislocationSlipBeyerleinUpdate for comparison with your modified Kalidindi model

[sapitts]

The crystal plasticity models in MOOSE are unfortunately slow. In addition to the solver tolerances there are constitutive model tolerances one can adjust. See the discussion here for more information: #28534

[p-antonioni]

Hi @sapitts,
Thanks for answering. I was looking at the CrystalPlasticityHCPDislocationSlipBayerleinUpdate, with the modifications I did, they look quite similar in some parts. Could I ask you why the equivalent slip increment is calculated in this way? It is like assuming that the slip system behavior does not depend on the current hardening state ?

  for (const auto i : make_range(_number_slip_systems))
  {
    if (MooseUtils::absoluteFuzzyEqual(_tau[_qp][i], 0.0))
      dslip_dtau[i] = 0.0;
    else
      dslip_dtau[i] = _slip_increment[_qp][i] /
                      (_rate_sensitivity_exponent * std::abs(_tau[_qp][i])) * _substep_dt;
  }
}

I would expect the derivative with respect to tau of the slip increment, something like this:

  for (const auto i : make_range(_number_slip_systems))
  {
    if (MooseUtils::absoluteFuzzyEqual(_tau[_qp][i], 0.0))
      dslip_dtau[i] = 0.0;
    else
      dslip_dtau[i] = _ao / _xm *
                      std::pow(std::abs(_tau[_qp][i] / _slip_resistance[_qp][i]), 1.0 / _xm - 1.0) /
                      _slip_resistance[_qp][i];
  }

In general, there are some preconditioning/executioner settings that you found helpful when running on many dofs?

[sapitts]

Hi @p-antonioni,
In case you haven't seen this, the HCP version of the derivative simply uses the definition of the slip increment to reduce the computational load. I've written out the intermediate steps in the image below, starting with the derivative as you copied from the Kalidindi model and ending with the HCP model, excluding the time step at the end

I've had the most luck with using the SMP preconditioner and asm petsc option as in the regression tests. I've heard others have had better luck with using hyper boomerang. Generally I find I get more performance improvements by changing the constitutive model tolerances, since most of the compute time for my runs is in the material class.

Hope this helps,

[p-antonioni]

Hi @sapitts,
Thanks a lot for the hints and the explanation.
In my case, using this approximation and same settings slows down the convergence when plasticity arises (I have two different phases).
I will try with Hyper boomerang and see how it goes.