solvers: select the integration scheme by name (explicit vs implicit from the deck)

include/numsim-materials/core/material_interface.h

@@ -70,19 +70,39 @@ class material_interface {
   /// Collects ALL missing materials and reports them in one error.
   void wire_materials(material_handler& handler) {
     std::vector<std::string> missing;
+    std::vector<std::string> wrong_type;
     for (auto& ref : m_material_refs) {
+      // Look-up and type-check are separate failures and must stay separate.
+      // A single catch(...) around both reported a wrongly-typed material as
+      // MISSING, so swapping a solver for one of another type told the user
+      // their solver did not exist.
+      material_interface* target = nullptr;
       try {
         auto& any_ref = handler.get(ref->target_name());
-        auto& mat = std::any_cast<
+        target = &std::any_cast<
             std::reference_wrapper<material_interface> const&>(any_ref).get();
-        ref->wire(mat);
       } catch (...) {
         missing.push_back(ref->target_name());
+        continue;
+      }
+      try {
+        ref->wire(*target);
+      } catch (...) {
+        wrong_type.push_back(ref->target_name());
       }
     }
-    if (!missing.empty()) {
-      std::string msg = "wire_materials(): material '" + m_name + "' references missing materials:";
-      for (auto& name : missing) msg += " '" + name + "'";
+    if (!missing.empty() || !wrong_type.empty()) {
+      std::string msg = "wire_materials(): material '" + m_name + "'";
+      if (!missing.empty()) {
+        msg += " references materials that do not exist:";
+        for (auto& name : missing) msg += " '" + name + "'";
+      }
+      if (!wrong_type.empty()) {
+        if (!missing.empty()) msg += ";";
+        msg += " references materials of the wrong type:";
+        for (auto& name : wrong_type) msg += " '" + name + "'";
+        msg += " (they exist, but are not the type this material requires)";
+      }
       throw std::runtime_error(msg);
     }
   }

include/numsim-materials/materials/j2_rk_plasticity.h

@@ -2,6 +2,7 @@
 #define NUMSIM_MATERIALS_J2_RK_PLASTICITY_H
 
 #include <cmath>
+#include <stdexcept>
 #include <vector>
 #include <tmech/tmech.h>
 #include "numsim-materials/core/material_base.h"
@@ -42,7 +43,8 @@ class j2_rk_plasticity final
         m_sigma_0(base::template get_parameter<value_type>("sigma_0")),
         m_tol(base::template get_parameter<value_type>("tolerance")),
         m_max_iter(base::template get_parameter<int>("max_iter")),
-        m_tableau(base::template get_parameter<const butcher_tableau*>("tableau")),
+        m_tableau(tableau_by_name(
+            base::template get_parameter<std::string>("tableau"))),
         m_hardening_source(base::template get_parameter<std::string>("hardening_source")),
         m_strain_source(base::template get_parameter<std::string>("strain_source")),
         m_strain(base::template add_input<tensor2>(
@@ -57,13 +59,30 @@ class j2_rk_plasticity final
     m_C_e = value_type{3} * m_K * IIvol +
             value_type{2} * m_G * plasticity_detail::make_IIdev<value_type, Dim>();
 
-    const int s = m_tableau->stages();
+    // The stage loop below accumulates only a(i,j) for j < i, plus the
+    // diagonal. A tableau with a(i,j) != 0 for j > i would have those terms
+    // silently DROPPED -- integrating with a method that is not the one named.
+    // Measured with gauss_legendre_4 before this guard: equivalent plastic
+    // strain of -8.8 (negative) and a yield residual of +10880.
+    //
+    // rk_integrator dispatches that case to a fully-implicit solve; this class
+    // has no such path, so it refuses rather than pretending.
+    if (!m_tableau.is_dirk())
+      throw std::invalid_argument(
+          "j2_rk_plasticity: the scheme '" +
+          base::template get_parameter<std::string>("tableau") +
+          "' is fully implicit (it has coupling above the diagonal), which this "
+          "return map cannot integrate -- its stage loop sums only j < i. Use a "
+          "DIRK or explicit scheme: forward_euler, explicit_midpoint, rk4, "
+          "implicit_euler, implicit_midpoint, crank_nicolson, sdirk3");
+
+    const int s = m_tableau.stages();
     m_dlambda.resize(s, value_type{0});
     m_N_stage.resize(s);
     m_is_implicit.resize(s);
     m_diag.resize(s);
     for (int i = 0; i < s; ++i) {
-      m_diag[i] = m_tableau->a(i, i);
+      m_diag[i] = m_tableau.a(i, i);
       m_is_implicit[i] = std::abs(m_diag[i]) >= 1e-30;
     }
   }
@@ -75,6 +94,12 @@ class j2_rk_plasticity final
     para.template insert<value_type>("K").template add<is_required>();
     para.template insert<value_type>("G").template add<is_required>();
     para.template insert<value_type>("sigma_0").template add<is_required>();
+    // The scheme by NAME. This parameter was never in the schema at all --
+    // the constructor read it while parameters() never declared it, so it
+    // could only ever be supplied from C++, where insert() bypasses the
+    // schema. That is why the explicit/implicit choice, which is the whole
+    // point of a tableau, was unreachable from a document.
+    para.template insert<std::string>("tableau").template add<is_required>();
     para.template insert<value_type>("tolerance")
         .template add<set_default>(value_type{1e-12});
     para.template insert<int>("max_iter")
@@ -102,7 +127,7 @@ class j2_rk_plasticity final
     }
 
     // Multi-stage return mapping
-    const auto& tab = *m_tableau;
+    const auto& tab = m_tableau;
     const int s = tab.stages();
 
     for (int i = 0; i < s; ++i)
@@ -189,7 +214,7 @@ class j2_rk_plasticity final
   const value_type& m_sigma_0;
   const value_type& m_tol;
   const int& m_max_iter;
-  const butcher_tableau* m_tableau;
+  const butcher_tableau m_tableau;
   const std::string& m_hardening_source;
   const std::string& m_strain_source;
 

include/numsim-materials/solvers/butcher_tableau.h

@@ -1,6 +1,8 @@
 #ifndef NUMSIM_MATERIALS_BUTCHER_TABLEAU_H
 #define NUMSIM_MATERIALS_BUTCHER_TABLEAU_H
 
+#include <stdexcept>
+#include <string>
 #include <Eigen/Dense>
 
 namespace numsim::materials {
@@ -118,6 +120,28 @@ inline butcher_tableau gauss_legendre_4() {
   return t;
 }
 
+/// Look a tableau up by name, so the integrator can be chosen from a document
+/// rather than by passing a pointer from C++.
+///
+/// The scheme IS the choice between explicit and implicit time integration --
+/// forward_euler and rk4 are explicit, sdirk3 and gauss_legendre_4 implicit --
+/// and that choice belongs in the deck, not in a recompile.
+inline butcher_tableau tableau_by_name(const std::string& name) {
+  if (name == "forward_euler")     return forward_euler();
+  if (name == "explicit_midpoint") return explicit_midpoint();
+  if (name == "rk4")               return rk4();
+  if (name == "implicit_euler")    return implicit_euler();
+  if (name == "implicit_midpoint") return implicit_midpoint();
+  if (name == "crank_nicolson")    return crank_nicolson();
+  if (name == "sdirk3")            return sdirk3();
+  if (name == "gauss_legendre_4")  return gauss_legendre_4();
+  throw std::invalid_argument(
+      "butcher_tableau: unknown scheme '" + name +
+      "' -- expected one of: forward_euler, explicit_midpoint, rk4, "
+      "implicit_euler, implicit_midpoint, crank_nicolson, sdirk3, "
+      "gauss_legendre_4");
+}
+
 } // namespace numsim::materials
 
 #endif // NUMSIM_MATERIALS_BUTCHER_TABLEAU_H

include/numsim-materials/solvers/rk_integrator.h

@@ -37,28 +37,29 @@ class rk_integrator final
         m_h(base::template get_parameter<value_type>("step_size")),
         m_tol(base::template get_parameter<value_type>("tolerance")),
         m_max_iter(base::template get_parameter<int>("max_iter")),
-        m_tableau(base::template get_parameter<const butcher_tableau*>("tableau")),
+        m_tableau(tableau_by_name(
+            base::template get_parameter<std::string>("tableau"))),
         m_func_name(base::template get_parameter<std::string>("function")),
         m_rate(base::template add_input<value_type>(
             m_func_name, "rate", EdgeKind::Local)),
         // rate_derivative only needed for implicit stages — not created for
         // explicit tableaux (the rate function may not provide it).
         // Safe: compute_explicit() never dereferences m_drate.
-        m_drate(m_tableau->is_explicit()
+        m_drate(m_tableau.is_explicit()
             ? nullptr
             : &base::template add_input<value_type>(
                 m_func_name, "rate_derivative", EdgeKind::Local)),
-        m_k(Eigen::VectorXd::Zero(m_tableau->stages()))
+        m_k(Eigen::VectorXd::Zero(m_tableau.stages()))
   {
-    const int s = m_tableau->stages();
-    m_is_explicit = m_tableau->is_explicit();
-    m_is_dirk = m_tableau->is_dirk();
+    const int s = m_tableau.stages();
+    m_is_explicit = m_tableau.is_explicit();
+    m_is_dirk = m_tableau.is_dirk();
 
     // Pre-compute diagonal properties for DIRK
     m_diag.resize(s);
     m_stage_implicit.resize(s);
     for (int i = 0; i < s; ++i) {
-      m_diag[i] = m_tableau->a(i, i);
+      m_diag[i] = m_tableau.a(i, i);
       m_stage_implicit[i] = std::abs(m_diag[i]) >= 1e-30;
     }
 
@@ -78,6 +79,12 @@ class rk_integrator final
         .template add<set_default>(value_type{1e-12});
     para.template insert<int>("max_iter")
         .template add<set_default>(int{50});
+        // The scheme by NAME. This parameter was never in the schema at all --
+    // the constructor read it while parameters() never declared it, so it
+    // could only ever be supplied from C++, where insert() bypasses the
+    // schema. That is why the explicit/implicit choice, which is the whole
+    // point of a tableau, was unreachable from a document.
+    para.template insert<std::string>("tableau").template add<is_required>();
     return para;
   }
 
@@ -89,7 +96,7 @@ class rk_integrator final
 
 private:
   void compute_explicit() {
-    const auto& tab = *m_tableau;
+    const auto& tab = m_tableau;
     const int s = tab.stages();
     const auto y_n = m_state.old_value();
     m_k.setZero();
@@ -105,7 +112,7 @@ class rk_integrator final
   }
 
   void compute_dirk() {
-    const auto& tab = *m_tableau;
+    const auto& tab = m_tableau;
     const int s = tab.stages();
     const auto y_n = m_state.old_value();
     m_k.setZero();
@@ -137,7 +144,7 @@ class rk_integrator final
   }
 
   void compute_fully_implicit() {
-    const auto& tab = *m_tableau;
+    const auto& tab = m_tableau;
     const int s = tab.stages();
     const auto y_n = m_state.old_value();
     m_k.setZero();
@@ -163,7 +170,7 @@ class rk_integrator final
   const value_type& m_h;
   const value_type& m_tol;
   const int& m_max_iter;
-  const butcher_tableau* m_tableau;
+  const butcher_tableau m_tableau;
   const std::string& m_func_name;
   const input_property<value_type, property_traits>& m_rate;
   const input_property<value_type, property_traits>* m_drate;

tests/test_j2_plasticity.cpp

@@ -233,8 +233,8 @@ TEST_F(J2TangentTest, ConsistentTangentAllSteps) {
 
 class RKPlasticityTest : public ::testing::Test {
 protected:
-  void setup_with_tableau(const numsim::materials::butcher_tableau& tab) {
-    m_tab = tab;
+  void setup_with_tableau(const std::string& name) {
+    m_tableau_name = name;
     param_type p;
 
     p.clear();
@@ -263,7 +263,7 @@ class RKPlasticityTest : public ::testing::Test {
     p.insert<T>("K", T{166.67});
     p.insert<T>("G", T{76.92});
     p.insert<T>("sigma_0", T{50.0});
-    p.insert<const numsim::materials::butcher_tableau*>("tableau", &m_tab);
+    p.insert<std::string>("tableau", m_tableau_name);
     ctx.create<numsim::materials::j2_rk_plasticity<policy>>(p);
 
     p.clear();
@@ -281,11 +281,11 @@ class RKPlasticityTest : public ::testing::Test {
   }
 
   ctx_type ctx;
-  numsim::materials::butcher_tableau m_tab;
+  std::string m_tableau_name;
 };
 
 TEST_F(RKPlasticityTest, ImplicitEulerMatchesMonolithic) {
-  setup_with_tableau(numsim::materials::implicit_euler());
+  setup_with_tableau("implicit_euler");
   T max_rel_error = 0;
   for (int i = 0; i < 20; ++i) {
     ctx.update();
@@ -300,7 +300,7 @@ TEST_F(RKPlasticityTest, ImplicitEulerMatchesMonolithic) {
 }
 
 TEST_F(RKPlasticityTest, SDIRK3TangentCheck) {
-  setup_with_tableau(numsim::materials::sdirk3());
+  setup_with_tableau("sdirk3");
   T max_rel_error = 0;
   for (int i = 0; i < 20; ++i) {
     ctx.update();
@@ -315,7 +315,7 @@ TEST_F(RKPlasticityTest, SDIRK3TangentCheck) {
 }
 
 TEST_F(RKPlasticityTest, PlasticStrainAccumulates) {
-  setup_with_tableau(numsim::materials::implicit_euler());
+  setup_with_tableau("implicit_euler");
   T prev_alpha = 0;
   for (int i = 0; i < 20; ++i) {
     ctx.update();
@@ -328,3 +328,68 @@ TEST_F(RKPlasticityTest, PlasticStrainAccumulates) {
 }
 
 } // namespace
+
+namespace {
+/// A fully implicit tableau must be refused, not silently mis-integrated.
+///
+/// The stage loop sums only a(i,j) for j < i plus the diagonal, so coupling
+/// above the diagonal is dropped. gauss_legendre_4 has a(0,1) != 0; before the
+/// guard it produced an equivalent plastic strain of -8.8 -- NEGATIVE -- and a
+/// yield residual of +10880, with no error. It became reachable from a deck
+/// when the tableau turned into a named parameter.
+TEST(J2RKScheme, RefusesAFullyImplicitTableau) {
+  using policy = numsim::materials::material_policy_default;
+  numsim::materials::material_context<policy> ctx;
+  policy::ParameterHandler p;
+  p.insert<std::string>("name", "stepper");
+  p.insert<T>("increment", T{0.01});
+  p.insert<std::vector<std::size_t>>("indices", {0, 0});
+  ctx.create<numsim::materials::tensor_component_stepper<2, policy>>(p);
+  p.clear();
+  p.insert<std::string>("name", "hardening");
+  p.insert<std::string>("source", "m");
+  p.insert<T>("K", T{1000.0});
+  ctx.create<numsim::materials::linear_isotropic_hardening<policy>>(p);
+  p.clear();
+  p.insert<std::string>("name", "m");
+  p.insert<std::string>("hardening_source", "hardening");
+  p.insert<std::string>("strain_source", "stepper");
+  p.insert<T>("K", T{166.67});
+  p.insert<T>("G", T{76.92});
+  p.insert<T>("sigma_0", T{50.0});
+  p.insert<std::string>("tableau", std::string("gauss_legendre_4"));
+  EXPECT_THROW(ctx.create<numsim::materials::j2_rk_plasticity<policy>>(p),
+               std::invalid_argument);
+}
+
+/// The DIRK and explicit schemes it CAN integrate must still be accepted, so
+/// the guard cannot pass by refusing everything.
+TEST(J2RKScheme, AcceptsEveryDirkAndExplicitScheme) {
+  using policy = numsim::materials::material_policy_default;
+  for (const char* scheme : {"forward_euler", "explicit_midpoint", "rk4",
+                             "implicit_euler", "implicit_midpoint",
+                             "crank_nicolson", "sdirk3"}) {
+    numsim::materials::material_context<policy> ctx;
+    policy::ParameterHandler p;
+    p.insert<std::string>("name", "stepper");
+    p.insert<T>("increment", T{0.01});
+    p.insert<std::vector<std::size_t>>("indices", {0, 0});
+    ctx.create<numsim::materials::tensor_component_stepper<2, policy>>(p);
+    p.clear();
+    p.insert<std::string>("name", "hardening");
+    p.insert<std::string>("source", "m");
+    p.insert<T>("K", T{1000.0});
+    ctx.create<numsim::materials::linear_isotropic_hardening<policy>>(p);
+    p.clear();
+    p.insert<std::string>("name", "m");
+    p.insert<std::string>("hardening_source", "hardening");
+    p.insert<std::string>("strain_source", "stepper");
+    p.insert<T>("K", T{166.67});
+    p.insert<T>("G", T{76.92});
+    p.insert<T>("sigma_0", T{50.0});
+    p.insert<std::string>("tableau", std::string(scheme));
+    EXPECT_NO_THROW(ctx.create<numsim::materials::j2_rk_plasticity<policy>>(p))
+        << scheme;
+  }
+}
+}  // namespace