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();