From 6979e0f62c47ef363e2505e9d7ef089f4c65a5c5 Mon Sep 17 00:00:00 2001 From: petlenz Date: Sun, 6 Sep 2026 14:00:17 +0200 Subject: [PATCH 1/3] solvers: select the integration scheme by name, and declare it in the schema "tableau" was never declared in parameters() at all. The constructor read get_parameter("tableau") while the schema never inserted it, so it could only be supplied from C++, where parameter_handler's insert() bypasses the schema. From a document it failed with "Key tableau not found". That mattered more than a missing parameter usually would: the tableau IS the explicit-versus-implicit choice. j2_rk_plasticity already handles explicit stages (m_is_implicit[i] = |a_ii| >= 1e-30), so forward_euler was always an explicit J2 integrator -- the capability existed and was unreachable from a deck. The parameter is now a scheme NAME resolved by tableau_by_name(), and it is declared. A raw pointer cannot be expressed in a document, which is the same shape as Drucker-Prager's yield_function object (#33) and vector_newton's zero_blocks (#17): a parameter the code reads but the schema does not declare is invisible to both the document layer and to validation. rk_integrator had the identical defect and is fixed the same way; leaving one of the two would have been worse than not starting. An unknown name now throws and lists the valid schemes, rather than leaving the parameter absent. Two tests: every published scheme resolves and reports the right explicit/implicit character, and a typo names itself and the alternatives. Verified from JSON: j2_rk_plasticity (sdirk3, implicit) OK j2_rk_plasticity EXPLICIT (forward_euler) OK a typo in the scheme name unknown scheme 'sdirk4' -- expected one of: ... --- .../materials/j2_rk_plasticity.h | 17 +++-- .../solvers/butcher_tableau.h | 24 +++++++ .../numsim-materials/solvers/rk_integrator.h | 29 +++++--- tests/test_j2_plasticity.cpp | 14 ++-- tests/test_rk_integrator.cpp | 67 ++++++++++++++----- 5 files changed, 112 insertions(+), 39 deletions(-) diff --git a/include/numsim-materials/materials/j2_rk_plasticity.h b/include/numsim-materials/materials/j2_rk_plasticity.h index 0a7c57a..1811bfc 100644 --- a/include/numsim-materials/materials/j2_rk_plasticity.h +++ b/include/numsim-materials/materials/j2_rk_plasticity.h @@ -42,7 +42,8 @@ class j2_rk_plasticity final m_sigma_0(base::template get_parameter("sigma_0")), m_tol(base::template get_parameter("tolerance")), m_max_iter(base::template get_parameter("max_iter")), - m_tableau(base::template get_parameter("tableau")), + m_tableau(tableau_by_name( + base::template get_parameter("tableau"))), m_hardening_source(base::template get_parameter("hardening_source")), m_strain_source(base::template get_parameter("strain_source")), m_strain(base::template add_input( @@ -57,13 +58,13 @@ class j2_rk_plasticity final m_C_e = value_type{3} * m_K * IIvol + value_type{2} * m_G * plasticity_detail::make_IIdev(); - const int s = m_tableau->stages(); + 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 +76,12 @@ class j2_rk_plasticity final para.template insert("K").template add(); para.template insert("G").template add(); para.template insert("sigma_0").template add(); + // 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("tableau").template add(); para.template insert("tolerance") .template add(value_type{1e-12}); para.template insert("max_iter") @@ -102,7 +109,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 +196,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; diff --git a/include/numsim-materials/solvers/butcher_tableau.h b/include/numsim-materials/solvers/butcher_tableau.h index 3d17385..a2bb220 100644 --- a/include/numsim-materials/solvers/butcher_tableau.h +++ b/include/numsim-materials/solvers/butcher_tableau.h @@ -1,6 +1,8 @@ #ifndef NUMSIM_MATERIALS_BUTCHER_TABLEAU_H #define NUMSIM_MATERIALS_BUTCHER_TABLEAU_H +#include +#include #include 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 diff --git a/include/numsim-materials/solvers/rk_integrator.h b/include/numsim-materials/solvers/rk_integrator.h index 163b116..cef9ea7 100644 --- a/include/numsim-materials/solvers/rk_integrator.h +++ b/include/numsim-materials/solvers/rk_integrator.h @@ -37,28 +37,29 @@ class rk_integrator final m_h(base::template get_parameter("step_size")), m_tol(base::template get_parameter("tolerance")), m_max_iter(base::template get_parameter("max_iter")), - m_tableau(base::template get_parameter("tableau")), + m_tableau(tableau_by_name( + base::template get_parameter("tableau"))), m_func_name(base::template get_parameter("function")), m_rate(base::template add_input( 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( 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(value_type{1e-12}); para.template insert("max_iter") .template add(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("tableau").template add(); 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& m_rate; const input_property* m_drate; diff --git a/tests/test_j2_plasticity.cpp b/tests/test_j2_plasticity.cpp index b24f045..b3a356a 100644 --- a/tests/test_j2_plasticity.cpp +++ b/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("K", T{166.67}); p.insert("G", T{76.92}); p.insert("sigma_0", T{50.0}); - p.insert("tableau", &m_tab); + p.insert("tableau", m_tableau_name); ctx.create>(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(); diff --git a/tests/test_rk_integrator.cpp b/tests/test_rk_integrator.cpp index 2146537..c1f912b 100644 --- a/tests/test_rk_integrator.cpp +++ b/tests/test_rk_integrator.cpp @@ -60,7 +60,7 @@ class exponential_decay final /// Run exponential decay with a given integrator type and tableau. /// Returns y at t=1.0 with N steps of size h=1/N. template -T run_decay(int N, const numsim::materials::butcher_tableau& tab, T lambda = 1.0) { +T run_decay(int N, const std::string& tableau, T lambda = 1.0) { ctx_type ctx; param_type p; @@ -68,7 +68,7 @@ T run_decay(int N, const numsim::materials::butcher_tableau& tab, T lambda = 1.0 p.insert("name", "integrator"); p.insert("function", "decay"); p.insert("step_size", T{1.0} / T(N)); - p.insert("tableau", &tab); + p.insert("tableau", tableau); auto& integ = ctx.create(p); p.clear(); @@ -100,13 +100,13 @@ const T exact = std::exp(-1.0); // y(1) = e^(-1) ≈ 0.367879... using RK = numsim::materials::rk_integrator; TEST(ExplicitRK, ForwardEulerConverges) { - auto tab = numsim::materials::forward_euler(); + const std::string tab = "forward_euler"; auto y = run_decay(100, tab); EXPECT_NEAR(y, exact, 0.01) << "Forward Euler with 100 steps should be close"; } TEST(ExplicitRK, ForwardEulerOrder1) { - auto tab = numsim::materials::forward_euler(); + const std::string tab = "forward_euler"; auto err_10 = std::abs(run_decay(10, tab) - exact); auto err_20 = std::abs(run_decay(20, tab) - exact); auto ratio = err_10 / err_20; @@ -116,7 +116,7 @@ TEST(ExplicitRK, ForwardEulerOrder1) { } TEST(ExplicitRK, RK4Order4) { - auto tab = numsim::materials::rk4(); + const std::string tab = "rk4"; auto err_10 = std::abs(run_decay(10, tab) - exact); auto err_20 = std::abs(run_decay(20, tab) - exact); auto ratio = err_10 / err_20; @@ -126,7 +126,7 @@ TEST(ExplicitRK, RK4Order4) { } TEST(ExplicitRK, RK4HighAccuracy) { - auto tab = numsim::materials::rk4(); + const std::string tab = "rk4"; auto y = run_decay(100, tab); EXPECT_NEAR(y, exact, 1e-10) << "RK4 with 100 steps should be very accurate"; } @@ -135,13 +135,13 @@ TEST(ExplicitRK, RK4HighAccuracy) { TEST(DIRK, ImplicitEulerConverges) { - auto tab = numsim::materials::implicit_euler(); + const std::string tab = "implicit_euler"; auto y = run_decay(100, tab); EXPECT_NEAR(y, exact, 0.01) << "Implicit Euler with 100 steps"; } TEST(DIRK, ImplicitMidpointOrder2) { - auto tab = numsim::materials::implicit_midpoint(); + const std::string tab = "implicit_midpoint"; auto err_10 = std::abs(run_decay(10, tab) - exact); auto err_20 = std::abs(run_decay(20, tab) - exact); auto ratio = err_10 / err_20; @@ -151,7 +151,7 @@ TEST(DIRK, ImplicitMidpointOrder2) { } TEST(DIRK, CrankNicolsonOrder2) { - auto tab = numsim::materials::crank_nicolson(); + const std::string tab = "crank_nicolson"; auto err_10 = std::abs(run_decay(10, tab) - exact); auto err_20 = std::abs(run_decay(20, tab) - exact); auto ratio = err_10 / err_20; @@ -164,7 +164,7 @@ TEST(DIRK, CrankNicolsonOrder2) { TEST(ImplicitRK, GaussLegendreOrder4) { - auto tab = numsim::materials::gauss_legendre_4(); + const std::string tab = "gauss_legendre_4"; auto err_10 = std::abs(run_decay(10, tab) - exact); auto err_20 = std::abs(run_decay(20, tab) - exact); auto ratio = err_10 / err_20; @@ -174,7 +174,7 @@ TEST(ImplicitRK, GaussLegendreOrder4) { } TEST(ImplicitRK, GaussLegendreHighAccuracy) { - auto tab = numsim::materials::gauss_legendre_4(); + const std::string tab = "gauss_legendre_4"; auto y = run_decay(50, tab); EXPECT_NEAR(y, exact, 1e-10) << "Gauss-Legendre with 50 steps"; } @@ -184,7 +184,7 @@ TEST(ImplicitRK, GaussLegendreHighAccuracy) { // All should converge to z ≈ 1 after enough steps at 80°C. template -T run_curing(int N, const numsim::materials::butcher_tableau& tab, T step_size = T{10}) { +T run_curing(int N, const std::string& tableau, T step_size = T{10}) { ctx_type ctx; param_type p; @@ -199,7 +199,7 @@ T run_curing(int N, const numsim::materials::butcher_tableau& tab, T step_size = p.insert("name", "integrator"); p.insert("function", "curing_rate"); p.insert("step_size", step_size); - p.insert("tableau", &tab); + p.insert("tableau", tableau); ctx.create(p); // Curing rate function — reads state from integrator @@ -236,14 +236,14 @@ T run_curing(int N, const numsim::materials::butcher_tableau& tab, T step_size = } TEST(CuringRK, ExplicitRK4ConvergesToFullCure) { - auto tab = numsim::materials::rk4(); + const std::string tab = "rk4"; auto z = run_curing(50, tab); std::println(" RK4 curing (500s): z = {:.6f}", z); EXPECT_GT(z, 0.90) << "RK4 should approach full cure"; } TEST(CuringRK, DIRKImplicitMidpointConverges) { - auto tab = numsim::materials::implicit_midpoint(); + const std::string tab = "implicit_midpoint"; // Smaller step for implicit — stiff initial phase needs h < 1/df_dy auto z = run_curing(500, tab, T{1}); // h=1, 500 steps std::println(" Implicit midpoint curing (500s, h=1): z = {:.6f}", z); @@ -251,10 +251,45 @@ TEST(CuringRK, DIRKImplicitMidpointConverges) { } TEST(CuringRK, FullyImplicitGaussLegendreConverges) { - auto tab = numsim::materials::gauss_legendre_4(); + const std::string tab = "gauss_legendre_4"; auto z = run_curing(500, tab, T{1}); // h=1, 500 steps std::println(" Gauss-Legendre curing (500s, h=1): z = {:.6f}", z); EXPECT_GT(z, 0.90) << "Gauss-Legendre should approach full cure"; } } // namespace + +namespace { +namespace nm_tb = numsim::materials; + +/// The scheme must be selectable by NAME, because that is the explicit-vs- +/// implicit choice and it belongs in the deck rather than in a recompile. +TEST(TableauByName, ResolvesEveryPublishedScheme) { + for (const char* n : {"forward_euler", "explicit_midpoint", "rk4", + "implicit_euler", "implicit_midpoint", + "crank_nicolson", "sdirk3", "gauss_legendre_4"}) { + const auto t = nm_tb::tableau_by_name(n); + EXPECT_GT(t.stages(), 0) << n; + } + // and the explicit/implicit split really is what the name selects + EXPECT_TRUE(nm_tb::tableau_by_name("forward_euler").is_explicit()); + EXPECT_TRUE(nm_tb::tableau_by_name("rk4").is_explicit()); + EXPECT_FALSE(nm_tb::tableau_by_name("sdirk3").is_explicit()); + EXPECT_FALSE(nm_tb::tableau_by_name("gauss_legendre_4").is_explicit()); +} + +/// A typo must name itself and list the alternatives, not fall back to a +/// default scheme -- silently integrating with the wrong method would change +/// results without changing anything visible. +TEST(TableauByName, RejectsAnUnknownSchemeAndListsTheValidOnes) { + try { + nm_tb::tableau_by_name("sdirk4"); + FAIL() << "an unknown scheme must throw"; + } catch (const std::invalid_argument& e) { + const std::string msg = e.what(); + EXPECT_NE(msg.find("sdirk4"), std::string::npos) << msg; + EXPECT_NE(msg.find("sdirk3"), std::string::npos) + << "the message should list what IS valid: " << msg; + } +} +} // namespace From 41dc5e6179265b44f6ab73ed10c3769ac83ccd92 Mon Sep 17 00:00:00 2001 From: petlenz Date: Sun, 6 Sep 2026 14:16:21 +0200 Subject: [PATCH 2/3] core: report a wrongly-typed material reference as wrong, not as missing wire_materials() wrapped the look-up and the type-check in one catch(...), so material_ref::wire()'s precise "is not of the requested type" was discarded and the reference reported as missing instead. The effect on a user: naming a backward_euler where a local_newton is required produced wire_materials(): material 'j2' references missing materials: 'solver' sending them to look for a material sitting right there in their document. That is what makes it worth fixing rather than merely tidying -- the message does not just omit information, it points the wrong way. Now: wire_materials(): material 'j2' references materials of the wrong type: 'solver' (they exist, but are not the type this material requires) The two failures are collected separately and both are reported, so a graph with one of each still names both in one error -- the collect-then-report behaviour was deliberate and is kept. Two tests, one per failure mode, and the wrong-type one is verified against the mutation: restoring the single catch(...) fails it. --- .../core/material_interface.h | 30 +++++- tests/test_materials.cpp | 96 +++++++++++++++++++ 2 files changed, 121 insertions(+), 5 deletions(-) diff --git a/include/numsim-materials/core/material_interface.h b/include/numsim-materials/core/material_interface.h index a14760b..107b433 100644 --- a/include/numsim-materials/core/material_interface.h +++ b/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 missing; + std::vector 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 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); } } diff --git a/tests/test_materials.cpp b/tests/test_materials.cpp index 74f1c2c..465078c 100644 --- a/tests/test_materials.cpp +++ b/tests/test_materials.cpp @@ -5,6 +5,9 @@ #include "numsim-materials/materials/linear_elasticity.h" #include "numsim-materials/materials/scalar_stepper.h" #include "numsim-materials/solvers/local_newton.h" +#include "numsim-materials/materials/j2_plasticity.h" +#include "numsim-materials/materials/linear_isotropic_hardening.h" +#include "numsim-materials/materials/tensor_component_stepper.h" #include "numsim-materials/materials/scalar_identity_weight.h" #include "numsim-materials/materials/autocatalytic_reaction.h" #include "numsim-materials/solvers/backward_euler.h" @@ -224,3 +227,96 @@ TEST(LocalNewtonSolver, SolvesAndReportsConvergence) { EXPECT_FALSE(bad.converged); } } // namespace + +namespace { +namespace nm_w = numsim::materials; + +/// A wrongly-typed reference must say so, not report the material as missing. +/// +/// wire_materials() wrapped look-up and type-check in one catch(...), so +/// swapping a solver for one of another type told the user their solver did +/// not exist -- sending them to look for a material that was sitting right +/// there in the document. +TEST(WireMaterials, WrongTypeIsNotReportedAsMissing) { + using policy = nm_w::material_policy_default; + using T2 = policy::value_type; + nm_w::material_context ctx; + policy::ParameterHandler p; + + p.insert("name", "stepper"); + p.insert("increment", T2{0.01}); + p.insert>("indices", {0, 0}); + ctx.create>(p); + + p.clear(); + p.insert("name", "hardening"); + p.insert("source", "j2"); + p.insert("K", T2{1000.0}); + ctx.create>(p); + + // A backward_euler where a local_newton is required: it EXISTS. + p.clear(); + p.insert("name", "solver"); + p.insert("function", "j2"); + ctx.create>(p); + + p.clear(); + p.insert("name", "j2"); + p.insert("hardening_source", "hardening"); + p.insert("strain_source", "stepper"); + p.insert("solver_source", "solver"); + p.insert("K", T2{166.67}); + p.insert("G", T2{76.92}); + p.insert("sigma_0", T2{50.0}); + ctx.create>(p); + + try { + ctx.finalize(); + FAIL() << "wiring a wrongly-typed solver must fail"; + } catch (const std::runtime_error& e) { + const std::string msg = e.what(); + EXPECT_NE(msg.find("wrong type"), std::string::npos) << msg; + EXPECT_NE(msg.find("'solver'"), std::string::npos) << msg; + EXPECT_EQ(msg.find("do not exist"), std::string::npos) + << "the material exists; saying otherwise sends the user hunting for " + "something that is right there: " << msg; + } +} + +/// The genuinely-absent case still reports absence. +TEST(WireMaterials, MissingIsStillReportedAsMissing) { + using policy = nm_w::material_policy_default; + using T2 = policy::value_type; + nm_w::material_context ctx; + policy::ParameterHandler p; + + p.insert("name", "stepper"); + p.insert("increment", T2{0.01}); + p.insert>("indices", {0, 0}); + ctx.create>(p); + p.clear(); + p.insert("name", "hardening"); + p.insert("source", "j2"); + p.insert("K", T2{1000.0}); + ctx.create>(p); + p.clear(); + p.insert("name", "j2"); + p.insert("hardening_source", "hardening"); + p.insert("strain_source", "stepper"); + p.insert("solver_source", "no_such_solver"); + p.insert("K", T2{166.67}); + p.insert("G", T2{76.92}); + p.insert("sigma_0", T2{50.0}); + ctx.create>(p); + + try { + ctx.finalize(); + FAIL() << "a missing solver must fail"; + } catch (const std::runtime_error& e) { + const std::string msg = e.what(); + EXPECT_NE(msg.find("do not exist"), std::string::npos) << msg; + EXPECT_NE(msg.find("no_such_solver"), std::string::npos) << msg; + EXPECT_EQ(msg.find("wrong type"), std::string::npos) << msg; + } +} +} // namespace From de18dd7331a6564e9e218cf8227e8b9aed163ace Mon Sep 17 00:00:00 2001 From: petlenz Date: Sun, 6 Sep 2026 17:22:09 +0200 Subject: [PATCH 3/3] materials: refuse a fully implicit tableau in j2_rk_plasticity Making the scheme a named parameter turned an unreachable wrong answer into a deck-selectable one, and the commit that did it advertised the new reachability without checking the numbers. The stage loop accumulates a(i,j) only for j < i, plus the diagonal. Any coupling above the diagonal is dropped silently. gauss_legendre_4 has a(0,1) = 0.25 - 1/(2*sqrt(3)), so it is integrated as something that is neither Gauss-Legendre nor fourth order. Measured against the monolithic backward-Euler reference on a 60-step uniaxial path: implicit_euler alpha=0.034372 ratio 1.000 F +0.0000 sdirk3 alpha=0.034372 ratio 1.000 F -0.0000 forward_euler alpha=0.034528 ratio 1.005 F -0.1923 rk4 alpha=0.033750 ratio 0.982 F +0.7663 gauss_legendre_4 alpha=-8.805943 ratio -256 F +10880.3 A NEGATIVE equivalent plastic strain, no error, no warning. The explicit schemes' small offsets are expected -- they do not enforce the yield condition exactly -- but the last row is not an accuracy question. rk_integrator dispatches this case to a fully-implicit solve. This class has no such path, so it now refuses rather than pretending, naming the schemes it can integrate. Adding that path would be the better fix and is a larger change than this PR should carry. Two tests: gauss_legendre_4 is rejected, and every DIRK or explicit scheme is still accepted, so the guard cannot pass by refusing everything. Found by an independent review of this PR. The same review claimed a 3x error from j2_rk_plasticity.h:135 passing G where effective_modulus(G) = 3G is expected. The argument IS inconsistent with the other two call sites, but the measured effect is not 3x -- hardening dominates the denominator (3G = 231 against H' = 1000) and the incremental process self-corrects. Left alone here and reported separately rather than fixed on a claim I could not reproduce. --- .../materials/j2_rk_plasticity.h | 18 +++++ tests/test_j2_plasticity.cpp | 65 +++++++++++++++++++ 2 files changed, 83 insertions(+) diff --git a/include/numsim-materials/materials/j2_rk_plasticity.h b/include/numsim-materials/materials/j2_rk_plasticity.h index 1811bfc..a0257ad 100644 --- a/include/numsim-materials/materials/j2_rk_plasticity.h +++ b/include/numsim-materials/materials/j2_rk_plasticity.h @@ -2,6 +2,7 @@ #define NUMSIM_MATERIALS_J2_RK_PLASTICITY_H #include +#include #include #include #include "numsim-materials/core/material_base.h" @@ -58,6 +59,23 @@ class j2_rk_plasticity final m_C_e = value_type{3} * m_K * IIvol + value_type{2} * m_G * plasticity_detail::make_IIdev(); + // 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("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); diff --git a/tests/test_j2_plasticity.cpp b/tests/test_j2_plasticity.cpp index b3a356a..5bead51 100644 --- a/tests/test_j2_plasticity.cpp +++ b/tests/test_j2_plasticity.cpp @@ -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 ctx; + policy::ParameterHandler p; + p.insert("name", "stepper"); + p.insert("increment", T{0.01}); + p.insert>("indices", {0, 0}); + ctx.create>(p); + p.clear(); + p.insert("name", "hardening"); + p.insert("source", "m"); + p.insert("K", T{1000.0}); + ctx.create>(p); + p.clear(); + p.insert("name", "m"); + p.insert("hardening_source", "hardening"); + p.insert("strain_source", "stepper"); + p.insert("K", T{166.67}); + p.insert("G", T{76.92}); + p.insert("sigma_0", T{50.0}); + p.insert("tableau", std::string("gauss_legendre_4")); + EXPECT_THROW(ctx.create>(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 ctx; + policy::ParameterHandler p; + p.insert("name", "stepper"); + p.insert("increment", T{0.01}); + p.insert>("indices", {0, 0}); + ctx.create>(p); + p.clear(); + p.insert("name", "hardening"); + p.insert("source", "m"); + p.insert("K", T{1000.0}); + ctx.create>(p); + p.clear(); + p.insert("name", "m"); + p.insert("hardening_source", "hardening"); + p.insert("strain_source", "stepper"); + p.insert("K", T{166.67}); + p.insert("G", T{76.92}); + p.insert("sigma_0", T{50.0}); + p.insert("tableau", std::string(scheme)); + EXPECT_NO_THROW(ctx.create>(p)) + << scheme; + } +} +} // namespace