diff --git a/modules/solid_mechanics/doc/content/source/materials/ComputeMultipleInelasticStressSmallStrain.md b/modules/solid_mechanics/doc/content/source/materials/ComputeMultipleInelasticStressSmallStrain.md
new file mode 100644
index 000000000000..2f18d4af0638
--- /dev/null
+++ b/modules/solid_mechanics/doc/content/source/materials/ComputeMultipleInelasticStressSmallStrain.md
@@ -0,0 +1,72 @@
+# ComputeMultipleInelasticStressSmallStrain
+
+!syntax description /Materials/ComputeMultipleInelasticStressSmallStrain
+
+## Description
+
+`ComputeMultipleInelasticStressSmallStrain` computes the stress and decomposition of strain into elastic and inelastic components for materials with multiple inelastic models (e.g., plasticity, creep) using a **small strain formulation**. This class is analogous to [ComputeMultipleInelasticStress](ComputeMultipleInelasticStress.md) but uses total strains instead of strain increments and does not include finite strain rotation effects.
+
+### Small Strain vs Finite Strain
+
+The key differences between this class and `ComputeMultipleInelasticStress` are:
+
+- **Strain Formulation**: Uses total small strains $\boldsymbol{\varepsilon} = \frac{1}{2}(\nabla \mathbf{u} + \nabla \mathbf{u}^T)$ instead of logarithmic strain increments
+- **No Rotations**: Does not perform finite strain rotations or use rotation increment tensors
+- **Parent Class**: Inherits from [ComputeStressBase](ComputeStressBase.md) instead of [ComputeFiniteStrainElasticStress](ComputeFiniteStrainElasticStress.md)
+- **Total Formulation**: Works with total strains rather than incremental formulation
+
+### Iterative Solution Procedure
+
+The material uses an iterative return mapping algorithm to find the admissible stress state that satisfies all inelastic models simultaneously:
+
+1. For each iteration:
+   - Loop over all inelastic models
+   - For each model, compute the current elastic strain by subtracting inelastic strains from other models
+   - Form trial stress: $\boldsymbol{\sigma}^{trial} = \mathbf{C} : \boldsymbol{\varepsilon}^{elastic}$
+   - Call the model's `updateState` method to compute admissible stress and inelastic strain increment
+   - Track min/max stress for convergence checking
+
+2. Check convergence based on the L2 norm of stress differences between models
+3. If not converged and iterations remain, repeat
+
+4. Once converged:
+   - Combine inelastic strains from all models using specified weights
+   - Compute final elastic strain: $\boldsymbol{\varepsilon}^{elastic} = \boldsymbol{\varepsilon}^{total} - \boldsymbol{\varepsilon}^{inelastic}$
+   - Compute Jacobian (tangent operator) if needed
+
+### Strain Decomposition
+
+The mechanical strain is decomposed as:
+
+\begin{equation}
+\boldsymbol{\varepsilon}^{mechanical} = \boldsymbol{\varepsilon}^{elastic} + \boldsymbol{\varepsilon}^{inelastic}
+\end{equation}
+
+where the combined inelastic strain is a weighted sum:
+
+\begin{equation}
+\boldsymbol{\varepsilon}^{inelastic} = \sum_i w_i \boldsymbol{\varepsilon}^{inelastic}_i
+\end{equation}
+
+## Usage Notes
+
+- This material should be paired with a small strain calculator such as [ComputeSmallStrain](ComputeSmallStrain.md)
+- Do NOT use with incremental strain formulations (e.g., `ComputeFiniteStrain`)
+- Inelastic models should be listed with creep models first, plasticity models last
+- The `perform_finite_strain_rotations` parameter from the finite strain version is not needed here
+
+## Example Input File Syntax
+
+!listing modules/solid_mechanics/test/tests/multiple_inelastic_stress_small_strain/two_models.i block=Materials/stress
+
+!syntax parameters /Materials/ComputeMultipleInelasticStressSmallStrain
+
+!syntax inputs /Materials/ComputeMultipleInelasticStressSmallStrain
+
+!syntax children /Materials/ComputeMultipleInelasticStressSmallStrain
+
+## See Also
+
+- [ComputeMultipleInelasticStress](ComputeMultipleInelasticStress.md) - Finite strain version
+- [ComputeSmallStrain](ComputeSmallStrain.md) - Small strain calculator to pair with this class
+- [StressUpdateBase](StressUpdateBase.md) - Base class for inelastic models
diff --git a/modules/solid_mechanics/include/materials/ComputeMultipleInelasticStressSmallStrain.h b/modules/solid_mechanics/include/materials/ComputeMultipleInelasticStressSmallStrain.h
new file mode 100644
index 000000000000..379edd69d060
--- /dev/null
+++ b/modules/solid_mechanics/include/materials/ComputeMultipleInelasticStressSmallStrain.h
@@ -0,0 +1,42 @@
+//* This file is part of the MOOSE framework
+//* https://mooseframework.inl.gov
+//*
+//* 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 "ComputeMultipleInelasticStressSmallStrainBase.h"
+
+/**
+ * ComputeMultipleInelasticStressSmallStrain computes the stress, the consistent tangent
+ * operator (or an approximation to it), and a decomposition of the strain
+ * into elastic and inelastic parts using small strain formulation.
+ *
+ * The elastic strain is calculated by subtracting the computed inelastic strain
+ * from the mechanical strain tensor. Mechanical strain is considered as the sum
+ * of the elastic and inelastic (plastic, creep, etc) strains.
+ *
+ * This material is used to call the recompute iterative materials of a number
+ * of specified inelastic models that inherit from StressUpdateBase. It iterates
+ * over the specified inelastic models until the change in stress is within
+ * a user-specified tolerance, in order to produce the stress, the consistent
+ * tangent operator and the elastic and inelastic strains for the time increment.
+ */
+
+class ComputeMultipleInelasticStressSmallStrain : public ComputeMultipleInelasticStressSmallStrainBase
+{
+public:
+  static InputParameters validParams();
+
+  ComputeMultipleInelasticStressSmallStrain(const InputParameters & parameters);
+
+protected:
+  virtual std::vector<MaterialName> getInelasticModelNames() override;
+
+  virtual void updateQpState(RankTwoTensor & elastic_strain,
+                             RankTwoTensor & inelastic_strain) override;
+};
diff --git a/modules/solid_mechanics/include/materials/ComputeMultipleInelasticStressSmallStrainBase.h b/modules/solid_mechanics/include/materials/ComputeMultipleInelasticStressSmallStrainBase.h
new file mode 100644
index 000000000000..9ee290e4666a
--- /dev/null
+++ b/modules/solid_mechanics/include/materials/ComputeMultipleInelasticStressSmallStrainBase.h
@@ -0,0 +1,171 @@
+//* This file is part of the MOOSE framework
+//* https://mooseframework.inl.gov
+//*
+//* 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 "ComputeStressBase.h"
+#include "GuaranteeConsumer.h"
+#include "DamageBase.h"
+#include "StressUpdateBase.h"
+#include "MultipleInelasticStressHelper.h"
+
+/**
+ * ComputeMultipleInelasticStressSmallStrainBase computes the stress, the consistent tangent
+ * operator (or an approximation to it), and a decomposition of the strain
+ * into elastic and inelastic parts using small strain formulation.
+ *
+ * The elastic strain is calculated by subtracting the computed inelastic strain
+ * from the mechanical strain tensor. Mechanical strain is considered as the sum
+ * of the elastic and inelastic (plastic, creep, etc) strains.
+ *
+ * This material is used to call the recompute iterative materials of a number
+ * of specified inelastic models that inherit from StressUpdateBase. It iterates
+ * over the specified inelastic models until the change in stress is within
+ * a user-specified tolerance, in order to produce the stress, the consistent
+ * tangent operator and the elastic and inelastic strains for the time increment.
+ *
+ * This class uses small strain (total) formulation, in contrast to
+ * ComputeMultipleInelasticStressBase which uses finite strain (incremental) formulation.
+ */
+
+class ComputeMultipleInelasticStressSmallStrainBase : public ComputeStressBase,
+                                                      public GuaranteeConsumer
+{
+public:
+  static InputParameters validParams();
+
+  ComputeMultipleInelasticStressSmallStrainBase(const InputParameters & parameters);
+
+  virtual void initialSetup() override;
+
+protected:
+  virtual std::vector<MaterialName> getInelasticModelNames() = 0;
+
+  virtual void initQpStatefulProperties() override;
+
+  virtual void computeQpStress() override;
+
+  /**
+   * Given the current strain, iterate over all of the user-specified
+   * recompute materials in order to find an admissible stress (which is placed
+   * into _stress[_qp]) and set of inelastic strains, as well as the tangent operator
+   * (which is placed into _Jacobian_mult[_qp]).
+   * @param elastic_strain The elastic strain after the iterative process has converged
+   * @param inelastic_strain The inelastic strain after the iterative process has converged
+   */
+  virtual void updateQpState(RankTwoTensor & elastic_strain, RankTwoTensor & inelastic_strain) = 0;
+
+  /**
+   * An optimised version of updateQpState that gets used when the number
+   * of plastic models is unity, or when we're cycling through models
+   * Given the current strain, find an admissible stress (which is
+   * put into _stress[_qp]) and inelastic strain, as well as the tangent operator
+   * (which is placed into _Jacobian_mult[_qp])
+   * @param model_number Use this model number
+   * @param elastic_strain The elastic strain
+   * @param inelastic_strain The inelastic strain
+   */
+  virtual void updateQpStateSingleModel(unsigned model_number,
+                                        RankTwoTensor & elastic_strain,
+                                        RankTwoTensor & inelastic_strain);
+
+  /**
+   * Using _elasticity_tensor[_qp] and the consistent tangent operators,
+   * _consistent_tangent_operator[...] computed by the inelastic models,
+   * compute _Jacobian_mult[_qp]
+   */
+  virtual void computeQpJacobianMult();
+
+  /**
+   * Given a trial stress (_stress[_qp]) and the current mechanical strain
+   * let the model_number model produce an admissible stress (gets placed back
+   * in _stress[_qp]), and compute the elastic and inelastic strains,
+   * as well as the consistent_tangent_operator
+   * @param model_number The inelastic model to use
+   * @param elastic_strain The elastic strain (computed)
+   * @param inelastic_strain The inelastic strain (computed)
+   * @param consistent_tangent_operator The consistent tangent operator
+   */
+  virtual void computeAdmissibleState(unsigned model_number,
+                                      RankTwoTensor & elastic_strain,
+                                      RankTwoTensor & inelastic_strain,
+                                      RankFourTensor & consistent_tangent_operator);
+
+  ///@{Input parameters associated with the recompute iteration to return the stress state to the yield surface
+  const unsigned int _max_iterations;
+  const Real _relative_tolerance;
+  const Real _absolute_tolerance;
+  const bool _internal_solve_full_iteration_history;
+  ///@}
+
+  /// Name of the elasticity tensor material property
+  const std::string _elasticity_tensor_name;
+  /// Elasticity tensor material property
+  const MaterialProperty<RankFourTensor> & _elasticity_tensor;
+
+  /// Old value of mechanical strain
+  const MaterialProperty<RankTwoTensor> & _mechanical_strain_old;
+
+  /// The sum of the inelastic strains that come from the plastic models
+  MaterialProperty<RankTwoTensor> & _inelastic_strain;
+
+  /// old value of inelastic strain
+  const MaterialProperty<RankTwoTensor> & _inelastic_strain_old;
+
+  /// Old value of stress
+  const MaterialProperty<RankTwoTensor> & _stress_old;
+
+  /// what sort of Tangent operator to calculate
+  const enum class TangentOperatorEnum { elastic, nonlinear } _tangent_operator_type;
+
+  /// number of plastic models
+  unsigned _num_models;
+
+  /// Flags to compute tangent during updateState call
+  std::vector<bool> _tangent_computation_flag;
+
+  /// Calculation method for the tangent modulus
+  TangentCalculationMethod _tangent_calculation_method;
+
+  /// _inelastic_strain = sum_i (_inelastic_weights_i * inelastic_strain_from_model_i)
+  std::vector<Real> _inelastic_weights;
+
+  /// the consistent tangent operators computed by each plastic model
+  std::vector<RankFourTensor> _consistent_tangent_operator;
+
+  /// whether to cycle through the models, using only one model per timestep
+  const bool _cycle_models;
+
+  MaterialProperty<Real> & _material_timestep_limit;
+
+  /**
+   * Rank four symmetric identity tensor
+   */
+  const RankFourTensor _identity_symmetric_four;
+
+  /**
+   * The user supplied list of inelastic models to use in the simulation
+   *
+   * Users should take care to list creep models first and plasticity
+   * models last to allow for the case when a creep model relaxes the stress state
+   * inside of the yield surface in an iteration.
+   */
+  std::vector<StressUpdateBase *> _models;
+
+  /// is the elasticity tensor guaranteed to be isotropic?
+  bool _is_elasticity_tensor_guaranteed_isotropic;
+
+  /// are all inelastic models inherently isotropic? (not the case for e.g. weak plane plasticity models)
+  bool _all_models_isotropic;
+
+  /// Pointer to the damage model
+  DamageBaseTempl<false> * _damage_model;
+
+  RankTwoTensor _undamaged_stress_old;
+};
diff --git a/modules/solid_mechanics/include/materials/MultipleInelasticStressHelper.h b/modules/solid_mechanics/include/materials/MultipleInelasticStressHelper.h
new file mode 100644
index 000000000000..86281bce93db
--- /dev/null
+++ b/modules/solid_mechanics/include/materials/MultipleInelasticStressHelper.h
@@ -0,0 +1,101 @@
+//* This file is part of the MOOSE framework
+//* https://mooseframework.inl.gov
+//*
+//* 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 "RankTwoTensor.h"
+#include "RankFourTensor.h"
+#include "StressUpdateBase.h"
+
+/**
+ * Helper class containing shared logic for multiple inelastic stress calculations.
+ * This allows code reuse between finite strain and small strain formulations.
+ */
+class MultipleInelasticStressHelper
+{
+public:
+  /**
+   * Compute the Jacobian multiplier using the elasticity tensor and consistent tangent operators
+   * @param tangent_calculation_method The method for computing the tangent
+   * @param elasticity_tensor The elasticity tensor
+   * @param consistent_tangent_operator Vector of consistent tangent operators from each model
+   * @param num_models Number of inelastic models
+   * @param identity_symmetric_four The rank-4 symmetric identity tensor
+   * @return The computed Jacobian multiplier
+   */
+  static RankFourTensor computeJacobianMult(
+      TangentCalculationMethod tangent_calculation_method,
+      const RankFourTensor & elasticity_tensor,
+      const std::vector<RankFourTensor> & consistent_tangent_operator,
+      unsigned int num_models,
+      const RankFourTensor & identity_symmetric_four);
+
+  /**
+   * Compute the combined inelastic strain increment from individual model contributions
+   * @param inelastic_strain_increment Vector of inelastic strain increments from each model
+   * @param inelastic_weights Weights for combining the inelastic strains
+   * @param num_models Number of inelastic models
+   * @return The weighted sum of inelastic strain increments
+   */
+  static RankTwoTensor computeCombinedInelasticStrainIncrement(
+      const std::vector<RankTwoTensor> & inelastic_strain_increment,
+      const std::vector<Real> & inelastic_weights,
+      unsigned int num_models);
+
+  /**
+   * Compute material timestep limit from all models
+   * @param models Vector of stress update models
+   * @param num_models Number of inelastic models
+   * @return The computed timestep limit
+   */
+  static Real computeMaterialTimestepLimit(const std::vector<StressUpdateBase *> & models,
+                                           unsigned int num_models);
+
+  /**
+   * Check convergence of the stress iteration
+   * @param stress_max Maximum stress from current iteration
+   * @param stress_min Minimum stress from current iteration
+   * @param counter Current iteration number
+   * @param max_iterations Maximum allowed iterations
+   * @param relative_tolerance Relative convergence tolerance
+   * @param absolute_tolerance Absolute convergence tolerance
+   * @param first_l2norm_delta_stress L2 norm from first iteration (for relative check)
+   * @param num_models Number of models (if 1, always converged)
+   * @return true if converged, false otherwise
+   */
+  static bool checkConvergence(const RankTwoTensor & stress_max,
+                               const RankTwoTensor & stress_min,
+                               unsigned int counter,
+                               unsigned int max_iterations,
+                               Real relative_tolerance,
+                               Real absolute_tolerance,
+                               Real first_l2norm_delta_stress,
+                               unsigned int num_models);
+
+  /**
+   * Update stress min/max for convergence checking
+   * @param stress Current stress state
+   * @param stress_max Maximum stress (updated in place)
+   * @param stress_min Minimum stress (updated in place)
+   * @param is_first_model True if this is the first model in the iteration
+   */
+  static void updateStressMinMax(const RankTwoTensor & stress,
+                                 RankTwoTensor & stress_max,
+                                 RankTwoTensor & stress_min,
+                                 bool is_first_model);
+
+  /**
+   * Compute the L2 norm of the stress difference for convergence checking
+   * @param stress_max Maximum stress
+   * @param stress_min Minimum stress
+   * @return L2 norm of the difference
+   */
+  static Real computeStressDifferenceNorm(const RankTwoTensor & stress_max,
+                                          const RankTwoTensor & stress_min);
+};
diff --git a/modules/solid_mechanics/src/materials/ComputeMultipleInelasticStress.C b/modules/solid_mechanics/src/materials/ComputeMultipleInelasticStress.C
index eff1144bf1c0..b4843f9fe75e 100644
--- a/modules/solid_mechanics/src/materials/ComputeMultipleInelasticStress.C
+++ b/modules/solid_mechanics/src/materials/ComputeMultipleInelasticStress.C
@@ -12,6 +12,7 @@
 #include "StressUpdateBase.h"
 #include "MooseException.h"
 #include "DamageBase.h"
+#include "MultipleInelasticStressHelper.h"
 #include "libmesh/int_range.h"
 
 registerMooseObject("SolidMechanicsApp", ComputeMultipleInelasticStress);
@@ -95,26 +96,16 @@ ComputeMultipleInelasticStress::updateQpState(RankTwoTensor & elastic_strain_inc
                              inelastic_strain_increment[i_rmm],
                              _consistent_tangent_operator[i_rmm]);
 
-      if (i_rmm == 0)
-      {
-        stress_max = _stress[_qp];
-        stress_min = _stress[_qp];
-      }
-      else
-      {
-        for (const auto i : make_range(Moose::dim))
-          for (const auto j : make_range(Moose::dim))
-            if (_stress[_qp](i, j) > stress_max(i, j))
-              stress_max(i, j) = _stress[_qp](i, j);
-            else if (stress_min(i, j) > _stress[_qp](i, j))
-              stress_min(i, j) = _stress[_qp](i, j);
-      }
+      // Update stress min/max for convergence checking
+      MultipleInelasticStressHelper::updateStressMinMax(
+          _stress[_qp], stress_max, stress_min, i_rmm == 0);
     }
 
     // now check convergence in the stress:
     // once the change in stress is within tolerance after each recompute material
     // consider the stress to be converged
-    l2norm_delta_stress = (stress_max - stress_min).L2norm();
+    l2norm_delta_stress =
+        MultipleInelasticStressHelper::computeStressDifferenceNorm(stress_max, stress_min);
     if (counter == 0 && l2norm_delta_stress > 0.0)
       first_l2norm_delta_stress = l2norm_delta_stress;
 
@@ -138,20 +129,13 @@ ComputeMultipleInelasticStress::updateQpState(RankTwoTensor & elastic_strain_inc
       (l2norm_delta_stress / first_l2norm_delta_stress) > _relative_tolerance)
     throw MooseException("Max stress iteration hit during ComputeMultipleInelasticStress solve!");
 
-  combined_inelastic_strain_increment.zero();
-  for (const auto i_rmm : make_range(_num_models))
-    combined_inelastic_strain_increment +=
-        _inelastic_weights[i_rmm] * inelastic_strain_increment[i_rmm];
+  combined_inelastic_strain_increment =
+      MultipleInelasticStressHelper::computeCombinedInelasticStrainIncrement(
+          inelastic_strain_increment, _inelastic_weights, _num_models);
 
   if (_fe_problem.currentlyComputingJacobian())
     computeQpJacobianMult();
 
-  _material_timestep_limit[_qp] = 0.0;
-  for (const auto i_rmm : make_range(_num_models))
-    _material_timestep_limit[_qp] += 1.0 / _models[i_rmm]->computeTimeStepLimit();
-
-  if (MooseUtils::absoluteFuzzyEqual(_material_timestep_limit[_qp], 0.0))
-    _material_timestep_limit[_qp] = std::numeric_limits<Real>::max();
-  else
-    _material_timestep_limit[_qp] = 1.0 / _material_timestep_limit[_qp];
+  _material_timestep_limit[_qp] =
+      MultipleInelasticStressHelper::computeMaterialTimestepLimit(_models, _num_models);
 }
diff --git a/modules/solid_mechanics/src/materials/ComputeMultipleInelasticStressBase.C b/modules/solid_mechanics/src/materials/ComputeMultipleInelasticStressBase.C
index 1e91d74d4fb8..6a1b883c623a 100644
--- a/modules/solid_mechanics/src/materials/ComputeMultipleInelasticStressBase.C
+++ b/modules/solid_mechanics/src/materials/ComputeMultipleInelasticStressBase.C
@@ -12,6 +12,7 @@
 #include "StressUpdateBase.h"
 #include "MooseException.h"
 #include "DamageBase.h"
+#include "MultipleInelasticStressHelper.h"
 #include "libmesh/int_range.h"
 
 InputParameters
@@ -282,23 +283,12 @@ ComputeMultipleInelasticStressBase::finiteStrainRotation(const bool force_elasti
 void
 ComputeMultipleInelasticStressBase::computeQpJacobianMult()
 {
-  if (_tangent_calculation_method == TangentCalculationMethod::ELASTIC)
-    _Jacobian_mult[_qp] = _elasticity_tensor[_qp];
-  else if (_tangent_calculation_method == TangentCalculationMethod::PARTIAL)
-  {
-    RankFourTensor A = _identity_symmetric_four;
-    for (const auto i_rmm : make_range(_num_models))
-      A += _consistent_tangent_operator[i_rmm];
-    mooseAssert(A.isSymmetric(), "Tangent operator isn't symmetric");
-    _Jacobian_mult[_qp] = A.invSymm() * _elasticity_tensor[_qp];
-  }
-  else
-  {
-    const RankFourTensor E_inv = _elasticity_tensor[_qp].invSymm();
-    _Jacobian_mult[_qp] = _consistent_tangent_operator[0];
-    for (const auto i_rmm : make_range(1u, _num_models))
-      _Jacobian_mult[_qp] = _consistent_tangent_operator[i_rmm] * E_inv * _Jacobian_mult[_qp];
-  }
+  _Jacobian_mult[_qp] = MultipleInelasticStressHelper::computeJacobianMult(
+      _tangent_calculation_method,
+      _elasticity_tensor[_qp],
+      _consistent_tangent_operator,
+      _num_models,
+      _identity_symmetric_four);
 }
 
 void
diff --git a/modules/solid_mechanics/src/materials/ComputeMultipleInelasticStressSmallStrain.C b/modules/solid_mechanics/src/materials/ComputeMultipleInelasticStressSmallStrain.C
new file mode 100644
index 000000000000..a43f8419f1da
--- /dev/null
+++ b/modules/solid_mechanics/src/materials/ComputeMultipleInelasticStressSmallStrain.C
@@ -0,0 +1,146 @@
+//* This file is part of the MOOSE framework
+//* https://mooseframework.inl.gov
+//*
+//* 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 "ComputeMultipleInelasticStressSmallStrain.h"
+
+#include "StressUpdateBase.h"
+#include "MooseException.h"
+#include "MultipleInelasticStressHelper.h"
+#include "libmesh/int_range.h"
+
+registerMooseObject("SolidMechanicsApp", ComputeMultipleInelasticStressSmallStrain);
+
+InputParameters
+ComputeMultipleInelasticStressSmallStrain::validParams()
+{
+  InputParameters params = ComputeMultipleInelasticStressSmallStrainBase::validParams();
+  params.addRequiredParam<std::vector<MaterialName>>(
+      "inelastic_models",
+      "The material objects to use to calculate stress and inelastic strains. "
+      "Note: specify creep models first and plasticity models second.");
+  params.addClassDescription("Compute state (stress and internal parameters such as plastic "
+                             "strains and internal parameters) using an iterative process with "
+                             "small strain formulation. Combinations of creep models and plastic "
+                             "models may be used.");
+  return params;
+}
+
+ComputeMultipleInelasticStressSmallStrain::ComputeMultipleInelasticStressSmallStrain(
+    const InputParameters & parameters)
+  : ComputeMultipleInelasticStressSmallStrainBase(parameters)
+{
+}
+
+std::vector<MaterialName>
+ComputeMultipleInelasticStressSmallStrain::getInelasticModelNames()
+{
+  return getParam<std::vector<MaterialName>>("inelastic_models");
+}
+
+void
+ComputeMultipleInelasticStressSmallStrain::updateQpState(RankTwoTensor & elastic_strain,
+                                                         RankTwoTensor & inelastic_strain)
+{
+  if (_internal_solve_full_iteration_history == true)
+  {
+    _console << std::endl
+             << "iteration output for ComputeMultipleInelasticStressSmallStrain solve:"
+             << " time=" << _t << " int_pt=" << _qp << std::endl;
+  }
+
+  Real l2norm_delta_stress;
+  Real first_l2norm_delta_stress = 1.0;
+  unsigned int counter = 0;
+
+  std::vector<RankTwoTensor> inelastic_strain_increment;
+  inelastic_strain_increment.resize(_num_models);
+
+  for (const auto i_rmm : index_range(_models))
+    inelastic_strain_increment[i_rmm].zero();
+
+  RankTwoTensor stress_max, stress_min;
+
+  do
+  {
+    for (const auto i_rmm : index_range(_models))
+    {
+      _models[i_rmm]->setQp(_qp);
+
+      // Compute current elastic strain by subtracting all other inelastic strains
+      // Start with total mechanical strain
+      RankTwoTensor current_elastic_strain = _mechanical_strain[_qp] - _inelastic_strain_old[_qp];
+
+      // Subtract all inelastic strain increments calculated so far except the current one
+      for (const auto j_rmm : make_range(_num_models))
+        if (i_rmm != j_rmm)
+          current_elastic_strain -= inelastic_strain_increment[j_rmm];
+
+      // Form the trial stress
+      _stress[_qp] = _elasticity_tensor[_qp] * current_elastic_strain;
+
+      // Compute strain increment for this iteration
+      RankTwoTensor strain_increment = current_elastic_strain - _elastic_strain[_qp];
+
+      // Given a trial stress and strain increment, let the model produce an admissible stress
+      // and decompose the strain into elastic and inelastic parts
+      computeAdmissibleState(i_rmm,
+                             elastic_strain,
+                             inelastic_strain_increment[i_rmm],
+                             _consistent_tangent_operator[i_rmm]);
+
+      // Update stress min/max for convergence checking
+      MultipleInelasticStressHelper::updateStressMinMax(
+          _stress[_qp], stress_max, stress_min, i_rmm == 0);
+    }
+
+    // Check convergence: once the change in stress is within tolerance after each model
+    // consider the stress to be converged
+    l2norm_delta_stress =
+        MultipleInelasticStressHelper::computeStressDifferenceNorm(stress_max, stress_min);
+
+    if (counter == 0 && l2norm_delta_stress > 0.0)
+      first_l2norm_delta_stress = l2norm_delta_stress;
+
+    if (_internal_solve_full_iteration_history == true)
+    {
+      _console << "stress iteration number = " << counter << "\n"
+               << " relative l2 norm delta stress = "
+               << (0 == first_l2norm_delta_stress ? 0
+                                                  : l2norm_delta_stress / first_l2norm_delta_stress)
+               << "\n"
+               << " stress convergence relative tolerance = " << _relative_tolerance << "\n"
+               << " absolute l2 norm delta stress = " << l2norm_delta_stress << "\n"
+               << " stress convergence absolute tolerance = " << _absolute_tolerance << std::endl;
+    }
+    ++counter;
+  } while (counter < _max_iterations && l2norm_delta_stress > _absolute_tolerance &&
+           (l2norm_delta_stress / first_l2norm_delta_stress) > _relative_tolerance &&
+           _num_models != 1);
+
+  if (counter == _max_iterations && l2norm_delta_stress > _absolute_tolerance &&
+      (l2norm_delta_stress / first_l2norm_delta_stress) > _relative_tolerance)
+    throw MooseException(
+        "Max stress iteration hit during ComputeMultipleInelasticStressSmallStrain solve!");
+
+  // Combine inelastic strains from all models
+  inelastic_strain = MultipleInelasticStressHelper::computeCombinedInelasticStrainIncrement(
+      inelastic_strain_increment, _inelastic_weights, _num_models);
+
+  // Add to old inelastic strain
+  inelastic_strain += _inelastic_strain_old[_qp];
+
+  // Compute final elastic strain
+  elastic_strain = _mechanical_strain[_qp] - inelastic_strain;
+
+  if (_fe_problem.currentlyComputingJacobian())
+    computeQpJacobianMult();
+
+  _material_timestep_limit[_qp] =
+      MultipleInelasticStressHelper::computeMaterialTimestepLimit(_models, _num_models);
+}
diff --git a/modules/solid_mechanics/src/materials/ComputeMultipleInelasticStressSmallStrainBase.C b/modules/solid_mechanics/src/materials/ComputeMultipleInelasticStressSmallStrainBase.C
new file mode 100644
index 000000000000..18fc9d35869f
--- /dev/null
+++ b/modules/solid_mechanics/src/materials/ComputeMultipleInelasticStressSmallStrainBase.C
@@ -0,0 +1,350 @@
+//* This file is part of the MOOSE framework
+//* https://mooseframework.inl.gov
+//*
+//* 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 "ComputeMultipleInelasticStressSmallStrainBase.h"
+
+#include "StressUpdateBase.h"
+#include "MooseException.h"
+#include "DamageBase.h"
+#include "libmesh/int_range.h"
+
+InputParameters
+ComputeMultipleInelasticStressSmallStrainBase::validParams()
+{
+  InputParameters params = ComputeStressBase::validParams();
+  params.addClassDescription("Compute state (stress and internal parameters such as plastic "
+                             "strains and internal parameters) using an iterative process with "
+                             "small strain formulation. Combinations of creep models and plastic "
+                             "models may be used.");
+  params.addParam<unsigned int>("max_iterations",
+                                30,
+                                "Maximum number of the stress update "
+                                "iterations over the stress change after all "
+                                "update materials are called");
+  params.addParam<Real>("relative_tolerance",
+                        1e-5,
+                        "Relative convergence tolerance for the stress "
+                        "update iterations over the stress change "
+                        "after all update materials are called");
+  params.addParam<Real>("absolute_tolerance",
+                        1e-5,
+                        "Absolute convergence tolerance for the stress "
+                        "update iterations over the stress change "
+                        "after all update materials are called");
+  params.addParam<bool>(
+      "internal_solve_full_iteration_history",
+      false,
+      "Set to true to output stress update iteration information over the stress change");
+  MooseEnum tangent_operator("elastic nonlinear", "nonlinear");
+  params.addParam<MooseEnum>(
+      "tangent_operator",
+      tangent_operator,
+      "Type of tangent operator to return.  'elastic': return the "
+      "elasticity tensor.  'nonlinear': return the full, general consistent tangent "
+      "operator.");
+  params.addParam<std::vector<Real>>("combined_inelastic_strain_weights",
+                                     "The combined_inelastic_strain Material Property is a "
+                                     "weighted sum of the model inelastic strains.  This parameter "
+                                     "is a vector of weights, of the same length as "
+                                     "inelastic_models.  Default = '1 1 ... 1'.  This "
+                                     "parameter is set to 1 if the number of models = 1");
+  params.addParam<bool>(
+      "cycle_models", false, "At timestep N use only inelastic model N % num_models.");
+  params.addParam<MaterialName>("damage_model", "Name of the damage model");
+
+  return params;
+}
+
+ComputeMultipleInelasticStressSmallStrainBase::ComputeMultipleInelasticStressSmallStrainBase(
+    const InputParameters & parameters)
+  : ComputeStressBase(parameters),
+    GuaranteeConsumer(this),
+    _max_iterations(parameters.get<unsigned int>("max_iterations")),
+    _relative_tolerance(parameters.get<Real>("relative_tolerance")),
+    _absolute_tolerance(parameters.get<Real>("absolute_tolerance")),
+    _internal_solve_full_iteration_history(getParam<bool>("internal_solve_full_iteration_history")),
+    _elasticity_tensor_name(_base_name + "elasticity_tensor"),
+    _elasticity_tensor(getMaterialPropertyByName<RankFourTensor>(_elasticity_tensor_name)),
+    _mechanical_strain_old(getMaterialPropertyOld<RankTwoTensor>(_base_name + "mechanical_strain")),
+    _inelastic_strain(declareProperty<RankTwoTensor>(_base_name + "combined_inelastic_strain")),
+    _inelastic_strain_old(
+        getMaterialPropertyOld<RankTwoTensor>(_base_name + "combined_inelastic_strain")),
+    _stress_old(getMaterialPropertyOld<RankTwoTensor>(_base_name + "stress")),
+    _tangent_operator_type(getParam<MooseEnum>("tangent_operator").getEnum<TangentOperatorEnum>()),
+    _tangent_calculation_method(TangentCalculationMethod::ELASTIC),
+    _cycle_models(getParam<bool>("cycle_models")),
+    _material_timestep_limit(declareProperty<Real>(_base_name + "material_timestep_limit")),
+    _identity_symmetric_four(RankFourTensor::initIdentitySymmetricFour),
+    _all_models_isotropic(true)
+{
+}
+
+void
+ComputeMultipleInelasticStressSmallStrainBase::initQpStatefulProperties()
+{
+  ComputeStressBase::initQpStatefulProperties();
+  _inelastic_strain[_qp].zero();
+}
+
+void
+ComputeMultipleInelasticStressSmallStrainBase::initialSetup()
+{
+  _damage_model = isParamValid("damage_model")
+                      ? dynamic_cast<DamageBaseTempl<false> *>(&getMaterial("damage_model"))
+                      : nullptr;
+
+  _is_elasticity_tensor_guaranteed_isotropic =
+      hasGuaranteedMaterialProperty(_elasticity_tensor_name, Guarantee::ISOTROPIC);
+
+  std::vector<MaterialName> models = getInelasticModelNames();
+
+  _num_models = models.size();
+  _tangent_computation_flag.resize(_num_models, false);
+  _consistent_tangent_operator.resize(_num_models);
+
+  _inelastic_weights = isParamValid("combined_inelastic_strain_weights")
+                           ? getParam<std::vector<Real>>("combined_inelastic_strain_weights")
+                           : std::vector<Real>(_num_models, 1.0);
+
+  if (_inelastic_weights.size() != _num_models)
+    mooseError("ComputeMultipleInelasticStressSmallStrainBase: combined_inelastic_strain_weights "
+               "must contain the same number of entries as inelastic_models ",
+               _inelastic_weights.size(),
+               " ",
+               _num_models);
+
+  for (const auto i : make_range(_num_models))
+  {
+    StressUpdateBase * rrr = dynamic_cast<StressUpdateBase *>(&getMaterialByName(models[i]));
+
+    if (rrr)
+    {
+      _models.push_back(rrr);
+      _all_models_isotropic = _all_models_isotropic && rrr->isIsotropic();
+      if (rrr->requiresIsotropicTensor() && !_is_elasticity_tensor_guaranteed_isotropic)
+        mooseError("Model " + models[i] +
+                   " requires an isotropic elasticity tensor, but the one supplied is not "
+                   "guaranteed isotropic");
+    }
+    else
+      mooseError("Model " + models[i] +
+                 " is not compatible with ComputeMultipleInelasticStressSmallStrainBase");
+  }
+
+  // Check if tangent calculation methods are consistent
+  if (_tangent_operator_type != TangentOperatorEnum::elastic)
+  {
+    bool full_present = false;
+    bool partial_present = false;
+    for (const auto i : make_range(_num_models))
+    {
+      if (_models[i]->getTangentCalculationMethod() == TangentCalculationMethod::FULL)
+      {
+        full_present = true;
+        _tangent_computation_flag[i] = true;
+        _tangent_calculation_method = TangentCalculationMethod::FULL;
+      }
+      else if (_models[i]->getTangentCalculationMethod() == TangentCalculationMethod::PARTIAL)
+      {
+        partial_present = true;
+        _tangent_computation_flag[i] = true;
+        _tangent_calculation_method = TangentCalculationMethod::PARTIAL;
+      }
+    }
+    if (full_present && partial_present)
+      mooseError("In ",
+                 _name,
+                 ": Models that calculate the full tangent operator and the partial tangent "
+                 "operator are being combined. Either set tangent_operator to elastic, implement "
+                 "the correct tangent formulations, or use different models.");
+  }
+
+  if (isParamValid("damage_model") && !_damage_model)
+    paramError("damage_model",
+               "Damage Model " + _damage_model->name() +
+                   " is not compatible with ComputeMultipleInelasticStressSmallStrainBase");
+
+  // Check for substepping capability
+  for (const auto model_number : index_range(_models))
+  {
+    const bool use_substep = _models[model_number]->substeppingCapabilityRequested();
+    if (use_substep && !_models[model_number]->substeppingCapabilityEnabled())
+    {
+      std::stringstream error_message;
+      error_message << "Usage of substepping has been requested, but the inelastic model "
+                    << _models[model_number]->name() << " does not implement substepping yet.";
+      mooseError(error_message.str());
+    }
+  }
+}
+
+void
+ComputeMultipleInelasticStressSmallStrainBase::computeQpStress()
+{
+  if (_damage_model)
+  {
+    _undamaged_stress_old = _stress_old[_qp];
+    _damage_model->setQp(_qp);
+    _damage_model->computeUndamagedOldStress(_undamaged_stress_old);
+  }
+
+  RankTwoTensor elastic_strain;
+  RankTwoTensor inelastic_strain;
+
+  if (_num_models == 0)
+  {
+    // No inelastic models: purely elastic response
+    elastic_strain = _mechanical_strain[_qp];
+    _elastic_strain[_qp] = elastic_strain;
+    _stress[_qp] = _elasticity_tensor[_qp] * elastic_strain;
+
+    if (_fe_problem.currentlyComputingJacobian())
+      _Jacobian_mult[_qp] = _elasticity_tensor[_qp];
+
+    _material_timestep_limit[_qp] = std::numeric_limits<Real>::max();
+    inelastic_strain.zero();
+  }
+  else
+  {
+    if (_num_models == 1 || _cycle_models)
+      updateQpStateSingleModel((_t_step - 1) % _num_models, elastic_strain, inelastic_strain);
+    else
+      updateQpState(elastic_strain, inelastic_strain);
+
+    _elastic_strain[_qp] = elastic_strain;
+  }
+
+  _inelastic_strain[_qp] = inelastic_strain;
+
+  if (_damage_model)
+  {
+    _damage_model->setQp(_qp);
+    _damage_model->updateDamage();
+    _damage_model->updateStressForDamage(_stress[_qp]);
+    _damage_model->updateJacobianMultForDamage(_Jacobian_mult[_qp]);
+
+    const Real damage_timestep_limit = _damage_model->computeTimeStepLimit();
+    if (_material_timestep_limit[_qp] > damage_timestep_limit)
+      _material_timestep_limit[_qp] = damage_timestep_limit;
+  }
+}
+
+void
+ComputeMultipleInelasticStressSmallStrainBase::computeQpJacobianMult()
+{
+  _Jacobian_mult[_qp] = MultipleInelasticStressHelper::computeJacobianMult(
+      _tangent_calculation_method,
+      _elasticity_tensor[_qp],
+      _consistent_tangent_operator,
+      _num_models,
+      _identity_symmetric_four);
+}
+
+void
+ComputeMultipleInelasticStressSmallStrainBase::updateQpStateSingleModel(
+    unsigned model_number,
+    RankTwoTensor & elastic_strain,
+    RankTwoTensor & inelastic_strain)
+{
+  for (auto model : _models)
+    model->setQp(_qp);
+
+  // Compute trial stress assuming elastic response
+  _stress[_qp] = _elasticity_tensor[_qp] * _mechanical_strain[_qp];
+
+  // Let the model compute admissible stress and strain decomposition
+  computeAdmissibleState(
+      model_number, elastic_strain, inelastic_strain, _consistent_tangent_operator[0]);
+
+  if (_fe_problem.currentlyComputingJacobian())
+  {
+    if (_tangent_calculation_method == TangentCalculationMethod::ELASTIC)
+      _Jacobian_mult[_qp] = _elasticity_tensor[_qp];
+    else if (_tangent_calculation_method == TangentCalculationMethod::PARTIAL)
+    {
+      RankFourTensor A = _identity_symmetric_four + _consistent_tangent_operator[0];
+      mooseAssert(A.isSymmetric(), "Tangent operator isn't symmetric");
+      _Jacobian_mult[_qp] = A.invSymm() * _elasticity_tensor[_qp];
+    }
+    else
+      _Jacobian_mult[_qp] = _consistent_tangent_operator[0];
+  }
+
+  _material_timestep_limit[_qp] = _models[0]->computeTimeStepLimit();
+
+  /* propagate internal variables, etc, to this timestep for those inelastic models where
+   * "updateState" is not called */
+  for (const auto i_rmm : make_range(_num_models))
+    if (i_rmm != model_number)
+      _models[i_rmm]->propagateQpStatefulProperties();
+}
+
+void
+ComputeMultipleInelasticStressSmallStrainBase::computeAdmissibleState(
+    unsigned model_number,
+    RankTwoTensor & elastic_strain,
+    RankTwoTensor & inelastic_strain,
+    RankFourTensor & consistent_tangent_operator)
+{
+  // Reset properties to the beginning of the time step (necessary if substepping is employed).
+  _models[model_number]->resetIncrementalMaterialProperties();
+
+  // For small strain, we work with total strains
+  // Compute strain increment for the stress update models
+  RankTwoTensor strain_increment = _mechanical_strain[_qp] - _mechanical_strain_old[_qp];
+
+  const bool jac = _fe_problem.currentlyComputingJacobian();
+
+  // Rotation increment is identity for small strain
+  RankTwoTensor rotation_increment;
+  rotation_increment.setToIdentity();
+
+  if (_damage_model)
+    _models[model_number]->updateState(strain_increment,
+                                       inelastic_strain,
+                                       rotation_increment,
+                                       _stress[_qp],
+                                       _undamaged_stress_old,
+                                       _elasticity_tensor[_qp],
+                                       _elastic_strain[_qp],
+                                       (jac && _tangent_computation_flag[model_number]),
+                                       consistent_tangent_operator);
+  else if (_models[model_number]->substeppingCapabilityEnabled() &&
+           _is_elasticity_tensor_guaranteed_isotropic)
+    _models[model_number]->updateStateSubstep(strain_increment,
+                                              inelastic_strain,
+                                              rotation_increment,
+                                              _stress[_qp],
+                                              _stress_old[_qp],
+                                              _elasticity_tensor[_qp],
+                                              _elastic_strain[_qp],
+                                              (jac && _tangent_computation_flag[model_number]),
+                                              consistent_tangent_operator);
+  else
+    _models[model_number]->updateState(strain_increment,
+                                       inelastic_strain,
+                                       rotation_increment,
+                                       _stress[_qp],
+                                       _stress_old[_qp],
+                                       _elasticity_tensor[_qp],
+                                       _elastic_strain[_qp],
+                                       (jac && _tangent_computation_flag[model_number]),
+                                       consistent_tangent_operator);
+
+  // Compute elastic strain from mechanical strain and inelastic strain
+  elastic_strain = _mechanical_strain[_qp] - _inelastic_strain_old[_qp] - inelastic_strain;
+
+  if (jac && !_tangent_computation_flag[model_number])
+  {
+    if (_tangent_calculation_method == TangentCalculationMethod::PARTIAL)
+      consistent_tangent_operator.zero();
+    else
+      consistent_tangent_operator = _elasticity_tensor[_qp];
+  }
+}
diff --git a/modules/solid_mechanics/src/materials/MultipleInelasticStressHelper.C b/modules/solid_mechanics/src/materials/MultipleInelasticStressHelper.C
new file mode 100644
index 000000000000..31ebfb0bd352
--- /dev/null
+++ b/modules/solid_mechanics/src/materials/MultipleInelasticStressHelper.C
@@ -0,0 +1,128 @@
+//* This file is part of the MOOSE framework
+//* https://mooseframework.inl.gov
+//*
+//* 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 "MultipleInelasticStressHelper.h"
+#include "MooseUtils.h"
+#include "libmesh/int_range.h"
+
+RankFourTensor
+MultipleInelasticStressHelper::computeJacobianMult(
+    TangentCalculationMethod tangent_calculation_method,
+    const RankFourTensor & elasticity_tensor,
+    const std::vector<RankFourTensor> & consistent_tangent_operator,
+    unsigned int num_models,
+    const RankFourTensor & identity_symmetric_four)
+{
+  if (tangent_calculation_method == TangentCalculationMethod::ELASTIC)
+    return elasticity_tensor;
+  else if (tangent_calculation_method == TangentCalculationMethod::PARTIAL)
+  {
+    RankFourTensor A = identity_symmetric_four;
+    for (const auto i_rmm : make_range(num_models))
+      A += consistent_tangent_operator[i_rmm];
+    mooseAssert(A.isSymmetric(), "Tangent operator isn't symmetric");
+    return A.invSymm() * elasticity_tensor;
+  }
+  else // FULL
+  {
+    const RankFourTensor E_inv = elasticity_tensor.invSymm();
+    RankFourTensor result = consistent_tangent_operator[0];
+    for (const auto i_rmm : make_range(1u, num_models))
+      result = consistent_tangent_operator[i_rmm] * E_inv * result;
+    return result;
+  }
+}
+
+RankTwoTensor
+MultipleInelasticStressHelper::computeCombinedInelasticStrainIncrement(
+    const std::vector<RankTwoTensor> & inelastic_strain_increment,
+    const std::vector<Real> & inelastic_weights,
+    unsigned int num_models)
+{
+  RankTwoTensor combined;
+  combined.zero();
+  for (const auto i_rmm : make_range(num_models))
+    combined += inelastic_weights[i_rmm] * inelastic_strain_increment[i_rmm];
+  return combined;
+}
+
+Real
+MultipleInelasticStressHelper::computeMaterialTimestepLimit(
+    const std::vector<StressUpdateBase *> & models, unsigned int num_models)
+{
+  Real limit = 0.0;
+  for (const auto i_rmm : make_range(num_models))
+    limit += 1.0 / models[i_rmm]->computeTimeStepLimit();
+
+  if (MooseUtils::absoluteFuzzyEqual(limit, 0.0))
+    return std::numeric_limits<Real>::max();
+  else
+    return 1.0 / limit;
+}
+
+bool
+MultipleInelasticStressHelper::checkConvergence(const RankTwoTensor & stress_max,
+                                                const RankTwoTensor & stress_min,
+                                                unsigned int counter,
+                                                unsigned int max_iterations,
+                                                Real relative_tolerance,
+                                                Real absolute_tolerance,
+                                                Real first_l2norm_delta_stress,
+                                                unsigned int num_models)
+{
+  // Single model is always converged
+  if (num_models == 1)
+    return true;
+
+  Real l2norm_delta_stress = computeStressDifferenceNorm(stress_max, stress_min);
+
+  // Check convergence criteria
+  if (counter >= max_iterations)
+    return false;
+
+  if (l2norm_delta_stress <= absolute_tolerance)
+    return true;
+
+  if (first_l2norm_delta_stress > 0.0 &&
+      (l2norm_delta_stress / first_l2norm_delta_stress) <= relative_tolerance)
+    return true;
+
+  return false;
+}
+
+void
+MultipleInelasticStressHelper::updateStressMinMax(const RankTwoTensor & stress,
+                                                  RankTwoTensor & stress_max,
+                                                  RankTwoTensor & stress_min,
+                                                  bool is_first_model)
+{
+  if (is_first_model)
+  {
+    stress_max = stress;
+    stress_min = stress;
+  }
+  else
+  {
+    for (const auto i : make_range(Moose::dim))
+      for (const auto j : make_range(Moose::dim))
+      {
+        if (stress(i, j) > stress_max(i, j))
+          stress_max(i, j) = stress(i, j);
+        else if (stress_min(i, j) > stress(i, j))
+          stress_min(i, j) = stress(i, j);
+      }
+  }
+}
+
+Real
+MultipleInelasticStressHelper::computeStressDifferenceNorm(const RankTwoTensor & stress_max,
+                                                           const RankTwoTensor & stress_min)
+{
+  return (stress_max - stress_min).L2norm();
+}
diff --git a/modules/solid_mechanics/test/tests/multiple_inelastic_stress_small_strain/gold/isotropic_plasticity_small_strain_out.e b/modules/solid_mechanics/test/tests/multiple_inelastic_stress_small_strain/gold/isotropic_plasticity_small_strain_out.e
new file mode 100644
index 000000000000..ca44678ad9b1
Binary files /dev/null and b/modules/solid_mechanics/test/tests/multiple_inelastic_stress_small_strain/gold/isotropic_plasticity_small_strain_out.e differ
diff --git a/modules/solid_mechanics/test/tests/multiple_inelastic_stress_small_strain/isotropic_plasticity_small_strain.i b/modules/solid_mechanics/test/tests/multiple_inelastic_stress_small_strain/isotropic_plasticity_small_strain.i
new file mode 100644
index 000000000000..093504f132ab
--- /dev/null
+++ b/modules/solid_mechanics/test/tests/multiple_inelastic_stress_small_strain/isotropic_plasticity_small_strain.i
@@ -0,0 +1,126 @@
+# Test for ComputeMultipleInelasticStressSmallStrain with isotropic plasticity
+# This test uses the small strain (total) formulation with J2 plasticity
+#
+# A single element is loaded in tension in the y direction while fixing
+# displacement in x and z directions. This test verifies that the small strain
+# multiple inelastic stress calculator produces correct stress and plastic strain.
+#
+# This is similar to isotropic_plasticity_incremental_strain.i but uses
+# ComputeMultipleInelasticStressSmallStrain instead of ComputeMultipleInelasticStress
+
+[Mesh]
+  type = GeneratedMesh
+  dim = 3
+  nx = 1
+  ny = 1
+  nz = 1
+[]
+
+[GlobalParams]
+  displacements = 'disp_x disp_y disp_z'
+[]
+
+[Functions]
+  [top_pull]
+    type = ParsedFunction
+    expression = 't * 0.01'
+  []
+[]
+
+[Physics/SolidMechanics/QuasiStatic]
+  [all]
+    strain = SMALL
+    incremental = false
+    add_variables = true
+    generate_output = 'stress_yy elastic_strain_yy'
+  []
+[]
+
+[BCs]
+  [y_pull_function]
+    type = FunctionDirichletBC
+    variable = disp_y
+    boundary = top
+    function = top_pull
+  []
+  [x_sides]
+    type = DirichletBC
+    variable = disp_x
+    boundary = 'left right'
+    value = 0.0
+  []
+  [y_bot]
+    type = DirichletBC
+    variable = disp_y
+    boundary = bottom
+    value = 0.0
+  []
+  [z_sides]
+    type = DirichletBC
+    variable = disp_z
+    boundary = 'back front'
+    value = 0.0
+  []
+[]
+
+[AuxVariables]
+  [plastic_strain_yy]
+    order = CONSTANT
+    family = MONOMIAL
+  []
+[]
+
+[AuxKernels]
+  [plastic_strain_yy]
+    type = RankTwoAux
+    rank_two_tensor = combined_inelastic_strain
+    variable = plastic_strain_yy
+    index_i = 1
+    index_j = 1
+    execute_on = timestep_end
+  []
+[]
+
+[Materials]
+  [elasticity_tensor]
+    type = ComputeIsotropicElasticityTensor
+    youngs_modulus = 2.5e5
+    poissons_ratio = 0.0
+  []
+  [isotropic_plasticity]
+    type = IsotropicPlasticityStressUpdate
+    yield_stress = 25.0
+    hardening_constant = 1000.0
+  []
+  [stress]
+    type = ComputeMultipleInelasticStressSmallStrain
+    inelastic_models = 'isotropic_plasticity'
+    tangent_operator = elastic
+  []
+[]
+
+[Executioner]
+  type = Transient
+  solve_type = 'PJFNK'
+
+  petsc_options = '-snes_ksp_ew'
+  petsc_options_iname = '-ksp_gmres_restart'
+  petsc_options_value = '101'
+
+  line_search = 'none'
+
+  l_max_its = 20
+  nl_max_its = 20
+  nl_rel_tol = 1e-10
+  nl_abs_tol = 1e-12
+  l_tol = 1e-9
+
+  start_time = 0.0
+  end_time = 0.01875
+  dt = 0.00125
+  dtmin = 0.0001
+[]
+
+[Outputs]
+  exodus = true
+[]
diff --git a/modules/solid_mechanics/test/tests/multiple_inelastic_stress_small_strain/tests b/modules/solid_mechanics/test/tests/multiple_inelastic_stress_small_strain/tests
new file mode 100644
index 000000000000..a021b4525496
--- /dev/null
+++ b/modules/solid_mechanics/test/tests/multiple_inelastic_stress_small_strain/tests
@@ -0,0 +1,13 @@
+[Tests]
+  design = 'ComputeMultipleInelasticStressSmallStrain.md'
+  issues = '#TBD'
+
+  [isotropic_plasticity_small_strain]
+    type = Exodiff
+    input = 'isotropic_plasticity_small_strain.i'
+    exodiff = 'isotropic_plasticity_small_strain_out.e'
+    requirement = 'The system shall compute the J2 isotropic plasticity stress and plastic strain '
+                  'response under tensile loading using the small strain (total) formulation with '
+                  'ComputeMultipleInelasticStressSmallStrain.'
+  []
+[]
diff --git a/modules/solid_mechanics/test/tests/multiple_inelastic_stress_small_strain/two_models_small_strain.i b/modules/solid_mechanics/test/tests/multiple_inelastic_stress_small_strain/two_models_small_strain.i
new file mode 100644
index 000000000000..16036dcc2cd2
--- /dev/null
+++ b/modules/solid_mechanics/test/tests/multiple_inelastic_stress_small_strain/two_models_small_strain.i
@@ -0,0 +1,149 @@
+# Test for ComputeMultipleInelasticStressSmallStrain with two inelastic models
+# This test combines power law creep and isotropic plasticity using small strain formulation
+#
+# A single element is loaded in tension. This verifies that the iterative
+# algorithm in ComputeMultipleInelasticStressSmallStrain correctly handles
+# multiple inelastic models and converges to the correct combined inelastic strain.
+
+[Mesh]
+  type = GeneratedMesh
+  dim = 3
+  nx = 1
+  ny = 1
+  nz = 1
+[]
+
+[GlobalParams]
+  displacements = 'disp_x disp_y disp_z'
+[]
+
+[Functions]
+  [top_pull]
+    type = ParsedFunction
+    expression = 't * 0.005'
+  []
+[]
+
+[Physics/SolidMechanics/QuasiStatic]
+  [all]
+    strain = SMALL
+    incremental = false
+    add_variables = true
+    generate_output = 'stress_yy'
+  []
+[]
+
+[BCs]
+  [y_pull_function]
+    type = FunctionDirichletBC
+    variable = disp_y
+    boundary = top
+    function = top_pull
+  []
+  [x_sides]
+    type = DirichletBC
+    variable = disp_x
+    boundary = 'left right'
+    value = 0.0
+  []
+  [y_bot]
+    type = DirichletBC
+    variable = disp_y
+    boundary = bottom
+    value = 0.0
+  []
+  [z_sides]
+    type = DirichletBC
+    variable = disp_z
+    boundary = 'back front'
+    value = 0.0
+  []
+[]
+
+[AuxVariables]
+  [combined_inelastic_strain_yy]
+    order = CONSTANT
+    family = MONOMIAL
+  []
+[]
+
+[AuxKernels]
+  [combined_inelastic_strain_yy]
+    type = RankTwoAux
+    rank_two_tensor = combined_inelastic_strain
+    variable = combined_inelastic_strain_yy
+    index_i = 1
+    index_j = 1
+    execute_on = timestep_end
+  []
+[]
+
+[Materials]
+  [elasticity_tensor]
+    type = ComputeIsotropicElasticityTensor
+    youngs_modulus = 2.0e5
+    poissons_ratio = 0.3
+  []
+
+  [creep]
+    type = PowerLawCreepStressUpdate
+    coefficient = 1.0e-15
+    n_exponent = 5.0
+    m_exponent = 0.0
+    activation_energy = 0.0
+  []
+
+  [plasticity]
+    type = IsotropicPlasticityStressUpdate
+    yield_stress = 150.0
+    hardening_constant = 500.0
+  []
+
+  [stress]
+    type = ComputeMultipleInelasticStressSmallStrain
+    inelastic_models = 'creep plasticity'
+    max_iterations = 50
+    relative_tolerance = 1e-8
+    absolute_tolerance = 1e-10
+  []
+[]
+
+[Preconditioning]
+  [smp]
+    type = SMP
+    full = true
+  []
+[]
+
+[Executioner]
+  type = Transient
+  solve_type = 'NEWTON'
+
+  petsc_options_iname = '-pc_type'
+  petsc_options_value = 'lu'
+
+  l_max_its = 20
+  nl_max_its = 20
+  nl_rel_tol = 1e-8
+  nl_abs_tol = 1e-10
+
+  start_time = 0.0
+  end_time = 2.0
+  dt = 0.5
+[]
+
+[Outputs]
+  exodus = true
+  csv = true
+[]
+
+[Postprocessors]
+  [stress_yy]
+    type = ElementAverageValue
+    variable = stress_yy
+  []
+  [combined_strain_yy]
+    type = ElementAverageValue
+    variable = combined_inelastic_strain_yy
+  []
+[]