Skip to content
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
31 changes: 30 additions & 1 deletion include/numsim-materials/solvers/vector_newton.h
Original file line number Diff line number Diff line change
Expand Up @@ -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.
Expand Down
67 changes: 67 additions & 0 deletions tests/test_vector_newton.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -180,6 +180,73 @@ TEST(VectorNewton, NonlinearSystemConvergesToAnalyticRoot) {
EXPECT_NEAR(f.ctx.get<T>("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<T, T>{f.ctx.get<T>("solver", "x"),
f.ctx.get<T>("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<T>("solver", "x"), first.first);

sys->template set_parameter<T>("a", 3.0);
sys->template set_parameter<T>("b", 5.0);
shared.solver->solve();
ASSERT_TRUE(shared.solver->converged());
EXPECT_DOUBLE_EQ(shared.ctx.get<T>("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<T>("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<T>("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<T>("a", 3.0);
sys->template set_parameter<T>("b", 5.0);
f.solver->solve();
ASSERT_TRUE(f.solver->converged());
EXPECT_DOUBLE_EQ(f.ctx.get<T>("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();
Expand Down
Loading