diff --git a/framework/include/materials/GenericConstantRealVectorValue.h b/framework/include/materials/GenericConstantRealVectorValue.h
new file mode 100644
index 000000000000..97accfb9d4ac
--- /dev/null
+++ b/framework/include/materials/GenericConstantRealVectorValue.h
@@ -0,0 +1,34 @@
+//* 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 "Material.h"
+
+/**
+ * Declares a constant material property of type RealVectorValue.
+ */
+template <bool is_ad>
+class GenericConstantRealVectorValueTempl : public Material
+{
+public:
+  static InputParameters validParams();
+
+  GenericConstantRealVectorValueTempl(const InputParameters & parameters);
+
+protected:
+  virtual void initQpStatefulProperties() override;
+  virtual void computeQpProperties() override;
+
+  const RealVectorValue _vector;
+  GenericMaterialProperty<RealVectorValue, is_ad> & _prop;
+};
+
+typedef GenericConstantRealVectorValueTempl<false> GenericConstantRealVectorValue;
+typedef GenericConstantRealVectorValueTempl<true> ADGenericConstantRealVectorValue;
diff --git a/framework/include/utils/MooseArray.h b/framework/include/utils/MooseArray.h
index 1d4941628bfa..ebde89e93773 100644
--- a/framework/include/utils/MooseArray.h
+++ b/framework/include/utils/MooseArray.h
@@ -173,10 +173,13 @@ class MooseArray
    */
   std::vector<T> stdVector() const;
 
+  ///@{
   /**
    * Reference to first element of array
    */
   const T * data() const { return _data; }
+  T * data() { return _data; }
+  ///@}
 
 private:
   /// Smart pointer storage
diff --git a/framework/include/utils/RankFourTensor.h b/framework/include/utils/RankFourTensor.h
index b930eae69105..350b730bf031 100644
--- a/framework/include/utils/RankFourTensor.h
+++ b/framework/include/utils/RankFourTensor.h
@@ -201,7 +201,7 @@ class RankFourTensorTempl
   void printReal(std::ostream & stm = Moose::out) const;
 
   /// copies values from a into this tensor
-  RankFourTensorTempl<T> & operator=(const RankFourTensorTempl<T> & a);
+  RankFourTensorTempl<T> & operator=(const RankFourTensorTempl<T> & a) = default;
 
   /**
    * Assignment-from-scalar operator.  Used only to zero out the tensor.
diff --git a/framework/include/utils/RankFourTensorImplementation.h b/framework/include/utils/RankFourTensorImplementation.h
index a5cd05461db0..06bb326ba52e 100644
--- a/framework/include/utils/RankFourTensorImplementation.h
+++ b/framework/include/utils/RankFourTensorImplementation.h
@@ -120,15 +120,6 @@ RankFourTensorTempl<T>::zero()
     _vals[i] = 0.0;
 }
 
-template <typename T>
-RankFourTensorTempl<T> &
-RankFourTensorTempl<T>::operator=(const RankFourTensorTempl<T> & a)
-{
-  for (auto i : make_range(N4))
-    _vals[i] = a._vals[i];
-  return *this;
-}
-
 template <typename T>
 template <template <typename> class Tensor, typename T2>
 auto
diff --git a/framework/include/utils/RankTwoTensor.h b/framework/include/utils/RankTwoTensor.h
index 3c0caf5afbb6..66e5e6f961eb 100644
--- a/framework/include/utils/RankTwoTensor.h
+++ b/framework/include/utils/RankTwoTensor.h
@@ -497,7 +497,7 @@ class RankTwoTensorTempl : public libMesh::TensorValue<T>
    * //       7 4 1 ]
    * ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    */
-  RankTwoTensorTempl<T> & operator=(const RankTwoTensorTempl<T> & a);
+  RankTwoTensorTempl<T> & operator=(const RankTwoTensorTempl<T> & a) = default;
 
   /**
    * @brief Assignment operator (from a ColumnMajorMatrixTempl<T>)
diff --git a/framework/include/utils/RankTwoTensorImplementation.h b/framework/include/utils/RankTwoTensorImplementation.h
index 27563a0e9106..b69ac378a50e 100644
--- a/framework/include/utils/RankTwoTensorImplementation.h
+++ b/framework/include/utils/RankTwoTensorImplementation.h
@@ -369,14 +369,6 @@ RankTwoTensorTempl<T>::transpose() const
   return libMesh::TensorValue<T>::transpose();
 }
 
-template <typename T>
-RankTwoTensorTempl<T> &
-RankTwoTensorTempl<T>::operator=(const RankTwoTensorTempl<T> & a)
-{
-  libMesh::TensorValue<T>::operator=(a);
-  return *this;
-}
-
 template <typename T>
 RankTwoTensorTempl<T> &
 RankTwoTensorTempl<T>::operator+=(const RankTwoTensorTempl<T> & a)
diff --git a/framework/include/utils/SymmetricRankFourTensor.h b/framework/include/utils/SymmetricRankFourTensor.h
index 7878acd7014e..c17ddd0c96c2 100644
--- a/framework/include/utils/SymmetricRankFourTensor.h
+++ b/framework/include/utils/SymmetricRankFourTensor.h
@@ -189,7 +189,7 @@ class SymmetricRankFourTensorTempl
   void printReal(std::ostream & stm = Moose::out) const;
 
   /// copies values from a into this tensor
-  SymmetricRankFourTensorTempl<T> & operator=(const SymmetricRankFourTensorTempl<T> & a);
+  SymmetricRankFourTensorTempl<T> & operator=(const SymmetricRankFourTensorTempl<T> & a) = default;
 
   /**
    * Assignment-from-scalar operator.  Used only to zero out the tensor.
diff --git a/framework/include/utils/SymmetricRankFourTensorImplementation.h b/framework/include/utils/SymmetricRankFourTensorImplementation.h
index 7a47d67790d4..e0004be3b9f3 100644
--- a/framework/include/utils/SymmetricRankFourTensorImplementation.h
+++ b/framework/include/utils/SymmetricRankFourTensorImplementation.h
@@ -159,14 +159,6 @@ SymmetricRankFourTensorTempl<T>::rotate(const TypeTensor<T> & R)
   (*this) = M * (*this) * M.transposeMajor();
 }
 
-template <typename T>
-SymmetricRankFourTensorTempl<T> &
-SymmetricRankFourTensorTempl<T>::operator=(const SymmetricRankFourTensorTempl<T> & a)
-{
-  std::copy(a._vals.begin(), a._vals.end(), _vals.begin());
-  return *this;
-}
-
 template <typename T>
 SymmetricRankFourTensorTempl<T> &
 SymmetricRankFourTensorTempl<T>::operator*=(const T & a)
diff --git a/framework/src/materials/GenericConstantRealVectorValue.C b/framework/src/materials/GenericConstantRealVectorValue.C
new file mode 100644
index 000000000000..41cc7667dccf
--- /dev/null
+++ b/framework/src/materials/GenericConstantRealVectorValue.C
@@ -0,0 +1,53 @@
+//* 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 "GenericConstantRealVectorValue.h"
+
+registerMooseObject("MooseApp", GenericConstantRealVectorValue);
+registerMooseObject("MooseApp", ADGenericConstantRealVectorValue);
+
+template <bool is_ad>
+InputParameters
+GenericConstantRealVectorValueTempl<is_ad>::validParams()
+{
+  InputParameters params = Material::validParams();
+  params.addClassDescription("Object for declaring a constant 3-vector as a material property.");
+  params.addRequiredParam<RealVectorValue>("vector_values", "Values defining the constant vector");
+  params.addRequiredParam<MaterialPropertyName>(
+      "vector_name", "Name of the vector material property to be created");
+  params.set<MooseEnum>("constant_on") = "SUBDOMAIN";
+  return params;
+}
+
+template <bool is_ad>
+GenericConstantRealVectorValueTempl<is_ad>::GenericConstantRealVectorValueTempl(
+    const InputParameters & parameters)
+  : Material(parameters),
+    _vector(getParam<RealVectorValue>("vector_values")),
+    _prop(declareGenericProperty<RealVectorValue, is_ad>(
+        getParam<MaterialPropertyName>("vector_name")))
+{
+}
+
+template <bool is_ad>
+void
+GenericConstantRealVectorValueTempl<is_ad>::initQpStatefulProperties()
+{
+  GenericConstantRealVectorValueTempl<is_ad>::computeQpProperties();
+}
+
+template <bool is_ad>
+void
+GenericConstantRealVectorValueTempl<is_ad>::computeQpProperties()
+{
+  _prop[_qp] = _vector;
+}
+
+template class GenericConstantRealVectorValueTempl<false>;
+template class GenericConstantRealVectorValueTempl<true>;
diff --git a/modules/combined/test/tests/optimization/invOpt_elasticity_modular/elasticity.i b/modules/combined/test/tests/optimization/invOpt_elasticity_modular/elasticity.i
index bfc34feabece..c40702f52d9a 100644
--- a/modules/combined/test/tests/optimization/invOpt_elasticity_modular/elasticity.i
+++ b/modules/combined/test/tests/optimization/invOpt_elasticity_modular/elasticity.i
@@ -1,14 +1,14 @@
 [Models]
   [adjoint_elasticity_model]
     type = LinearIsotropicElasticity
-    youngs_modulus = 5.0
-    poisson_ratio = 0.3
+    coefficients = '5.0 0.3'
+    coefficient_types = 'YOUNGS_MODULUS POISSONS_RATIO'
     strain = 'forces/E'
   []
   [forward_elasticity_model]
     type = LinearIsotropicElasticity
-    youngs_modulus = 5.0
-    poisson_ratio = 0.3
+    coefficients = '5.0 0.3'
+    coefficient_types = 'YOUNGS_MODULUS POISSONS_RATIO'
     strain = 'forces/E'
   []
 []
diff --git a/modules/combined/test/tests/optimization/invOpt_elasticity_modular/grad.i b/modules/combined/test/tests/optimization/invOpt_elasticity_modular/grad.i
index b02824f7b7e8..2f9b40927c5d 100644
--- a/modules/combined/test/tests/optimization/invOpt_elasticity_modular/grad.i
+++ b/modules/combined/test/tests/optimization/invOpt_elasticity_modular/grad.i
@@ -69,7 +69,6 @@
   verbose = true
   device = 'cpu'
   [forward]
-    enable_AD = true
     model = 'forward_elasticity_model'
 
     moose_input_types = 'MATERIAL'
diff --git a/modules/solid_mechanics/contrib/neml2 b/modules/solid_mechanics/contrib/neml2
index 91c735e9d92e..f0891dce8016 160000
--- a/modules/solid_mechanics/contrib/neml2
+++ b/modules/solid_mechanics/contrib/neml2
@@ -1 +1 @@
-Subproject commit 91c735e9d92e1c101bb4af2ec89fed5666402eb7
+Subproject commit f0891dce80164f55cbb38a640bd85527599430c9
diff --git a/modules/solid_mechanics/include/materials/lagrangian/ComputeLagrangianCauchyCustomStress.h b/modules/solid_mechanics/include/materials/lagrangian/ComputeLagrangianCauchyCustomStress.h
new file mode 100644
index 000000000000..37a8b935595f
--- /dev/null
+++ b/modules/solid_mechanics/include/materials/lagrangian/ComputeLagrangianCauchyCustomStress.h
@@ -0,0 +1,30 @@
+//* 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 "ComputeLagrangianStressCauchy.h"
+#include "DerivativeMaterialPropertyNameInterface.h"
+
+/// Provide the Cauchy stress and jacobian directly
+///
+class ComputeLagrangianCauchyCustomStress : public ComputeLagrangianStressCauchy,
+                                            public DerivativeMaterialPropertyNameInterface
+{
+public:
+  static InputParameters validParams();
+  ComputeLagrangianCauchyCustomStress(const InputParameters & parameters);
+
+protected:
+  /// Implement the copy
+  virtual void computeQpCauchyStress() override;
+
+  const MaterialProperty<RankTwoTensor> & _custom_stress;
+  const MaterialProperty<RankFourTensor> & _custom_jacobian;
+};
diff --git a/modules/solid_mechanics/include/materials/lagrangian/ComputeLagrangianStrainBase.h b/modules/solid_mechanics/include/materials/lagrangian/ComputeLagrangianStrainBase.h
index 11e4df4bd258..a3fd4a6514f0 100644
--- a/modules/solid_mechanics/include/materials/lagrangian/ComputeLagrangianStrainBase.h
+++ b/modules/solid_mechanics/include/materials/lagrangian/ComputeLagrangianStrainBase.h
@@ -88,6 +88,9 @@ class ComputeLagrangianStrainBase : public Material, public G
   /// Strain increment
   MaterialProperty<RankTwoTensor> & _strain_increment;
 
+  /// Spatial velocity gradient increment
+  MaterialProperty<RankTwoTensor> & _spatial_velocity_increment;
+
   /// Vorticity increment
   MaterialProperty<RankTwoTensor> & _vorticity_increment;
 
diff --git a/modules/solid_mechanics/include/neml2/actions/NEML2Action.h b/modules/solid_mechanics/include/neml2/actions/NEML2Action.h
index db156e989ab8..edaf5a25969f 100644
--- a/modules/solid_mechanics/include/neml2/actions/NEML2Action.h
+++ b/modules/solid_mechanics/include/neml2/actions/NEML2Action.h
@@ -15,6 +15,8 @@
 
 #include "Action.h"
 
+class NEML2ActionCommon;
+
 /**
  * Action to set up NEML2 objects.
  */
@@ -28,6 +30,8 @@ class NEML2Action : public Action
   virtual void act() override;
 
 protected:
+  const NEML2ActionCommon & getCommonAction() const;
+
 #ifdef NEML2_ENABLED
 
   enum class MOOSEIOType
@@ -130,6 +134,35 @@ class NEML2Action : public Action
   /// Blocks this sub-block action applies to
   const std::vector<SubdomainName> _block;
 
+  /// Material property initial conditions
+  std::map<MaterialPropertyName, MaterialPropertyName> _initialize_output_values;
+
   /// Material property additional outputs
   std::map<MaterialPropertyName, std::vector<OutputName>> _export_output_targets;
+
+private:
+  /// Get parameter lists for mapping between MOOSE and NEML2 quantities
+  template <typename EnumType, typename T1, typename T2>
+  std::tuple<std::vector<EnumType>, std::vector<T1>, std::vector<T2>>
+  getInputParameterMapping(const std::string & moose_type_opt,
+                           const std::string & moose_name_opt,
+                           const std::string & neml2_name_opt) const
+  {
+    const auto moose_types = getParam<MultiMooseEnum>(moose_type_opt).getSetValueIDs<EnumType>();
+    const auto moose_names = getParam<std::vector<T1>>(moose_name_opt);
+    const auto neml2_names = getParam<std::vector<T2>>(neml2_name_opt);
+
+    if (moose_types.size() != moose_names.size())
+      paramError(moose_name_opt, moose_name_opt, " must have the same length as ", moose_type_opt);
+    if (moose_names.size() != neml2_names.size())
+      paramError(moose_name_opt, moose_name_opt, " must have the same length as ", neml2_name_opt);
+
+    return {moose_types, moose_names, neml2_names};
+  }
+
+  /// Print a summary of the NEML2 model
+  void printSummary(const neml2::Model &) const;
+
+  /// Get the maximum length of all MOOSE names (for printing purposes)
+  std::size_t getLongestMOOSEName() const;
 };
diff --git a/modules/solid_mechanics/include/neml2/actions/NEML2ActionCommon.h b/modules/solid_mechanics/include/neml2/actions/NEML2ActionCommon.h
index acc6840c91e7..67d198a59913 100644
--- a/modules/solid_mechanics/include/neml2/actions/NEML2ActionCommon.h
+++ b/modules/solid_mechanics/include/neml2/actions/NEML2ActionCommon.h
@@ -47,13 +47,19 @@ class NEML2ActionCommon : public Action
   /// Parameters that can be specified EITHER under the common area OR under sub-blocks
   static InputParameters commonParams();
 
+  /// Parameters that can ONLY be specified under the common area
   static InputParameters validParams();
 
   NEML2ActionCommon(const InputParameters &);
 
   virtual void act() override;
 
+  const FileName & fname() const { return _fname; }
+
 protected:
   /// Name of the NEML2 input file
   const FileName _fname;
+
+  /// List of cli-args
+  const std::vector<std::string> _cli_args;
 };
diff --git a/modules/solid_mechanics/include/neml2/interfaces/MOOSEToNEML2.h b/modules/solid_mechanics/include/neml2/interfaces/MOOSEToNEML2.h
index 1ce75f9a2e56..6ca535548a0c 100644
--- a/modules/solid_mechanics/include/neml2/interfaces/MOOSEToNEML2.h
+++ b/modules/solid_mechanics/include/neml2/interfaces/MOOSEToNEML2.h
@@ -13,7 +13,6 @@
 #include "InputParameters.h"
 
 #ifdef NEML2_ENABLED
-#include "neml2/tensors/LabeledVector.h"
 #include "neml2/models/Model.h"
 #endif
 
@@ -70,7 +69,7 @@ class MOOSEToNEML2
   virtual neml2::Tensor gatheredData() const = 0;
 
   /// Insert the gathered data into the NEML2 material model
-  void insertInto(neml2::Model &) const;
+  void insertInto(neml2::ValueMap &, std::map<std::string, neml2::Tensor> &) const;
 
 protected:
   /// Whether we should insert into NEML2 input variable or NEML2 model parameter
diff --git a/modules/solid_mechanics/include/neml2/interfaces/NEML2ModelInterface.h b/modules/solid_mechanics/include/neml2/interfaces/NEML2ModelInterface.h
index ff7ed5c37204..7116ce7fe559 100644
--- a/modules/solid_mechanics/include/neml2/interfaces/NEML2ModelInterface.h
+++ b/modules/solid_mechanics/include/neml2/interfaces/NEML2ModelInterface.h
@@ -9,6 +9,8 @@
 
 #pragma once
 
+#include "NEML2Utils.h"
+
 #ifdef NEML2_ENABLED
 #include "neml2/models/Model.h"
 #include "neml2/misc/parser_utils.h"
@@ -16,8 +18,6 @@
 #include "UserObject.h"
 #endif
 
-#include "NEML2Utils.h"
-
 /**
  * Interface class to provide common input parameters, members, and methods for MOOSEObjects that
  * use NEML2 models.
@@ -35,29 +35,17 @@ class NEML2ModelInterface : public T
 
 protected:
   /**
-   * Validate the NEML2 material model. This method should throw a moose error with the first
-   * encountered problem. Note that the developer is responsible for calling this method at the
-   * appropriate times, for example, at initialSetup().
+   * Validate the NEML2 material model. Note that the developer is responsible for calling this
+   * method at the appropriate times, for example, at initialSetup().
    */
   virtual void validateModel() const;
 
-  /// Initialize the model with a batch shape
-  void initModel(neml2::TensorShapeRef batch_shape);
-
   /// Get the NEML2 model
   neml2::Model & model() const { return _model; }
 
   /// Get the target compute device
   const torch::Device & device() const { return _device; }
 
-  /**
-   * @brief Convert a raw string to a neml2::VariableName
-   *
-   * @param raw_str
-   * @return neml2::VariableName
-   */
-  neml2::VariableName getNEML2VariableName(const std::string & raw_str) const;
-
 private:
   /// The NEML2 material model
   neml2::Model & _model;
@@ -86,11 +74,6 @@ NEML2ModelInterface<T>::validParams()
       ":<device-index> optionally specifies a device index. For example, device='cpu' sets the "
       "target compute device to be CPU, and device='cuda:1' sets the target compute device to be "
       "CUDA with device ID 1.");
-  params.addParam<bool>("enable_AD",
-                        false,
-                        "Set to true to enable PyTorch AD. When set to false (default), no "
-                        "function graph or gradient is computed, which speeds up model "
-                        "evaluation.");
   return params;
 }
 
@@ -109,7 +92,7 @@ template <class T>
 template <typename... P>
 NEML2ModelInterface<T>::NEML2ModelInterface(const InputParameters & params, P &&... args)
   : T(params, args...),
-    _model(neml2::get_model(params.get<std::string>("model"), params.get<bool>("enable_AD"))),
+    _model(neml2::get_model(params.get<std::string>("model"))),
     _device(params.get<std::string>("device"))
 {
 }
@@ -121,18 +104,4 @@ NEML2ModelInterface<T>::validateModel() const
   neml2::diagnose(_model);
 }
 
-template <class T>
-void
-NEML2ModelInterface<T>::initModel(neml2::TensorShapeRef batch_shape)
-{
-  _model.reinit(batch_shape, /*deriv_order=*/1, _device, /*dtype=*/torch::kFloat64);
-}
-
-template <class T>
-neml2::VariableName
-NEML2ModelInterface<T>::getNEML2VariableName(const std::string & raw_str) const
-{
-  return neml2::utils::parse<neml2::VariableName>(raw_str);
-}
-
 #endif // NEML2_ENABLED
diff --git a/modules/solid_mechanics/include/neml2/materials/NEML2ToMOOSEMaterialProperty.h b/modules/solid_mechanics/include/neml2/materials/NEML2ToMOOSEMaterialProperty.h
index 98421b5ae770..d3c92d528d99 100644
--- a/modules/solid_mechanics/include/neml2/materials/NEML2ToMOOSEMaterialProperty.h
+++ b/modules/solid_mechanics/include/neml2/materials/NEML2ToMOOSEMaterialProperty.h
@@ -32,7 +32,7 @@ class NEML2ToMOOSEMaterialProperty : public Material
   void computeProperties() override;
 
 protected:
-  void initQpStatefulProperties() override;
+  void initQpStatefulProperties() override {}
 
   /// User object managing the execution of the NEML2 model
   const NEML2ModelExecutor & _execute_neml2_model;
@@ -43,8 +43,8 @@ class NEML2ToMOOSEMaterialProperty : public Material
   /// Initial condition
   const MaterialProperty<T> * _prop0;
 
-  /// labled view to the requested output (or its derivative)
-  const neml2::Tensor & _output_view;
+  /// Reference to the requested output (or its derivative) value
+  const neml2::Tensor & _value;
 #endif
 };
 
@@ -54,4 +54,6 @@ class NEML2ToMOOSEMaterialProperty : public Material
 DefineNEML2ToMOOSEMaterialPropertyAlias(Real, Real);
 DefineNEML2ToMOOSEMaterialPropertyAlias(SymmetricRankTwoTensor, SymmetricRankTwoTensor);
 DefineNEML2ToMOOSEMaterialPropertyAlias(SymmetricRankFourTensor, SymmetricRankFourTensor);
-DefineNEML2ToMOOSEMaterialPropertyAlias(std::vector<Real>, StdVector);
+DefineNEML2ToMOOSEMaterialPropertyAlias(RealVectorValue, RealVectorValue);
+DefineNEML2ToMOOSEMaterialPropertyAlias(RankTwoTensor, RankTwoTensor);
+DefineNEML2ToMOOSEMaterialPropertyAlias(RankFourTensor, RankFourTensor);
diff --git a/modules/solid_mechanics/include/neml2/userobjects/MOOSEMaterialPropertyToNEML2.h b/modules/solid_mechanics/include/neml2/userobjects/MOOSEMaterialPropertyToNEML2.h
index ad515c1e07d6..ede0196eb02d 100644
--- a/modules/solid_mechanics/include/neml2/userobjects/MOOSEMaterialPropertyToNEML2.h
+++ b/modules/solid_mechanics/include/neml2/userobjects/MOOSEMaterialPropertyToNEML2.h
@@ -18,7 +18,7 @@
  * Gather a MOOSE material property for insertion into the NEML2 model.
  */
 template <typename T, unsigned int state>
-class MOOSEMaterialPropertyToNEML2 : public MOOSEToNEML2Batched
+class MOOSEMaterialPropertyToNEML2 : public MOOSEToNEML2Batched<T>
 {
 public:
   static InputParameters validParams();
@@ -27,7 +27,7 @@ class MOOSEMaterialPropertyToNEML2 : public MOOSEToNEML2Batched
 
 #ifdef NEML2_ENABLED
 protected:
-  torch::Tensor convertQpMOOSEData() const override;
+  const MooseArray<T> & elemMOOSEData() const override { return _mat_prop.get(); }
 
   /// MOOSE material property to read data from
   const MaterialProperty<T> & _mat_prop;
@@ -41,4 +41,4 @@ class MOOSEMaterialPropertyToNEML2 : public MOOSEToNEML2Batched
 DefineMOOSEMaterialPropertyToNEML2Alias(Real, Real);
 DefineMOOSEMaterialPropertyToNEML2Alias(RankTwoTensor, RankTwoTensor);
 DefineMOOSEMaterialPropertyToNEML2Alias(SymmetricRankTwoTensor, SymmetricRankTwoTensor);
-DefineMOOSEMaterialPropertyToNEML2Alias(std::vector<Real>, StdVector);
+DefineMOOSEMaterialPropertyToNEML2Alias(RealVectorValue, RealVectorValue);
diff --git a/modules/solid_mechanics/include/neml2/userobjects/MOOSEToNEML2Batched.h b/modules/solid_mechanics/include/neml2/userobjects/MOOSEToNEML2Batched.h
index 496b788cd6cd..4485f274297c 100644
--- a/modules/solid_mechanics/include/neml2/userobjects/MOOSEToNEML2Batched.h
+++ b/modules/solid_mechanics/include/neml2/userobjects/MOOSEToNEML2Batched.h
@@ -12,6 +12,20 @@
 #include "MOOSEToNEML2.h"
 #include "ElementUserObject.h"
 
+/**
+ * @brief Generic gatherer for collecting "batched" MOOSE data for NEML2
+ *
+ * It is generic in the sense that it can be used for most MOOSE data types that take the form of
+ * MooseArray<T>.
+ *
+ * It is not so generic in the sense that the collected data is always a std::vector of
+ * MooseArray<T>, where the vector size is generally the number of elements this ElementUserObject
+ * operates on, and the MooseArray<T> size is generally the number of quadrature points in each
+ * element.
+ *
+ * @tparam T Type of the underlying MOOSE data, e.g., Real, SymmetricRankTwoTensor, etc.
+ */
+template <typename T>
 class MOOSEToNEML2Batched : public MOOSEToNEML2, public ElementUserObject
 {
 public:
@@ -32,17 +46,71 @@ class MOOSEToNEML2Batched : public MOOSEToNEML2, public ElementUserObject
 
   neml2::Tensor gatheredData() const override;
 
-  // The number of gathered data items (for setting the model batch size)
+  // The number of batches
   std::size_t size() const { return _buffer.size(); }
 
 protected:
-  /// Convert the underlying MOOSE data to a torch::Tensor
-  virtual torch::Tensor convertQpMOOSEData() const = 0;
+  /// MOOSE data for the current element
+  virtual const MooseArray<T> & elemMOOSEData() const = 0;
 
   /// Intermediate data buffer, filled during the element loop
-  std::vector<torch::Tensor> _buffer;
-
-  /// Current element's quadrature point indexing
-  unsigned int _qp;
+  std::vector<T> _buffer;
 #endif
 };
+
+template <typename T>
+InputParameters
+MOOSEToNEML2Batched<T>::validParams()
+{
+  auto params = MOOSEToNEML2::validParams();
+  params += ElementUserObject::validParams();
+
+  // Since we use the NEML2 model to evaluate the residual AND the Jacobian at the same time, we
+  // want to execute this user object only at execute_on = LINEAR (i.e. during residual evaluation).
+  // The NONLINEAR exec flag below is for computing Jacobian during automatic scaling.
+  ExecFlagEnum execute_options = MooseUtils::getDefaultExecFlagEnum();
+  execute_options = {EXEC_INITIAL, EXEC_LINEAR, EXEC_NONLINEAR};
+  params.set<ExecFlagEnum>("execute_on") = execute_options;
+
+  return params;
+}
+
+template <typename T>
+MOOSEToNEML2Batched<T>::MOOSEToNEML2Batched(const InputParameters & params)
+  : MOOSEToNEML2(params), ElementUserObject(params)
+{
+}
+
+#ifdef NEML2_ENABLED
+template <typename T>
+void
+MOOSEToNEML2Batched<T>::initialize()
+{
+  _buffer.clear();
+}
+
+template <typename T>
+void
+MOOSEToNEML2Batched<T>::execute()
+{
+  const auto & elem_data = this->elemMOOSEData();
+  for (auto i : index_range(elem_data))
+    _buffer.push_back(elem_data[i]);
+}
+
+template <typename T>
+void
+MOOSEToNEML2Batched<T>::threadJoin(const UserObject & uo)
+{
+  // append vectors
+  const auto & m2n = static_cast<const MOOSEToNEML2Batched<T> &>(uo);
+  _buffer.insert(_buffer.end(), m2n._buffer.begin(), m2n._buffer.end());
+}
+
+template <typename T>
+neml2::Tensor
+MOOSEToNEML2Batched<T>::gatheredData() const
+{
+  return NEML2Utils::fromBlob(_buffer);
+}
+#endif
diff --git a/modules/solid_mechanics/include/neml2/userobjects/MOOSEVariableToNEML2.h b/modules/solid_mechanics/include/neml2/userobjects/MOOSEVariableToNEML2.h
index 1a9f3ff614bf..c46bde8b4a64 100644
--- a/modules/solid_mechanics/include/neml2/userobjects/MOOSEVariableToNEML2.h
+++ b/modules/solid_mechanics/include/neml2/userobjects/MOOSEVariableToNEML2.h
@@ -15,7 +15,7 @@
  * Gather a MOOSE variable for insertion into the NEML2 model.
  */
 template <unsigned int state>
-class MOOSEVariableToNEML2Templ : public MOOSEToNEML2Batched
+class MOOSEVariableToNEML2Templ : public MOOSEToNEML2Batched<Real>
 {
 public:
   static InputParameters validParams();
@@ -24,7 +24,7 @@ class MOOSEVariableToNEML2Templ : public MOOSEToNEML2Batched
 
 #ifdef NEML2_ENABLED
 protected:
-  torch::Tensor convertQpMOOSEData() const override;
+  const MooseArray<Real> & elemMOOSEData() const override { return _moose_variable; }
 
   /// Coupled MOOSE variable to read data from
   const VariableValue & _moose_variable;
diff --git a/modules/solid_mechanics/include/neml2/userobjects/NEML2ModelExecutor.h b/modules/solid_mechanics/include/neml2/userobjects/NEML2ModelExecutor.h
index 2d7a6408ca44..4780682e1028 100644
--- a/modules/solid_mechanics/include/neml2/userobjects/NEML2ModelExecutor.h
+++ b/modules/solid_mechanics/include/neml2/userobjects/NEML2ModelExecutor.h
@@ -13,12 +13,6 @@
 #include "GeneralUserObject.h"
 #include "NEML2BatchIndexGenerator.h"
 
-#ifdef NEML2_ENABLED
-#include "neml2/tensors/LabeledVector.h"
-#include "neml2/tensors/LabeledMatrix.h"
-#include "neml2/models/Model.h"
-#endif
-
 class MOOSEToNEML2;
 
 /**
@@ -51,15 +45,15 @@ class NEML2ModelExecutor : public NEML2ModelInterface<GeneralUserObject>
   std::size_t getBatchIndex(dof_id_type elem_id) const;
 
   /// Get a reference(!) to the requested output view
-  const neml2::Tensor & getOutputView(const neml2::VariableName & output_name) const;
+  const neml2::Tensor & getOutput(const neml2::VariableName & output_name) const;
 
   /// Get a reference(!) to the requested output derivative view
-  const neml2::Tensor & getOutputDerivativeView(const neml2::VariableName & output_name,
-                                                const neml2::VariableName & input_name) const;
+  const neml2::Tensor & getOutputDerivative(const neml2::VariableName & output_name,
+                                            const neml2::VariableName & input_name) const;
 
   /// Get a reference(!) to the requested output parameter derivative view
-  const neml2::Tensor & getOutputParameterDerivativeView(const neml2::VariableName & output_name,
-                                                         const std::string & parameter_name) const;
+  const neml2::Tensor & getOutputParameterDerivative(const neml2::VariableName & output_name,
+                                                     const std::string & parameter_name) const;
 
   /// check if the output is fully computed and ready to be fetched
   bool outputReady() const { return _output_ready; }
@@ -92,14 +86,17 @@ class NEML2ModelExecutor : public NEML2ModelInterface<GeneralUserObject>
   /// flag that indicates if output data has been fully computed
   bool _output_ready;
 
-  /// The input vector of the material model
-  const neml2::LabeledVector & _in;
+  /// The model parameters to update (gathered from MOOSE)
+  std::map<std::string, neml2::Tensor> _model_params;
+
+  /// The input variables of the material model
+  neml2::ValueMap _in;
 
-  /// The output vector of the material model
-  const neml2::LabeledVector & _out;
+  /// The output variables of the material model
+  neml2::ValueMap _out;
 
-  /// The derivative of the output vector w.r.t. the input vector
-  const neml2::LabeledMatrix & _dout_din;
+  /// The derivative of the output variables w.r.t. the input variables
+  neml2::DerivMap _dout_din;
 
   // set of variables to skip
   std::set<neml2::VariableName> _skip_vars;
@@ -114,14 +111,13 @@ class NEML2ModelExecutor : public NEML2ModelInterface<GeneralUserObject>
   std::vector<const MOOSEToNEML2 *> _gatherers;
 
   /// set of output variables that were retrieved (by other objects)
-  mutable std::map<neml2::VariableName, neml2::Tensor> _retrieved_outputs;
+  mutable neml2::ValueMap _retrieved_outputs;
 
   /// set of derivatives that were retrieved (by other objects)
-  mutable std::map<std::pair<neml2::VariableName, neml2::VariableName>, neml2::Tensor>
-      _retrieved_derivatives;
+  mutable neml2::DerivMap _retrieved_derivatives;
 
   /// set of parameter derivatives that were retrieved (by other objects)
-  mutable std::map<std::pair<neml2::VariableName, std::string>, neml2::Tensor>
+  mutable std::map<neml2::VariableName, std::map<std::string, neml2::Tensor>>
       _retrieved_parameter_derivatives;
 #endif
 };
diff --git a/modules/solid_mechanics/include/neml2/utils/NEML2Utils.h b/modules/solid_mechanics/include/neml2/utils/NEML2Utils.h
index 19df14c1d6c5..c712929d474e 100644
--- a/modules/solid_mechanics/include/neml2/utils/NEML2Utils.h
+++ b/modules/solid_mechanics/include/neml2/utils/NEML2Utils.h
@@ -11,8 +11,9 @@
 
 #ifdef NEML2_ENABLED
 
-#include "neml2/models/Model.h"
 #include "neml2/misc/parser_utils.h"
+#include "neml2/tensors/tensors.h"
+#include "neml2/models/LabeledAxisAccessor.h"
 #include "RankTwoTensor.h"
 #include "RankFourTensor.h"
 #include "SymmetricRankTwoTensor.h"
@@ -22,19 +23,13 @@
 #endif
 
 #include "InputParameters.h"
+#include "MooseArray.h"
+#include <type_traits>
 
 class MooseObject;
 class Action;
 class SubProblem;
 
-#ifdef NEML2_ENABLED
-namespace neml2
-{
-/// Pretty print the summary of a NEML2 model
-std::ostream & operator<<(std::ostream &, const Model &);
-}
-#endif
-
 namespace NEML2Utils
 {
 
@@ -45,76 +40,101 @@ void assertVariable(const neml2::VariableName &);
 /// Assert that the NEML2 variable name sits on either the old_forces or the old_state subaxis
 void assertOldVariable(const neml2::VariableName &);
 
-/// Convert a MOOSE data structure to its NEML2 counterpart
-template <typename T>
-neml2::Tensor toNEML2(const T &);
+/// Parse a raw string into NEML2 variable name
+neml2::VariableName parseVariableName(const std::string &);
 
-/**
- * Convert a wrapped (batched) MOOSE data structure to its NEML2 counterpart
- * The wrapper should implement size()
- */
 template <typename T>
-neml2::Tensor toNEML2Batched(const T & data);
-
-// @{ Template specializations
-template <>
-neml2::Tensor toNEML2(const Real & v);
+struct Layout
+{
+};
 template <>
-neml2::Tensor toNEML2(const RankTwoTensor & r2t);
+struct Layout<Real>
+{
+  static constexpr std::array<neml2::Size, 0> shape{};
+  static constexpr std::array<neml2::Size, 1> strides{1};
+};
 template <>
-neml2::Tensor toNEML2(const std::vector<Real> & v);
-// @}
-
-/**
- * Decompose a `std::vector<std::tuple<Args...>>` into a `std::tuple<std::vector<Args>...>`
- * The original data structure (the batched tuple) is what BatchMaterial gathered for us, but it
- * allows for inhomogeneous batches, i.e., the batch data structure could change from batch to
- * batch. The batched tuple is efficient in terms of data gathering, but it is disadvantageous when
- * we want to convert and set all batches at once. This method essentially "homogenizes" the batched
- * tuple to a tuple of batched data structures.
- */
-template <typename... Args>
-std::tuple<std::vector<Args>...>
-homogenizeBatchedTuple(const std::vector<std::tuple<Args...>> & from);
-
-template <size_t I = 0, typename... Args>
-void homogenizeBatchedTupleInner(const std::vector<std::tuple<Args...>> & from,
-                                 std::tuple<std::vector<Args>...> & to);
-
-/// Convert a NEML2 data structure to its MOOSE counterpart
-template <typename T>
-T toMOOSE(const neml2::Tensor &);
-
-// @{ Template specializations
+struct Layout<RealVectorValue>
+{
+  static constexpr std::array<neml2::Size, 1> shape{3};
+  static constexpr std::array<neml2::Size, 2> strides{3, 1};
+};
 template <>
-Real toMOOSE(const neml2::Tensor & t);
+struct Layout<RankTwoTensor>
+{
+  static constexpr std::array<neml2::Size, 2> shape{3, 3};
+  static constexpr std::array<neml2::Size, 3> strides{9, 3, 1};
+};
 template <>
-SymmetricRankTwoTensor toMOOSE(const neml2::Tensor & t);
+struct Layout<SymmetricRankTwoTensor>
+{
+  static constexpr std::array<neml2::Size, 1> shape{6};
+  static constexpr std::array<neml2::Size, 2> strides{6, 1};
+};
 template <>
-std::vector<Real> toMOOSE(const neml2::Tensor & t);
+struct Layout<RankFourTensor>
+{
+  static constexpr std::array<neml2::Size, 4> shape{3, 3, 3, 3};
+  static constexpr std::array<neml2::Size, 5> strides{81, 27, 9, 3, 1};
+};
 template <>
-SymmetricRankFourTensor toMOOSE(const neml2::Tensor & t);
-// @}
+struct Layout<SymmetricRankFourTensor>
+{
+  static constexpr std::array<neml2::Size, 2> shape{6, 6};
+  static constexpr std::array<neml2::Size, 3> strides{36, 6, 1};
+};
 
 /**
- * Convert a MOOSE data structure to its NEML2 counterpart and copy the values into a NEML2
- * LabeledVector
+ * @brief Map from MooseArray to neml2::Tensor without copying the data
+ *
+ * This method is used in gatherers which gather data from MOOSE as input variables to the NEML2
+ * material model. So in theory, we only need to overload MOOSE types that can potentially be used
+ * as input variables.
  */
-template <size_t I = 0, typename T, typename... Ts>
-void set(neml2::LabeledVector & v,
-         const std::vector<neml2::VariableName> & indices,
-         const T * t0,
-         const Ts *... t);
+template <typename T>
+neml2::Tensor
+fromBlob(const std::vector<T> & data)
+{
+  // The const_cast is fine because torch works with non-const ptr so that it can optionally handle
+  // deallocation. But we are not going to let torch do that.
+  const auto torch_tensor =
+      torch::from_blob(const_cast<T *>(data.data()),
+                       neml2::utils::add_shapes(data.size(), Layout<T>::shape),
+                       torch::TensorOptions().dtype(torch::kFloat64));
+  return neml2::Tensor(torch_tensor, 1);
+}
 
 /**
- * Convert a wrapped (batched) MOOSE data structure to its NEML2 counterpart and copy the values
- * into a NEML2 LabeledVector
+ * @brief Directly copy a contiguous chunk of memory of a torch::Tensor to a MOOSE data of type T
+ *
+ * This assumes the torch::Tensor and T have the same layout, for example both row-major
+ * with T = RankTwoTensor. If the layouts are different, we may need to reshape/reorder/transpose
+ * the torch::Tensor before copying.
+ *
+ * For this method to work,
+ * 1. the address of \p dest must align with the first element of the data,
+ * 2. the number of elements in \p dest must match the number of elements in \p src,
+ * 3. the \p src tensor must be of type torch::kFloat64, and
+ * 4. data in \p dest must be reinterpretable as Real.
  */
-template <size_t I = 0, typename T, typename... Ts>
-void setBatched(neml2::LabeledVector & v,
-                const std::vector<neml2::VariableName> & indices,
-                const T * t0,
-                const Ts *... t);
+template <typename T>
+void
+copyTensorToMOOSEData(const torch::Tensor & src, T & dest)
+{
+  if (src.dtype() != torch::kFloat64)
+    mooseError(
+        "Cannot copy torch::Tensor with dtype ", src.dtype(), " into ", demangle(typeid(T).name()));
+  if (src.numel() != Layout<T>::strides[0])
+    mooseError("Cannot copy torch::Tensor with shape ",
+               src.sizes(),
+               " into ",
+               demangle(typeid(T).name()),
+               " with different number of elements.");
+  auto dest_tensor = torch::from_blob(reinterpret_cast<Real *>(&dest),
+                                      Layout<T>::shape,
+                                      torch::TensorOptions().dtype(torch::kFloat64));
+  dest_tensor.copy_(src.reshape(Layout<T>::shape));
+}
 
 static std::string NEML2_help_message = R""""(
 ==============================================================================
@@ -130,78 +150,6 @@ To debug NEML2 related issues:
 ==============================================================================
 )"""";
 
-////////////////////////////////////////////////////////////////////////////////
-// Implementations
-////////////////////////////////////////////////////////////////////////////////
-template <typename T>
-neml2::Tensor
-toNEML2Batched(const T & data)
-{
-  std::vector<torch::Tensor> res(data.size());
-  for (const auto i : index_range(data))
-    res[i] = toNEML2<typename T::value_type>(data[i]);
-  return neml2::Tensor(torch::stack(res, 0), 1);
-}
-
-template <typename... Args>
-std::tuple<std::vector<Args>...>
-homogenizeBatchedTuple(const std::vector<std::tuple<Args...>> & from)
-{
-  std::tuple<std::vector<Args>...> to;
-  homogenizeBatchedTupleInner(from, to);
-  return to;
-}
-
-template <size_t I, typename... Args>
-void
-homogenizeBatchedTupleInner(const std::vector<std::tuple<Args...>> & from,
-                            std::tuple<std::vector<Args>...> & to)
-{
-  typedef typename std::tuple_element<I, std::tuple<Args...>>::type Arg;
-  std::vector<Arg> toi;
-  std::transform(from.cbegin(),
-                 from.cend(),
-                 std::back_inserter(toi),
-                 [](const auto & b) { return std::get<I>(b); });
-  std::get<I>(to) = toi;
-
-  // recursively act on the rest of the args
-  if constexpr ((I + 1) < sizeof...(Args))
-    homogenizeBatchedTupleInner<I + 1>(from, to);
-}
-
-template <size_t I, typename T, typename... Ts>
-void
-set(neml2::LabeledVector & v,
-    const std::vector<neml2::VariableName> & indices,
-    const T * t0,
-    const Ts *... t)
-{
-  if (t0)
-    v.base_index_put_(indices[I], toNEML2(*t0));
-
-  // Recursively act on the rest of the data
-  // The compiler should be able to easily deduce the rest of the template parameters...
-  if constexpr (sizeof...(Ts) > 0)
-    set<I + 1>(v, indices, t...);
-}
-
-template <size_t I, typename T, typename... Ts>
-void
-setBatched(neml2::LabeledVector & v,
-           const std::vector<neml2::VariableName> & indices,
-           const T * t0,
-           const Ts *... t)
-{
-  if (t0)
-    v.base_index_put_(indices[I], toNEML2Batched(*t0));
-
-  // recursively act on the rest of the data
-  // The compiler should be able to easily deduce the rest of the template parameters...
-  if constexpr (sizeof...(Ts) > 0)
-    setBatched<I + 1>(v, indices, t...);
-}
-
 #endif // NEML2_ENABLED
 
 /// Determine whether the NEML2 material model should be evaluated
diff --git a/modules/solid_mechanics/src/materials/lagrangian/ComputeLagrangianCauchyCustomStress.C b/modules/solid_mechanics/src/materials/lagrangian/ComputeLagrangianCauchyCustomStress.C
new file mode 100644
index 000000000000..6bae30765224
--- /dev/null
+++ b/modules/solid_mechanics/src/materials/lagrangian/ComputeLagrangianCauchyCustomStress.C
@@ -0,0 +1,42 @@
+//* 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 "ComputeLagrangianCauchyCustomStress.h"
+
+registerMooseObject("SolidMechanicsApp", ComputeLagrangianCauchyCustomStress);
+
+InputParameters
+ComputeLagrangianCauchyCustomStress::validParams()
+{
+  InputParameters params = ComputeLagrangianStressCauchy::validParams();
+
+  params.addClassDescription("Custom stress update providing the Cauchy stress");
+
+  params.addRequiredParam<MaterialPropertyName>("custom_cauchy_stress",
+                                                "The name of the Cauchy stress.");
+  params.addRequiredParam<MaterialPropertyName>("custom_cauchy_jacobian",
+                                                "The name of the Cauchy stress Jacobian.");
+
+  return params;
+}
+
+ComputeLagrangianCauchyCustomStress::ComputeLagrangianCauchyCustomStress(
+    const InputParameters & parameters)
+  : ComputeLagrangianStressCauchy(parameters),
+    _custom_stress(getMaterialProperty<RankTwoTensor>("custom_cauchy_stress")),
+    _custom_jacobian(getMaterialProperty<RankFourTensor>("custom_cauchy_jacobian"))
+{
+}
+
+void
+ComputeLagrangianCauchyCustomStress::computeQpCauchyStress()
+{
+  _cauchy_stress[_qp] = _custom_stress[_qp];
+  _cauchy_jacobian[_qp] = _custom_jacobian[_qp];
+}
diff --git a/modules/solid_mechanics/src/materials/lagrangian/ComputeLagrangianStrainBase.C b/modules/solid_mechanics/src/materials/lagrangian/ComputeLagrangianStrainBase.C
index 57204926a86c..9445be440d19 100644
--- a/modules/solid_mechanics/src/materials/lagrangian/ComputeLagrangianStrainBase.C
+++ b/modules/solid_mechanics/src/materials/lagrangian/ComputeLagrangianStrainBase.C
@@ -51,6 +51,8 @@ ComputeLagrangianStrainBase<G>::ComputeLagrangianStrainBase(const InputParameter
     _mechanical_strain(declareProperty<RankTwoTensor>(_base_name + "mechanical_strain")),
     _mechanical_strain_old(getMaterialPropertyOld<RankTwoTensor>(_base_name + "mechanical_strain")),
     _strain_increment(declareProperty<RankTwoTensor>(_base_name + "strain_increment")),
+    _spatial_velocity_increment(
+        declareProperty<RankTwoTensor>(_base_name + "spatial_velocity_increment")),
     _vorticity_increment(declareProperty<RankTwoTensor>(_base_name + "vorticity_increment")),
     _F_ust(declareProperty<RankTwoTensor>(_base_name + "unstabilized_deformation_gradient")),
     _F_avg(declareProperty<RankTwoTensor>(_base_name + "average_deformation_gradient")),
@@ -148,6 +150,9 @@ ComputeLagrangianStrainBase<G>::computeQpIncrementalStrains(const RankTwoTensor
   // Increment the mechanical strain
   _mechanical_strain[_qp] = _mechanical_strain_old[_qp] + _strain_increment[_qp];
 
+  // Yes, this does make sense to do it here
+  _spatial_velocity_increment[_qp] = _vorticity_increment[_qp] + _strain_increment[_qp];
+
   // Faked rotation increment for ComputeStressBase materials
   _rotation_increment[_qp] = RankTwoTensor::Identity();
 }
diff --git a/modules/solid_mechanics/src/neml2/actions/NEML2Action.C b/modules/solid_mechanics/src/neml2/actions/NEML2Action.C
index 1c83a1dfba68..48eb0cfa8573 100644
--- a/modules/solid_mechanics/src/neml2/actions/NEML2Action.C
+++ b/modules/solid_mechanics/src/neml2/actions/NEML2Action.C
@@ -29,7 +29,7 @@ const std::map<neml2::TensorType, std::string> tensor_type_map = {
     {neml2::TensorType::kR2, "RankTwoTensor"},
     {neml2::TensorType::kSSR4, "SymmetricRankFourTensor"},
     {neml2::TensorType::kR4, "RankFourTensor"},
-};
+    {neml2::TensorType::kRot, "RealVectorValue"}};
 // NEML2 (output, input) type --> NEML2 derivative type
 const std::map<std::pair<neml2::TensorType, neml2::TensorType>, neml2::TensorType> deriv_type_map =
     {
@@ -81,13 +81,21 @@ NEML2Action::NEML2Action(const InputParameters & params)
   // Apply parameters under the common area, i.e., under [NEML2]
   const auto & all_params = _app.getInputParameterWarehouse().getInputParameters();
   auto & sub_block_params = *(all_params.find(uniqueActionName())->second.get());
-  auto common_block = _awh.getActions<NEML2ActionCommon>();
-  mooseAssert(common_block.size() == 1, "There must exist one and only one common NEML2 action.");
-  sub_block_params.applyParameters(common_block[0]->parameters());
+  const auto & common_action = getCommonAction();
+  sub_block_params.applyParameters(common_action.parameters());
 
   // verbosity
   _verbose = getParam<bool>("verbose");
 
+  // Set up optional output variable initialization
+  auto init_vars = getParam<std::vector<MaterialPropertyName>>("initialize_outputs");
+  auto init_vals = getParam<std::vector<MaterialPropertyName>>("initialize_output_values");
+  if (init_vars.size() != init_vals.size())
+    paramError("initialize_outputs",
+               "initialize_outputs should have the same length as initialize_output_values");
+  for (auto i : index_range(init_vars))
+    _initialize_output_values[init_vars[i]] = init_vals[i];
+
   // Set up additional outputs for each requested material property
   auto outputs = getParam<std::vector<MaterialPropertyName>>("export_outputs");
   auto output_targets = getParam<std::vector<std::vector<OutputName>>>("export_output_targets");
@@ -98,6 +106,14 @@ NEML2Action::NEML2Action(const InputParameters & params)
     _export_output_targets[outputs[i]] = output_targets[i];
 }
 
+const NEML2ActionCommon &
+NEML2Action::getCommonAction() const
+{
+  auto common_block = _awh.getActions<NEML2ActionCommon>();
+  mooseAssert(common_block.size() == 1, "There must exist one and only one common NEML2 action.");
+  return *common_block[0];
+}
+
 #ifndef NEML2_ENABLED
 
 void
@@ -113,24 +129,7 @@ NEML2Action::act()
   if (_current_task == "add_user_object")
   {
     // Get the NEML2 model so that we can introspect variable tensor types
-    auto & model = neml2::get_model(getParam<std::string>("model"));
-
-    // List inputs, outputs, and parameters of the model
-    if (_verbose)
-    {
-      _console << COLOR_YELLOW << "*** NEML2 MATERIAL MODEL SUMMARY BEGIN ***" << COLOR_DEFAULT
-               << std::endl;
-      _console << "Input file location:       " << parameters().blockLocation() << std::endl;
-      _console << "Input file path:           " << parameters().blockFullpath() << std::endl;
-      _console << "Material model Name:       " << getParam<std::string>("model") << std::endl;
-      _console << "Evaluation device:         " << getParam<std::string>("device") << std::endl;
-      _console << "Automatic differentiation: "
-               << (getParam<bool>("enable_AD") ? "Enabled" : "Disabled") << std::endl;
-      _console << "Material model structure:  " << std::endl;
-      _console << model;
-      _console << COLOR_YELLOW << "*** NEML2 MATERIAL MODEL SUMMARY END ***\n"
-               << COLOR_DEFAULT << std::endl;
-    }
+    const auto & model = neml2::get_model(getParam<std::string>("model"));
 
     setupInputMappings(model);
     setupParameterMappings(model);
@@ -138,6 +137,9 @@ NEML2Action::act()
     setupDerivativeMappings(model);
     setupParameterDerivativeMappings(model);
 
+    if (_verbose)
+      printSummary(model);
+
     // MOOSEToNEML2 input gatherers
     std::vector<UserObjectName> gatherers;
     for (const auto & input : _inputs)
@@ -146,8 +148,10 @@ NEML2Action::act()
       {
         auto obj_name = "__moose(" + input.moose.name + ")->neml2(" +
                         neml2::utils::stringify(input.neml2.name) + ")_" + name() + "__";
+        if (!tensor_type_map.count(input.neml2.type))
+          mooseError("NEML2 type ", input.neml2.type, " not yet mapped to MOOSE");
         auto obj_moose_type = tensor_type_map.at(input.neml2.type) + "MaterialProperty";
-        if (input.neml2.name.start_with("old_forces") || input.neml2.name.start_with("old_state"))
+        if (input.neml2.name.is_old_force() || input.neml2.name.is_old_state())
           obj_moose_type = "Old" + obj_moose_type;
         auto obj_type = "MOOSE" + obj_moose_type + "ToNEML2";
         auto obj_params = _factory.getValidParams(obj_type);
@@ -162,7 +166,7 @@ NEML2Action::act()
         auto obj_name = "__moose(" + input.moose.name + ")->neml2(" +
                         neml2::utils::stringify(input.neml2.name) + ")_" + name() + "__";
         std::string obj_moose_type = "Variable";
-        if (input.neml2.name.start_with("old_forces") || input.neml2.name.start_with("old_state"))
+        if (input.neml2.name.is_old_force() || input.neml2.name.is_old_state())
           obj_moose_type = "Old" + obj_moose_type;
         auto obj_type = "MOOSE" + obj_moose_type + "ToNEML2";
         auto obj_params = _factory.getValidParams(obj_type);
@@ -177,7 +181,7 @@ NEML2Action::act()
         auto obj_name = "__moose(" + input.moose.name + ")->neml2(" +
                         neml2::utils::stringify(input.neml2.name) + ")" + name() + "__";
         auto obj_moose_type = std::string("Postprocessor");
-        if (input.neml2.name.start_with("old_forces") || input.neml2.name.start_with("old_state"))
+        if (input.neml2.name.is_old_force() || input.neml2.name.is_old_state())
           obj_moose_type = "Old" + obj_moose_type;
         auto obj_type = "MOOSE" + obj_moose_type + "ToNEML2";
         auto obj_params = _factory.getValidParams(obj_type);
@@ -200,6 +204,8 @@ NEML2Action::act()
       {
         auto obj_name =
             "__moose(" + param.moose.name + ")->neml2(" + param.neml2.name + ")_" + name() + "__";
+        if (!tensor_type_map.count(param.neml2.type))
+          mooseError("NEML2 type ", param.neml2.type, " not yet mapped to MOOSE");
         auto obj_moose_type = tensor_type_map.at(param.neml2.type);
         auto obj_type = "MOOSE" + obj_moose_type + "MaterialPropertyToNEML2";
         auto obj_params = _factory.getValidParams(obj_type);
@@ -268,12 +274,17 @@ NEML2Action::act()
       {
         auto obj_name = "__neml2(" + neml2::utils::stringify(output.neml2.name) + ")->moose(" +
                         output.moose.name + ")_" + name() + "__";
+        if (!tensor_type_map.count(output.neml2.type))
+          mooseError("NEML2 type ", output.neml2.type, " not yet mapped to MOOSE");
         auto obj_type = "NEML2ToMOOSE" + tensor_type_map.at(output.neml2.type) + "MaterialProperty";
         auto obj_params = _factory.getValidParams(obj_type);
         obj_params.set<UserObjectName>("neml2_executor") = _executor_name;
         obj_params.set<MaterialPropertyName>("to_moose") = output.moose.name;
         obj_params.set<std::string>("from_neml2") = neml2::utils::stringify(output.neml2.name);
         obj_params.set<std::vector<SubdomainName>>("block") = _block;
+        if (_initialize_output_values.count(output.moose.name))
+          obj_params.set<MaterialPropertyName>("moose_material_property_init") =
+              _initialize_output_values[output.moose.name];
         if (_export_output_targets.count(output.moose.name))
           obj_params.set<std::vector<OutputName>>("outputs") =
               _export_output_targets[output.moose.name];
@@ -293,10 +304,16 @@ NEML2Action::act()
         auto obj_name = "__neml2(d(" + neml2::utils::stringify(deriv.neml2.y.name) + ")/d(" +
                         neml2::utils::stringify(deriv.neml2.x.name) + "))->moose(" +
                         deriv.moose.name + ")_" + name() + "__";
-        auto obj_type =
-            "NEML2ToMOOSE" +
-            tensor_type_map.at(deriv_type_map.at({deriv.neml2.y.type, deriv.neml2.x.type})) +
-            "MaterialProperty";
+        if (!deriv_type_map.count({deriv.neml2.y.type, deriv.neml2.x.type}))
+          mooseError("NEML2 derivative type for d(",
+                     deriv.neml2.y.type,
+                     ")/d(",
+                     deriv.neml2.x.type,
+                     ") not yet mapped to MOOSE");
+        auto deriv_type = deriv_type_map.at({deriv.neml2.y.type, deriv.neml2.x.type});
+        if (!tensor_type_map.count(deriv_type))
+          mooseError("NEML2 type ", deriv_type, " not yet mapped to MOOSE");
+        auto obj_type = "NEML2ToMOOSE" + tensor_type_map.at(deriv_type) + "MaterialProperty";
         auto obj_params = _factory.getValidParams(obj_type);
         obj_params.set<UserObjectName>("neml2_executor") = _executor_name;
         obj_params.set<MaterialPropertyName>("to_moose") = deriv.moose.name;
@@ -323,10 +340,16 @@ NEML2Action::act()
         auto obj_name = "__neml2(d(" + neml2::utils::stringify(param_deriv.neml2.y.name) + ")/d(" +
                         param_deriv.neml2.x.name + "))->moose(" + param_deriv.moose.name + ")_" +
                         name() + "__";
-        auto obj_type = "NEML2ToMOOSE" +
-                        tensor_type_map.at(deriv_type_map.at(
-                            {param_deriv.neml2.y.type, param_deriv.neml2.x.type})) +
-                        "MaterialProperty";
+        if (!deriv_type_map.count({param_deriv.neml2.y.type, param_deriv.neml2.x.type}))
+          mooseError("NEML2 derivative type for d(",
+                     param_deriv.neml2.y.type,
+                     ")/d(",
+                     param_deriv.neml2.x.type,
+                     ") not yet mapped to MOOSE");
+        auto deriv_type = deriv_type_map.at({param_deriv.neml2.y.type, param_deriv.neml2.x.type});
+        if (!tensor_type_map.count(deriv_type))
+          mooseError("NEML2 type ", deriv_type, " not yet mapped to MOOSE");
+        auto obj_type = "NEML2ToMOOSE" + tensor_type_map.at(deriv_type) + "MaterialProperty";
         auto obj_params = _factory.getValidParams(obj_type);
         obj_params.set<UserObjectName>("neml2_executor") = _executor_name;
         obj_params.set<MaterialPropertyName>("to_moose") = param_deriv.moose.name;
@@ -350,22 +373,16 @@ NEML2Action::act()
 void
 NEML2Action::setupInputMappings(const neml2::Model & model)
 {
-  const auto moose_input_types =
-      getParam<MultiMooseEnum>("moose_input_types").getSetValueIDs<MOOSEIOType>();
-  const auto moose_inputs = getParam<std::vector<std::string>>("moose_inputs");
-  const auto neml2_inputs = getParam<std::vector<std::string>>("neml2_inputs");
-
-  if (moose_inputs.size() != moose_input_types.size())
-    paramError("moose_inputs", "moose_inputs must have the same length as moose_input_types");
-  if (moose_inputs.size() != neml2_inputs.size())
-    paramError("moose_inputs", "moose_inputs must have the same length as neml2_inputs");
+  const auto [moose_input_types, moose_inputs, neml2_inputs] =
+      getInputParameterMapping<MOOSEIOType, std::string, std::string>(
+          "moose_input_types", "moose_inputs", "neml2_inputs");
 
   for (auto i : index_range(moose_inputs))
   {
-    auto neml2_input = neml2::utils::parse<neml2::VariableName>(neml2_inputs[i]);
+    auto neml2_input = NEML2Utils::parseVariableName(neml2_inputs[i]);
     _inputs.push_back({
         {moose_inputs[i], moose_input_types[i]},
-        {neml2_input, model.input_type(neml2_input)},
+        {neml2_input, model.input_variable(neml2_input).type()},
     });
   }
 }
@@ -373,17 +390,9 @@ NEML2Action::setupInputMappings(const neml2::Model & model)
 void
 NEML2Action::setupParameterMappings(const neml2::Model & model)
 {
-  const auto moose_param_types =
-      getParam<MultiMooseEnum>("moose_parameter_types").getSetValueIDs<MOOSEIOType>();
-  const auto moose_params = getParam<std::vector<std::string>>("moose_parameters");
-  const auto neml2_params = getParam<std::vector<std::string>>("neml2_parameters");
-
-  if (moose_params.size() != moose_param_types.size())
-    paramError("moose_parameters",
-               "moose_params must have the same length as moose_parameter_types");
-  if (moose_params.size() != neml2_params.size())
-    paramError("moose_parameters",
-               "moose_parameters must have the same length as neml2_parameters");
+  const auto [moose_param_types, moose_params, neml2_params] =
+      getInputParameterMapping<MOOSEIOType, std::string, std::string>(
+          "moose_parameter_types", "moose_parameters", "neml2_parameters");
 
   for (auto i : index_range(moose_params))
     _params.push_back({{moose_params[i], moose_param_types[i]},
@@ -393,22 +402,16 @@ NEML2Action::setupParameterMappings(const neml2::Model & model)
 void
 NEML2Action::setupOutputMappings(const neml2::Model & model)
 {
-  const auto moose_output_types =
-      getParam<MultiMooseEnum>("moose_output_types").getSetValueIDs<MOOSEIOType>();
-  const auto moose_outputs = getParam<std::vector<std::string>>("moose_outputs");
-  const auto neml2_outputs = getParam<std::vector<std::string>>("neml2_outputs");
-
-  if (moose_outputs.size() != moose_output_types.size())
-    paramError("moose_outputs", "moose_outputs must have the same length as moose_output_types");
-  if (moose_outputs.size() != neml2_outputs.size())
-    paramError("moose_outputs", "moose_outputs must have the same length as neml2_outputs");
+  const auto [moose_output_types, moose_outputs, neml2_outputs] =
+      getInputParameterMapping<MOOSEIOType, std::string, std::string>(
+          "moose_output_types", "moose_outputs", "neml2_outputs");
 
   for (auto i : index_range(moose_outputs))
   {
-    auto neml2_output = neml2::utils::parse<neml2::VariableName>(neml2_outputs[i]);
+    auto neml2_output = NEML2Utils::parseVariableName(neml2_outputs[i]);
     _outputs.push_back({
         {moose_outputs[i], moose_output_types[i]},
-        {neml2_output, model.output_type(neml2_output)},
+        {neml2_output, model.output_variable(neml2_output).type()},
     });
   }
 }
@@ -416,28 +419,21 @@ NEML2Action::setupOutputMappings(const neml2::Model & model)
 void
 NEML2Action::setupDerivativeMappings(const neml2::Model & model)
 {
-  const auto moose_deriv_types =
-      getParam<MultiMooseEnum>("moose_derivative_types").getSetValueIDs<MOOSEIOType>();
-  const auto moose_derivs = getParam<std::vector<std::string>>("moose_derivatives");
-  const auto neml2_derivs = getParam<std::vector<std::vector<std::string>>>("neml2_derivatives");
-
-  if (moose_derivs.size() != moose_deriv_types.size())
-    paramError("moose_derivatives",
-               "moose_derivatives must have the same length as moose_derivative_types");
-  if (moose_derivs.size() != neml2_derivs.size())
-    paramError("moose_derivatives",
-               "moose_derivatives must have the same length as neml2_derivatives");
+  const auto [moose_deriv_types, moose_derivs, neml2_derivs] =
+      getInputParameterMapping<MOOSEIOType, std::string, std::vector<std::string>>(
+          "moose_derivative_types", "moose_derivatives", "neml2_derivatives");
 
   for (auto i : index_range(moose_derivs))
   {
     if (neml2_derivs[i].size() != 2)
       paramError("neml2_derivatives", "The length of each pair in neml2_derivatives must be 2.");
 
-    auto neml2_y = neml2::utils::parse<neml2::VariableName>(neml2_derivs[i][0]);
-    auto neml2_x = neml2::utils::parse<neml2::VariableName>(neml2_derivs[i][1]);
+    auto neml2_y = NEML2Utils::parseVariableName(neml2_derivs[i][0]);
+    auto neml2_x = NEML2Utils::parseVariableName(neml2_derivs[i][1]);
     _derivs.push_back({
         {moose_derivs[i], moose_deriv_types[i]},
-        {{neml2_y, model.output_type(neml2_y)}, {neml2_x, model.input_type(neml2_x)}},
+        {{neml2_y, model.output_variable(neml2_y).type()},
+         {neml2_x, model.input_variable(neml2_x).type()}},
     });
   }
 }
@@ -445,20 +441,11 @@ NEML2Action::setupDerivativeMappings(const neml2::Model & model)
 void
 NEML2Action::setupParameterDerivativeMappings(const neml2::Model & model)
 {
-  const auto moose_param_deriv_types =
-      getParam<MultiMooseEnum>("moose_parameter_derivative_types").getSetValueIDs<MOOSEIOType>();
-  const auto moose_param_derivs = getParam<std::vector<std::string>>("moose_parameter_derivatives");
-  const auto neml2_param_derivs =
-      getParam<std::vector<std::vector<std::string>>>("neml2_parameter_derivatives");
-
-  if (moose_param_derivs.size() != moose_param_deriv_types.size())
-    paramError("moose_parameter_derivatives",
-               "moose_parameter_derivatives must have the same length as "
-               "moose_parameter_derivative_types");
-  if (moose_param_derivs.size() != neml2_param_derivs.size())
-    paramError(
-        "moose_parameter_derivatives",
-        "moose_parameter_derivatives must have the same length as neml2_parameter_derivatives");
+  const auto [moose_param_deriv_types, moose_param_derivs, neml2_param_derivs] =
+      getInputParameterMapping<MOOSEIOType, std::string, std::vector<std::string>>(
+          "moose_parameter_derivative_types",
+          "moose_parameter_derivatives",
+          "neml2_parameter_derivatives");
 
   for (auto i : index_range(moose_param_derivs))
   {
@@ -466,12 +453,104 @@ NEML2Action::setupParameterDerivativeMappings(const neml2::Model & model)
       paramError("neml2_parameter_derivatives",
                  "The length of each pair in neml2_parameter_derivatives must be 2.");
 
-    auto neml2_y = neml2::utils::parse<neml2::VariableName>(neml2_param_derivs[i][0]);
+    auto neml2_y = NEML2Utils::parseVariableName(neml2_param_derivs[i][0]);
     auto neml2_x = neml2_param_derivs[i][1];
     _param_derivs.push_back({
         {moose_param_derivs[i], moose_param_deriv_types[i]},
-        {{neml2_y, model.output_type(neml2_y)}, {neml2_x, model.get_parameter(neml2_x).type()}},
+        {{neml2_y, model.output_variable(neml2_y).type()},
+         {neml2_x, model.get_parameter(neml2_x).type()}},
     });
   }
 }
+
+void
+NEML2Action::printSummary(const neml2::Model & model) const
+{
+  // Save formatting of the output stream so that we can restore it later
+  const auto flags = _console.flags();
+
+  // Default width for the summary
+  const int width = 79;
+
+  _console << std::endl;
+  _console << COLOR_CYAN << std::setw(width) << std::setfill('*') << std::left
+           << "NEML2 MATERIAL MODEL SUMMARY BEGIN " << std::setfill(' ') << COLOR_DEFAULT
+           << std::endl;
+
+  // Metadata
+  _console << "NEML2 input file location: " << getCommonAction().fname() << std::endl;
+  _console << "NEML2 action path:         " << parameters().blockFullpath() << std::endl;
+
+  // List inputs, outputs, and parameters of the model
+  _console << COLOR_CYAN << std::setw(width) << std::setfill('-') << std::left
+           << "Material model structure " << std::setfill(' ') << COLOR_DEFAULT << std::endl;
+  _console << model;
+
+  // List transfer between MOOSE and NEML2
+  _console << COLOR_CYAN << std::setw(width) << std::setfill('-') << std::left
+           << "Transfer between MOOSE and NEML2 " << std::setfill(' ') << COLOR_DEFAULT
+           << std::endl;
+
+  // Figure out the longest name length so that we could align the arrows
+  const auto max_moose_name_length = getLongestMOOSEName();
+
+  // List input transfer, MOOSE -> NEML2
+  for (const auto & input : _inputs)
+  {
+    _console << std::setw(max_moose_name_length) << std::right
+             << (input.neml2.name.is_old_force() || input.neml2.name.is_old_state()
+                     ? ("(old) " + input.moose.name)
+                     : input.moose.name)
+             << " --> " << input.neml2.name << std::endl;
+  }
+
+  // List parameter transfer, MOOSE -> NEML2
+  for (const auto & param : _params)
+    _console << std::setw(max_moose_name_length) << std::right << param.moose.name << " --> "
+             << param.neml2.name << std::endl;
+
+  // List output transfer, NEML2 -> MOOSE
+  for (const auto & output : _outputs)
+    _console << std::setw(max_moose_name_length) << std::right << output.moose.name << " <-- "
+             << output.neml2.name << std::endl;
+
+  // List derivative transfer, NEML2 -> MOOSE
+  for (const auto & deriv : _derivs)
+    _console << std::setw(max_moose_name_length) << std::right << deriv.moose.name << " <-- d("
+             << deriv.neml2.y.name << ")/d(" << deriv.neml2.x.name << ")" << std::endl;
+
+  // List parameter derivative transfer, NEML2 -> MOOSE
+  for (const auto & param_deriv : _param_derivs)
+    _console << std::setw(max_moose_name_length) << std::right << param_deriv.moose.name
+             << " <-- d(" << param_deriv.neml2.y.name << ")/d(" << param_deriv.neml2.x.name << ")"
+             << std::endl;
+
+  _console << COLOR_CYAN << std::setw(width) << std::setfill('*') << std::left
+           << "NEML2 MATERIAL MODEL SUMMARY END " << std::setfill(' ') << COLOR_DEFAULT << std::endl
+           << std::endl;
+
+  // Restore the formatting of the output stream
+  _console.flags(flags);
+}
+
+std::size_t
+NEML2Action::getLongestMOOSEName() const
+{
+  std::size_t max_moose_name_length = 0;
+  for (const auto & input : _inputs)
+    max_moose_name_length =
+        std::max(max_moose_name_length,
+                 input.neml2.name.is_old_force() || input.neml2.name.is_old_state()
+                     ? input.moose.name.size() + 6
+                     : input.moose.name.size()); // 6 is the length of "(old) "
+  for (const auto & param : _params)
+    max_moose_name_length = std::max(max_moose_name_length, param.moose.name.size());
+  for (const auto & output : _outputs)
+    max_moose_name_length = std::max(max_moose_name_length, output.moose.name.size());
+  for (const auto & deriv : _derivs)
+    max_moose_name_length = std::max(max_moose_name_length, deriv.moose.name.size());
+  for (const auto & param_deriv : _param_derivs)
+    max_moose_name_length = std::max(max_moose_name_length, param_deriv.moose.name.size());
+  return max_moose_name_length;
+}
 #endif // NEML2_ENABLED
diff --git a/modules/solid_mechanics/src/neml2/actions/NEML2ActionCommon.C b/modules/solid_mechanics/src/neml2/actions/NEML2ActionCommon.C
index b919a130abc3..37e60985b7d9 100644
--- a/modules/solid_mechanics/src/neml2/actions/NEML2ActionCommon.C
+++ b/modules/solid_mechanics/src/neml2/actions/NEML2ActionCommon.C
@@ -14,6 +14,7 @@
 
 #ifdef NEML2_ENABLED
 #include "neml2/base/Factory.h"
+#include "neml2/misc/parser_utils.h"
 #endif
 
 registerMooseAction("SolidMechanicsApp", NEML2ActionCommon, "parse_neml2");
@@ -111,6 +112,21 @@ NEML2ActionCommon::commonParams()
                         true,
                         NEML2Utils::docstring("Whether to print additional information about the "
                                               "NEML2 model at the beginning of the simulation"));
+
+  params.addParam<std::vector<MaterialPropertyName>>(
+      "initialize_outputs",
+      {},
+      "List of MOOSE material properties to be initialized. Each these properties must correspond "
+      "to a stateful NEML2 variable (which appears on both the input old state sub-axis and the "
+      "output state sub-axis). These MOOSE material properties will be initialized with the values "
+      "of properties specified in the initialize_output_values list.");
+  params.addParam<std::vector<MaterialPropertyName>>(
+      "initialize_output_values",
+      {},
+      "List of MOOSE material properties whose initial values (evaluated at the 0th time step) "
+      "will be used to initialize stateful properties. See the description of initialize_outputs "
+      "for more details.");
+
   params.addParam<std::vector<MaterialPropertyName>>(
       "export_outputs",
       {},
@@ -122,6 +138,7 @@ NEML2ActionCommon::commonParams()
       {},
       "The export targets corresponding to each MOOSE material property specified in "
       "export_outputs.");
+
   return params;
 }
 
@@ -133,28 +150,26 @@ NEML2ActionCommon::validParams()
   params.addRequiredParam<DataFileName>(
       "input",
       NEML2Utils::docstring("Path to the NEML2 input file containing the NEML2 model(s)."));
+  params.addParam<std::vector<std::string>>(
+      "cli_args",
+      {},
+      "Additional command line arguments to use when parsing the NEML2 input file.");
   return params;
 }
 
 NEML2ActionCommon::NEML2ActionCommon(const InputParameters & params)
-  : Action(params), _fname(getDataFileName("input"))
+  : Action(params),
+    _fname(getParam<DataFileName>("input")),
+    _cli_args(getParam<std::vector<std::string>>("cli_args"))
 {
   NEML2Utils::assertNEML2Enabled();
 }
 
-#ifndef NEML2_ENABLED
-
-void
-NEML2ActionCommon::act()
-{
-}
-
-#else
-
 void
 NEML2ActionCommon::act()
 {
+#ifdef NEML2_ENABLED
   if (_current_task == "parse_neml2")
-    neml2::load_input(std::string(_fname));
-}
+    neml2::load_input(std::string(_fname), neml2::utils::join(_cli_args, " "));
 #endif
+}
diff --git a/modules/solid_mechanics/src/neml2/interfaces/MOOSEToNEML2.C b/modules/solid_mechanics/src/neml2/interfaces/MOOSEToNEML2.C
index 05050c9561ce..ec002cc33a3f 100644
--- a/modules/solid_mechanics/src/neml2/interfaces/MOOSEToNEML2.C
+++ b/modules/solid_mechanics/src/neml2/interfaces/MOOSEToNEML2.C
@@ -36,7 +36,7 @@ MOOSEToNEML2::setMode(MOOSEToNEML2::Mode m) const
   _mode = m;
 
   if (_mode == Mode::VARIABLE || _mode == Mode::OLD_VARIABLE)
-    _neml2_variable = neml2::utils::parse<neml2::VariableName>(_raw_name);
+    _neml2_variable = NEML2Utils::parseVariableName(_raw_name);
   else if (_mode == Mode::PARAMETER)
     _neml2_parameter = _raw_name;
   else
@@ -70,12 +70,13 @@ MOOSEToNEML2::NEML2ParameterName() const
 }
 
 void
-MOOSEToNEML2::insertInto(neml2::Model & model) const
+MOOSEToNEML2::insertInto(neml2::ValueMap & input,
+                         std::map<std::string, neml2::Tensor> & params) const
 {
   if (_mode == Mode::VARIABLE || _mode == Mode::OLD_VARIABLE)
-    model.input_storage().base_index_put_(_neml2_variable, gatheredData());
+    input[_neml2_variable] = gatheredData();
   else if (_mode == Mode::PARAMETER)
-    model.set_parameter(_neml2_parameter, gatheredData());
+    params[_neml2_parameter] = gatheredData();
   else
     mooseError("Encountered invalid Mode in MOOSEToNEML2::insertInto");
 }
diff --git a/modules/solid_mechanics/src/neml2/materials/NEML2ToMOOSEMaterialProperty.C b/modules/solid_mechanics/src/neml2/materials/NEML2ToMOOSEMaterialProperty.C
index a3725490eea2..ce169955e562 100644
--- a/modules/solid_mechanics/src/neml2/materials/NEML2ToMOOSEMaterialProperty.C
+++ b/modules/solid_mechanics/src/neml2/materials/NEML2ToMOOSEMaterialProperty.C
@@ -16,7 +16,9 @@
 RegisterNEML2ToMOOSEMaterialProperty(Real);
 RegisterNEML2ToMOOSEMaterialProperty(SymmetricRankTwoTensor);
 RegisterNEML2ToMOOSEMaterialProperty(SymmetricRankFourTensor);
-RegisterNEML2ToMOOSEMaterialProperty(StdVector);
+RegisterNEML2ToMOOSEMaterialProperty(RealVectorValue);
+RegisterNEML2ToMOOSEMaterialProperty(RankTwoTensor);
+RegisterNEML2ToMOOSEMaterialProperty(RankFourTensor);
 
 template <typename T>
 InputParameters
@@ -65,45 +67,42 @@ NEML2ToMOOSEMaterialProperty<T>::NEML2ToMOOSEMaterialProperty(const InputParamet
     _prop0(isParamValid("moose_material_property_init")
                ? &getMaterialProperty<T>("moose_material_property_init")
                : nullptr),
-    _output_view(
+    _value(
         !isParamValid("neml2_input_derivative")
             ? (!isParamValid("neml2_parameter_derivative")
-                   ? _execute_neml2_model.getOutputView(neml2::utils::parse<neml2::VariableName>(
-                         getParam<std::string>("from_neml2")))
-                   : _execute_neml2_model.getOutputParameterDerivativeView(
-                         neml2::utils::parse<neml2::VariableName>(
-                             getParam<std::string>("from_neml2")),
+                   ? _execute_neml2_model.getOutput(
+                         NEML2Utils::parseVariableName(getParam<std::string>("from_neml2")))
+                   : _execute_neml2_model.getOutputParameterDerivative(
+                         NEML2Utils::parseVariableName(getParam<std::string>("from_neml2")),
                          getParam<std::string>("neml2_parameter_derivative")))
-            : _execute_neml2_model.getOutputDerivativeView(
-                  neml2::utils::parse<neml2::VariableName>(getParam<std::string>("from_neml2")),
-                  neml2::utils::parse<neml2::VariableName>(
-                      getParam<std::string>("neml2_input_derivative"))))
+            : _execute_neml2_model.getOutputDerivative(
+                  NEML2Utils::parseVariableName(getParam<std::string>("from_neml2")),
+                  NEML2Utils::parseVariableName(getParam<std::string>("neml2_input_derivative"))))
 #endif
 {
   NEML2Utils::assertNEML2Enabled();
 }
 
 #ifdef NEML2_ENABLED
-template <typename T>
-void
-NEML2ToMOOSEMaterialProperty<T>::initQpStatefulProperties()
-{
-  if (_prop0)
-    _prop[_qp] = (*_prop0)[_qp];
-}
-
 template <typename T>
 void
 NEML2ToMOOSEMaterialProperty<T>::computeProperties()
 {
+  // See issue #28971: Using _prop0 to set initial condition for this possibly stateful property may
+  // not work. As a workaround, we set the initial condition here when _t_step == 0.
+  if (_t_step == 0 && _prop0)
+  {
+    _prop.set() = _prop0->get();
+    return;
+  }
+
   if (!_execute_neml2_model.outputReady())
     return;
 
   // look up start index for current element
   const auto i = _execute_neml2_model.getBatchIndex(_current_elem->id());
-
-  for (std::size_t qp = 0; qp < _qrule->n_points(); qp++)
-    _prop[qp] = NEML2Utils::toMOOSE<T>(_output_view.batch_index({neml2::Size(i + qp)}));
+  for (_qp = 0; _qp < _qrule->n_points(); ++_qp)
+    NEML2Utils::copyTensorToMOOSEData(_value.batch_index({neml2::Size(i + _qp)}), _prop[_qp]);
 }
 #endif
 
@@ -112,4 +111,6 @@ NEML2ToMOOSEMaterialProperty<T>::computeProperties()
 InstantiateNEML2ToMOOSEMaterialProperty(Real);
 InstantiateNEML2ToMOOSEMaterialProperty(SymmetricRankTwoTensor);
 InstantiateNEML2ToMOOSEMaterialProperty(SymmetricRankFourTensor);
-InstantiateNEML2ToMOOSEMaterialProperty(std::vector<Real>);
+InstantiateNEML2ToMOOSEMaterialProperty(RealVectorValue);
+InstantiateNEML2ToMOOSEMaterialProperty(RankTwoTensor);
+InstantiateNEML2ToMOOSEMaterialProperty(RankFourTensor);
diff --git a/modules/solid_mechanics/src/neml2/userobjects/MOOSEMaterialPropertyToNEML2.C b/modules/solid_mechanics/src/neml2/userobjects/MOOSEMaterialPropertyToNEML2.C
index 722fa8d29eab..37b92da7788a 100644
--- a/modules/solid_mechanics/src/neml2/userobjects/MOOSEMaterialPropertyToNEML2.C
+++ b/modules/solid_mechanics/src/neml2/userobjects/MOOSEMaterialPropertyToNEML2.C
@@ -16,40 +16,31 @@
 RegisterMOOSEMaterialPropertyToNEML2(Real);
 RegisterMOOSEMaterialPropertyToNEML2(RankTwoTensor);
 RegisterMOOSEMaterialPropertyToNEML2(SymmetricRankTwoTensor);
-RegisterMOOSEMaterialPropertyToNEML2(StdVector);
+RegisterMOOSEMaterialPropertyToNEML2(RealVectorValue);
 
 template <typename T, unsigned int state>
 InputParameters
 MOOSEMaterialPropertyToNEML2<T, state>::validParams()
 {
-  auto params = MOOSEToNEML2Batched::validParams();
+  auto params = MOOSEToNEML2Batched<T>::validParams();
   params.addClassDescription(NEML2Utils::docstring(
       "Gather a MOOSE material property of type " + demangle(typeid(T).name()) +
       " for insertion into the specified input or model parameter of a NEML2 model."));
-  params.addRequiredParam<MaterialPropertyName>(
+  params.template addRequiredParam<MaterialPropertyName>(
       "from_moose", NEML2Utils::docstring("MOOSE material property to read from"));
   return params;
 }
 
 template <typename T, unsigned int state>
 MOOSEMaterialPropertyToNEML2<T, state>::MOOSEMaterialPropertyToNEML2(const InputParameters & params)
-  : MOOSEToNEML2Batched(params)
+  : MOOSEToNEML2Batched<T>(params)
 #ifdef NEML2_ENABLED
     ,
-    _mat_prop(getGenericMaterialProperty<T, false>("from_moose", state))
+    _mat_prop(this->template getGenericMaterialProperty<T, false>("from_moose", state))
 #endif
 {
 }
 
-#ifdef NEML2_ENABLED
-template <typename T, unsigned int state>
-torch::Tensor
-MOOSEMaterialPropertyToNEML2<T, state>::convertQpMOOSEData() const
-{
-  return NEML2Utils::toNEML2<T>(_mat_prop[_qp]);
-}
-#endif
-
 #define InstantiateMOOSEMaterialPropertyToNEML2(T)                                                 \
   template class MOOSEMaterialPropertyToNEML2<T, 0>;                                               \
   template class MOOSEMaterialPropertyToNEML2<T, 1>
@@ -57,4 +48,4 @@ MOOSEMaterialPropertyToNEML2<T, state>::convertQpMOOSEData() const
 InstantiateMOOSEMaterialPropertyToNEML2(Real);
 InstantiateMOOSEMaterialPropertyToNEML2(RankTwoTensor);
 InstantiateMOOSEMaterialPropertyToNEML2(SymmetricRankTwoTensor);
-InstantiateMOOSEMaterialPropertyToNEML2(std::vector<Real>);
+InstantiateMOOSEMaterialPropertyToNEML2(RealVectorValue);
diff --git a/modules/solid_mechanics/src/neml2/userobjects/MOOSEPostprocessorToNEML2.C b/modules/solid_mechanics/src/neml2/userobjects/MOOSEPostprocessorToNEML2.C
index 0bd60de3f34c..7407de9cc62b 100644
--- a/modules/solid_mechanics/src/neml2/userobjects/MOOSEPostprocessorToNEML2.C
+++ b/modules/solid_mechanics/src/neml2/userobjects/MOOSEPostprocessorToNEML2.C
@@ -50,7 +50,7 @@ template <unsigned int state>
 neml2::Tensor
 MOOSEPostprocessorToNEML2Templ<state>::gatheredData() const
 {
-  return NEML2Utils::toNEML2<Real>(_moose_pp);
+  return neml2::Scalar::full(_moose_pp);
 }
 #endif
 
diff --git a/modules/solid_mechanics/src/neml2/userobjects/MOOSEToNEML2Batched.C b/modules/solid_mechanics/src/neml2/userobjects/MOOSEToNEML2Batched.C
deleted file mode 100644
index 48261a6ad13a..000000000000
--- a/modules/solid_mechanics/src/neml2/userobjects/MOOSEToNEML2Batched.C
+++ /dev/null
@@ -1,60 +0,0 @@
-//* 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 "MOOSEToNEML2Batched.h"
-
-InputParameters
-MOOSEToNEML2Batched::validParams()
-{
-  auto params = MOOSEToNEML2::validParams();
-  params += ElementUserObject::validParams();
-
-  // Since we use the NEML2 model to evaluate the residual AND the Jacobian at the same time, we
-  // want to execute this user object only at execute_on = LINEAR (i.e. during residual evaluation).
-  // The NONLINEAR exec flag below is for computing Jacobian during automatic scaling.
-  ExecFlagEnum execute_options = MooseUtils::getDefaultExecFlagEnum();
-  execute_options = {EXEC_INITIAL, EXEC_LINEAR, EXEC_NONLINEAR};
-  params.set<ExecFlagEnum>("execute_on") = execute_options;
-
-  return params;
-}
-
-MOOSEToNEML2Batched::MOOSEToNEML2Batched(const InputParameters & params)
-  : MOOSEToNEML2(params), ElementUserObject(params)
-{
-}
-
-#ifdef NEML2_ENABLED
-void
-MOOSEToNEML2Batched::initialize()
-{
-  _buffer.clear();
-}
-
-void
-MOOSEToNEML2Batched::execute()
-{
-  for (_qp = 0; _qp < _qrule->n_points(); _qp++)
-    _buffer.push_back(convertQpMOOSEData());
-}
-
-void
-MOOSEToNEML2Batched::threadJoin(const UserObject & uo)
-{
-  // append vectors
-  const auto & m2n = static_cast<const MOOSEToNEML2Batched &>(uo);
-  _buffer.insert(_buffer.end(), m2n._buffer.begin(), m2n._buffer.end());
-}
-
-neml2::Tensor
-MOOSEToNEML2Batched::gatheredData() const
-{
-  return neml2::Tensor(torch::stack(_buffer, 0), 1);
-}
-#endif
diff --git a/modules/solid_mechanics/src/neml2/userobjects/MOOSEVariableToNEML2.C b/modules/solid_mechanics/src/neml2/userobjects/MOOSEVariableToNEML2.C
index f88b7c92719c..bcb8f1621425 100644
--- a/modules/solid_mechanics/src/neml2/userobjects/MOOSEVariableToNEML2.C
+++ b/modules/solid_mechanics/src/neml2/userobjects/MOOSEVariableToNEML2.C
@@ -44,14 +44,5 @@ MOOSEVariableToNEML2Templ<1>::MOOSEVariableToNEML2Templ(const InputParameters &
 {
 }
 
-#ifdef NEML2_ENABLED
-template <unsigned int state>
-torch::Tensor
-MOOSEVariableToNEML2Templ<state>::convertQpMOOSEData() const
-{
-  return NEML2Utils::toNEML2<Real>(_moose_variable[_qp]);
-}
-#endif
-
 template class MOOSEVariableToNEML2Templ<0>;
 template class MOOSEVariableToNEML2Templ<1>;
diff --git a/modules/solid_mechanics/src/neml2/userobjects/NEML2ModelExecutor.C b/modules/solid_mechanics/src/neml2/userobjects/NEML2ModelExecutor.C
index b432fa0149fa..3e66ebcc05e4 100644
--- a/modules/solid_mechanics/src/neml2/userobjects/NEML2ModelExecutor.C
+++ b/modules/solid_mechanics/src/neml2/userobjects/NEML2ModelExecutor.C
@@ -25,9 +25,10 @@ NEML2ModelExecutor::actionParams()
   params.addParam<std::vector<std::string>>(
       "skip_inputs",
       {},
-      NEML2Utils::docstring("List of NEML2 variables to skip when setting up the model input. If "
-                            "an input variable is skipped, its value will stay zero. If a required "
-                            "input variable is not skipped, an error will be raised."));
+      NEML2Utils::docstring(
+          "List of NEML2 variables to skip error checking when setting up the model input. If an "
+          "input variable is skipped, its value will stay zero. If a required input variable is "
+          "not skipped, an error will be raised."));
   return params;
 }
 
@@ -66,23 +67,21 @@ NEML2ModelExecutor::NEML2ModelExecutor(const InputParameters & params)
 #ifdef NEML2_ENABLED
     ,
     _batch_index_generator(getUserObject<NEML2BatchIndexGenerator>("batch_index_generator")),
-    _output_ready(false),
-    _in(model().input_storage()),
-    _out(model().output_storage()),
-    _dout_din(model().derivative_storage())
+    _output_ready(false)
 #endif
 {
 #ifdef NEML2_ENABLED
+  validateModel();
+
   // add user object dependencies by name (the UOs do not need to exist yet for this)
   for (const auto & gatherer_name : getParam<std::vector<UserObjectName>>("gatherers"))
     _depend_uo.insert(gatherer_name);
   for (const auto & gatherer_name : getParam<std::vector<UserObjectName>>("param_gatherers"))
     _depend_uo.insert(gatherer_name);
 
+  // variables to skip error checking (converting vector to set to prevent duplicate checks)
   for (const auto & var_name : getParam<std::vector<std::string>>("skip_inputs"))
-    _skip_vars.insert(neml2::utils::parse<neml2::VariableName>(var_name));
-
-  validateModel();
+    _skip_vars.insert(NEML2Utils::parseVariableName(var_name));
 #endif
 }
 
@@ -95,26 +94,24 @@ NEML2ModelExecutor::initialSetup()
   {
     // gather coupled user objects late to ensure they are constructed. Do not add them as
     // dependencies (that's already done in the constructor).
-    const auto & uo = getUserObjectByName<MOOSEToNEML2>(gatherer_name, /*is_dependency = */ false);
+    const auto & uo = getUserObjectByName<MOOSEToNEML2>(gatherer_name, /*is_dependency=*/false);
 
-    // introspect the NEML2 model to figure out if the gatherer UO is gathering for a NEML2 input
-    // variable or for a NEML2 model parameter
-    if (model().input_axis().has_variable(getNEML2VariableName(uo.NEML2Name())))
-    {
-      auto varname = getNEML2VariableName(uo.NEML2Name());
-      if (varname.start_with("old_forces") || varname.start_with("old_state"))
-        uo.setMode(MOOSEToNEML2::Mode::OLD_VARIABLE);
-      else
-        uo.setMode(MOOSEToNEML2::Mode::VARIABLE);
-      addGatheredVariable(gatherer_name, uo.NEML2VariableName());
-    }
-    else
+    // the target neml2 variable must exist on the input axis
+    if (!model().input_axis().has_variable(NEML2Utils::parseVariableName(uo.NEML2Name())))
       mooseError("The MOOSEToNEML2 gatherer named '",
                  gatherer_name,
                  "' is gathering MOOSE data for a non-existent NEML2 input variable named '",
                  uo.NEML2Name(),
                  "'.");
 
+    // tell the gatherer to gather for a model input variable
+    const auto varname = NEML2Utils::parseVariableName(uo.NEML2Name());
+    if (varname.is_old_force() || varname.is_old_state())
+      uo.setMode(MOOSEToNEML2::Mode::OLD_VARIABLE);
+    else
+      uo.setMode(MOOSEToNEML2::Mode::VARIABLE);
+
+    addGatheredVariable(gatherer_name, uo.NEML2VariableName());
     _gatherers.push_back(&uo);
   }
 
@@ -123,32 +120,32 @@ NEML2ModelExecutor::initialSetup()
   {
     // gather coupled user objects late to ensure they are constructed. Do not add them as
     // dependencies (that's already done in the constructor).
-    const auto & uo = getUserObjectByName<MOOSEToNEML2>(gatherer_name, /*is_dependency = */ false);
+    const auto & uo = getUserObjectByName<MOOSEToNEML2>(gatherer_name, /*is_dependency=*/false);
 
     // introspect the NEML2 model to figure out if the gatherer UO is gathering for a NEML2 input
     // variable or for a NEML2 model parameter
-    if (model().named_parameters().has_key(uo.NEML2Name()))
-    {
-      uo.setMode(MOOSEToNEML2::Mode::PARAMETER);
-      addGatheredParameter(gatherer_name, uo.NEML2ParameterName());
-    }
-    else
+    if (!model().named_parameters().has_key(uo.NEML2Name()))
       mooseError("The MOOSEToNEML2 gatherer named '",
                  gatherer_name,
                  "' is gathering MOOSE data for a non-existent NEML2 model parameter named '",
                  uo.NEML2Name(),
                  "'.");
 
+    // tell the gatherer to gather for a model parameter
+    uo.setMode(MOOSEToNEML2::Mode::PARAMETER);
+
+    addGatheredParameter(gatherer_name, uo.NEML2ParameterName());
     _gatherers.push_back(&uo);
   }
 
   // iterate over set of required inputs and error out if we find one that is not provided
-  std::set<neml2::VariableName> required_inputs = model().input_axis().variable_names();
+  std::vector<neml2::VariableName> required_inputs = model().input_axis().variable_names();
   for (const auto & input : required_inputs)
   {
     if (_skip_vars.count(input))
       continue;
-    if (input.start_with("state"))
+    // skip input state variables because they are "initial guesses" to the nonlinear system
+    if (input.is_state())
       continue;
     if (!_gathered_variable_names.count(input))
       paramError("gatherers", "The required model input `", input, "` is not gathered");
@@ -158,7 +155,7 @@ NEML2ModelExecutor::initialSetup()
   // sufficient for stateful data management, but that's the best we can do here without being
   // overly restrictive.
   for (const auto & input : required_inputs)
-    if (input.start_with("old_state") && !_retrieved_outputs.count(input.remount("state")))
+    if (input.is_old_state() && !_retrieved_outputs.count(input.current()))
       mooseError(
           "The NEML2 model requires a stateful input variable `",
           input,
@@ -219,12 +216,6 @@ NEML2ModelExecutor::execute()
 
   try
   {
-    auto batch_shape = neml2::TensorShape{neml2::Size(_batch_index_generator.getBatchIndex())};
-
-    // Reallocate the variable storage only when the batch shape has changed
-    if (neml2::TensorShapeRef(batch_shape) != model().batch_sizes())
-      initModel(batch_shape);
-
     fillInputs();
 
     if (_t_step > 0)
@@ -255,50 +246,75 @@ void
 NEML2ModelExecutor::fillInputs()
 {
   for (const auto & uo : _gatherers)
-    uo->insertInto(model());
+    uo->insertInto(_in, _model_params);
+
+  // Update model parameters
+  model().set_parameters(_model_params);
+  _model_params.clear();
 
   // Request gradient for the model parameters that we request AD for
-  for (auto & [deriv, tensor] : _retrieved_parameter_derivatives)
-    model().get_parameter(deriv.second).requires_grad_(true);
+  for (const auto & [y, dy] : _retrieved_parameter_derivatives)
+    for (const auto & [p, tensor] : dy)
+      model().get_parameter(p).requires_grad_(true);
 }
 
 void
 NEML2ModelExecutor::applyPredictor()
 {
+  if (!model().input_axis().has_state())
+    return;
+
   // Set trial state variables (i.e., initial guesses).
   // Right now we hard-code to use the old state as the trial state.
   // TODO: implement other predictors
-  if (model().input_axis().has_subaxis("state") && model().input_axis().has_subaxis("old_state"))
-    _in.slice("state").fill(_in.slice("old_state"));
+  const auto input_state = model().input_axis().subaxis(neml2::STATE);
+  const auto input_old_state = model().input_axis().subaxis(neml2::OLD_STATE);
+  for (const auto & var : input_state.variable_names())
+    if (input_old_state.has_variable(var))
+      _in[var.prepend(neml2::STATE)] = _in[var.prepend(neml2::OLD_STATE)];
 }
 
 void
 NEML2ModelExecutor::solve()
 {
   // Evaluate the NEML2 material model
-  model().value_and_dvalue(_in);
+  std::tie(_out, _dout_din) = model().value_and_dvalue(_in);
+  _in.clear();
 }
 
 void
 NEML2ModelExecutor::extractOutputs()
 {
+  const auto N = _batch_index_generator.getBatchIndex();
+
   // retrieve outputs
-  for (auto & [y, tensor] : _retrieved_outputs)
-    tensor = _out.base_index(y).to(torch::kCPU);
+  for (auto & [y, target] : _retrieved_outputs)
+    target = _out[y].to(torch::kCPU);
+
+  // retrieve parameter derivatives
+  for (auto & [y, dy] : _retrieved_parameter_derivatives)
+    for (auto & [p, target] : dy)
+      target = neml2::math::jacrev(_out[y],
+                                   model().get_parameter(p),
+                                   /*retain_graph=*/true,
+                                   /*create_graph=*/false,
+                                   /*allow_unused=*/false)
+                   .to(torch::kCPU);
+
+  // clear output
+  _out.clear();
 
   // retrieve derivatives
-  for (auto & [deriv, tensor] : _retrieved_derivatives)
-  {
-    auto & [y, x] = deriv;
-    tensor = _dout_din.base_index({y, x}).to(torch::kCPU);
-  }
+  for (auto & [y, dy] : _retrieved_derivatives)
+    for (auto & [x, target] : dy)
+    {
+      const auto & source = _dout_din[y][x];
+      if (source.defined())
+        target = source.to(torch::kCPU).batch_expand({neml2::Size(N)});
+    }
 
-  // retrieve parameter derivatives
-  for (auto & [deriv, tensor] : _retrieved_parameter_derivatives)
-  {
-    auto & [y, p] = deriv;
-    tensor = neml2::math::jacrev(_out.base_index(y), model().get_parameter(p)).to(torch::kCPU);
-  }
+  // clear derivatives
+  _dout_din.clear();
 }
 
 void
@@ -310,7 +326,7 @@ NEML2ModelExecutor::checkExecutionStage() const
 }
 
 const neml2::Tensor &
-NEML2ModelExecutor::getOutputView(const neml2::VariableName & output_name) const
+NEML2ModelExecutor::getOutput(const neml2::VariableName & output_name) const
 {
   checkExecutionStage();
 
@@ -321,8 +337,8 @@ NEML2ModelExecutor::getOutputView(const neml2::VariableName & output_name) const
 }
 
 const neml2::Tensor &
-NEML2ModelExecutor::getOutputDerivativeView(const neml2::VariableName & output_name,
-                                            const neml2::VariableName & input_name) const
+NEML2ModelExecutor::getOutputDerivative(const neml2::VariableName & output_name,
+                                        const neml2::VariableName & input_name) const
 {
   checkExecutionStage();
 
@@ -340,12 +356,12 @@ NEML2ModelExecutor::getOutputDerivativeView(const neml2::VariableName & output_n
                input_name,
                "', but the NEML2 input variable does not exist.");
 
-  return _retrieved_derivatives[{output_name, input_name}];
+  return _retrieved_derivatives[output_name][input_name];
 }
 
 const neml2::Tensor &
-NEML2ModelExecutor::getOutputParameterDerivativeView(const neml2::VariableName & output_name,
-                                                     const std::string & parameter_name) const
+NEML2ModelExecutor::getOutputParameterDerivative(const neml2::VariableName & output_name,
+                                                 const std::string & parameter_name) const
 {
   checkExecutionStage();
 
@@ -363,7 +379,7 @@ NEML2ModelExecutor::getOutputParameterDerivativeView(const neml2::VariableName &
                parameter_name,
                "', but the NEML2 model parameter does not exist.");
 
-  return _retrieved_parameter_derivatives[{output_name, parameter_name}];
+  return _retrieved_parameter_derivatives[output_name][parameter_name];
 }
 
 #endif
diff --git a/modules/solid_mechanics/src/neml2/utils/NEML2Utils.C b/modules/solid_mechanics/src/neml2/utils/NEML2Utils.C
index bbccc4a6b60a..63fd01fa75ed 100644
--- a/modules/solid_mechanics/src/neml2/utils/NEML2Utils.C
+++ b/modules/solid_mechanics/src/neml2/utils/NEML2Utils.C
@@ -10,43 +10,8 @@
 #include "NEML2Utils.h"
 #include "SubProblem.h"
 
-#ifdef NEML2_ENABLED
-
-#include "VariadicTable.h"
-
-namespace neml2
-{
-
-std::ostream &
-operator<<(std::ostream & os, const Model & model)
-{
-  VariadicTable<std::string, std::string, std::string, Size> table(
-      {"Type", "Name", "Tensor type", "Storage size"});
-
-  for (const auto && [name, var] : model.input_variables())
-    table.addRow("input", utils::stringify(name), utils::stringify(var.type()), var.base_storage());
-
-  for (const auto && [name, var] : model.output_variables())
-    table.addRow(
-        "output", utils::stringify(name), utils::stringify(var.type()), var.base_storage());
-
-  for (auto && [name, param] : model.named_parameters())
-    table.addRow("parameter", name, utils::stringify(param.type()), Tensor(param).base_storage());
-
-  for (auto && [name, buffer] : model.named_buffers())
-    table.addRow("buffer", name, utils::stringify(buffer.type()), Tensor(buffer).base_storage());
-
-  table.print(os);
-
-  return os;
-}
-} // namespace neml2
-
-#endif // NEML2_ENABLED
-
 namespace NEML2Utils
 {
-
 #ifdef NEML2_ENABLED
 void
 assertVariable(const neml2::VariableName & v)
@@ -54,12 +19,8 @@ assertVariable(const neml2::VariableName & v)
   if (v.empty())
     mooseError("Empty NEML2 variable");
 
-  if (!v.start_with("forces") && !v.start_with("state"))
-    mooseError("The NEML2 variable '",
-               v,
-               "' should be defined on the forces or the state sub-axis, got ",
-               v.slice(0, 1),
-               " instead");
+  if (!v.is_force() && !v.is_state())
+    mooseError("The NEML2 variable '", v, "' is on the wrong subaxis.");
 }
 
 void
@@ -68,86 +29,15 @@ assertOldVariable(const neml2::VariableName & v)
   if (v.empty())
     mooseError("Empty NEML2 variable");
 
-  if (!v.start_with("old_forces") && !v.start_with("old_state"))
-    mooseError("The NEML2 variable '",
-               v,
-               "' should be defined on the old_forces or the old_state sub-axis, got ",
-               v.slice(0, 1),
-               " instead");
-}
-
-template <>
-neml2::Tensor
-toNEML2(const Real & v)
-{
-  return neml2::Scalar::full(v);
-}
-
-// FIXME: This is an unfortunately specialization because the models I included for testing use
-// symmetric tensors everywhere. Once I tested all the models with full tensors (i.e. not in Mandel
-// notation), I should be able to "fix" this specialization.
-template <>
-neml2::Tensor
-toNEML2(const RankTwoTensor & r2t)
-{
-  return neml2::SR2::fill(r2t(0, 0), r2t(1, 1), r2t(2, 2), r2t(1, 2), r2t(0, 2), r2t(0, 1));
-}
-
-template <>
-neml2::Tensor
-toNEML2(const SymmetricRankTwoTensor & r2t)
-{
-  return neml2::SR2::fill(r2t(0, 0), r2t(1, 1), r2t(2, 2), r2t(1, 2), r2t(0, 2), r2t(0, 1));
-}
-
-template <>
-neml2::Tensor
-toNEML2(const std::vector<Real> & v)
-{
-  return neml2::Tensor(torch::tensor(v, neml2::default_tensor_options()), 0);
-}
-
-template <>
-Real
-toMOOSE(const neml2::Tensor & t)
-{
-  return t.item<Real>();
-}
-
-template <>
-SymmetricRankTwoTensor
-toMOOSE(const neml2::Tensor & t)
-{
-  using symr2t = SymmetricRankTwoTensor;
-  return symr2t(t.base_index({0}).item<neml2::Real>() / symr2t::mandelFactor(0),
-                t.base_index({1}).item<neml2::Real>() / symr2t::mandelFactor(1),
-                t.base_index({2}).item<neml2::Real>() / symr2t::mandelFactor(2),
-                t.base_index({3}).item<neml2::Real>() / symr2t::mandelFactor(3),
-                t.base_index({4}).item<neml2::Real>() / symr2t::mandelFactor(4),
-                t.base_index({5}).item<neml2::Real>() / symr2t::mandelFactor(5));
-}
-
-template <>
-std::vector<Real>
-toMOOSE(const neml2::Tensor & t)
-{
-  auto tc = t.contiguous();
-  return std::vector<Real>(tc.data_ptr<neml2::Real>(), tc.data_ptr<neml2::Real>() + tc.numel());
+  if (!v.is_old_force() && !v.is_old_state())
+    mooseError("The NEML2 variable '", v, "' is on the wrong subaxis.");
 }
 
-template <>
-SymmetricRankFourTensor
-toMOOSE(const neml2::Tensor & t)
+neml2::VariableName
+parseVariableName(const std::string & s)
 {
-  // Well I don't see a good constructor for this, so let me fill out all the components.
-  SymmetricRankFourTensor symsymr4t;
-  for (const auto a : make_range(6))
-    for (const auto b : make_range(6))
-      symsymr4t(a, b) = t.base_index({a, b}).item<neml2::Real>();
-
-  return symsymr4t;
+  return neml2::utils::parse<neml2::VariableName>(s);
 }
-
 #endif // NEML2_ENABLED
 
 static const std::string missing_neml2 = "The `NEML2` library is required but not enabled. Refer "
diff --git a/modules/solid_mechanics/test/neml2/include/NEML2TestModel.h b/modules/solid_mechanics/test/neml2/include/NEML2TestModel.h
index 7de9326b8ed3..561f1a6d6a43 100644
--- a/modules/solid_mechanics/test/neml2/include/NEML2TestModel.h
+++ b/modules/solid_mechanics/test/neml2/include/NEML2TestModel.h
@@ -21,6 +21,8 @@ class NEML2TestModel : public Model
 
   NEML2TestModel(const OptionSet & options);
 
+  void request_AD() override;
+
 protected:
   // input variables
   const Variable<Scalar> & _input_a;
diff --git a/modules/solid_mechanics/test/neml2/src/NEML2TestModel.cxx b/modules/solid_mechanics/test/neml2/src/NEML2TestModel.cxx
index 85e32088c883..54892ec53275 100644
--- a/modules/solid_mechanics/test/neml2/src/NEML2TestModel.cxx
+++ b/modules/solid_mechanics/test/neml2/src/NEML2TestModel.cxx
@@ -22,9 +22,6 @@ NEML2TestModel::expected_options()
   options.set<VariableName>("B") = VariableName("forces", "B");
   options.set<VariableName>("sum") = VariableName("state", "internal", "sum");
   options.set<VariableName>("product") = VariableName("state", "internal", "product");
-  // Use AD
-  options.set<bool>("_use_AD_first_derivative") = true;
-  options.set<bool>("_use_AD_second_derivative") = true;
   return options;
 }
 
@@ -39,6 +36,16 @@ NEML2TestModel::NEML2TestModel(const OptionSet & options)
 {
 }
 
+void
+NEML2TestModel::request_AD()
+{
+  std::vector<const VariableBase *> inputs = {&_input_a, &_input_b};
+
+  // First derivatives
+  _sum.request_AD(inputs);
+  _product.request_AD(inputs);
+}
+
 void
 NEML2TestModel::set_value(bool out, bool dout_din, bool d2out_din2)
 {
diff --git a/modules/solid_mechanics/test/tests/neml2/crystal_plasticity.i b/modules/solid_mechanics/test/tests/neml2/crystal_plasticity.i
new file mode 100644
index 000000000000..0bd3f7da5ef3
--- /dev/null
+++ b/modules/solid_mechanics/test/tests/neml2/crystal_plasticity.i
@@ -0,0 +1,118 @@
+N = 2
+
+[GlobalParams]
+  displacements = 'disp_x disp_y disp_z'
+[]
+
+[Mesh]
+  [gmg]
+    type = GeneratedMeshGenerator
+    dim = 3
+    nx = ${N}
+    ny = ${N}
+    nz = ${N}
+  []
+[]
+
+[Physics]
+  [SolidMechanics]
+    [QuasiStatic]
+      [all]
+        strain = FINITE
+        new_system = true
+        add_variables = true
+        formulation = TOTAL
+        volumetric_locking_correction = true
+      []
+    []
+  []
+[]
+
+[NEML2]
+  input = 'models/crystal_plasticity.i'
+  [all]
+    model = 'model'
+    verbose = true
+    device = 'cpu'
+
+    moose_input_types = 'MATERIAL                           POSTPROCESSOR POSTPROCESSOR MATERIAL                  MATERIAL                  MATERIAL'
+    moose_inputs = '     spatial_velocity_increment         time          time          elastic_strain            orientation               slip_hardening'
+    neml2_inputs = '     forces/spatial_velocity_increment  forces/t      old_forces/t  old_state/elastic_strain  old_state/orientation     old_state/internal/slip_hardening'
+
+    moose_output_types = 'MATERIAL                            MATERIAL                  MATERIAL                  MATERIAL'
+    moose_outputs = '     neml2_cauchy_stress                 elastic_strain            orientation               slip_hardening'
+    neml2_outputs = '     state/internal/full_cauchy_stress   state/elastic_strain      state/orientation         state/internal/slip_hardening'
+
+    moose_derivative_types = 'MATERIAL'
+    moose_derivatives = '     neml2_cauchy_jacobian'
+    neml2_derivatives = '     state/internal/full_cauchy_stress forces/spatial_velocity_increment'
+
+    initialize_outputs = '      orientation'
+    initialize_output_values = 'initial_orientation'
+  []
+[]
+
+[Postprocessors]
+  [time]
+    type = TimePostprocessor
+    execute_on = 'INITIAL TIMESTEP_BEGIN'
+  []
+[]
+
+[Materials]
+  [copy]
+    type = ComputeLagrangianCauchyCustomStress
+    custom_cauchy_stress = 'neml2_cauchy_stress'
+    custom_cauchy_jacobian = 'neml2_cauchy_jacobian'
+  []
+  [initial_orientation]
+    type = GenericConstantRealVectorValue
+    vector_name = 'initial_orientation'
+    vector_values = '-0.54412095 -0.34931944 0.12600655'
+  []
+[]
+
+[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
+  []
+  [xdisp]
+    type = FunctionDirichletBC
+    variable = disp_x
+    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 = 5e-3
+  dtmin = 1e-3
+  num_steps = 5
+  residual_and_jacobian_together = true
+[]
+
+[Outputs]
+  file_base = 'crystal_plasticity'
+  exodus = true
+[]
diff --git a/modules/solid_mechanics/test/tests/neml2/custom_model.i b/modules/solid_mechanics/test/tests/neml2/custom_model.i
index 312311f23bb9..a7a85e4b2496 100644
--- a/modules/solid_mechanics/test/tests/neml2/custom_model.i
+++ b/modules/solid_mechanics/test/tests/neml2/custom_model.i
@@ -13,7 +13,6 @@
     model = 'model'
     verbose = true
     device = 'cpu'
-    enable_AD = true
 
     moose_input_types = 'VARIABLE MATERIAL'
     moose_inputs = '     a        b'
diff --git a/modules/solid_mechanics/test/tests/neml2/gold/crystal_plasticity.e b/modules/solid_mechanics/test/tests/neml2/gold/crystal_plasticity.e
new file mode 100644
index 000000000000..88f76a60081b
Binary files /dev/null and b/modules/solid_mechanics/test/tests/neml2/gold/crystal_plasticity.e differ
diff --git a/modules/solid_mechanics/test/tests/neml2/models/crystal_plasticity.i b/modules/solid_mechanics/test/tests/neml2/models/crystal_plasticity.i
new file mode 100644
index 000000000000..2459ca64e6b6
--- /dev/null
+++ b/modules/solid_mechanics/test/tests/neml2/models/crystal_plasticity.i
@@ -0,0 +1,128 @@
+[Tensors]
+  [a]
+    type = Scalar
+    values = '1.0'
+  []
+  [sdirs]
+    type = FillMillerIndex
+    values = '1 1 0'
+  []
+  [splanes]
+    type = FillMillerIndex
+    values = '1 1 1'
+  []
+[]
+
+[Solvers]
+  [newton]
+    type = NewtonWithLineSearch
+    max_linesearch_iterations = 5
+  []
+[]
+
+[Data]
+  [crystal_geometry]
+    type = CubicCrystal
+    lattice_parameter = "a"
+    slip_directions = "sdirs"
+    slip_planes = "splanes"
+  []
+[]
+
+[Models]
+  [spatial_velocity_gradient]
+    type = R2IncrementalRate
+    variable = 'forces/spatial_velocity_increment'
+    time = 'forces/t'
+    rate = 'forces/spatial_velocity_gradient'
+  []
+  [split_to_deformation_rate]
+    type = R2toSR2
+    input = 'forces/spatial_velocity_gradient'
+    output = 'forces/deformation_rate'
+  []
+  [split_to_vorticity]
+    type = R2toWR2
+    input = 'forces/spatial_velocity_gradient'
+    output = 'forces/vorticity'
+  []
+  [euler_rodrigues]
+    type = RotationMatrix
+    from = 'state/orientation'
+    to = 'state/orientation_matrix'
+  []
+  [elasticity]
+    type = LinearIsotropicElasticity
+    coefficient_types = 'YOUNGS_MODULUS POISSONS_RATIO'
+    coefficients = '1e5 0.25'
+    strain = "state/elastic_strain"
+    stress = "state/internal/cauchy_stress"
+  []
+  [resolved_shear]
+    type = ResolvedShear
+  []
+  [elastic_stretch]
+    type = ElasticStrainRate
+  []
+  [plastic_spin]
+    type = PlasticVorticity
+  []
+  [plastic_deformation_rate]
+    type = PlasticDeformationRate
+  []
+  [orientation_rate]
+    type = OrientationRate
+  []
+  [sum_slip_rates]
+    type = SumSlipRates
+  []
+  [slip_rule]
+    type = PowerLawSlipRule
+    n = 8.0
+    gamma0 = 2.0e-1
+  []
+  [slip_strength]
+    type = SingleSlipStrengthMap
+    constant_strength = 50.0
+  []
+  [voce_hardening]
+    type = VoceSingleSlipHardeningRule
+    initial_slope = 500.0
+    saturated_hardening = 50.0
+  []
+  [integrate_slip_hardening]
+    type = ScalarBackwardEulerTimeIntegration
+    variable = 'state/internal/slip_hardening'
+  []
+  [integrate_elastic_strain]
+    type = SR2BackwardEulerTimeIntegration
+    variable = 'state/elastic_strain'
+  []
+  [integrate_orientation]
+    type = WR2ImplicitExponentialTimeIntegration
+    variable = 'state/orientation'
+  []
+
+  [implicit_rate]
+    type = ComposedModel
+    models = "spatial_velocity_gradient split_to_deformation_rate split_to_vorticity euler_rodrigues elasticity orientation_rate resolved_shear
+              elastic_stretch plastic_deformation_rate plastic_spin
+              sum_slip_rates slip_rule slip_strength voce_hardening
+              integrate_slip_hardening integrate_elastic_strain integrate_orientation"
+  []
+  [model_without_stress]
+    type = ImplicitUpdate
+    implicit_model = 'implicit_rate'
+    solver = 'newton'
+  []
+  [full_stress]
+    type = SR2toR2
+    input = 'state/internal/cauchy_stress'
+    output = 'state/internal/full_cauchy_stress'
+  []
+  [model]
+    type = ComposedModel
+    models = 'model_without_stress elasticity full_stress'
+    additional_outputs = 'state/elastic_strain'
+  []
+[]
diff --git a/modules/solid_mechanics/test/tests/neml2/models/elasticity.i b/modules/solid_mechanics/test/tests/neml2/models/elasticity.i
index c2a54c6ed04a..f5e9f15b7013 100644
--- a/modules/solid_mechanics/test/tests/neml2/models/elasticity.i
+++ b/modules/solid_mechanics/test/tests/neml2/models/elasticity.i
@@ -1,8 +1,8 @@
 [Models]
   [model]
     type = LinearIsotropicElasticity
-    youngs_modulus = 1
-    poisson_ratio = 0.3
+    coefficients = '1 0.3'
+    coefficient_types = 'YOUNGS_MODULUS POISSONS_RATIO'
     strain = 'forces/E'
   []
 []
diff --git a/modules/solid_mechanics/test/tests/neml2/models/elasticity_2.i b/modules/solid_mechanics/test/tests/neml2/models/elasticity_2.i
index fc11758a371b..a06662ac9108 100644
--- a/modules/solid_mechanics/test/tests/neml2/models/elasticity_2.i
+++ b/modules/solid_mechanics/test/tests/neml2/models/elasticity_2.i
@@ -1,14 +1,14 @@
 [Models]
   [model1]
     type = LinearIsotropicElasticity
-    youngs_modulus = 1
-    poisson_ratio = 0.3
+    coefficients = '1 0.3'
+    coefficient_types = 'YOUNGS_MODULUS POISSONS_RATIO'
     strain = 'forces/E'
   []
   [model2]
     type = LinearIsotropicElasticity
-    youngs_modulus = 2
-    poisson_ratio = 0.3
+    coefficients = '2 0.3'
+    coefficient_types = 'YOUNGS_MODULUS POISSONS_RATIO'
     strain = 'forces/E'
   []
 []
diff --git a/modules/solid_mechanics/test/tests/neml2/models/error.i b/modules/solid_mechanics/test/tests/neml2/models/error.i
index d1a3b30f24b7..a35727029365 100644
--- a/modules/solid_mechanics/test/tests/neml2/models/error.i
+++ b/modules/solid_mechanics/test/tests/neml2/models/error.i
@@ -50,8 +50,8 @@
   []
   [elasticity]
     type = LinearIsotropicElasticity
-    youngs_modulus = 1e5
-    poisson_ratio = 0.3
+    coefficients = '1e5 0.3'
+    coefficient_types = 'YOUNGS_MODULUS POISSONS_RATIO'
     rate_form = true
   []
   [integrate_stress]
diff --git a/modules/solid_mechanics/test/tests/neml2/models/rate_independent_plasticity_isoharden.i b/modules/solid_mechanics/test/tests/neml2/models/rate_independent_plasticity_isoharden.i
index ab20052a632f..4893a9e0e42b 100644
--- a/modules/solid_mechanics/test/tests/neml2/models/rate_independent_plasticity_isoharden.i
+++ b/modules/solid_mechanics/test/tests/neml2/models/rate_independent_plasticity_isoharden.i
@@ -18,8 +18,8 @@
   []
   [elasticity]
     type = LinearIsotropicElasticity
-    youngs_modulus = 1e5
-    poisson_ratio = 0.3
+    coefficients = '1e5 0.3'
+    coefficient_types = 'YOUNGS_MODULUS POISSONS_RATIO'
   []
   [mandel_stress]
     type = IsotropicMandelStress
diff --git a/modules/solid_mechanics/test/tests/neml2/models/rate_independent_plasticity_isokinharden.i b/modules/solid_mechanics/test/tests/neml2/models/rate_independent_plasticity_isokinharden.i
index 82fb3d591ac9..18a611ea53bd 100644
--- a/modules/solid_mechanics/test/tests/neml2/models/rate_independent_plasticity_isokinharden.i
+++ b/modules/solid_mechanics/test/tests/neml2/models/rate_independent_plasticity_isokinharden.i
@@ -22,8 +22,8 @@
   []
   [elasticity]
     type = LinearIsotropicElasticity
-    youngs_modulus = 1e5
-    poisson_ratio = 0.3
+    coefficients = '1e5 0.3'
+    coefficient_types = 'YOUNGS_MODULUS POISSONS_RATIO'
   []
   [mandel_stress]
     type = IsotropicMandelStress
diff --git a/modules/solid_mechanics/test/tests/neml2/models/rate_independent_plasticity_kinharden.i b/modules/solid_mechanics/test/tests/neml2/models/rate_independent_plasticity_kinharden.i
index 7c6b149a198d..dc2adb553e5b 100644
--- a/modules/solid_mechanics/test/tests/neml2/models/rate_independent_plasticity_kinharden.i
+++ b/modules/solid_mechanics/test/tests/neml2/models/rate_independent_plasticity_kinharden.i
@@ -17,8 +17,8 @@
   []
   [elasticity]
     type = LinearIsotropicElasticity
-    youngs_modulus = 1e5
-    poisson_ratio = 0.3
+    coefficients = '1e5 0.3'
+    coefficient_types = 'YOUNGS_MODULUS POISSONS_RATIO'
   []
   [mandel_stress]
     type = IsotropicMandelStress
diff --git a/modules/solid_mechanics/test/tests/neml2/models/rate_independent_plasticity_perfect.i b/modules/solid_mechanics/test/tests/neml2/models/rate_independent_plasticity_perfect.i
index 852c3b626605..0da308dd49c6 100644
--- a/modules/solid_mechanics/test/tests/neml2/models/rate_independent_plasticity_perfect.i
+++ b/modules/solid_mechanics/test/tests/neml2/models/rate_independent_plasticity_perfect.i
@@ -13,8 +13,8 @@
   []
   [elasticity]
     type = LinearIsotropicElasticity
-    youngs_modulus = 1e5
-    poisson_ratio = 0.3
+    coefficients = '1e5 0.3'
+    coefficient_types = 'YOUNGS_MODULUS POISSONS_RATIO'
   []
   [mandel_stress]
     type = IsotropicMandelStress
diff --git a/modules/solid_mechanics/test/tests/neml2/models/viscoplasticity_chaboche.i b/modules/solid_mechanics/test/tests/neml2/models/viscoplasticity_chaboche.i
index b7e55f8f12fe..16ebe0a9b506 100644
--- a/modules/solid_mechanics/test/tests/neml2/models/viscoplasticity_chaboche.i
+++ b/modules/solid_mechanics/test/tests/neml2/models/viscoplasticity_chaboche.i
@@ -85,8 +85,8 @@
   []
   [elasticity]
     type = LinearIsotropicElasticity
-    youngs_modulus = 1e5
-    poisson_ratio = 0.3
+    coefficients = '1e5 0.3'
+    coefficient_types = 'YOUNGS_MODULUS POISSONS_RATIO'
     rate_form = true
   []
   [integrate_ep]
diff --git a/modules/solid_mechanics/test/tests/neml2/models/viscoplasticity_isoharden.i b/modules/solid_mechanics/test/tests/neml2/models/viscoplasticity_isoharden.i
index f56c3fa1de1e..2652dd7d89e9 100644
--- a/modules/solid_mechanics/test/tests/neml2/models/viscoplasticity_isoharden.i
+++ b/modules/solid_mechanics/test/tests/neml2/models/viscoplasticity_isoharden.i
@@ -57,8 +57,8 @@
   []
   [elasticity]
     type = LinearIsotropicElasticity
-    youngs_modulus = 1e5
-    poisson_ratio = 0.3
+    coefficients = '1e5 0.3'
+    coefficient_types = 'YOUNGS_MODULUS POISSONS_RATIO'
     rate_form = true
   []
   [integrate_stress]
diff --git a/modules/solid_mechanics/test/tests/neml2/models/viscoplasticity_isokinharden.i b/modules/solid_mechanics/test/tests/neml2/models/viscoplasticity_isokinharden.i
index 54305d68ff0f..d8e40c076e60 100644
--- a/modules/solid_mechanics/test/tests/neml2/models/viscoplasticity_isokinharden.i
+++ b/modules/solid_mechanics/test/tests/neml2/models/viscoplasticity_isokinharden.i
@@ -70,8 +70,8 @@
   []
   [elasticity]
     type = LinearIsotropicElasticity
-    youngs_modulus = 1e5
-    poisson_ratio = 0.3
+    coefficients = '1e5 0.3'
+    coefficient_types = 'YOUNGS_MODULUS POISSONS_RATIO'
     rate_form = true
   []
   [integrate_ep]
diff --git a/modules/solid_mechanics/test/tests/neml2/models/viscoplasticity_kinharden.i b/modules/solid_mechanics/test/tests/neml2/models/viscoplasticity_kinharden.i
index 77a928dd3a36..aca3cc984184 100644
--- a/modules/solid_mechanics/test/tests/neml2/models/viscoplasticity_kinharden.i
+++ b/modules/solid_mechanics/test/tests/neml2/models/viscoplasticity_kinharden.i
@@ -62,8 +62,8 @@
   []
   [elasticity]
     type = LinearIsotropicElasticity
-    youngs_modulus = 1e5
-    poisson_ratio = 0.3
+    coefficients = '1e5 0.3'
+    coefficient_types = 'YOUNGS_MODULUS POISSONS_RATIO'
     rate_form = true
   []
   [integrate_Kp]
diff --git a/modules/solid_mechanics/test/tests/neml2/models/viscoplasticity_perfect.i b/modules/solid_mechanics/test/tests/neml2/models/viscoplasticity_perfect.i
index b276b92bddeb..f88733dc3995 100644
--- a/modules/solid_mechanics/test/tests/neml2/models/viscoplasticity_perfect.i
+++ b/modules/solid_mechanics/test/tests/neml2/models/viscoplasticity_perfect.i
@@ -49,8 +49,8 @@
   []
   [elasticity]
     type = LinearIsotropicElasticity
-    youngs_modulus = 1e5
-    poisson_ratio = 0.3
+    coefficients = '1e5 0.3'
+    coefficient_types = 'YOUNGS_MODULUS POISSONS_RATIO'
     rate_form = true
   []
   [integrate_stress]
diff --git a/modules/solid_mechanics/test/tests/neml2/models/viscoplasticity_radial_return.i b/modules/solid_mechanics/test/tests/neml2/models/viscoplasticity_radial_return.i
index 2aef40ce234c..37ff9d878164 100644
--- a/modules/solid_mechanics/test/tests/neml2/models/viscoplasticity_radial_return.i
+++ b/modules/solid_mechanics/test/tests/neml2/models/viscoplasticity_radial_return.i
@@ -17,8 +17,8 @@
   []
   [trial_cauchy_stress]
     type = LinearIsotropicElasticity
-    youngs_modulus = 1e5
-    poisson_ratio = 0.3
+    coefficients = '1e5 0.3'
+    coefficient_types = 'YOUNGS_MODULUS POISSONS_RATIO'
     strain = 'forces/Ee_trial'
     stress = 'forces/S_trial'
   []
@@ -79,8 +79,8 @@
   []
   [cauchy_stress]
     type = LinearIsotropicElasticity
-    youngs_modulus = 1e5
-    poisson_ratio = 0.3
+    coefficients = '1e5 0.3'
+    coefficient_types = 'YOUNGS_MODULUS POISSONS_RATIO'
   []
   [mandel_stress]
     type = IsotropicMandelStress
diff --git a/modules/solid_mechanics/test/tests/neml2/tests b/modules/solid_mechanics/test/tests/neml2/tests
index 018fc91ee2fc..9046793f9bae 100644
--- a/modules/solid_mechanics/test/tests/neml2/tests
+++ b/modules/solid_mechanics/test/tests/neml2/tests
@@ -69,6 +69,16 @@
     issues = '#27493'
     collections = 'SYSTEM'
     requirement = 'The system shall be capable of modularly interfacing with the NEML2 library to perform constitutive update '
+    [crystal_plasticity]
+      detail = 'in a crystal plasticity problem.'
+      type = Exodiff
+      input = 'crystal_plasticity.i'
+      exodiff = 'crystal_plasticity.e'
+      required_submodule = 'contrib/neml2'
+      libtorch = true
+      recover = false
+      valgrind = heavy
+    []
     [elasticity]
       detail = 'to solve an elasticity problem;'
       type = Exodiff
diff --git a/modules/solid_mechanics/unit/src/NEML2UtilsTest.C b/modules/solid_mechanics/unit/src/NEML2UtilsTest.C
new file mode 100644
index 000000000000..49c4fef37dc6
--- /dev/null
+++ b/modules/solid_mechanics/unit/src/NEML2UtilsTest.C
@@ -0,0 +1,174 @@
+//* 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 "gtest/gtest.h"
+
+#include "NEML2Utils.h"
+
+#ifdef NEML2_ENABLED
+TEST(NEML2Utils, fromBlob_Real)
+{
+  std::vector<Real> data(3);
+  data[0] = 1.0;
+  data[1] = 2.0;
+  data[2] = 3.0;
+
+  auto tensor = NEML2Utils::fromBlob(data);
+  ASSERT_TRUE(tensor.defined());
+  ASSERT_TRUE(tensor.dim() == 1);
+  ASSERT_TRUE(tensor.batch_dim() == 1);
+  ASSERT_TRUE(tensor.base_dim() == 0);
+  ASSERT_TRUE(tensor.size(0) == 3);
+
+  for (neml2::Size n : index_range(data))
+    EXPECT_NEAR(tensor.index({n}).item<Real>(), data[n], 1e-12);
+}
+
+TEST(NEML2Utils, fromBlob_RealVectorValue)
+{
+  std::vector<RealVectorValue> data(3);
+  data[0] = RealVectorValue(1.0, 2.0, 3.0);
+  data[1] = RealVectorValue(4.0, 5.0, 6.0);
+  data[2] = RealVectorValue(7.0, 8.0, 9.0);
+
+  auto tensor = NEML2Utils::fromBlob(data);
+  ASSERT_TRUE(tensor.defined());
+  ASSERT_TRUE(tensor.dim() == 2);
+  ASSERT_TRUE(tensor.batch_dim() == 1);
+  ASSERT_TRUE(tensor.base_dim() == 1);
+  ASSERT_TRUE(tensor.size(0) == 3);
+  ASSERT_TRUE(tensor.size(1) == 3);
+
+  for (neml2::Size n : index_range(data))
+    for (neml2::Size i : make_range(3))
+      EXPECT_NEAR(tensor.index({n, i}).item<Real>(), data[n](i), 1e-12);
+}
+
+TEST(NEML2Utils, fromBlob_RankTwoTensor)
+{
+  std::vector<RankTwoTensor> data(2);
+  data[0] = RankTwoTensor(1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0);
+  data[1] = RankTwoTensor(-1.0, -2.0, -3.0, -4.0, -5.0, -6.0, -7.0, -8.0, -9.0);
+
+  auto tensor = NEML2Utils::fromBlob(data);
+  ASSERT_TRUE(tensor.defined());
+  ASSERT_TRUE(tensor.dim() == 3);
+  ASSERT_TRUE(tensor.batch_dim() == 1);
+  ASSERT_TRUE(tensor.base_dim() == 2);
+  ASSERT_TRUE(tensor.size(0) == 2);
+  ASSERT_TRUE(tensor.size(1) == 3);
+  ASSERT_TRUE(tensor.size(2) == 3);
+
+  for (neml2::Size n : index_range(data))
+    for (neml2::Size i : make_range(3))
+      for (neml2::Size j : make_range(3))
+        EXPECT_NEAR(tensor.index({n, i, j}).item<Real>(), data[n](i, j), 1e-12);
+}
+
+TEST(NEML2Utils, fromBlob_SymmetricRankTwoTensor)
+{
+  std::vector<SymmetricRankTwoTensor> data(2);
+  data[0] = SymmetricRankTwoTensor(1.0, 2.0, 3.0, 4.0, 5.0, 6.0);
+  data[1] = SymmetricRankTwoTensor(-1.0, -2.0, -3.0, -4.0, -5.0, -6.0);
+
+  auto tensor = NEML2Utils::fromBlob(data);
+  ASSERT_TRUE(tensor.defined());
+  ASSERT_TRUE(tensor.dim() == 2);
+  ASSERT_TRUE(tensor.batch_dim() == 1);
+  ASSERT_TRUE(tensor.base_dim() == 1);
+  ASSERT_TRUE(tensor.size(0) == 2);
+  ASSERT_TRUE(tensor.size(1) == 6);
+
+  for (neml2::Size n : index_range(data))
+    for (neml2::Size i : make_range(6))
+      EXPECT_NEAR(tensor.index({n, i}).item<Real>(), data[n](i), 1e-12);
+}
+
+TEST(NEML2Utils, copyTensorToMOOSEData_Real)
+{
+  Real data;
+  const auto tensor = torch::tensor(1.0, torch::TensorOptions().dtype(torch::kFloat64));
+  NEML2Utils::copyTensorToMOOSEData(tensor, data);
+
+  EXPECT_NEAR(tensor.item<Real>(), data, 1e-12);
+}
+
+TEST(NEML2Utils, copyTensorToMOOSEData_RealVectorValue)
+{
+  RealVectorValue data;
+  const auto tensor = torch::tensor({1.0, 2.0, 3.0}, torch::TensorOptions().dtype(torch::kFloat64));
+  NEML2Utils::copyTensorToMOOSEData(tensor, data);
+
+  for (neml2::Size i : make_range(3))
+    EXPECT_NEAR(tensor.index({i}).item<Real>(), data(i), 1e-12);
+}
+
+TEST(NEML2Utils, copyTensorToMOOSEData_RankTwoTensor)
+{
+  RankTwoTensor data;
+  const auto tensor = torch::tensor({{1.0, 2.0, 3.0}, {4.0, 5.0, 6.0}, {7.0, 8.0, 9.0}},
+                                    torch::TensorOptions().dtype(torch::kFloat64));
+  NEML2Utils::copyTensorToMOOSEData(tensor, data);
+
+  for (neml2::Size i : make_range(3))
+    for (neml2::Size j : make_range(3))
+      EXPECT_NEAR(tensor.index({i, j}).item<Real>(), data(i, j), 1e-12);
+}
+
+TEST(NEML2Utils, copyTensorToMOOSEData_SymmetricRankTwoTensor)
+{
+  SymmetricRankTwoTensor data;
+  const auto tensor =
+      torch::tensor({1.0, 2.0, 3.0, 4.0, 5.0, 6.0}, torch::TensorOptions().dtype(torch::kFloat64));
+  NEML2Utils::copyTensorToMOOSEData(tensor, data);
+
+  for (neml2::Size i : make_range(6))
+    EXPECT_NEAR(tensor.index({i}).item<Real>(), data(i), 1e-12);
+}
+
+TEST(NEML2Utils, copyTensorToMOOSEData_RankFourTensor)
+{
+  RankFourTensor data;
+  const auto tensor =
+      torch::tensor({{{{1.0, 2.0, 3.0}, {4.0, 5.0, 6.0}, {7.0, 8.0, 9.0}},
+                      {{10.0, 11.0, 12.0}, {13.0, 14.0, 15.0}, {16.0, 17.0, 18.0}},
+                      {{19.0, 20.0, 21.0}, {22.0, 23.0, 24.0}, {25.0, 26.0, 27.0}}},
+                     {{{28.0, 29.0, 30.0}, {31.0, 32.0, 33.0}, {34.0, 35.0, 36.0}},
+                      {{37.0, 38.0, 39.0}, {40.0, 41.0, 42.0}, {43.0, 44.0, 45.0}},
+                      {{46.0, 47.0, 48.0}, {49.0, 50.0, 51.0}, {52.0, 53.0, 54.0}}},
+                     {{{55.0, 56.0, 57.0}, {58.0, 59.0, 60.0}, {61.0, 62.0, 63.0}},
+                      {{64.0, 65.0, 66.0}, {67.0, 68.0, 69.0}, {70.0, 71.0, 72.0}},
+                      {{73.0, 74.0, 75.0}, {76.0, 77.0, 78.0}, {79.0, 80.0, 81.0}}}},
+                    torch::TensorOptions().dtype(torch::kFloat64));
+  NEML2Utils::copyTensorToMOOSEData(tensor, data);
+
+  for (neml2::Size i : make_range(3))
+    for (neml2::Size j : make_range(3))
+      for (neml2::Size k : make_range(3))
+        for (neml2::Size l : make_range(3))
+          EXPECT_NEAR(tensor.index({i, j, k, l}).item<Real>(), data(i, j, k, l), 1e-12);
+}
+
+TEST(NEML2Utils, copyTensorToMOOSEData_SymmetricRankFourTensor)
+{
+  SymmetricRankFourTensor data;
+  const auto tensor = torch::tensor({{1.0, 2.0, 3.0, 4.0, 5.0, 6.0},
+                                     {2.0, 3.0, 4.0, 5.0, 6.0, 7.0},
+                                     {3.0, 4.0, 5.0, 6.0, 7.0, 8.0},
+                                     {4.0, 5.0, 6.0, 7.0, 8.0, 9.0},
+                                     {5.0, 6.0, 7.0, 8.0, 9.0, 10.0},
+                                     {6.0, 7.0, 8.0, 9.0, 10.0, 11.0}},
+                                    torch::TensorOptions().dtype(torch::kFloat64));
+  NEML2Utils::copyTensorToMOOSEData(tensor, data);
+
+  for (neml2::Size i : make_range(6))
+    for (neml2::Size j : make_range(6))
+      EXPECT_NEAR(tensor.index({i, j}).item<Real>(), data(i, j), 1e-12);
+}
+#endif