From dd4c3b394549b7921f0a246076d9eb262fb6d7a8 Mon Sep 17 00:00:00 2001 From: petlenz Date: Sat, 5 Sep 2026 11:47:52 +0200 Subject: [PATCH 1/6] solvers: split the scalar Newton into an algorithm and two interfaces backward_euler was two solvers in one class, selected by whether a "function" parameter happened to be non-empty: function set graph-driven: wires residual/jacobian inputs, engine calls update(), result leaves through the "delta" property function empty callback-driven: no inputs, caller passes a lambda to solve(eval), result is returned Almost every oddity in that file traced back to the split. The inputs existed in one mode only, so they were raw pointers with a null guard rather than references -- input_property already tracks wiring through is_wired(), but an input that is never created cannot be checked. update() had no way to report failure, so m_converged was never set on the graph path at all: converging, stalling on a singular jacobian, and exhausting the iteration budget were indistinguishable from outside. And m_converged started true, so a solver that had never run reported success. Worst of it: "function" defaulted to empty, so OMITTING it silently chose callback mode. No input was registered, so wire_inputs() had nothing to validate; update() was never bound; "delta" stayed zero. A graph-driven consumer then read zero forever. Measured with autocatalytic_reaction: a cure that reaches 1.000000 with the parameter set sat at 0.010000 -- its start value -- for 30 steps without an error anywhere. Now three pieces: newton_scalar the algorithm. No properties, no graph, no material_base. Returns {x, converged, iterations}. backward_euler the graph-driven material. "function" is REQUIRED, so the silent case cannot be expressed. Reports convergence. local_newton the material-driven one. No inputs, no function; exposes solve(eval). Referenced by the return maps, which cannot be graph-driven: they solve twice per update with different residuals and pick the branch on the first solve's convergence, which one edge carrying one number cannot do. Convergence now travels WITH the result instead of being queried from the solver afterwards. Drucker-Prager relied on that side channel between its smooth and apex solves, where it went stale by construction. The clamps move to the callers that own them. newton_scalar does not clamp at all: std::max(x, 0) is a statement about a plastic multiplier and abs(x) about a curing degree, neither about Newton's method, and a general solver enforcing one silently is wrong for every other caller. That is #13, resolved by relocation rather than by argument. Verified: J2 and Drucker-Prager reproduce their pre-refactor trajectories to ~1 ULP, unchanged from before this commit. The graph-driven path keeps its own tests. Three new tests cover the setup that used to be silent, a function material lacking residual/jacobian, and local_newton reporting a root it cannot find. --- .../materials/drucker_prager_plasticity.h | 38 +++-- .../materials/j2_plasticity.h | 12 +- .../numsim-materials/solvers/backward_euler.h | 130 ++++++++---------- .../numsim-materials/solvers/local_newton.h | 69 ++++++++++ .../numsim-materials/solvers/newton_scalar.h | 73 ++++++++++ tests/plot_data.cpp | 6 +- tests/test_drucker_prager.cpp | 14 +- tests/test_j2_plasticity.cpp | 6 +- tests/test_materials.cpp | 65 +++++++++ 9 files changed, 309 insertions(+), 104 deletions(-) create mode 100644 include/numsim-materials/solvers/local_newton.h create mode 100644 include/numsim-materials/solvers/newton_scalar.h diff --git a/include/numsim-materials/materials/drucker_prager_plasticity.h b/include/numsim-materials/materials/drucker_prager_plasticity.h index 67aed65..bdc87dd 100644 --- a/include/numsim-materials/materials/drucker_prager_plasticity.h +++ b/include/numsim-materials/materials/drucker_prager_plasticity.h @@ -3,6 +3,7 @@ #include #include +#include #include #include #include @@ -10,7 +11,7 @@ #include "numsim-materials/core/material_ref.h" #include "numsim-materials/materials/drucker_prager_yield_function.h" #include "numsim-materials/materials/plasticity_utils.h" -#include "numsim-materials/solvers/backward_euler.h" +#include "numsim-materials/solvers/local_newton.h" namespace numsim::materials { @@ -29,7 +30,7 @@ class drucker_prager_plasticity final using tensor2 = tmech::tensor; using tensor4 = tmech::tensor; using yield_fn = drucker_prager_yield_function; - using solver_type = backward_euler; + using solver_type = local_newton; template explicit drucker_prager_plasticity(Args&&... args) @@ -106,18 +107,20 @@ class drucker_prager_plasticity final return; } - const auto dlambda = solve_smooth_newton(ts.eval.modified_sig_eq, kappa_n); + const auto smooth = solve_smooth_newton(ts.eval.modified_sig_eq, kappa_n); - // If smooth Newton fails and apex is available, try apex as fallback. - if (!m_solver.get().converged()) { - do_apex_return(ts.eval.sig, C_e, kappa_n); - if (!m_solver.get().converged()) + // If the smooth return fails, the apex is the remaining branch. + if (!smooth.converged) { + if (!do_apex_return(ts.eval.sig, C_e, kappa_n)) throw std::runtime_error( "drucker_prager_plasticity: both smooth and apex Newton failed"); return; } - do_smooth_return(ts.eval, C_e, kappa_n, dlambda); + // dlambda >= 0 is a statement about the plastic multiplier, enforced here + // rather than inside a general scalar solver (see #13). + do_smooth_return(ts.eval, C_e, kappa_n, + std::max(smooth.x, value_type{0})); } private: @@ -131,8 +134,12 @@ class drucker_prager_plasticity final /// Scalar Newton solve: r(Δλ) = phi - G_eff·Δλ - Y0 - H(κ_n + Δλ) = 0. /// Used for both the smooth and apex returns with different (phi, G_eff). - value_type solve_scalar_return(value_type phi_trial, value_type G_eff, - value_type kappa_n) { + /// Returns the solver's result, not a bare number: convergence travels WITH + /// the value instead of being queried from the solver afterwards, where it + /// went stale between the smooth and apex solves. + typename solver_type::result solve_scalar_return(value_type phi_trial, + value_type G_eff, + value_type kappa_n) { auto eval = [&](value_type dl) -> std::pair { m_kappa.new_value() = kappa_n + dl; m_H.update_source(); @@ -143,7 +150,8 @@ class drucker_prager_plasticity final } /// Smooth-cone return Newton: phi = modified_sig_eq, G_eff from yield function. - value_type solve_smooth_newton(value_type phi_trial, value_type kappa_n) { + typename solver_type::result solve_smooth_newton(value_type phi_trial, + value_type kappa_n) { return solve_scalar_return(phi_trial, m_yf.effective_modulus(m_G), kappa_n); } @@ -172,11 +180,14 @@ class drucker_prager_plasticity final /// Apex return: deviatoric stress vanishes, only volumetric Newton. /// dev(ε_p) = dev(ε), tr(ε_p) += β·Δκ. Tangent is rank-1 volumetric. - void do_apex_return(const tensor2& sig_trial, const tensor4& C_e, + /// @return whether the apex Newton converged. + bool do_apex_return(const tensor2& sig_trial, const tensor4& C_e, value_type kappa_n) { const auto phi_apex = m_yf.apex_modified_sig_eq(sig_trial); const auto G_eff_apex = m_yf.apex_effective_modulus(); - const auto dkappa = solve_scalar_return(phi_apex, G_eff_apex, kappa_n); + const auto sol = solve_scalar_return(phi_apex, G_eff_apex, kappa_n); + if (!sol.converged) return false; + const auto dkappa = std::max(sol.x, value_type{0}); m_eps_p.new_value() = m_yf.apex_plastic_strain( m_strain.get(), m_eps_p.old_value(), dkappa); @@ -187,6 +198,7 @@ class drucker_prager_plasticity final // remain on the active apex branch (q stays at 0). m_H.update_source(); m_tangent = m_yf.apex_tangent(m_dH.get()); + return true; } private: diff --git a/include/numsim-materials/materials/j2_plasticity.h b/include/numsim-materials/materials/j2_plasticity.h index dee151c..359d8a1 100644 --- a/include/numsim-materials/materials/j2_plasticity.h +++ b/include/numsim-materials/materials/j2_plasticity.h @@ -2,13 +2,14 @@ #define NUMSIM_MATERIALS_J2_PLASTICITY_H #include +#include #include #include #include #include "numsim-materials/core/material_base.h" #include "numsim-materials/core/material_ref.h" #include "numsim-materials/materials/plasticity_utils.h" -#include "numsim-materials/solvers/backward_euler.h" +#include "numsim-materials/solvers/local_newton.h" namespace numsim::materials { @@ -59,7 +60,7 @@ class j2_plasticity final using base::Dim; using tensor2 = tmech::tensor; using tensor4 = tmech::tensor; - using solver_type = backward_euler; + using solver_type = local_newton; template explicit j2_plasticity(Args&&... args) @@ -150,10 +151,13 @@ class j2_plasticity final return {sig_eq - G_eff * dl - m_sigma_0 - m_H.get(), -G_eff - m_dH.get()}; }; - const auto dlambda = m_solver.get().solve(eval); - if (!m_solver.get().converged()) + const auto sol = m_solver.get().solve(eval); + if (!sol.converged) throw std::runtime_error( "j2_plasticity: return-mapping Newton failed to converge"); + // dlambda >= 0 is a statement about the plastic multiplier, so it is + // enforced here rather than inside a general scalar solver (see #13). + const auto dlambda = std::max(sol.x, value_type{0}); m_eps_p.new_value() = m_eps_p.old_value() + dlambda * N; m_kappa.new_value() = kappa_n + dlambda; diff --git a/include/numsim-materials/solvers/backward_euler.h b/include/numsim-materials/solvers/backward_euler.h index 75dff6b..9bf11ad 100644 --- a/include/numsim-materials/solvers/backward_euler.h +++ b/include/numsim-materials/solvers/backward_euler.h @@ -2,20 +2,32 @@ #define NUMSIM_MATERIALS_BACKWARD_EULER_H #include +#include #include "numsim-materials/core/material_base.h" +#include "numsim-materials/solvers/newton_scalar.h" namespace numsim::materials { -/// Backward Euler solver as a material. +/// Backward-Euler update solved by a scalar Newton the PROPERTY GRAPH drives. /// -/// Consumes: function_name::residual, function_name::jacobian (Local edges) -/// Produces: "delta" (converged increment) +/// Reads "residual" and "jacobian" from the material named by "function", and +/// publishes the increment as "delta". The function material re-evaluates its +/// residual against the current delta through update_source(), which is the +/// circularity that makes this work at all. /// -/// The Newton iteration calls update_source() on the residual/jacobian inputs -/// to re-evaluate them at each trial increment. +/// "function" is REQUIRED. It used to default to an empty string, which +/// silently selected a second, callback-driven mode inside this same class: no +/// inputs were created, update() was never bound, and "delta" stayed at zero. +/// A consumer reading it -- autocatalytic_reaction does -- then froze at its +/// start value for the whole analysis with no error anywhere. Measured: a cure +/// that should reach 1.0 sat at 0.01 for 30 steps. That mode is now +/// local_newton, chosen by naming a different type rather than by omitting a +/// parameter. +/// +/// Parameters: +/// "name", "function", "tolerance", "max_iter" template -class backward_euler final - : public material_base, Traits> { +class backward_euler final : public material_base, Traits> { public: using base = material_base, Traits>; using value_type = typename base::value_type; @@ -26,89 +38,59 @@ class backward_euler final : base(std::forward(args)...), m_delta(base::template add_output("delta")), m_func_name(base::template get_parameter("function")), - m_tol(base::template get_parameter("tolerance")), - m_max_iter(base::template get_parameter("max_iter")) + m_solver(base::template get_parameter("tolerance"), + base::template get_parameter("max_iter")), + m_residual(base::template add_input( + m_func_name, "residual", EdgeKind::Local)), + m_jacobian(base::template add_input( + m_func_name, "jacobian", EdgeKind::Local)) { - // If a function name is provided, set up graph-driven iteration - if (!m_func_name.empty()) { - m_residual = &base::template add_input( - m_func_name, "residual", EdgeKind::Local); - m_jacobian = &base::template add_input( - m_func_name, "jacobian", EdgeKind::Local); - // Bind update callback for graph-driven mode - if (auto p = base::m_property_handler.find(base::m_name, "delta")) - (*p)->traits().update = [this]() { this->update(); }; - } + if (auto p = base::m_property_handler.find(base::m_name, "delta")) + (*p)->traits().update = [this]() { this->update(); }; } static input_parameter_controller parameters() { input_parameter_controller para{base::parameters()}; - para.template insert("function") - .template add(std::string{}); + para.template insert("function").template add(); para.template insert("tolerance") - .template add(value_type{5e-12}); - para.template insert("max_iter") - .template add(100); + .template add(value_type{1e-10}); + para.template insert("max_iter").template add(50); return para; } - /// Property-graph driven iteration: reads residual/jacobian via update_source. - /// Used when the solver is in the graph (e.g., curing reaction). + /// Whether the last update() converged. + /// + /// Previously never set on this path at all: the loop broke on tolerance, on + /// a singular jacobian and on exhausting its budget, and all three looked + /// identical from outside. + [[nodiscard]] bool converged() const noexcept { return m_converged; } + void update() override { - if (!m_residual || !m_jacobian) return; - m_delta = value_type{5e-12}; - for (int i = 0; i < m_max_iter; ++i) { - m_residual->update_source(); - const auto& r = m_residual->get(); - if (std::abs(r) <= m_tol) break; - m_jacobian->update_source(); - const auto& j = m_jacobian->get(); - if (std::abs(j) < value_type{1e-30}) break; - auto step = r / j; - // Damped Newton: halve step if it produces NaN - for (int k = 0; k < 5; ++k) { - auto candidate = m_delta - step; - m_delta = candidate; - m_residual->update_source(); - auto r_new = m_residual->get(); - if (!std::isnan(r_new) && !std::isinf(r_new)) break; - m_delta = candidate + step; // restore - step *= value_type{0.5}; - } - } - // Ensure positive increment (curing degree can only increase) - m_delta = std::abs(m_delta); - } + // The seed is nonzero because a jacobian evaluated at exactly zero is + // singular for the rate laws this drives. + const auto r = m_solver.solve( + [this](value_type x) -> std::pair { + m_delta = x; + m_residual.update_source(); + const auto res = m_residual.get(); + m_jacobian.update_source(); + return {res, m_jacobian.get()}; + }, + value_type{5e-12}); - /// Direct call: another material provides eval(x) → {residual, jacobian}. - /// Used when the caller drives the iteration (e.g., plasticity return mapping). - /// Sets m_converged to indicate whether the iteration converged. - /// The returned value is clamped to be non-negative — for plasticity, a - /// negative plastic-multiplier increment is unphysical (backward plastic flow). - template - value_type solve(Eval&& eval, value_type x0 = value_type{0}) { - auto x = x0; - for (int i = 0; i < m_max_iter; ++i) { - auto [r, dr] = eval(x); - if (std::abs(r) < m_tol) { m_converged = true; return std::max(x, value_type{0}); } - if (std::abs(dr) < value_type{1e-30}) { m_converged = false; return std::max(x, value_type{0}); } - x -= r / dr; - } - m_converged = false; - return std::max(x, value_type{0}); + m_converged = r.converged; + // Sign convention owned here rather than by the solver: the quantities this + // integrates -- degree of cure, and similar -- only increase. + m_delta = std::abs(r.x); } - /// Whether the last solve() call converged. - bool converged() const { return m_converged; } - private: value_type& m_delta; const std::string& m_func_name; - const value_type& m_tol; - const int& m_max_iter; - const input_property* m_residual{nullptr}; - const input_property* m_jacobian{nullptr}; - bool m_converged{true}; + newton_scalar m_solver; + const input_property& m_residual; + const input_property& m_jacobian; + bool m_converged{false}; }; } // namespace numsim::materials diff --git a/include/numsim-materials/solvers/local_newton.h b/include/numsim-materials/solvers/local_newton.h new file mode 100644 index 0000000..9b9bd48 --- /dev/null +++ b/include/numsim-materials/solvers/local_newton.h @@ -0,0 +1,69 @@ +#ifndef NUMSIM_MATERIALS_LOCAL_NEWTON_H +#define NUMSIM_MATERIALS_LOCAL_NEWTON_H + +#include +#include "numsim-materials/core/material_base.h" +#include "numsim-materials/solvers/newton_scalar.h" + +namespace numsim::materials { + +/// A scalar Newton another MATERIAL drives, by calling solve() with its own +/// residual. +/// +/// The counterpart to backward_euler, which the property graph drives. The two +/// used to be one class selected by whether a "function" parameter happened to +/// be set, which meant an empty "function" silently chose this behaviour -- +/// leaving a graph-driven consumer reading a delta of 0 forever, with no error +/// anywhere. Splitting them makes the choice a type rather than a defaulted +/// string. +/// +/// Publishes no properties and takes no inputs. It is a material only so its +/// tolerance and iteration budget are configurable from the same document as +/// everything else, and so one instance can serve several materials. +/// +/// Used by the return maps, which cannot be graph-driven: they solve twice per +/// update with different residuals (smooth cone, then apex) and choose the +/// branch on the first solve's convergence. A graph edge carries one number, +/// not that. +/// +/// Parameters: +/// "name", "tolerance", "max_iter" +template +class local_newton final : public material_base, Traits> { +public: + using base = material_base, Traits>; + using value_type = typename base::value_type; + using input_parameter_controller = typename base::input_parameter_controller; + using result = typename newton_scalar::result; + + template + explicit local_newton(Args&&... args) + : base(std::forward(args)...), + m_solver(base::template get_parameter("tolerance"), + base::template get_parameter("max_iter")) {} + + static input_parameter_controller parameters() { + input_parameter_controller para{base::parameters()}; + para.template insert("tolerance") + .template add(value_type{1e-10}); + para.template insert("max_iter").template add(50); + return para; + } + + /// Solve with the caller's residual. @p eval maps x -> {residual, jacobian}. + /// + /// The result carries its own convergence flag, so it cannot be read without + /// being available -- unlike a converged() queried separately, which a caller + /// can forget and which goes stale between solves. + template + result solve(Eval&& eval, value_type x0 = value_type{}) const { + return m_solver.solve(std::forward(eval), x0); + } + +private: + newton_scalar m_solver; +}; + +} // namespace numsim::materials + +#endif // NUMSIM_MATERIALS_LOCAL_NEWTON_H diff --git a/include/numsim-materials/solvers/newton_scalar.h b/include/numsim-materials/solvers/newton_scalar.h new file mode 100644 index 0000000..8e97a95 --- /dev/null +++ b/include/numsim-materials/solvers/newton_scalar.h @@ -0,0 +1,73 @@ +#ifndef NUMSIM_MATERIALS_NEWTON_SCALAR_H +#define NUMSIM_MATERIALS_NEWTON_SCALAR_H + +#include +#include + +namespace numsim::materials { + +/// Scalar Newton iteration. Plain algorithm: no properties, no graph presence, +/// no material_base. +/// +/// Split out of backward_euler, which had grown two Newton loops -- one driven +/// by graph properties, one by a caller's lambda -- with separate damping and +/// separate clamping. The algorithm is the same in both cases; only where the +/// residual comes from differs, and that belongs to the caller. +/// +/// Deliberately does NOT clamp its result. backward_euler forced x >= 0 for +/// plasticity's dlambda and |x| for a curing degree; both are statements about +/// a particular unknown, not about Newton's method, and a general solver that +/// silently enforces one is wrong for every other caller (see #13). Callers +/// clamp what they own. +template +class newton_scalar { +public: + struct result { + T x{}; ///< the iterate, converged or not + bool converged{false}; + int iterations{0}; + }; + + newton_scalar(T tolerance, int max_iterations) noexcept + : m_tol(tolerance), m_max_iter(max_iterations) {} + + /// @param eval x -> {residual, jacobian} + /// @param x0 initial iterate + /// + /// Reports convergence rather than returning a bare number: a caller that + /// does not check gets a non-converged iterate, which is the same failure + /// backward_euler's graph path had -- there it was unreportable, here it is + /// merely unchecked. + template + result solve(Eval&& eval, T x0 = T{}) const { + result r{x0, false, 0}; + for (int i = 0; i < m_max_iter; ++i) { + const auto [residual, jacobian] = eval(r.x); + ++r.iterations; + if (std::abs(residual) < m_tol) { + r.converged = true; + return r; + } + // A vanishing jacobian is a stall, not a convergence. + if (std::abs(jacobian) < T{1e-30}) return r; + r.x -= residual / jacobian; + } + // Exhausted the budget: the last step may still have landed on the root, + // so the residual is checked once more rather than assumed bad. + const auto [residual, jacobian] = eval(r.x); + (void)jacobian; + r.converged = std::abs(residual) < m_tol; + return r; + } + + [[nodiscard]] T tolerance() const noexcept { return m_tol; } + [[nodiscard]] int max_iterations() const noexcept { return m_max_iter; } + +private: + T m_tol; + int m_max_iter; +}; + +} // namespace numsim::materials + +#endif // NUMSIM_MATERIALS_NEWTON_SCALAR_H diff --git a/tests/plot_data.cpp b/tests/plot_data.cpp index b32db4b..9b9984e 100644 --- a/tests/plot_data.cpp +++ b/tests/plot_data.cpp @@ -8,7 +8,7 @@ #include "numsim-materials/materials/drucker_prager_yield_function.h" #include "numsim-materials/materials/drucker_prager_plasticity.h" #include "numsim-materials/materials/j2_plasticity.h" -#include "numsim-materials/solvers/backward_euler.h" +#include "numsim-materials/solvers/local_newton.h" #include "numsim-materials/postprocessing/numerical_diff_checker.h" using policy = numsim::materials::material_policy_default; @@ -68,7 +68,7 @@ run_result run_j2(T increment, int steps, p.clear(); p.insert("name", "solver"); - ctx.create>(p); + ctx.create>(p); p.clear(); p.insert("name", "hardening"); @@ -133,7 +133,7 @@ run_result run_dp(T increment, int steps, p.clear(); p.insert("name", "solver"); - ctx.create>(p); + ctx.create>(p); p.clear(); p.insert("name", "hardening"); diff --git a/tests/test_drucker_prager.cpp b/tests/test_drucker_prager.cpp index c495c32..af7269e 100644 --- a/tests/test_drucker_prager.cpp +++ b/tests/test_drucker_prager.cpp @@ -10,7 +10,7 @@ #include "numsim-materials/materials/drucker_prager_yield_function.h" #include "numsim-materials/materials/drucker_prager_plasticity.h" #include "numsim-materials/materials/rk_plasticity.h" -#include "numsim-materials/solvers/backward_euler.h" +#include "numsim-materials/solvers/local_newton.h" #include "numsim-materials/solvers/butcher_tableau.h" #include "numsim-materials/postprocessing/numerical_diff_checker.h" @@ -104,7 +104,7 @@ class DruckerPragerTest : public ::testing::Test { p.clear(); p.insert("name", "solver"); - ctx.create>(p); + ctx.create>(p); p.clear(); p.insert("name", "hardening"); @@ -209,7 +209,7 @@ class DPTangentTest : public ::testing::Test { p.clear(); p.insert("name", "solver"); - ctx.create>(p); + ctx.create>(p); p.clear(); p.insert("name", "hardening"); @@ -282,7 +282,7 @@ T run_dp_max_tangent_error(T increment, int steps) { p.clear(); p.insert("name", "solver"); - ctx.create>(p); + ctx.create>(p); p.clear(); p.insert("name", "hardening"); @@ -370,7 +370,7 @@ T max_tangent_error(std::vector direction, T increment, int steps) { p.clear(); p.insert("name", "solver"); - ctx.create>(p); + ctx.create>(p); p.clear(); p.insert("name", "hardening"); @@ -467,7 +467,7 @@ TEST(DruckerPragerApex, HydrostaticTensionReachesTheApex) { p.clear(); p.insert("name", "solver"); - ctx.create>(p); + ctx.create>(p); p.clear(); p.insert("name", "hardening"); @@ -527,7 +527,7 @@ TEST(DruckerPragerApex, ApexStateIsAdmissible) { ctx.create>(p); p.clear(); p.insert("name", "solver"); - ctx.create>(p); + ctx.create>(p); p.clear(); p.insert("name", "hardening"); p.insert("source", "dp"); diff --git a/tests/test_j2_plasticity.cpp b/tests/test_j2_plasticity.cpp index a7b3429..71f7ea0 100644 --- a/tests/test_j2_plasticity.cpp +++ b/tests/test_j2_plasticity.cpp @@ -9,7 +9,7 @@ #include "numsim-materials/materials/drucker_prager_plasticity.h" #include "numsim-materials/materials/j2_plasticity.h" #include "numsim-materials/materials/rk_plasticity.h" -#include "numsim-materials/solvers/backward_euler.h" +#include "numsim-materials/solvers/local_newton.h" #include "numsim-materials/solvers/butcher_tableau.h" #include "numsim-materials/postprocessing/numerical_diff_checker.h" @@ -44,7 +44,7 @@ class J2PlasticityTest : public ::testing::Test { // Newton-Raphson solver p.clear(); p.insert("name", "solver"); - ctx.create>(p); + ctx.create>(p); // Linear isotropic hardening (Local edge — called in inner loop) p.clear(); @@ -158,7 +158,7 @@ class J2TangentTest : public ::testing::Test { p.clear(); p.insert("name", "solver"); - ctx.create>(p); + ctx.create>(p); p.clear(); p.insert("name", "hardening"); diff --git a/tests/test_materials.cpp b/tests/test_materials.cpp index 44e09cb..74f1c2c 100644 --- a/tests/test_materials.cpp +++ b/tests/test_materials.cpp @@ -4,6 +4,7 @@ #include "numsim-materials/materials/tensor_component_stepper.h" #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/scalar_identity_weight.h" #include "numsim-materials/materials/autocatalytic_reaction.h" #include "numsim-materials/solvers/backward_euler.h" @@ -159,3 +160,67 @@ TEST(CuringSimulation, ConvergesToFullCure) { } } // namespace + +namespace { +namespace nm_be = numsim::materials; + +/// backward_euler is the GRAPH-driven solver, so "function" is required. +/// +/// It used to default to empty, which silently selected a second, +/// callback-driven mode inside the same class: no inputs were created, update() +/// was never bound, and "delta" stayed at zero. A consumer reading it froze at +/// its start value for the whole analysis with no error -- measured, a cure +/// that should reach 1.0 sat at 0.01 for 30 steps. That mode is local_newton +/// now, chosen by naming a type rather than by omitting a parameter. +TEST(BackwardEulerSetup, RequiresAFunctionToSolve) { + using policy = nm_be::material_policy_default; + nm_be::material_context ctx; + policy::ParameterHandler p; + p.insert("name", "solver"); + EXPECT_THROW(ctx.create>(p), + std::invalid_argument) + << "omitting \"function\" must fail at setup, not produce an inert " + "solver whose consumers silently never advance"; +} + +/// A function material that does not publish residual/jacobian is caught at +/// finalize, by name. +TEST(BackwardEulerSetup, RejectsAFunctionWithoutResidualOrJacobian) { + using policy = nm_be::material_policy_default; + using T2 = policy::value_type; + nm_be::material_context ctx; + policy::ParameterHandler p; + p.insert("name", "drv"); + p.insert("increment", T2{0.1}); + ctx.create>(p); + p.clear(); + p.insert("name", "solver"); + p.insert("function", std::string("drv")); + ctx.create>(p); + EXPECT_THROW(ctx.finalize(), std::runtime_error); +} + +/// local_newton is the material-driven one: no function, no graph inputs, and +/// its result carries convergence so a caller cannot read the number without +/// the flag being at hand. +TEST(LocalNewtonSolver, SolvesAndReportsConvergence) { + using policy = nm_be::material_policy_default; + using T2 = policy::value_type; + nm_be::material_context ctx; + policy::ParameterHandler p; + p.insert("name", "solver"); + auto& s = ctx.create>(p); + ctx.finalize(); + + // x^2 - 4 = 0 from x0 = 3 + const auto ok = s.solve( + [](T2 x) { return std::pair{x * x - T2{4}, T2{2} * x}; }, T2{3}); + EXPECT_TRUE(ok.converged); + EXPECT_NEAR(ok.x, 2.0, 1e-10); + + // No root: x^2 + 1 = 0. Must report failure rather than a plausible number. + const auto bad = s.solve( + [](T2 x) { return std::pair{x * x + T2{1}, T2{2} * x}; }, T2{1}); + EXPECT_FALSE(bad.converged); +} +} // namespace From d952cf38ad528a8d6c7ead45b43c66d13ae097ac Mon Sep 17 00:00:00 2001 From: petlenz Date: Sat, 5 Sep 2026 12:05:30 +0200 Subject: [PATCH 2/6] docs: describe the plasticity materials after the split How the three plasticity models are wired now that small_strain_plasticity is gone, what each requires, and which assumptions carry weight. Written to record the reasoning that is not visible in the code: - isotropic C_e is REQUIRED, not preferred, which is why each material builds its own rather than reading a rank-4 tangent from an elastic material - linear_elasticity's stress is wrong inside a plasticity graph, with the measured divergence, and why the material is still valid elsewhere - the return map cannot be graph-driven, which is why local_newton exists alongside backward_euler rather than replacing it - the apex tangent is a BRANCH tangent, valid only on the branch Includes the two coverage lessons this work turned up: a single load path proves one path (the apex return was executed by no test at all, because it is unreachable from a uniaxial path), and a tolerance set by the worst step licenses errors in every other one. Performance quoted as paired interleaved medians with ranges, since the same binary varies 1171-2253 ns on this machine and single-run comparisons are worth about one significant digit. --- docs/plasticity.md | 286 +++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 286 insertions(+) create mode 100644 docs/plasticity.md diff --git a/docs/plasticity.md b/docs/plasticity.md new file mode 100644 index 0000000..a47e714 --- /dev/null +++ b/docs/plasticity.md @@ -0,0 +1,286 @@ +# Plasticity materials + +How the plasticity models are put together after the split of +`small_strain_plasticity`, what each one requires, and which assumptions are +load-bearing. + +Current as of `refactor/scalar-newton-split` (PRs #39 and #40). + +--- + +## 1. The layout + +``` + ┌──────────────────────┐ + │ newton_scalar │ plain algorithm + │ no graph, no props │ returns {x, converged, iterations} + └──────────┬───────────┘ + used by │ + ┌────────────────┴────────────────┐ + │ │ + ┌──────────▼──────────┐ ┌───────────▼───────────┐ + │ backward_euler │ │ local_newton │ + │ MATERIAL, graph- │ │ MATERIAL, driven by │ + │ driven │ │ another material │ + │ "function" REQUIRED│ │ exposes solve(eval) │ + └──────────┬──────────┘ └───────────┬───────────┘ + │ │ material_ref + ┌──────────▼──────────┐ ┌───────────▼───────────┐ + │ autocatalytic_ │ │ j2_plasticity │ + │ reaction, curing │ │ drucker_prager_ │ + │ │ │ plasticity │ + └─────────────────────┘ └───────────────────────┘ + + rk_plasticity iterates its own Butcher tableau and uses neither. +``` + +Three plasticity materials, none of them templated on a yield function any +more: + +| material | integrator | yield surface | flow | +|---|---|---|---| +| `j2_plasticity` | backward Euler (radial return) | von Mises cylinder | associative | +| `drucker_prager_plasticity` | backward Euler + apex branch | DP cone | **non**-associative (β ≠ η) | +| `rk_plasticity` | Runge–Kutta, any Butcher tableau | J2 | associative | + +--- + +## 2. Interfaces + +### `j2_plasticity` + +| parameter | meaning | +|---|---| +| `hardening_source` | material publishing `hardening_stress` and `hardening_modulus` | +| `strain_source` | material publishing `strain` | +| `solver_source` | a `local_newton` | +| `K`, `G` | bulk and shear moduli — the elastic stiffness is built from these | +| `sigma_0` | initial yield stress | + +Outputs: `stress`, `tangent`, and the history pair `plastic_strain`, +`equivalent_plastic_strain`. + +### `drucker_prager_plasticity` + +Same, minus `K`, plus the cone: + +| parameter | meaning | +|---|---| +| `eta` | friction (pressure coefficient) | +| `beta` | dilatancy — `beta != eta` is what makes the flow non-associative | +| `K_bulk` | bulk modulus, used both for the cone apex and for the elastic stiffness | + +These used to arrive inside a `yield_function` **C++ object**, which the JSON +reader cannot convert — so the material could not be configured from a document +at all. As plain scalars it can (issue #33). + +### `rk_plasticity` + +Takes `K`, `G`, `sigma_0`, the two sources, plus `tolerance`, `max_iter` and a +`tableau` pointer. Available tableaus: `forward_euler`, `explicit_midpoint`, +`rk4`, `implicit_euler`, `implicit_midpoint`, `crank_nicolson`, `sdirk3`, +`gauss_legendre_4`. + +### Hardening + +Both publish `hardening_stress` (H) and `hardening_modulus` (dH/dκ), and read +`equivalent_plastic_strain` from the plasticity material — a `Local` edge, +because H depends on κ, which is the unknown being solved for. + +| material | law | parameters | +|---|---|---| +| `linear_isotropic_hardening` | `H = K κ` | `source`, `K` | +| `exponential_isotropic_hardening` | `H = K_inf (1 − e^{−δκ})` | `source`, `K_inf`, `delta` | + +--- + +## 3. What is assumed, and where it bites + +### Isotropic elasticity is required, not preferred + +Every closed form below rests on `C_e` being isotropic: + +``` +C_e : N = 2G dev(N) + K tr(N) I +N : C_e : N = 3G (J2, N deviatoric) +C_e : X : C_e = 4G² X (X deviatoric in both index pairs) +``` + +This is why the plasticity materials **build their own** `C_e` from `K` and `G` +rather than reading a rank-4 tangent from an elastic material. Accepting an +arbitrary `C_e` advertised a generality none of them can honour — a caller +wiring an anisotropic tangent would have got silently wrong answers. + +It also keeps `linear_elasticity` out of plasticity graphs, which matters for a +second reason: its `stress` output is `C : ε` with `ε_p` ignored, so inside a +plasticity graph it is not the stress of anything and diverges as plastic strain +accumulates. + +| step | κ | `stress` (real) | `elastic::stress` | error | +|---|---|---|---|---| +| 40 | 0.0094 | 106.25 | 107.69 | +1.4% | +| 120 | 0.1094 | 306.25 | 323.08 | +5.5% | +| 200 | 0.2094 | 506.25 | 538.46 | **+6.4%** | + +`linear_elasticity` remains a valid material — for a genuinely elastic model +`C : ε` *is* the answer, and `isotropic_damage` consumes both its stress and its +tangent legitimately. It is simply the wrong dependency for plasticity. + +### Clamps belong to the caller + +`newton_scalar` does not clamp its result. `max(x, 0)` is a statement about a +plastic multiplier and `abs(x)` about a curing degree; neither is about Newton's +method. Each lives in the material that owns the assumption (issue #13). + +### The return map cannot be graph-driven + +`local_newton` exists rather than everything using `backward_euler` because a +return map solves **twice per update** — smooth cone, then apex — with different +residuals, and chooses the branch on the first solve's convergence. A property +edge carries one number; it cannot carry "and it failed, so take the other +branch". + +Convergence therefore travels *with* the result (`{x, converged, iterations}`) +rather than being queried from the solver afterwards, where it went stale +between the two solves. + +--- + +## 4. The consistent tangent + +For the smooth return, with `M` the yield normal and `N` the flow normal +(`M == N` when associative): + +``` +A = C_e − Δλ (C_e : dN/dσ : C_e) +dλ/dε = (M : C_e) / (M : C_e : N + H′) +C_consistent = A − (C_e : N) ⊗ dλ/dε +``` + +`M ≠ N` under non-associative flow makes the tangent **major-asymmetric**, +measured at 0.6–0.8 % of peak magnitude on multiaxial DP paths. + +J2 collapses this to the standard closed form: + +``` +C = C_e − (6G²Δλ/σ_eq) IIdev + (4G²Δλ/σ_eq − 4G²/(3G+H′)) N ⊗ N +``` + +which is algebraically identical, not an approximation — confirmed by +bit-identical stress against the general path. + +### Apex return (Drucker-Prager only) + +When the deviatoric correction would overshoot the cone tip +(`G Δλ ≥ √J₂`), the return goes to the apex instead: deviatoric stress to zero, +pressure pinned at `(k + H)/η`, volumetric plastic flow from `β`. Its tangent is +a **branch** tangent — valid only for perturbations staying on the apex — because +the return map is non-smooth there. + +--- + +## 5. Verification + +`tangent_checker` compares the analytical tangent against a central difference +through the graph. Coverage: + +| path | purpose | +|---|---| +| uniaxial | the historical case | +| pure shear, biaxial, triaxial, mixed dev+shear | exercise genuine non-associativity | +| hydrostatic | the **only** path that reaches the apex branch | + +All agree to ~1e-10, the apex to 8e-9. + +Two coverage lessons worth keeping: + +- **A single load path proves one path.** The apex return was executed by *no* + test for its entire existence — instrumenting the predicate gave + `APEX_HITS=0` across every binary — because every path was uniaxial and the + apex sits on the hydrostatic axis. It was unreachable by construction, not by + oversight. +- **A tolerance set by the worst step licenses errors in all the others.** + `J2TangentTest` bounded a whole run at `0.1` because the elastic→plastic + transition step is genuinely inexact. A 0.5 % error injected into the tangent + landed at 2.4e-4 and passed, while real plastic steps sit at 4.4e-10. The two + regimes are now bounded separately. + +--- + +## 6. Performance + +Two binaries built from the actual commits, run interleaved so each pair sees +the same machine state; 15 pairs; speedup computed per pair. + +| | old (median) | new (median) | paired speedup | range | +|---|---|---|---|---| +| J2 | 1466.2 ns/step | **211.5** | **7.00×** | 5.54–8.84× | +| Drucker-Prager | 1176.3 ns/step | **652.4** | **1.79×** | 1.58–1.90× | + +The range is the honest figure: the *same* binary measured 1171–2253 ns across +runs on this machine, so any single-run comparison is worth about one +significant digit. + +J2 gains more because it also sheds machinery it never used. DP keeps the +non-associative structure and the apex branch and gains only the arithmetic — +which is the correct outcome. **The generality DP pays for is generality DP +uses.** + +--- + +## 7. Building a model + +```cpp +// solver — one instance can serve several materials +p.insert("name", "solver"); +ctx.create>(p); + +// hardening — reads back from the plasticity material +p.clear(); +p.insert("name", "hardening"); +p.insert("source", "j2"); +p.insert("K", 1000.0); +ctx.create>(p); + +// plasticity — no elastic material anywhere +p.clear(); +p.insert("name", "j2"); +p.insert("hardening_source", "hardening"); +p.insert("strain_source", "stepper"); +p.insert("solver_source", "solver"); +p.insert("K", 166.67); +p.insert("G", 76.92); +p.insert("sigma_0", 50.0); +ctx.create>(p); +``` + +Note `hardening.source = "j2"` while `j2.hardening_source = "hardening"`. The +cycle is deliberate and is why those edges are `Local`: H depends on κ, which is +what the return map solves for, so the hardening material is re-evaluated inside +the Newton loop through `update_source()`. + +The same model in JSON, which is possible for Drucker-Prager only since its cone +parameters became plain scalars: + +```json +{"type": "drucker_prager_plasticity", "name": "dp", + "hardening_source": "hardening", "strain_source": "stepper", + "solver_source": "solver", + "G": 76.92, "sigma_0": 20.0, + "eta": 0.1, "beta": 0.05, "K_bulk": 166.67} +``` + +--- + +## 8. Known gaps + +- **`rk_plasticity` is still templated** on a yield function with one + instantiation (`j2_rk_plasticity`) — the same shape that was removed from + `small_strain_plasticity`. +- **`material_ref`** exists for exactly two call sites (the two backward-Euler + return maps) and costs ~99 lines of core machinery. It bypasses the + topological sort, so the plasticity↔solver ordering is not an edge the engine + knows about. Safe today because `local_newton` holds no per-solve state. +- **The apex tangent is a branch tangent.** Verified consistent *on* the branch; + a perturbation that leaves the apex back onto the smooth cone is not covered, + and cannot be by a central difference. From b67b84081dfc7d54b1a3b7492b71261eb81239f3 Mon Sep 17 00:00:00 2001 From: petlenz Date: Sat, 5 Sep 2026 16:08:02 +0200 Subject: [PATCH 3/6] materials: collapse rk_plasticity into a dedicated j2_rk_plasticity rk_plasticity had exactly one instantiation, the j2_rk_plasticity alias -- the same shape already removed from small_strain_plasticity. A template parameter with one argument is indirection, not generality. The yield function becomes a fixed member (j2_yield_function) as it did for Drucker-Prager, and the file is renamed to match the class. 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 No plasticity class is templated on a yield function now. What remains generic is plasticity_utils -- compute_trial, evaluate_at_state, compute_tangent -- and that generality is genuine: drucker_prager_plasticity instantiates it with the cone, j2_rk_plasticity with the cylinder. Two callers, two yield functions, shared return-mapping algebra. Documented in docs/plasticity.md, since the difference between that and a one-argument template is the whole point. --- docs/plasticity.md | 24 +++++++++++++------ .../{rk_plasticity.h => j2_rk_plasticity.h} | 23 ++++++++---------- tests/test_drucker_prager.cpp | 2 +- tests/test_j2_plasticity.cpp | 2 +- 4 files changed, 29 insertions(+), 22 deletions(-) rename include/numsim-materials/materials/{rk_plasticity.h => j2_rk_plasticity.h} (92%) diff --git a/docs/plasticity.md b/docs/plasticity.md index a47e714..44aece9 100644 --- a/docs/plasticity.md +++ b/docs/plasticity.md @@ -31,7 +31,7 @@ Current as of `refactor/scalar-newton-split` (PRs #39 and #40). │ │ │ plasticity │ └─────────────────────┘ └───────────────────────┘ - rk_plasticity iterates its own Butcher tableau and uses neither. + j2_rk_plasticity iterates its own Butcher tableau and uses neither. ``` Three plasticity materials, none of them templated on a yield function any @@ -41,7 +41,7 @@ more: |---|---|---|---| | `j2_plasticity` | backward Euler (radial return) | von Mises cylinder | associative | | `drucker_prager_plasticity` | backward Euler + apex branch | DP cone | **non**-associative (β ≠ η) | -| `rk_plasticity` | Runge–Kutta, any Butcher tableau | J2 | associative | +| `j2_rk_plasticity` | Runge–Kutta, any Butcher tableau | von Mises cylinder | associative | --- @@ -74,7 +74,7 @@ These used to arrive inside a `yield_function` **C++ object**, which the JSON reader cannot convert — so the material could not be configured from a document at all. As plain scalars it can (issue #33). -### `rk_plasticity` +### `j2_rk_plasticity` Takes `K`, `G`, `sigma_0`, the two sources, plus `tolerance`, `max_iter` and a `tableau` pointer. Available tableaus: `forward_euler`, `explicit_midpoint`, @@ -272,11 +272,21 @@ parameters became plain scalars: --- -## 8. Known gaps +## 8. What is still generic, and why + +`plasticity_utils`' free functions — `compute_trial`, `evaluate_at_state`, +`compute_tangent` — remain templated on a yield function, and that generality is +real: `drucker_prager_plasticity` instantiates them with the DP cone and +`j2_rk_plasticity` with the von Mises cylinder. Two callers, two yield +functions, shared return-mapping algebra. + +That is the distinction worth holding onto. A template parameter with **one** +argument is indirection — it was removed from `small_strain_plasticity` and from +`rk_plasticity`. A template parameter with two genuinely different arguments is +what templates are for. + +## 9. Known gaps -- **`rk_plasticity` is still templated** on a yield function with one - instantiation (`j2_rk_plasticity`) — the same shape that was removed from - `small_strain_plasticity`. - **`material_ref`** exists for exactly two call sites (the two backward-Euler return maps) and costs ~99 lines of core machinery. It bypasses the topological sort, so the plasticity↔solver ordering is not an edge the engine diff --git a/include/numsim-materials/materials/rk_plasticity.h b/include/numsim-materials/materials/j2_rk_plasticity.h similarity index 92% rename from include/numsim-materials/materials/rk_plasticity.h rename to include/numsim-materials/materials/j2_rk_plasticity.h index 05d40bd..0a7c57a 100644 --- a/include/numsim-materials/materials/rk_plasticity.h +++ b/include/numsim-materials/materials/j2_rk_plasticity.h @@ -1,5 +1,5 @@ -#ifndef NUMSIM_MATERIALS_RK_PLASTICITY_H -#define NUMSIM_MATERIALS_RK_PLASTICITY_H +#ifndef NUMSIM_MATERIALS_J2_RK_PLASTICITY_H +#define NUMSIM_MATERIALS_J2_RK_PLASTICITY_H #include #include @@ -17,23 +17,23 @@ namespace numsim::materials { /// Each implicit stage solves F = 0 for Δλ_i. Explicit stages use /// the consistency condition. Shares trial/tangent code with /// small_strain_plasticity via plasticity_utils.h. -template -class rk_plasticity final - : public material_base, Traits> { +template +class j2_rk_plasticity final + : public material_base, Traits> { public: - using base = material_base, Traits>; + using base = material_base, Traits>; using value_type = typename base::value_type; using input_parameter_controller = typename base::input_parameter_controller; static constexpr auto Dim = base::Dim; using tensor2 = tmech::tensor; using tensor4 = tmech::tensor; - using yield_fn = YieldFunction; + using yield_fn = j2_yield_function; template - explicit rk_plasticity(Args&&... args) + explicit j2_rk_plasticity(Args&&... args) : base(std::forward(args)...), m_stress(base::template add_output( - "stress", &rk_plasticity::compute)), + "stress", &j2_rk_plasticity::compute)), 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")), @@ -211,10 +211,7 @@ class rk_plasticity final std::vector m_diag; }; -template -using j2_rk_plasticity = rk_plasticity>; } // namespace numsim::materials -#endif // NUMSIM_MATERIALS_RK_PLASTICITY_H +#endif // NUMSIM_MATERIALS_J2_RK_PLASTICITY_H diff --git a/tests/test_drucker_prager.cpp b/tests/test_drucker_prager.cpp index af7269e..2f03c7d 100644 --- a/tests/test_drucker_prager.cpp +++ b/tests/test_drucker_prager.cpp @@ -9,7 +9,7 @@ #include "numsim-materials/materials/linear_isotropic_hardening.h" #include "numsim-materials/materials/drucker_prager_yield_function.h" #include "numsim-materials/materials/drucker_prager_plasticity.h" -#include "numsim-materials/materials/rk_plasticity.h" +#include "numsim-materials/materials/j2_rk_plasticity.h" #include "numsim-materials/solvers/local_newton.h" #include "numsim-materials/solvers/butcher_tableau.h" #include "numsim-materials/postprocessing/numerical_diff_checker.h" diff --git a/tests/test_j2_plasticity.cpp b/tests/test_j2_plasticity.cpp index 71f7ea0..b24f045 100644 --- a/tests/test_j2_plasticity.cpp +++ b/tests/test_j2_plasticity.cpp @@ -8,7 +8,7 @@ #include "numsim-materials/materials/linear_isotropic_hardening.h" #include "numsim-materials/materials/drucker_prager_plasticity.h" #include "numsim-materials/materials/j2_plasticity.h" -#include "numsim-materials/materials/rk_plasticity.h" +#include "numsim-materials/materials/j2_rk_plasticity.h" #include "numsim-materials/solvers/local_newton.h" #include "numsim-materials/solvers/butcher_tableau.h" #include "numsim-materials/postprocessing/numerical_diff_checker.h" From a37020f5e00ed9231b848ae3aab651a254d1cfd4 Mon Sep 17 00:00:00 2001 From: petlenz Date: Sun, 6 Sep 2026 12:14:44 +0200 Subject: [PATCH 4/6] docs: correct why the return maps are not graph-driven The reason given was that they solve twice per update and branch on convergence. That is true of Drucker-Prager and NOT of J2, which solves once and has no branch -- so the justification did not cover the case it was written for. The real reason is common to both and simpler: a property's update callback runs whenever the graph updates, and plasticity only solves when the trial state exceeds yield. A conditional computation cannot be an unconditional callback. That is not merely tidier, it is what the cost is: elastic step ~92 ns plastic step ~168 ns Graph-driving would run a Newton at every elastic point -- roughly 80% more work, paid where a real analysis spends most of its time. And at dl = 0 on an elastic step the residual is negative, so Newton drives dl negative and needs a clamp to hold it at zero, which is exactly the max(x, 0) that #13 objects to. The clamp and the graph-driving are the same problem. Drucker-Prager's branching is now stated as a SECOND, independent reason rather than the primary one, and J2's position is stated plainly: it could be graph-driven, it would just cost more than it saves. --- docs/plasticity.md | 46 +++++++++++++++++++++++++++++++++++++++------- 1 file changed, 39 insertions(+), 7 deletions(-) diff --git a/docs/plasticity.md b/docs/plasticity.md index 44aece9..08a8d6f 100644 --- a/docs/plasticity.md +++ b/docs/plasticity.md @@ -132,17 +132,49 @@ tangent legitimately. It is simply the wrong dependency for plasticity. plastic multiplier and `abs(x)` about a curing degree; neither is about Newton's method. Each lives in the material that owns the assumption (issue #13). -### The return map cannot be graph-driven +### The return map solves conditionally; a graph property is evaluated unconditionally -`local_newton` exists rather than everything using `backward_euler` because a -return map solves **twice per update** — smooth cone, then apex — with different -residuals, and chooses the branch on the first solve's convergence. A property -edge carries one number; it cannot carry "and it failed, so take the other -branch". +This is why `local_newton` exists rather than everything using +`backward_euler`. + +A property's update callback runs whenever the graph updates. **Plasticity only +solves when the trial state exceeds yield** — and in a real analysis most +integration points are elastic most of the time. Both return maps short-circuit +before touching the solver: + +```cpp +if (sig_eq - sigma_0 - H <= 0) { /* elastic: stress = trial, tangent = C_e */ return; } +``` + +Making the solve a property callback would run a Newton at every elastic point: + +| | measured | +|---|---| +| elastic step | ~92 ns | +| plastic step | ~168 ns | + +Roughly 80 % more work, paid exactly where a real analysis spends most of its +time. + +It is worse than a cost. At `Δλ = 0` on an elastic step the residual is +negative, so Newton drives `Δλ` negative; holding it at zero needs a clamp — +which is precisely the `max(x, 0)` that issue #13 objects to. **The clamp and +the graph-driving are the same problem**: the graph mode has no way to express +*do not solve*. + +Drucker-Prager cannot be graph-driven for a second, independent reason: it +solves up to **twice** per update — an apex pre-check, a smooth-cone solve, then +an apex fallback if that fails to converge — and chooses the branch on the first +solve's convergence. A property edge carries one number, not "and it failed, so +take the other branch". + +J2 solves once and has no such branch, so *only* the conditional argument +applies to it. It could be graph-driven; it would simply cost more than it +saves, and would leave two patterns where there is now one. Convergence therefore travels *with* the result (`{x, converged, iterations}`) rather than being queried from the solver afterwards, where it went stale -between the two solves. +between Drucker-Prager's two solves. --- From b7db7b007b1933b06b39aa764a4057d87d485b51 Mon Sep 17 00:00:00 2001 From: petlenz Date: Sun, 6 Sep 2026 12:53:38 +0200 Subject: [PATCH 5/6] docs: record that material_ref is a decision, not an open gap It was listed under "known gaps", which reads as outstanding work and invites someone to reopen it. It is a choice, so the reasoning is written down instead. The alternative -- a "yielding" flag letting a graph-driven backward_euler know when to iterate -- would work, and the pattern is already native here (strain_threshold_yield publishes is_yielding, isotropic_damage consumes it). It would retire local_newton and material_ref together, about 99 lines of core machinery for two call sites. Rejected on what it costs to express: the material splits into two or three properties where it now has one, and a flag cannot carry Drucker-Prager's apex fallback, which depends on whether the smooth solve converged -- known only after it runs. That becomes two solver instances, three flags and four phases in place of one if-statement. The constraint the decision carries is recorded with it: material_ref bypasses the topological sort, so the plasticity/solver ordering is not an edge the engine knows about. That is safe only because local_newton holds no per-solve state -- solve() is const and returns everything it computes. Give it mutable state and the ordering becomes real and unenforced. --- docs/plasticity.md | 41 ++++++++++++++++++++++++++++++++++++----- 1 file changed, 36 insertions(+), 5 deletions(-) diff --git a/docs/plasticity.md b/docs/plasticity.md index 08a8d6f..482ca08 100644 --- a/docs/plasticity.md +++ b/docs/plasticity.md @@ -317,12 +317,43 @@ argument is indirection — it was removed from `small_strain_plasticity` and fr `rk_plasticity`. A template parameter with two genuinely different arguments is what templates are for. -## 9. Known gaps +## 9. Settled: the solver is reached by `material_ref`, not by the graph + +The return maps hold a `material_ref` and call `solve()`. The +alternative — publishing a `yielding` flag so a graph-driven `backward_euler` +knows when to iterate — was considered and **rejected**. + +The pattern would work: `strain_threshold_yield` already publishes +`is_yielding` and `isotropic_damage` consumes it, so a gating flag is native to +this codebase. It would retire `local_newton` and, with it, `material_ref` — +about 99 lines of core machinery whose only two call sites are these. + +It was rejected on what it would cost to express: + +- **The material splits into phases.** `compute()` currently does trial → + solve → update in one callback. Graph-driven it becomes publish + `yielding`/`residual`/`jacobian`, solver runs, read `delta`, update. Two or + three properties where there is one, and graph dispatch measures ~30 ns even + for a trivial graph. +- **A flag cannot carry Drucker-Prager's branch.** The apex fallback fires when + the *smooth solve fails*, which is not known until after it runs. Expressible + only as trial → smooth solver → `smooth_converged` → apex solver gated on + that flag → update: two solver instances, three flags, four phases, replacing + `if (!smooth.converged) do_apex_return(...)`. +- **The elastic saving is partial** anyway — the solver's callback still fires + and returns early on the flag. + +The cost of the machinery is real but bounded; the cost of removing it is a +graph harder to read than the code it replaces. + +**One constraint this decision carries.** `material_ref` bypasses the +topological sort, so the plasticity↔solver ordering is not an edge the engine +knows about. That is safe only because `local_newton` holds **no per-solve +state** — its `solve()` is `const` and returns everything it computes. Give it +mutable state and the ordering becomes real and unenforced. + +## 10. Known gaps -- **`material_ref`** exists for exactly two call sites (the two backward-Euler - return maps) and costs ~99 lines of core machinery. It bypasses the - topological sort, so the plasticity↔solver ordering is not an edge the engine - knows about. Safe today because `local_newton` holds no per-solve state. - **The apex tangent is a branch tangent.** Verified consistent *on* the branch; a perturbation that leaves the apex back onto the smooth cone is not covered, and cannot be by a central difference. From 9cd963a656631a313cfce72d14e900f9f5ca17a2 Mon Sep 17 00:00:00 2001 From: petlenz Date: Sun, 6 Sep 2026 20:43:21 +0200 Subject: [PATCH 6/6] headers: name each include guard after its own file The 3 header(s) this branch introduces follow the convention set on feature/drucker-prager: the guard is the file's own name, no NUMSIM_MATERIALS_ prefix. Checked against every dependency header and /usr/include for a prior #define of each new name -- none. --- include/numsim-materials/materials/j2_rk_plasticity.h | 6 +++--- include/numsim-materials/solvers/local_newton.h | 6 +++--- include/numsim-materials/solvers/newton_scalar.h | 6 +++--- 3 files changed, 9 insertions(+), 9 deletions(-) diff --git a/include/numsim-materials/materials/j2_rk_plasticity.h b/include/numsim-materials/materials/j2_rk_plasticity.h index 0a7c57a..e796536 100644 --- a/include/numsim-materials/materials/j2_rk_plasticity.h +++ b/include/numsim-materials/materials/j2_rk_plasticity.h @@ -1,5 +1,5 @@ -#ifndef NUMSIM_MATERIALS_J2_RK_PLASTICITY_H -#define NUMSIM_MATERIALS_J2_RK_PLASTICITY_H +#ifndef J2_RK_PLASTICITY_H +#define J2_RK_PLASTICITY_H #include #include @@ -214,4 +214,4 @@ class j2_rk_plasticity final } // namespace numsim::materials -#endif // NUMSIM_MATERIALS_J2_RK_PLASTICITY_H +#endif // J2_RK_PLASTICITY_H diff --git a/include/numsim-materials/solvers/local_newton.h b/include/numsim-materials/solvers/local_newton.h index 9b9bd48..0ae13d8 100644 --- a/include/numsim-materials/solvers/local_newton.h +++ b/include/numsim-materials/solvers/local_newton.h @@ -1,5 +1,5 @@ -#ifndef NUMSIM_MATERIALS_LOCAL_NEWTON_H -#define NUMSIM_MATERIALS_LOCAL_NEWTON_H +#ifndef LOCAL_NEWTON_H +#define LOCAL_NEWTON_H #include #include "numsim-materials/core/material_base.h" @@ -66,4 +66,4 @@ class local_newton final : public material_base, Traits> { } // namespace numsim::materials -#endif // NUMSIM_MATERIALS_LOCAL_NEWTON_H +#endif // LOCAL_NEWTON_H diff --git a/include/numsim-materials/solvers/newton_scalar.h b/include/numsim-materials/solvers/newton_scalar.h index 8e97a95..5cd4aa4 100644 --- a/include/numsim-materials/solvers/newton_scalar.h +++ b/include/numsim-materials/solvers/newton_scalar.h @@ -1,5 +1,5 @@ -#ifndef NUMSIM_MATERIALS_NEWTON_SCALAR_H -#define NUMSIM_MATERIALS_NEWTON_SCALAR_H +#ifndef NEWTON_SCALAR_H +#define NEWTON_SCALAR_H #include #include @@ -70,4 +70,4 @@ class newton_scalar { } // namespace numsim::materials -#endif // NUMSIM_MATERIALS_NEWTON_SCALAR_H +#endif // NEWTON_SCALAR_H