solvers: split the scalar Newton into an algorithm and two interfaces (closes #13) by petlenz · Pull Request #40 · NumSim-Stack/numsim-materials

petlenz · 2026-09-05T09:53:48Z

Stacked on #39, because it changes how the plasticity classes that PR introduces reach their solver.

backward_euler was two solvers in one class, selected by whether a "function" parameter happened to be non-empty:

`function`	mode
set	graph-driven — wires residual/jacobian inputs, engine calls `update()`, result leaves through the `delta` property
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 silent one

"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 against autocatalytic_reaction:

function SET (correct)     degree of cure after 30 steps = 1.000000
function OMITTED           degree of cure after 30 steps = 0.010000   <- start value

No error anywhere. After this change:

Parameter function is missing!

The mode is now a type, not a defaulted string, so the case cannot be expressed.

The others, which fell out rather than needing fixes

Raw pointer members. The inputs existed in one mode only, and a reference cannot be absent — so they were const input_property<...>* with a null guard. input_property already tracks wiring through is_wired(), but an input that is never created cannot be checked. They are references again.
m_converged was never set on the graph path at all. The loop broke on tolerance, on a singular jacobian, and on exhausting its budget — all three indistinguishable from outside. It is reported now.
m_converged{true} initially, so a solver that had never run reported success.
Convergence travels with the result. Drucker-Prager queried converged() between its smooth and apex solves, a side channel that goes stale by construction.

Three pieces

newton_scalar    the algorithm. No properties, no graph, no material_base.
                 Returns {x, converged, iterations}.

backward_euler   graph-driven material. "function" REQUIRED.

local_newton     material-driven. No inputs, no function; exposes solve(eval).

Why the return maps do not use the graph-driven one

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:

if (sig_eq - sigma_0 - H <= 0) { /* elastic */ 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. And at dl = 0 on an elastic step the residual is negative, so Newton drives dl negative; holding it at zero needs a clamp — which is precisely the max(x, 0) this PR removes from the shared solver. The clamp and the graph-driving are the same problem: the graph mode cannot express do not solve.

Drucker-Prager has a second, independent reason: it solves up to twice per update — apex pre-check, smooth cone, then an apex fallback if that fails — and picks the branch on the first solve's convergence. One edge carrying one number cannot express that.

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.

Closes #13 by relocation

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 is about Newton's method, and a general solver silently enforcing one is wrong for every other caller. Each moved to the material that owns the assumption.

Performance

Neutral by design — this splits interfaces, it does not change the arithmetic. Confirmed against
the same interleaved harness used for #39 (15 paired runs): the numbers are those of #39, since
the Newton loop itself is unchanged. J2 and Drucker-Prager reproduce their pre-split trajectories
to ~1 ULP, which is the stronger statement.

Verification

53/53 tests
J2 and Drucker-Prager reproduce their pre-refactor trajectories to ~1 ULP, unchanged from before this commit — so the Newton rewrite moved nothing
Three new tests: the setup that used to be silent, a function material lacking residual/jacobian, and local_newton correctly failing to find a root of x^2 + 1

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.

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.

rk_plasticity<Traits, YieldFunction> 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.

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.

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.

petlenz

Critical review, prompted by a question this PR does not answer: what if I want a different solver?

Probed rather than reasoned about. Three documents, one that works and two that a user would reasonably expect to:

j2 + local_newton (as designed)         OK
j2 + backward_euler (swap the solver)   wire_materials(): material 'j2' references missing materials...
j2_rk_plasticity (pick a tableau)       Key tableau not found

The solver TYPE is compile-time fixed

using solver_type = local_newton<Traits>;
m_solver(base::template add_material_ref<solver_type>(
    base::template get_parameter<std::string>("solver_source")))

solver_source names an INSTANCE. The type is baked in, and material_ref::wire() does a dynamic_cast<T*> that throws on anything else. So a document can say which local_newton, never which kind of solver.

The error is also misleading: naming a backward_euler reports the material as missing, not as the wrong type — because wire_materials() catches every exception and reports them all as absent. A user swapping solvers is told their solver does not exist.

This is not a regression from this PR — small_strain_plasticity hardcoded backward_euler the same way. But this PR is the one that makes the solver a first-class choice with two implementations, so it is the natural place for the limitation to become visible, and arguably the place it hardens: there are now two solver types and plasticity is bound to exactly one.

And `j2_rk_plasticity` cannot be configured from a document at all

m_tableau(base::template get_parameter<const butcher_tableau*>("tableau"))

A raw pointer — the JSON reader has no converter, so "tableau": "sdirk3" fails with Key tableau not found. This is exactly the shape of #33: a C++ object reaching a material through a parameter the document layer cannot express. I removed that pattern from Drucker-Prager in #43 (eta, beta, K_bulk as plain scalars) and left it standing here, in a material this stack renamed and de-templated.

So the honest state is: a document can choose J2 vs Drucker-Prager, and their moduli, but it cannot choose the integrator. That is a strange place to stop for a library whose stated direction is config-driven.

What would fix it, in increasing order of work

"tableau": "sdirk3" as a string. The tableaus are already named factory functions (forward_euler, implicit_euler, sdirk3, gauss_legendre_4, ...). A string-to-tableau converter in the JSON reader registry is small and self-contained, and it makes the explicit/implicit choice a document decision. This one I would do regardless of the rest.
Report the real error. wire_materials() should distinguish not found from wrong type; the information is already there in material_ref::wire() and is being discarded by a blanket catch (...).
A solver interface. For plasticity to accept any scalar solver, material_ref needs a base type — scalar_solver_base<Traits> with a virtual solve(eval). That is a real design change: it makes the solve virtual on a path that runs per integration point, and solve() is currently a template taking any callable, which a virtual cannot be without type-erasing the residual.

The third is the one worth arguing about. The others are gaps, not trade-offs.

What I checked and could not fault

The split itself is sound: newton_scalar has no graph presence, backward_euler requires function, local_newton exposes solve(), and the silent-freeze case is unreachable.
Convergence travelling with the result rather than being queried afterwards is a genuine improvement — Drucker-Prager relied on a side channel that went stale between its two solves.
J2 and Drucker-Prager reproduce their pre-split trajectories to ~1 ULP.
53/53, CI green.

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.

# Conflicts: # include/numsim-materials/materials/j2_rk_plasticity.h

The plasticity chain landed on main while this branch was open, and it moved three things this branch depends on. backward_euler.h conflicted. Both sides were fixing the same defect -- a std::max(x, 0) clamp buried inside a general scalar Newton, which turned a diverged solve into a plausible-looking non-negative answer. This branch split the method into solve() and solve_nonnegative(); main (#40) deleted the callback mode from backward_euler entirely, moving it to local_newton over newton_scalar, which never clamps and returns convergence WITH the value rather than by a separate accessor. main's version is kept whole. Nothing is lost: the non-negativity is a KKT statement and now lives at the three plasticity call sites, where #46 further replaced clamping with a throw for the ill-posed softening case. drucker_prager_plasticity.h did NOT conflict, which is the dangerous part. Git followed the small_strain_plasticity.h -> drucker_prager_plasticity.h rename and silently applied this branch's edit to the renamed file: - return m_solver.get().solve(eval); + return m_solver.get().solve_nonnegative(eval); That would have reintroduced clamping inside the solve -- undoing #46 -- and it was reported as a clean auto-merge. main's version is kept verbatim. Three of this branch's tests in test_materials.cpp exercised the deleted callback API. Two are superseded: NonnegativeVariantClampsConvergedRoot tests a method that no longer exists, and FailurePathIsNotClamped is covered by NewtonScalar.ABudgetExhaustedShortOfTheRootReportsFailure. The third asserted something main does NOT cover -- that a negative root survives the solver -- so it is ported to NewtonScalar.ANegativeRootSurvives rather than dropped. The umat fixtures needed three adaptations to the split plasticity API: - small_strain_plasticity.h supplied a j2_plasticity ALIAS; it is a real class now, so the five umat tests include j2_plasticity.h. They were not unused includes, as they first appeared -- the build caught that. - The fixtures built a backward_euler as the plasticity solver, relying on the callback mode #40 removed. That role is local_newton. The one solver that drives curing is genuinely graph-driven and stays backward_euler. - j2_plasticity builds its own elastic tangent from K and G (#44) and has no elastic_source, so the plasticity blocks pass K. 207/207 tests pass. Against main this branch is now purely additive: 7580 insertions, no deletions. Left alone deliberately: the fixtures still create a linear_elasticity that plasticity no longer consumes -- harmless, and removing it would move the statev layout these tests assert -- and ExternalScalarSource uses temp = 353.0 under the same degrees-C/K confusion fixed on main, which makes both sides of that comparison saturate instantly. Neither belongs in a merge resolution.

Drucker-Prager was held out of the factory with a stated condition: register it once the yield function is expressible from a document. #43 met that condition by making eta, beta and K_bulk plain required scalars instead of members of a C++ yield_function object the JSON reader could not convert. Registered now, and the test that pinned its absence is replaced by two that pin the reason the absence was needed: - a complete Drucker-Prager document builds and yields - a document missing "eta" throws and leaves no material behind, rather than silently getting eta = beta = k = 0 and running as elasticity local_newton is registered too. Without it #33 would still not be closed: every return map names its solver through "solver_source", so a deck could name j2_plasticity but not the solver it requires, and the model still could not be built from a document. The J2 document in this file also needed updating -- it named backward_euler as the plasticity solver, which is the callback mode #40 removed, and passed elastic_source, which #44 removed when plasticity took ownership of its own elastic tangent. One honest limitation recorded in the test rather than papered over: the exception for a missing parameter carries only a COUNT, "missing 1 required parameter(s)". numsim-core's input_parameter_controller prints the names to stdout and throws the count separately, so a deck typo is loud but not self-explanatory. Fixing that is a numsim-core change; what this test pins is that the cone cannot be built without its parameters. 258/258 tests pass on the merge result, verified locally.

TheOptOutListHasNoStaleEntries did exactly what it exists for: kNotForJson exempted small_strain_plasticity and rk_plasticity, and neither header exists any more -- #43 and #40 split them into j2_plasticity, drucker_prager_plasticity and j2_rk_plasticity, each its own header and, since #36, its own factory entry. The exemption for tensor_component_stepper stays: it is a template over Rank, registered as tensor_component_stepper_rank1 and _rank2. 263/263 tests pass on the merge result, verified locally.

petlenz mentioned this pull request Sep 5, 2026

backward_euler::solve() hardcodes a non-negativity clamp (std::max(x,0)) — plasticity-specific, not a general scalar Newton solver #13

Closed

petlenz added 3 commits September 5, 2026 16:42

petlenz changed the base branch from refactor/split-j2-plasticity to refactor/plasticity-owns-tangent September 5, 2026 14:45

petlenz force-pushed the refactor/scalar-newton-split branch from 8cc5866 to b67b840 Compare September 5, 2026 14:45

petlenz mentioned this pull request Sep 5, 2026

materials: split J2 from Drucker-Prager, and give plasticity its own elastic tangent #39

Closed

petlenz added 2 commits September 6, 2026 12:14

petlenz commented Sep 6, 2026

View reviewed changes

This was referenced Sep 6, 2026

materials: plasticity builds its own elastic tangent #44

Merged

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

Merged

petlenz changed the base branch from refactor/plasticity-owns-tangent to main September 6, 2026 19:24

Merge remote-tracking branch 'refs/remotes/tmp/mainx' into wk

bdcf7ed

# Conflicts: # include/numsim-materials/materials/j2_rk_plasticity.h

petlenz merged commit 4cb9b91 into main Sep 6, 2026
1 check passed

petlenz deleted the refactor/scalar-newton-split branch September 6, 2026 19:27

petlenz mentioned this pull request Sep 6, 2026

docs: retire small_strain_plasticity.md, keeping what it alone had #49

Merged

Navigation Menu

Uh oh!

solvers: split the scalar Newton into an algorithm and two interfaces (closes #13) - #40

solvers: split the scalar Newton into an algorithm and two interfaces (closes #13)#40
petlenz merged 7 commits into
mainfrom
refactor/scalar-newton-split

petlenz commented Sep 5, 2026 •

edited

Loading

Uh oh!

petlenz left a comment

Uh oh!

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

1 participant

Uh oh!

Conversation

petlenz commented Sep 5, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

The silent one

The others, which fell out rather than needing fixes

Three pieces

Why the return maps do not use the graph-driven one

Closes #13 by relocation

Performance

Verification

Uh oh!

petlenz left a comment

Choose a reason for hiding this comment

The solver TYPE is compile-time fixed

And j2_rk_plasticity cannot be configured from a document at all

What would fix it, in increasing order of work

What I checked and could not fault

Uh oh!

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

1 participant

petlenz commented Sep 5, 2026 •

edited

Loading

And `j2_rk_plasticity` cannot be configured from a document at all