Add Composite Power Law Creep for multiphases

 - Add PowerLawCreepStressUpdate that takes into account multiple different values based on weights so different regions in the mesh can have different plasticity values.
 - Designed it so it is similar to Composite Elasticity Tensor
 - Two tests file for CSVDIFF added to:
    1) results with the same materials but two phases give the same result as power_law_creep.i
    2) creep computation with two different material properties work
 - Three tests file for RunException to test the vector length of the input parameters.

closes idaholab#28368

modules/solid_mechanics/doc/content/source/materials/CompositePowerLawCreepStressUpdate.md

@@ -0,0 +1,47 @@
+# Composite Power Law Creep Stress Update
+
+!syntax description /Materials/CompositePowerLawCreepStressUpdate
+
+## Description
+
+!include modules/solid_mechanics/common/supplementalRadialReturnStressUpdate.md
+
+!include source/materials/PowerLawCreepStressUpdate.md
+
+The increment of inelastic strain is computed from the creep rate in this class.
+
+\begin{equation}
+  \label{eq:composite_power_law_creep}
+  \dot{\epsilon} = \sum_{i=1}^N h_i \left[ A_i \left( \sigma^{trial}_{effective} - 3 G \Delta p \right)^{n_i} \exp \left( \frac{-Q_i}{RT} \right) \right] \left(t - t_o \right)^m
+\end{equation}
+
+where subscript $i$ denotes phase specific material property, $h$ is material properties to interpolate different phases, $A_i$ is the power law creep coefficient, also known as Dorn's Constant, $\sigma^{trial}_{effective}$ is the scalar von Mises trial stress, $G$ is
+the isotropic shear modulus, $Q$ is the activation energy, $R$ is the universal
+gas constant, $T$ is the temperature, $t$ and $t_o$ are the current and initial
+times, respectively, and $n$ and $m$ are exponent values.
+
+This class calculates an effective trial stress, an effective creep strain rate
+increment and the derivative of the creep strain rate, and an effective scalar
+inelastic strain increment based on different contributions from different phases; these values are passed to the
+[ADRadialReturnStressUpdate](/ADRadialReturnStressUpdate.md) to compute the radial
+return stress increment. This isotropic plasticity class also computes the
+plastic strain as a stateful material property.
+
+This class is based on the implicit integration algorithm in
+[!cite](dunne2005introduction) pg. 146 - 149.
+
+`CompositePowerLawCreepStressUpdate` must be run in conjunction with an inelastic
+strain return mapping stress calculator such as
+[ADComputeMultipleInelasticStress](ADComputeMultipleInelasticStress.md)
+
+## Example Input File Syntax
+
+!listing modules/combined/test/tests/power_law_creep/composite_power_law_creep.i block=Materials/power_law_creep
+
+!syntax parameters /Materials/CompositePowerLawCreepStressUpdate
+
+!syntax inputs /Materials/CompositePowerLawCreepStressUpdate
+
+!syntax children /Materials/CompositePowerLawCreepStressUpdate
+
+!bibtex bibliography

modules/solid_mechanics/include/materials/CompositePowerLawCreepStressUpdate.h

@@ -0,0 +1,106 @@
+//* This file is part of the MOOSE framework
+//* https://www.mooseframework.org
+//*
+//* All rights reserved, see COPYRIGHT for full restrictions
+//* https://github.com/idaholab/moose/blob/master/COPYRIGHT
+//*
+//* Licensed under LGPL 2.1, please see LICENSE for details
+//* https://www.gnu.org/licenses/lgpl-2.1.html
+
+#pragma once
+
+#include "RadialReturnCreepStressUpdateBase.h"
+
+/**
+ * This class uses the stress update material in a radial return isotropic creep
+ * model.  This class is one of the basic radial return constitutive models; more complex
+ * constitutive models combine creep and plasticity.
+ *
+ * This class inherits from RadialReturnCreepStressUpdateBase and must be used
+ * in conjunction with ComputeMultipleInelasticStress.  This class calculates
+ * creep based on stress, temperature, and time effects.  This class also
+ * computes the creep strain as a stateful material property. This class extends
+ * from PowerLawCreepStressUpdate to include multiple different phases.
+ */
+template <bool is_ad>
+class CompositePowerLawCreepStressUpdateTempl : public RadialReturnCreepStressUpdateBaseTempl<is_ad>
+{
+public:
+  static InputParameters validParams();
+
+  CompositePowerLawCreepStressUpdateTempl(const InputParameters & parameters);
+
+  virtual Real computeStrainEnergyRateDensity(
+      const GenericMaterialProperty<RankTwoTensor, is_ad> & stress,
+      const GenericMaterialProperty<RankTwoTensor, is_ad> & strain_rate) override;
+
+  virtual bool substeppingCapabilityEnabled() override;
+
+  virtual void resetIncrementalMaterialProperties() override;
+
+  virtual void
+  computeStressInitialize(const GenericReal<is_ad> & effective_trial_stress,
+                          const GenericRankFourTensor<is_ad> & elasticity_tensor) override;
+  virtual GenericReal<is_ad> computeResidual(const GenericReal<is_ad> & effective_trial_stress,
+                                             const GenericReal<is_ad> & scalar) override
+  {
+    return computeResidualInternal<GenericReal<is_ad>>(effective_trial_stress, scalar);
+  }
+  virtual GenericReal<is_ad> computeDerivative(const GenericReal<is_ad> & effective_trial_stress,
+                                               const GenericReal<is_ad> & scalar) override;
+  virtual void
+  computeStressFinalize(const GenericRankTwoTensor<is_ad> & plastic_strain_increment) override;
+
+protected:
+  virtual GenericChainedReal<is_ad>
+  computeResidualAndDerivative(const GenericReal<is_ad> & effective_trial_stress,
+                               const GenericChainedReal<is_ad> & scalar) override
+  {
+    return computeResidualInternal<GenericChainedReal<is_ad>>(effective_trial_stress, scalar);
+  }
+
+  /// Temperature variable value
+  const GenericVariableValue<is_ad> * const _temperature;
+
+  /// Leading coefficient
+  std::vector<Real> _coefficient;
+
+  /// Exponent on the effective stress
+  std::vector<Real> _n_exponent;
+
+  /// Exponent on time
+  const Real _m_exponent;
+
+  /// Activation energy for exp term
+  std::vector<Real> _activation_energy;
+
+  /// Gas constant for exp term
+  const Real _gas_constant;
+
+  /// Simulation start time
+  const Real _start_time;
+
+  /// Exponential calculated from current time
+  Real _exp_time;
+
+  /// vector to keep the material property name for switching function material
+  const std::vector<MaterialPropertyName> _switchingFuncNames;
+  unsigned int _num_materials;
+
+  /// switching functions for each phase
+  std::vector<const GenericMaterialProperty<Real, is_ad> *> _switchingFunc;
+
+  usingTransientInterfaceMembers;
+  using RadialReturnCreepStressUpdateBaseTempl<is_ad>::_qp;
+  using RadialReturnCreepStressUpdateBaseTempl<is_ad>::_three_shear_modulus;
+  using RadialReturnCreepStressUpdateBaseTempl<is_ad>::_creep_strain;
+  using RadialReturnCreepStressUpdateBaseTempl<is_ad>::_creep_strain_old;
+
+private:
+  template <typename ScalarType>
+  ScalarType computeResidualInternal(const GenericReal<is_ad> & effective_trial_stress,
+                                     const ScalarType & scalar);
+};
+
+typedef CompositePowerLawCreepStressUpdateTempl<false> CompositePowerLawCreepStressUpdate;
+typedef CompositePowerLawCreepStressUpdateTempl<true> ADCompositePowerLawCreepStressUpdate;

modules/solid_mechanics/src/materials/CompositePowerLawCreepStressUpdate.C

@@ -0,0 +1,189 @@
+//* This file is part of the MOOSE framework
+//* https://www.mooseframework.org
+//*
+//* All rights reserved, see COPYRIGHT for full restrictions
+//* https://github.com/idaholab/moose/blob/master/COPYRIGHT
+//*
+//* Licensed under LGPL 2.1, please see LICENSE for details
+//* https://www.gnu.org/licenses/lgpl-2.1.html
+
+#include "CompositePowerLawCreepStressUpdate.h"
+#include <cmath>
+
+registerMooseObject("SolidMechanicsApp", CompositePowerLawCreepStressUpdate);
+registerMooseObject("SolidMechanicsApp", ADCompositePowerLawCreepStressUpdate);
+
+template <bool is_ad>
+InputParameters
+CompositePowerLawCreepStressUpdateTempl<is_ad>::validParams()
+{
+  InputParameters params = RadialReturnCreepStressUpdateBaseTempl<is_ad>::validParams();
+  params.addClassDescription(
+      "This class uses the stress update material in a radial return isotropic power law creep "
+      "model. This class can be used in conjunction with other creep and plasticity materials "
+      "for more complex simulations. This class is an extension to include multi-phase "
+      "capability.");
+
+  // Linear strain hardening parameters
+  params.addRequiredCoupledVar("temperature", "Coupled temperature");
+  params.addRequiredParam<std::vector<Real>>(
+      "coefficient",
+      "a vector of leading coefficient / Dorn Constant in power-law equation for each material.");
+  params.addRequiredParam<std::vector<Real>>(
+      "n_exponent",
+      "a vector of Exponent on effective stress in power-law equation for each material");
+  params.addParam<Real>("m_exponent", 0.0, "Exponent on time in power-law equation");
+  params.addRequiredParam<std::vector<Real>>("activation_energy",
+                                             "a vector of Activation energy for Arrhenius-type "
+                                             "equation of the Dorn Constant for each material");
+  params.addParam<Real>("gas_constant", 8.3143, "Universal gas constant");
+  params.addParam<Real>("start_time", 0.0, "Start time (if not zero)");
+  params.addRequiredParam<std::vector<MaterialPropertyName>>(
+      "switching_functions", "a vector of switching functions for each material");
+  return params;
+}
+
+template <bool is_ad>
+CompositePowerLawCreepStressUpdateTempl<is_ad>::CompositePowerLawCreepStressUpdateTempl(
+    const InputParameters & parameters)
+  : RadialReturnCreepStressUpdateBaseTempl<is_ad>(parameters),
+    _temperature(this->isParamValid("temperature")
+                     ? &this->template coupledGenericValue<is_ad>("temperature")
+                     : nullptr),
+    _coefficient(this->template getParam<std::vector<Real>>("coefficient")),
+    _n_exponent(this->template getParam<std::vector<Real>>("n_exponent")),
+    _m_exponent(this->template getParam<Real>("m_exponent")),
+    _activation_energy(this->template getParam<std::vector<Real>>("activation_energy")),
+    _gas_constant(this->template getParam<Real>("gas_constant")),
+    _start_time(this->template getParam<Real>("start_time")),
+    _switchingFuncNames(
+        this->template getParam<std::vector<MaterialPropertyName>>("switching_functions"))
+
+{
+  _num_materials = _switchingFuncNames.size();
+  if (_n_exponent.size() != _num_materials)
+    this->paramError("n_exponent", "n exponent must be equal to the number of switching functions");
+
+  if (_coefficient.size() != _num_materials)
+    this->paramError("coefficient",
+                     "number of Dorn constant must be equal to the number of switching functions");
+
+  if (_activation_energy.size() != _num_materials)
+    this->paramError("activation_energy",
+                     "activation energy must be equal to the number of swithing functions");
+
+  if (_start_time < this->_app.getStartTime() && (std::trunc(_m_exponent) != _m_exponent))
+    this->paramError("start_time",
+                     "Start time must be equal to or greater than the Executioner start_time if a "
+                     "non-integer m_exponent is used");
+  _switchingFunc.resize(_num_materials);
+  // set switching functions material properties for each phase
+  for (unsigned int i = 0; i < _num_materials; ++i)
+  {
+    _switchingFunc[i] =
+        &this->template getGenericMaterialProperty<Real, is_ad>(_switchingFuncNames[i]);
+  }
+}
+
+template <bool is_ad>
+void
+CompositePowerLawCreepStressUpdateTempl<is_ad>::computeStressInitialize(
+    const GenericReal<is_ad> & effective_trial_stress,
+    const GenericRankFourTensor<is_ad> & elasticity_tensor)
+{
+  RadialReturnStressUpdateTempl<is_ad>::computeStressInitialize(effective_trial_stress,
+                                                                elasticity_tensor);
+  _exp_time = std::pow(_t - _start_time, _m_exponent);
+}
+
+template <bool is_ad>
+template <typename ScalarType>
+ScalarType
+CompositePowerLawCreepStressUpdateTempl<is_ad>::computeResidualInternal(
+    const GenericReal<is_ad> & effective_trial_stress, const ScalarType & scalar)
+{
+  const ScalarType stress_delta = effective_trial_stress - _three_shear_modulus * scalar;
+  ScalarType creep_rate =
+      _coefficient[0] * std::pow(stress_delta, _n_exponent[0]) * (*_switchingFunc[0])[_qp] *
+      std::exp(-_activation_energy[0] / (_gas_constant * (*_temperature)[_qp])) * _exp_time;
+  for (unsigned int n = 1; n < _num_materials; ++n)
+  {
+    creep_rate +=
+        _coefficient[n] * std::pow(stress_delta, _n_exponent[n]) * (*_switchingFunc[n])[_qp] *
+        std::exp(-_activation_energy[n] / (_gas_constant * (*_temperature)[_qp])) * _exp_time;
+  }
+  return creep_rate * _dt - scalar;
+}
+
+template <bool is_ad>
+GenericReal<is_ad>
+CompositePowerLawCreepStressUpdateTempl<is_ad>::computeDerivative(
+    const GenericReal<is_ad> & effective_trial_stress, const GenericReal<is_ad> & scalar)
+{
+  const GenericReal<is_ad> stress_delta = effective_trial_stress - _three_shear_modulus * scalar;
+  GenericReal<is_ad> creep_rate_derivative =
+      -_coefficient[0] * _three_shear_modulus * _n_exponent[0] *
+      std::pow(stress_delta, _n_exponent[0] - 1.0) * (*_switchingFunc[0])[_qp] *
+      std::exp(-_activation_energy[0] / (_gas_constant * (*_temperature)[_qp])) * _exp_time;
+
+  for (unsigned int n = 1; n < _num_materials; ++n)
+  {
+    creep_rate_derivative +=
+        -_coefficient[n] * _three_shear_modulus * _n_exponent[n] *
+        std::pow(stress_delta, _n_exponent[n] - 1.0) * (*_switchingFunc[n])[_qp] *
+        std::exp(-_activation_energy[n] / (_gas_constant * (*_temperature)[_qp])) * _exp_time;
+  }
+  return creep_rate_derivative * _dt - 1.0;
+}
+
+template <bool is_ad>
+Real
+CompositePowerLawCreepStressUpdateTempl<is_ad>::computeStrainEnergyRateDensity(
+    const GenericMaterialProperty<RankTwoTensor, is_ad> & stress,
+    const GenericMaterialProperty<RankTwoTensor, is_ad> & strain_rate)
+{
+  GenericReal<is_ad> interpolated_exponent = 0.0;
+  for (unsigned int n = 0; n < _num_materials; ++n)
+  {
+    interpolated_exponent += (_n_exponent[n] / (_n_exponent[n] + 1)) * (*_switchingFunc[n])[_qp];
+  }
+  return MetaPhysicL::raw_value(interpolated_exponent *
+                                stress[_qp].doubleContraction((strain_rate)[_qp]));
+}
+
+template <bool is_ad>
+void
+CompositePowerLawCreepStressUpdateTempl<is_ad>::computeStressFinalize(
+    const GenericRankTwoTensor<is_ad> & plastic_strain_increment)
+{
+  _creep_strain[_qp] += plastic_strain_increment;
+}
+
+template <bool is_ad>
+void
+CompositePowerLawCreepStressUpdateTempl<is_ad>::resetIncrementalMaterialProperties()
+{
+  _creep_strain[_qp] = _creep_strain_old[_qp];
+}
+
+template <bool is_ad>
+bool
+CompositePowerLawCreepStressUpdateTempl<is_ad>::substeppingCapabilityEnabled()
+{
+  return this->_use_substepping != RadialReturnStressUpdateTempl<is_ad>::SubsteppingType::NONE;
+}
+
+template class CompositePowerLawCreepStressUpdateTempl<false>;
+template class CompositePowerLawCreepStressUpdateTempl<true>;
+template Real
+CompositePowerLawCreepStressUpdateTempl<false>::computeResidualInternal<Real>(const Real &,
+                                                                              const Real &);
+template ADReal
+CompositePowerLawCreepStressUpdateTempl<true>::computeResidualInternal<ADReal>(const ADReal &,
+                                                                               const ADReal &);
+template ChainedReal
+CompositePowerLawCreepStressUpdateTempl<false>::computeResidualInternal<ChainedReal>(
+    const Real &, const ChainedReal &);
+template ChainedADReal
+CompositePowerLawCreepStressUpdateTempl<true>::computeResidualInternal<ChainedADReal>(
+    const ADReal &, const ChainedADReal &);