From 1264719dc03f4c1cc56b873a9e9dab94bfc3714c Mon Sep 17 00:00:00 2001 From: petlenz Date: Sun, 13 Sep 2026 18:37:28 +0200 Subject: [PATCH] solvers: cold-start vector_newton instead of seeding from its own outputs solve() seeded m_x by reading back its own output properties. That looks like a warm start and is not one. The unknowns are created with add_output, not add_history_output, so statev_map never collects them -- it skips anything !is_history() -- and unpack() never restores them. One instance serves every integration point on a thread: umat_interface builds one material_context per thread and states why that is sufficient, "a consequence of the evaluator being stateless: all point state lives in the host's STATEV array". material_point_evaluator is blunter: "anything the evaluator remembered across calls would be state from a trial the host has since discarded". Seeding from the properties is exactly such a memory. Point N starts from point N-1's converged answer, and after a failure from a diverged iterate. For a system with several real roots that does not merely change the iteration count, it selects WHICH root is returned -- converged, plausible, and not the root this point is on. Cold start every solve, which makes the solver point-local by construction. That is what backward_euler already does: it passes a fixed literal seed rather than reading back its own "delta". Two tests, and the first is the one that matters. It drives two different problems through ONE instance, changing them in place with set_parameter, and requires each answer to equal a fresh instance's. x^2 + y = a; x + y^2 = b has several real roots, so the seed is observable in the result rather than only in the iteration count. The second does the same with a diverging problem in between, for the failure path. An earlier version of this change kept the property seeding and only held it back after a failure. That was the wrong level: it removed failure-dependence and left point-order-dependence, and added a second piece of hidden cross-call state to a class whose surrounding contract is that nothing is retained between calls. SolvesArePointLocal fails against that version too, which is the point of writing it this way. Warm-starting is still worth having -- plane_stress_evaluator keeps eps_33 in its own STATEV slot precisely because a cold start costs several extra graph evaluations per Gauss point -- but it has to be per point, which means the unknowns becoming history properties. Left to #21 with the zero-seed caveat for a model whose jacobian is singular at the origin. No production caller: vector_newton appears only in the factory registration and its own tests. 282/282 pass. --- .../numsim-materials/solvers/vector_newton.h | 31 ++++++++- tests/test_vector_newton.cpp | 67 +++++++++++++++++++ 2 files changed, 97 insertions(+), 1 deletion(-) diff --git a/include/numsim-materials/solvers/vector_newton.h b/include/numsim-materials/solvers/vector_newton.h index 95d1cd0..49584e9 100644 --- a/include/numsim-materials/solvers/vector_newton.h +++ b/include/numsim-materials/solvers/vector_newton.h @@ -165,7 +165,36 @@ class vector_newton final void solve() { m_converged = false; - gather_state(m_x); // seed from whatever the properties currently hold + + // Cold start, every solve. This class used to seed from its own output + // properties, which looks like a warm start and is not one. + // + // The unknowns are plain outputs, not history properties, so statev_map + // never collects them (it skips anything !is_history()) and unpack() never + // restores them. One instance serves every integration point on a thread -- + // umat_interface builds one material_context per thread and says why that + // is sufficient: "a consequence of the evaluator being stateless: all point + // state lives in the host's STATEV array". material_point_evaluator is + // blunter still: "anything the evaluator remembered across calls would be + // state from a trial the host has since discarded". + // + // Seeding from the properties is exactly such a memory. It seeds point N + // from point N-1's converged answer, and after a failure from a diverged + // iterate. For a system with several roots that silently selects a + // different one -- converged, plausible, and not the root this point is on. + // + // backward_euler is point-local by construction for the same reason: it + // passes a fixed literal seed rather than reading back its own "delta". + // + // Warm-starting is worth having, and plane_stress_evaluator shows the shape + // -- it keeps eps_33 in its own STATEV slot precisely because a cold start + // costs several extra graph evaluations per Gauss point. But it has to be + // PER POINT, which means the unknowns becoming history properties. Tracked + // in #21; until then, correct and slower beats fast and point-dependent. + // + // Zero is the seed. A model whose jacobian is singular there needs an + // explicit initial guess, which is the same #21 follow-up. + m_x.setZero(m_N); // Evaluate-first: the residual AFTER the final update is the one tested, so // a solve that converges on its last allowed update is reported converged. diff --git a/tests/test_vector_newton.cpp b/tests/test_vector_newton.cpp index eb9cefe..238b40f 100644 --- a/tests/test_vector_newton.cpp +++ b/tests/test_vector_newton.cpp @@ -180,6 +180,73 @@ TEST(VectorNewton, NonlinearSystemConvergesToAnalyticRoot) { EXPECT_NEAR(f.ctx.get("solver", "y"), 2.0, 1e-9); } +/// One solver instance, two different problems: each must land where a fresh +/// instance lands. +/// +/// This is the property the class needs and did not have. vector_newton's +/// unknowns are plain outputs, so statev_map never collects them and unpack() +/// never restores them; one instance serves every integration point on a +/// thread, driven by an evaluator whose stated design is that nothing is +/// retained between calls. solve() nevertheless seeded itself by reading those +/// properties back, so point N started from point N-1's answer. +/// +/// x^2 + y = a ; x + y^2 = b has several real roots, so a seed does not merely +/// change the iteration count -- it selects WHICH root is returned, converged +/// and plausible and wrong for this point. +TEST(VectorNewton, SolvesArePointLocal) { + auto fresh = [](T a, T b) { + scalar_fixture f(system_mode::nonlinear, a, b); + f.solver->solve(); + EXPECT_TRUE(f.solver->converged()); + return std::pair{f.ctx.get("solver", "x"), + f.ctx.get("solver", "y")}; + }; + const auto first = fresh(5.248, 27.4); + const auto second = fresh(3.0, 5.0); + + // Now drive BOTH through one instance, changing the problem in place. + scalar_fixture shared(system_mode::nonlinear, 5.248, 27.4); + auto* sys = shared.ctx.find("sys"); + ASSERT_NE(sys, nullptr); + shared.solver->solve(); + ASSERT_TRUE(shared.solver->converged()); + EXPECT_DOUBLE_EQ(shared.ctx.get("solver", "x"), first.first); + + sys->template set_parameter("a", 3.0); + sys->template set_parameter("b", 5.0); + shared.solver->solve(); + ASSERT_TRUE(shared.solver->converged()); + EXPECT_DOUBLE_EQ(shared.ctx.get("solver", "x"), second.first) + << "the second solve was seeded by the first: this instance returned a " + "different root than a fresh one does for the same problem"; + EXPECT_DOUBLE_EQ(shared.ctx.get("solver", "y"), second.second); +} + +/// And a FAILED solve in between must not move the one after it either. +/// +/// On failure the output properties deliberately keep the last raw iterate, so +/// that a failed solve cannot pass for a plausible answer. Seeding the next +/// solve from them turned that honesty into a hazard. +TEST(VectorNewton, AFailedSolveDoesNotMoveTheNext) { + scalar_fixture reference(system_mode::nonlinear, 3.0, 5.0); + reference.solver->solve(); + ASSERT_TRUE(reference.solver->converged()); + const auto x_ref = reference.ctx.get("solver", "x"); + + scalar_fixture f(system_mode::nonlinear, -10.0, -10.0); + f.solver->solve(); + ASSERT_FALSE(f.solver->converged()) << "the setup must actually fail"; + + auto* sys = f.ctx.find("sys"); + ASSERT_NE(sys, nullptr); + sys->template set_parameter("a", 3.0); + sys->template set_parameter("b", 5.0); + f.solver->solve(); + ASSERT_TRUE(f.solver->converged()); + EXPECT_DOUBLE_EQ(f.ctx.get("solver", "x"), x_ref) + << "the solve after a failure was seeded from the diverged iterate"; +} + TEST(VectorNewton, SingularJacobianReportsFailure) { scalar_fixture f(system_mode::singular, 4.0, 7.0); // inconsistent: a != b f.solver->solve();