From 6892d40f948cd8e4d5111feb7c5501d261976dc1 Mon Sep 17 00:00:00 2001 From: petlenz Date: Sat, 5 Sep 2026 10:38:00 +0200 Subject: [PATCH 1/3] materials: drop the redundant stress contraction in the return mapping Profiling the step showed one rank-4 : rank-2 contraction at ~33 ns against ~245 ns for a whole J2 step, and the return map did two of them: once for the trial stress, once for the returned stress. The second is redundant: sigma = C_e : (eps - eps_p_old - dl N) = sig_trial - dl (C_e : N) and for isotropic C_e, C_e : N = 2G dev(N) + K tr(N) I -- the same identity the tangent collapse already uses. J2's flow is deviatoric so tr(N) = 0 and it reduces to 2G N; Drucker-Prager's is not, so the volumetric term stays. J2 (same-process A/B): 245.2 -> 195.4 ns/step DP (interleaved A/B): ~8%, medians 1037 -> 953 ns/step NOT bit-identical this time, unlike the earlier steps: the result agrees to about 1 ULP (s00 108.7500531151484 vs ...843) because the arithmetic is algebraically equal but ordered differently. All 50 tests pass, including the FD tangent checks on six load paths. The apex return keeps its contraction: there eps_p comes from apex_plastic_strain rather than dl*N, so the identity does not apply, and the branch is rare. Measurement note: absolute ns figures in this branch's earlier commits were single runs taken at different times, and this machine's load moves them by 30%+. The J2 comparison above is a same-process A/B; the DP one is interleaved across alternating runs. Ratios are meaningful, absolutes are not comparable across commits. --- .../materials/drucker_prager_plasticity.h | 13 ++++++++++++- include/numsim-materials/materials/j2_plasticity.h | 9 +++++++-- 2 files changed, 19 insertions(+), 3 deletions(-) diff --git a/include/numsim-materials/materials/drucker_prager_plasticity.h b/include/numsim-materials/materials/drucker_prager_plasticity.h index 7797fd8..89ef2a8 100644 --- a/include/numsim-materials/materials/drucker_prager_plasticity.h +++ b/include/numsim-materials/materials/drucker_prager_plasticity.h @@ -41,6 +41,7 @@ class drucker_prager_plasticity final m_kappa(base::template add_history_output("equivalent_plastic_strain")), m_G(base::template get_parameter("G")), m_sigma_0(base::template get_parameter("sigma_0")), + m_K_bulk(base::template get_parameter("K_bulk")), m_solver(base::template add_material_ref( base::template get_parameter("solver_source"))), m_C_e(base::template add_input( @@ -151,7 +152,16 @@ class drucker_prager_plasticity final const tensor4& C_e, value_type kappa_n, value_type dlambda) { m_eps_p.new_value() = m_eps_p.old_value() + dlambda * ts.N; m_kappa.new_value() = kappa_n + dlambda; - m_stress = tmech::dcontract(C_e, m_strain.get() - m_eps_p.new_value()); + + // sigma = C_e : (eps - eps_p_new) = sig_trial - dl (C_e : N). + // For isotropic C_e that is 2G dev(N) + K tr(N) I -- no rank-4 : rank-2 + // contraction, which measures 33 ns against ~720 for the whole step. The + // flow is non-associative, so tr(N) = beta is generally nonzero and the + // volumetric term does not drop out as it does for J2. + const auto I = tmech::eye(); + const tensor2 Ce_N{value_type{2} * m_G * tmech::dev(ts.N) + + m_K_bulk * tmech::trace(ts.N) * I}; + m_stress = ts.sig - dlambda * Ce_N; // Tangent at trial state (return mapping uses N_trial). // For J2, trial = converged. For DP, they differ. @@ -187,6 +197,7 @@ class drucker_prager_plasticity final const value_type& m_G; const value_type& m_sigma_0; + const value_type& m_K_bulk; material_ref& m_solver; const input_property& m_C_e; diff --git a/include/numsim-materials/materials/j2_plasticity.h b/include/numsim-materials/materials/j2_plasticity.h index 7322fd4..07b352e 100644 --- a/include/numsim-materials/materials/j2_plasticity.h +++ b/include/numsim-materials/materials/j2_plasticity.h @@ -143,8 +143,13 @@ class j2_plasticity final m_eps_p.new_value() = m_eps_p.old_value() + dlambda * N; m_kappa.new_value() = kappa_n + dlambda; - m_stress = tmech::dcontract( - C_e, tensor2(m_strain.get() - m_eps_p.new_value())); + + // sigma = C_e : (eps - eps_p_new) + // = C_e : (eps - eps_p_old) - dl (C_e : N) + // = sig_trial - 2G dl N, since C_e : N = 2G N. + // The same identity the tangent uses. Avoids a second rank-4 : rank-2 + // contraction, which measures 33 ns against ~245 for the whole step. + m_stress = sig_trial - value_type{2} * m_G * dlambda * N; m_H.update_source(); const auto GG = m_G * m_G; From 8b2fc73c7470e4143ecc778fb3090d371706b1a8 Mon Sep 17 00:00:00 2001 From: petlenz Date: Sat, 5 Sep 2026 11:07:57 +0200 Subject: [PATCH 2/3] materials: plasticity builds its own elastic tangent Both plasticity classes read the elastic stiffness from an elastic_source's "tangent" property. That looked like flexibility and was not: the closed forms REQUIRE an isotropic C_e -- it is what makes C_e : N = 2G N and N : C_e : N = 3G true -- so accepting an arbitrary rank-4 tangent advertised a generality neither class can honour. It also dragged a linear_elasticity into every plasticity graph, and that material's own "stress" output is not merely unused there, it is WRONG: it is C : eps with eps_p ignored, so once yielding starts it is not the stress of anything. Measured on a 200-step path it over-predicts by a growing margin: step 40 alpha 0.0094 j2 106.25 elastic 107.69 +1.4% step 120 alpha 0.1094 j2 306.25 elastic 323.08 +5.5% step 200 alpha 0.2094 j2 506.25 elastic 538.46 +6.4% A postprocessor logging elastic::stress from a plasticity graph gets that, under a name that reads as authoritative. Both classes now build C_e from moduli they already hold or now take: j2_plasticity gains a "K" parameter; drucker_prager_plasticity needs nothing new, since K_bulk and G were already required for the cone. elastic_source is gone from both, and a plasticity graph no longer contains an elastic material at all. J2: 195.4 -> 168.1 ns/step DP: ~720 -> ~608 ns/step Values agree with the pre-refactor implementation to about 1 ULP on alpha, stress and tangent. linear_elasticity is untouched and keeps its users: for a genuinely elastic model C : eps IS the answer, and isotropic_damage consumes both its stress and its tangent legitimately. rk_plasticity still takes an elastic_source; it is the remaining templated class and is left alone here. An earlier draft of this reasoning proposed dead-property elimination in the property engine, with a declared-outputs mechanism to tell an unread property from one a host reads through ctx.get(). That was solving the symptom. The property should not be in the graph. --- .../materials/drucker_prager_plasticity.h | 20 ++++++---- .../materials/j2_plasticity.h | 37 +++++++++++++------ tests/plot_data.cpp | 3 +- tests/test_drucker_prager.cpp | 6 --- tests/test_j2_plasticity.cpp | 8 ++-- 5 files changed, 44 insertions(+), 30 deletions(-) diff --git a/include/numsim-materials/materials/drucker_prager_plasticity.h b/include/numsim-materials/materials/drucker_prager_plasticity.h index 89ef2a8..67aed65 100644 --- a/include/numsim-materials/materials/drucker_prager_plasticity.h +++ b/include/numsim-materials/materials/drucker_prager_plasticity.h @@ -44,9 +44,6 @@ class drucker_prager_plasticity final m_K_bulk(base::template get_parameter("K_bulk")), m_solver(base::template add_material_ref( base::template get_parameter("solver_source"))), - m_C_e(base::template add_input( - base::template get_parameter("elastic_source"), - "tangent", EdgeKind::Global)), m_strain(base::template add_input( base::template get_parameter("strain_source"), "strain", EdgeKind::Global)), @@ -59,12 +56,15 @@ class drucker_prager_plasticity final { m_yf = yield_fn(base::template get_parameter("eta"), base::template get_parameter("beta"), - base::template get_parameter("K_bulk")); + m_K_bulk); + const auto I = tmech::eye(); + const tensor4 IIvol{tmech::otimes(I, I) / value_type{Dim}}; + m_C_e = value_type{3} * m_K_bulk * IIvol + + value_type{2} * m_G * plasticity_detail::make_IIdev(); } static input_parameter_controller parameters() { input_parameter_controller para{base::parameters()}; - para.template insert("elastic_source").template add(); para.template insert("hardening_source").template add(); para.template insert("strain_source").template add(); para.template insert("solver_source").template add(); @@ -81,7 +81,7 @@ class drucker_prager_plasticity final } void compute() { - const auto& C_e = m_C_e.get(); + const auto& C_e = m_C_e; const auto kappa_n = m_kappa.old_value(); m_kappa.new_value() = kappa_n; @@ -200,11 +200,17 @@ class drucker_prager_plasticity final const value_type& m_K_bulk; material_ref& m_solver; - const input_property& m_C_e; const input_property& m_strain; const input_property& m_H; const input_property& m_dH; yield_fn m_yf{}; + + /// The elastic stiffness, built here rather than read from another material. + /// The tangent collapse requires an isotropic C_e -- C_e : X = 2G X for + /// deviatoric X -- so accepting an arbitrary rank-4 tangent advertised a + /// generality this material cannot honour. K_bulk and G are already + /// parameters, so nothing new is asked of the caller. + tensor4 m_C_e{}; }; diff --git a/include/numsim-materials/materials/j2_plasticity.h b/include/numsim-materials/materials/j2_plasticity.h index 07b352e..dee151c 100644 --- a/include/numsim-materials/materials/j2_plasticity.h +++ b/include/numsim-materials/materials/j2_plasticity.h @@ -45,9 +45,10 @@ namespace numsim::materials { /// may be nonlinear (see exponential_isotropic_hardening). /// /// Parameters: -/// "name", "elastic_source", "hardening_source", "strain_source", -/// "solver_source", "G", "sigma_0" -/// -- the same set small_strain_plasticity takes, so this is a drop-in. +/// "name", "hardening_source", "strain_source", "solver_source", +/// "K", "G", "sigma_0" +/// +/// No elastic_source: the stiffness is built from K and G here. template class j2_plasticity final : public material_base, Traits> { @@ -71,11 +72,9 @@ class j2_plasticity final "equivalent_plastic_strain")), m_G(base::template get_parameter("G")), m_sigma_0(base::template get_parameter("sigma_0")), + m_K(base::template get_parameter("K")), m_solver(base::template add_material_ref( base::template get_parameter("solver_source"))), - m_C_e(base::template add_input( - base::template get_parameter("elastic_source"), - "tangent", EdgeKind::Global)), m_strain(base::template add_input( base::template get_parameter("strain_source"), "strain", EdgeKind::Global)), @@ -85,26 +84,41 @@ class j2_plasticity final m_dH(base::template add_input( base::template get_parameter("hardening_source"), "hardening_modulus", EdgeKind::Local)), - m_IIdev(plasticity_detail::make_IIdev()) + m_IIdev(plasticity_detail::make_IIdev()), + m_C_e(build_elastic_tangent(m_K, m_G, m_IIdev)) {} + /// The elastic stiffness, built here rather than read from another material. + /// + /// The closed forms below REQUIRE an isotropic C_e -- that is what makes + /// C_e : N = 2G N and N : C_e : N = 3G true. Accepting an arbitrary rank-4 + /// tangent from an elastic_source advertised a generality this material + /// cannot honour, and dragged a linear_elasticity into every plasticity graph + /// whose own "stress" output (C : eps, ignoring eps_p) is meaningless once + /// yielding starts. + static tensor4 build_elastic_tangent(value_type K, value_type G, + const tensor4& IIdev) { + const auto I = tmech::eye(); + return value_type{3} * K * (tmech::otimes(I, I) / value_type{Dim}) + + value_type{2} * G * IIdev; + } + static input_parameter_controller parameters() { input_parameter_controller para{base::parameters()}; - para.template insert("elastic_source") - .template add(); para.template insert("hardening_source") .template add(); para.template insert("strain_source") .template add(); para.template insert("solver_source") .template add(); + para.template insert("K").template add(); para.template insert("G").template add(); para.template insert("sigma_0").template add(); return para; } void compute() { - const auto& C_e = m_C_e.get(); + const auto& C_e = m_C_e; const auto kappa_n = m_kappa.old_value(); m_kappa.new_value() = kappa_n; @@ -167,14 +181,15 @@ class j2_plasticity final const value_type& m_G; const value_type& m_sigma_0; + const value_type& m_K; material_ref& m_solver; - const input_property& m_C_e; const input_property& m_strain; const input_property& m_H; const input_property& m_dH; const tensor4 m_IIdev; + const tensor4 m_C_e; }; } // namespace numsim::materials diff --git a/tests/plot_data.cpp b/tests/plot_data.cpp index b32ed3f..b32db4b 100644 --- a/tests/plot_data.cpp +++ b/tests/plot_data.cpp @@ -78,10 +78,10 @@ run_result run_j2(T increment, int steps, p.clear(); p.insert("name", "j2"); - p.insert("elastic_source", "elastic"); p.insert("hardening_source", "hardening"); p.insert("strain_source", "stepper"); p.insert("solver_source", "solver"); + p.insert("K", K_val); p.insert("G", G_val); p.insert("sigma_0", sigma_0); ctx.create(p); @@ -144,7 +144,6 @@ run_result run_dp(T increment, int steps, p.clear(); p.insert("name", "dp"); - p.insert("elastic_source", "elastic"); p.insert("hardening_source", "hardening"); p.insert("strain_source", "stepper"); p.insert("solver_source", "solver"); diff --git a/tests/test_drucker_prager.cpp b/tests/test_drucker_prager.cpp index aa8f4fc..c495c32 100644 --- a/tests/test_drucker_prager.cpp +++ b/tests/test_drucker_prager.cpp @@ -116,7 +116,6 @@ class DruckerPragerTest : public ::testing::Test { p.clear(); p.insert("name", "dp"); - p.insert("elastic_source", "elastic"); p.insert("hardening_source", "hardening"); p.insert("strain_source", "stepper"); p.insert("solver_source", "solver"); @@ -221,7 +220,6 @@ class DPTangentTest : public ::testing::Test { p.clear(); p.insert("name", "dp"); - p.insert("elastic_source", "elastic"); p.insert("hardening_source", "hardening"); p.insert("strain_source", "stepper"); p.insert("solver_source", "solver"); @@ -295,7 +293,6 @@ T run_dp_max_tangent_error(T increment, int steps) { p.clear(); p.insert("name", "dp"); - p.insert("elastic_source", "elastic"); p.insert("hardening_source", "hardening"); p.insert("strain_source", "stepper"); p.insert("solver_source", "solver"); @@ -383,7 +380,6 @@ T max_tangent_error(std::vector direction, T increment, int steps) { p.clear(); p.insert("name", "dp"); - p.insert("elastic_source", "elastic"); p.insert("hardening_source", "hardening"); p.insert("strain_source", "stepper"); p.insert("solver_source", "solver"); @@ -481,7 +477,6 @@ TEST(DruckerPragerApex, HydrostaticTensionReachesTheApex) { p.clear(); p.insert("name", "dp"); - p.insert("elastic_source", "elastic"); p.insert("hardening_source", "hardening"); p.insert("strain_source", "stepper"); p.insert("solver_source", "solver"); @@ -540,7 +535,6 @@ TEST(DruckerPragerApex, ApexStateIsAdmissible) { ctx.create>(p); p.clear(); p.insert("name", "dp"); - p.insert("elastic_source", "elastic"); p.insert("hardening_source", "hardening"); p.insert("strain_source", "stepper"); p.insert("solver_source", "solver"); diff --git a/tests/test_j2_plasticity.cpp b/tests/test_j2_plasticity.cpp index 03ff086..bc88a7a 100644 --- a/tests/test_j2_plasticity.cpp +++ b/tests/test_j2_plasticity.cpp @@ -56,10 +56,10 @@ class J2PlasticityTest : public ::testing::Test { // J2 plasticity — solver passed as pointer p.clear(); p.insert("name", "j2"); - p.insert("elastic_source", "elastic"); - p.insert("hardening_source", "hardening"); + p.insert("hardening_source", "hardening"); p.insert("strain_source", "stepper"); p.insert("solver_source", "solver"); + p.insert("K", K); p.insert("G", G); p.insert("sigma_0", sigma_0); ctx.create>(p); @@ -168,10 +168,10 @@ class J2TangentTest : public ::testing::Test { p.clear(); p.insert("name", "j2"); - p.insert("elastic_source", "elastic"); - p.insert("hardening_source", "hardening"); + p.insert("hardening_source", "hardening"); p.insert("strain_source", "stepper"); p.insert("solver_source", "solver"); + p.insert("K", T{166.67}); p.insert("G", T{76.92}); p.insert("sigma_0", T{50.0}); ctx.create>(p); From 384f9310c09a8dfec9357b40706456c381060ae7 Mon Sep 17 00:00:00 2001 From: petlenz Date: Sat, 5 Sep 2026 11:13:58 +0200 Subject: [PATCH 3/3] materials: rk_plasticity builds its own elastic tangent too Same change as the other two plasticity classes, for the same reason. rk_plasticity read the stiffness from an elastic_source while assuming isotropy twice over: effective_modulus(G) = 3G in its stage residuals, and compute_tangent's C_e : X = 2G X collapse. An arbitrary rank-4 tangent could not have been honoured by either. Takes "K" and builds C_e in the constructor; elastic_source is gone. Bit- identical over a 40-step SDIRK3 path, all 17 digits: before: 0.01062302967272258 108.7500531151484 250.00062497968736 76.554019397938148 after : 0.01062302967272258 108.7500531151484 250.00062497968736 76.554019397938148 That comparison matters here because the suite's rk-vs-j2 equivalence test now has both sides changed; this one holds rk against its own pre-change output. SDIRK3TangentCheck, which compares the tangent against a numerical derivative, is the independent check and passes unchanged. No plasticity material takes an elastic_source now. The one remaining consumer is isotropic_damage, which reads elastic::stress legitimately -- damage scales an elastic stress, and with no plastic strain in that model C : eps IS the stress. --- .../materials/rk_plasticity.h | 23 +++++++++++++------ tests/test_j2_plasticity.cpp | 2 +- 2 files changed, 17 insertions(+), 8 deletions(-) diff --git a/include/numsim-materials/materials/rk_plasticity.h b/include/numsim-materials/materials/rk_plasticity.h index cf5b0b1..05d40bd 100644 --- a/include/numsim-materials/materials/rk_plasticity.h +++ b/include/numsim-materials/materials/rk_plasticity.h @@ -37,16 +37,14 @@ class rk_plasticity final m_tangent(base::template add_output("tangent")), m_eps_p(base::template add_history_output("plastic_strain")), m_kappa(base::template add_history_output("equivalent_plastic_strain")), + m_K(base::template get_parameter("K")), m_G(base::template get_parameter("G")), 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_elastic_source(base::template get_parameter("elastic_source")), m_hardening_source(base::template get_parameter("hardening_source")), m_strain_source(base::template get_parameter("strain_source")), - m_C_e(base::template add_input( - m_elastic_source, "tangent", EdgeKind::Global)), m_strain(base::template add_input( m_strain_source, "strain", EdgeKind::Global)), m_H(base::template add_input( @@ -54,6 +52,11 @@ class rk_plasticity final m_dH(base::template add_input( m_hardening_source, "hardening_modulus", EdgeKind::Local)) { + const auto I = tmech::eye(); + const tensor4 IIvol{tmech::otimes(I, I) / value_type{Dim}}; + m_C_e = value_type{3} * m_K * IIvol + + value_type{2} * m_G * plasticity_detail::make_IIdev(); + const int s = m_tableau->stages(); m_dlambda.resize(s, value_type{0}); m_N_stage.resize(s); @@ -67,9 +70,9 @@ class rk_plasticity final static input_parameter_controller parameters() { input_parameter_controller para{base::parameters()}; - para.template insert("elastic_source").template add(); para.template insert("hardening_source").template add(); para.template insert("strain_source").template add(); + para.template insert("K").template add(); para.template insert("G").template add(); para.template insert("sigma_0").template add(); para.template insert("tolerance") @@ -80,7 +83,7 @@ class rk_plasticity final } void compute() { - const auto& C_e = m_C_e.get(); + const auto& C_e = m_C_e; const auto& eps = m_strain.get(); const auto kappa_n = m_kappa.old_value(); const auto eps_p_n = m_eps_p.old_value(); @@ -181,20 +184,26 @@ class rk_plasticity final history_property& m_eps_p; history_property& m_kappa; + const value_type& m_K; const value_type& m_G; const value_type& m_sigma_0; const value_type& m_tol; const int& m_max_iter; const butcher_tableau* m_tableau; - const std::string& m_elastic_source; const std::string& m_hardening_source; const std::string& m_strain_source; - const input_property& m_C_e; const input_property& m_strain; const input_property& m_H; const input_property& m_dH; + /// Built here, not read from an elastic_source. compute_tangent's collapse + /// (C_e : X = 2G X) and effective_modulus(G) = 3G both require an isotropic + /// C_e, so taking an arbitrary rank-4 tangent promised more than this class + /// can deliver -- and pulled a linear_elasticity into the graph whose own + /// "stress" output is meaningless once eps_p is nonzero. + tensor4 m_C_e{}; + yield_fn m_yf{}; std::vector m_dlambda; std::vector m_N_stage; diff --git a/tests/test_j2_plasticity.cpp b/tests/test_j2_plasticity.cpp index bc88a7a..a7b3429 100644 --- a/tests/test_j2_plasticity.cpp +++ b/tests/test_j2_plasticity.cpp @@ -258,9 +258,9 @@ class RKPlasticityTest : public ::testing::Test { p.clear(); p.insert("name", "j2"); - p.insert("elastic_source", "elastic"); 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", &m_tab);