diff --git a/include/numsim_codegen/recipe.h b/include/numsim_codegen/recipe.h index 1ee63b0..d584f22 100644 --- a/include/numsim_codegen/recipe.h +++ b/include/numsim_codegen/recipe.h @@ -219,6 +219,23 @@ struct EvolutionEquation { std::string doc; }; +// Phase D strain-coupled (verified-reachable 2026-06-15): an IMPLICIT evolution. +// Instead of a rate `dx/dt = f(x)`, the state x is defined by a residual +// `R(x, inputs) = 0` solved by a Newton solver (numsim-materials' backward_euler +// on the graph-coupled path). Unlike a rate — which the rk_integrator contract +// forbids from referencing inputs (the integrator owns discretization) — a +// residual is EXPECTED to depend on a tensor strain input (that is the coupling), +// so its expression is `tensor_to_scalar`-typed: a scalar R built from the scalar +// state/params AND tensor inputs (e.g. `x - c*trace(eps)`). This unlocks real +// return-map / plasticity-class models and a strain-coupled consistent tangent +// dσ/dε = ∂σ/∂ε + ∂σ/∂x·(−∂R/∂ε / ∂R/∂x). A state variable carries EITHER a rate +// (EvolutionEquation) OR a residual (this) — never both. +struct ResidualEquation { + std::size_t state_variable_idx; + cas::expression_holder residual; + std::string doc; +}; + // Phase 3a-2 (issue #75): tuning for the in-function Newton solve emitted // by LocalNewtonLoweringPass. `tol` is the absolute residual threshold for // convergence; `max_iter` caps the iteration count (no line search / @@ -481,74 +498,11 @@ class ConstitutiveModel { ScalarStateVariableHandle const &state_var, cas::expression_holder rate, std::string doc = "") -> void { - // PR #69 round-1 #3: cross-recipe hijack defense. The handle's - // `model_token` MUST equal `this` — otherwise the user passed a - // handle that came from a different ConstitutiveModel (or - // manually-constructed). Pure name-based matching would silently - // bind to a coincidentally-named state variable on this model. - if (state_var.model_token == nullptr) { - throw std::runtime_error(std::format( - "ConstitutiveModel '{}': add_scalar_evolution_equation handle " - "has no model token. Handles must come from " - "add_scalar_state_variable on this model — manually-constructed " - "handles are not accepted.", - m_name)); - } - if (state_var.model_token != this) { - throw std::runtime_error(std::format( - "ConstitutiveModel '{}': add_scalar_evolution_equation handle " - "came from a different ConstitutiveModel (handle.model_token = " - "{}, this = {}). Each handle is bound to the recipe that " - "created it; cross-recipe use would silently bind by name and " - "discretise the wrong state variable.", - m_name, static_cast(state_var.model_token), - static_cast(this))); - } - - // PR #69 round-1 review (CRITICAL): the handle's `current` MUST be a - // bare scalar leaf (i.e. `cas::scalar`), not a compound expression. - // `expression_holder::get()` would use an - // `assert(dynamic_cast != nullptr)` path that's stripped under - // `NDEBUG`, then UB on the unchecked `static_cast`. Use an explicit - // runtime dynamic_cast that throws a clear diagnostic instead. - auto const *typed = dynamic_cast( - state_var.current.data().get()); - if (!typed) { - throw std::runtime_error(std::format( - "ConstitutiveModel '{}': add_scalar_evolution_equation handle's " - "`current` is not a bare scalar leaf symbol. Handles must come " - "from add_scalar_state_variable — you cannot synthesise one " - "from a compound expression like `K * alpha`.", - m_name)); - } - auto const &sv_name = typed->name(); - - std::size_t found_idx = m_state_variables.size(); - for (std::size_t i = 0; i < m_state_variables.size(); ++i) { - if (m_state_variables[i].kind == SymbolDecl::Kind::Scalar && - m_state_variables[i].name == sv_name) { - found_idx = i; - break; - } - } - if (found_idx == m_state_variables.size()) { - // List registered scalar state variables to help debug. - std::string registered; - for (auto const &sv : m_state_variables) { - if (sv.kind == SymbolDecl::Kind::Scalar) { - if (!registered.empty()) registered += ", "; - registered += sv.name; - } - } - throw std::runtime_error(std::format( - "ConstitutiveModel '{}': add_scalar_evolution_equation handle " - "names scalar state variable '{}' but no such state variable " - "is registered on this model. Did the handle come from a " - "different ConstitutiveModel, or was the state variable added " - "with add_tensor_state_variable instead? Registered scalar " - "state variables: [{}].", - m_name, sv_name, registered)); - } + auto const found_idx = resolve_scalar_state_var_index_( + state_var, "add_scalar_evolution_equation"); + auto const &sv_name = m_state_variables[found_idx].name; + assert_state_var_unbound_(found_idx, sv_name, + "add_scalar_evolution_equation"); // Validate that the rate expression's leaves are all declared // symbols on this model (PR #69 round-1 #4). Fail fast in the @@ -562,6 +516,36 @@ class ConstitutiveModel { m_evolution_equations.push_back(std::move(eq)); } + // Phase D strain-coupled: declare an IMPLICIT residual `R(x, inputs) = 0` for + // an already-added scalar state variable, solved by a Newton solver. Unlike a + // rate, the residual MAY (and typically does) reference tensor inputs (strain) + // — that is the coupling — hence the `tensor_to_scalar`-typed residual. A state + // variable carries EITHER a rate OR a residual, never both (enforced). + // + // Scope (current): the residual is `tensor_to_scalar`-typed, so a purely + // scalar implicit residual (no tensor dependence) is intentionally NOT + // expressible here — a strain-independent scalar evolution is the rate path's + // job (add_scalar_evolution_equation). The residual also cannot reference the + // framework time step `dt` (it is not auto-registered on this path): the + // first target is rate-INDEPENDENT (e.g. return-map plasticity); a + // rate-dependent (viscoplastic) residual needing `dt` is a follow-up. The + // residual's differentiability (∂R/∂x, ∂R/∂ε) is checked at EMIT time, where + // a non-differentiable t2s node (e.g. a piecewise if_then_else, cas#241) + // surfaces a clear cas error. + auto add_scalar_residual_equation( + ScalarStateVariableHandle const &state_var, + cas::expression_holder residual, + std::string doc = "") -> void { + auto const found_idx = resolve_scalar_state_var_index_( + state_var, "add_scalar_residual_equation"); + auto const &sv_name = m_state_variables[found_idx].name; + assert_state_var_unbound_(found_idx, sv_name, + "add_scalar_residual_equation"); + validate_residual_expression_leaves_(residual, sv_name); + ResidualEquation eq{found_idx, std::move(residual), std::move(doc)}; + m_residual_equations.push_back(std::move(eq)); + } + // ─── Output declarations ──────────────────────────────────────── void add_output(std::string name, @@ -769,6 +753,11 @@ class ConstitutiveModel { return m_evolution_equations; } + [[nodiscard]] auto residual_equations() const noexcept + -> std::span { + return m_residual_equations; + } + // Phase 3a-2 (issue #75): opt into in-function local Newton solving. // When enabled (and the recipe has evolution equations), // `emit_compute_function` registers `LocalNewtonLoweringPass` instead @@ -994,6 +983,126 @@ class ConstitutiveModel { m_parameters_cache.back().is_time_step = true; } + // Resolve a ScalarStateVariableHandle to its index in m_state_variables, + // with the cross-recipe-hijack + bare-leaf defenses (PR #69 round-1 #3). + // `caller` names the public method for the diagnostic. Shared by + // add_scalar_evolution_equation and add_scalar_residual_equation. + [[nodiscard]] auto resolve_scalar_state_var_index_( + ScalarStateVariableHandle const &state_var, std::string_view caller) const + -> std::size_t { + if (state_var.model_token == nullptr) { + throw std::runtime_error(std::format( + "ConstitutiveModel '{}': {} handle has no model token. Handles must " + "come from add_scalar_state_variable on this model — " + "manually-constructed handles are not accepted.", + m_name, caller)); + } + if (state_var.model_token != this) { + throw std::runtime_error(std::format( + "ConstitutiveModel '{}': {} handle came from a different " + "ConstitutiveModel (handle.model_token = {}, this = {}). Each handle " + "is bound to the recipe that created it; cross-recipe use would " + "silently bind by name and discretise the wrong state variable.", + m_name, caller, static_cast(state_var.model_token), + static_cast(this))); + } + auto const *typed = + dynamic_cast(state_var.current.data().get()); + if (!typed) { + throw std::runtime_error(std::format( + "ConstitutiveModel '{}': {} handle's `current` is not a bare scalar " + "leaf symbol. Handles must come from add_scalar_state_variable — you " + "cannot synthesise one from a compound expression like `K * alpha`.", + m_name, caller)); + } + auto const &sv_name = typed->name(); + for (std::size_t i = 0; i < m_state_variables.size(); ++i) { + if (m_state_variables[i].kind == SymbolDecl::Kind::Scalar && + m_state_variables[i].name == sv_name) { + return i; + } + } + std::string registered; + for (auto const &sv : m_state_variables) { + if (sv.kind == SymbolDecl::Kind::Scalar) { + if (!registered.empty()) registered += ", "; + registered += sv.name; + } + } + throw std::runtime_error(std::format( + "ConstitutiveModel '{}': {} handle names scalar state variable '{}' but " + "no such state variable is registered on this model. Did the handle " + "come from a different ConstitutiveModel, or was the state variable " + "added with add_tensor_state_variable instead? Registered scalar state " + "variables: [{}].", + m_name, caller, sv_name, registered)); + } + + // A scalar state variable carries EXACTLY ONE evolution mechanism — a rate + // (EvolutionEquation) or an implicit residual (ResidualEquation). Reject a + // second binding rather than emitting a contradictory material. + void assert_state_var_unbound_(std::size_t idx, std::string_view sv_name, + std::string_view caller) const { + for (auto const &e : m_evolution_equations) { + if (e.state_variable_idx == idx) { + throw std::runtime_error(std::format( + "ConstitutiveModel '{}': {} — state variable '{}' already has a " + "rate evolution equation; a state variable carries at most one " + "rate or residual.", + m_name, caller, sv_name)); + } + } + for (auto const &r : m_residual_equations) { + if (r.state_variable_idx == idx) { + throw std::runtime_error(std::format( + "ConstitutiveModel '{}': {} — state variable '{}' already has an " + "implicit residual equation; a state variable carries at most one " + "rate or residual.", + m_name, caller, sv_name)); + } + } + } + + // Like validate_rate_expression_leaves_ but for a t2s residual: every leaf + // (scalar AND tensor) must be a declared symbol, and the state must itself + // appear (else ∂R/∂x ≡ 0 — a singular Newton Jacobian). + void validate_residual_expression_leaves_( + cas::expression_holder const &residual, + std::string_view sv_name) const { + LeafCollector lc; + lc.collect_t2s(residual); + std::vector missing; + for (auto const &leaf_name : lc.scalar_names()) { + bool found = false; + for (auto const &[name, _] : m_scalar_symbols) + if (name == leaf_name) { found = true; break; } + if (!found) missing.push_back("scalar '" + leaf_name + "'"); + } + for (auto const &leaf_name : lc.tensor_names()) { + bool found = false; + for (auto const &[name, _] : m_tensor_symbols) + if (name == leaf_name) { found = true; break; } + if (!found) missing.push_back("tensor '" + leaf_name + "'"); + } + if (!missing.empty()) { + std::string msg = std::format( + "ConstitutiveModel '{}': add_scalar_residual_equation residual for " + "state variable '{}' references undeclared symbol(s):", + m_name, sv_name); + for (auto const &m : missing) msg += std::format("\n - {}", m); + msg += "\nCall add_scalar_input / add_tensor_input / add_parameter before " + "referencing a symbol in a residual."; + throw std::runtime_error(msg); + } + if (!lc.scalar_names().contains(std::string(sv_name))) { + throw std::runtime_error(std::format( + "ConstitutiveModel '{}': add_scalar_residual_equation residual for " + "state variable '{}' does not reference its own state — the residual " + "must depend on '{}' (else ∂R/∂{} ≡ 0, a singular Newton Jacobian).", + m_name, sv_name, sv_name, sv_name)); + } + } + // PR #69 round-1 #4: validate that the leaves of a rate expression // are all registered symbols on this model. Called from // `add_scalar_evolution_equation` to fail fast in the user's stack @@ -1043,6 +1152,7 @@ class ConstitutiveModel { std::vector m_outputs; std::vector m_state_variables; // Phase 2.1 std::vector m_evolution_equations; // Phase 2.2 + std::vector m_residual_equations; // Phase D strain-coupled bool m_local_newton = false; // Phase 3a-2 (issue #75) NewtonOptions m_newton_options{}; // Phase 3a-2 (issue #75) std::vector m_tangents; // Phase 3b-1 (issue #35) diff --git a/tests/RecipeTest.cpp b/tests/RecipeTest.cpp index da1ba37..ebd222e 100644 --- a/tests/RecipeTest.cpp +++ b/tests/RecipeTest.cpp @@ -94,4 +94,149 @@ TEST(Recipe, MultipleOutputsShareCseAcrossBody) { << src; } +// ─── Phase D: implicit residual equations (strain-coupled state) ───────────── + +// A scalar state defined by an implicit residual R(z, ε)=0. Unlike a rate, the +// residual MAY reference a tensor strain input — that is the coupling. +TEST(Recipe, AddScalarResidualEquationStoresStrainCoupledResidual) { + ConstitutiveModel m("ReturnMap"); + auto c = m.add_parameter("c", 2.0); + auto eps = m.add_tensor_input("strain", 3, 2, roles::Strain); + auto z = m.add_scalar_state_variable( + "z", cas::make_expression(0.0)); + m.add_scalar_residual_equation(z, z.current - c * trace(eps)); // R = z - c·tr(ε) + ASSERT_EQ(m.residual_equations().size(), 1u); + EXPECT_TRUE(m.residual_equations()[0].residual.is_valid()); + EXPECT_TRUE(m.evolution_equations().empty()); +} + +// The residual's leaves must all be declared symbols. +TEST(Recipe, ResidualRejectsUndeclaredLeaf) { + ConstitutiveModel m("Bad"); + auto eps = m.add_tensor_input("strain", 3, 2, roles::Strain); + auto z = m.add_scalar_state_variable( + "z", cas::make_expression(0.0)); + auto bogus = cas::make_expression("T_bogus"); + try { + m.add_scalar_residual_equation(z, z.current - bogus * trace(eps)); + FAIL() << "expected throw on undeclared leaf"; + } catch (std::exception const &e) { + EXPECT_NE(std::string(e.what()).find("T_bogus"), std::string::npos) + << e.what(); + } +} + +// The residual must depend on its own state, else ∂R/∂z ≡ 0 (singular Jacobian). +TEST(Recipe, ResidualMustReferenceOwnState) { + ConstitutiveModel m("NoState"); + auto c = m.add_parameter("c", 2.0); + auto eps = m.add_tensor_input("strain", 3, 2, roles::Strain); + auto z = m.add_scalar_state_variable( + "z", cas::make_expression(0.0)); + try { + m.add_scalar_residual_equation(z, c * trace(eps)); // no z + FAIL() << "expected throw on state-independent residual"; + } catch (std::exception const &e) { + EXPECT_NE(std::string(e.what()).find("does not reference its own state"), + std::string::npos) + << e.what(); + } +} + +// A state carries either a rate OR a residual — never both (rate first). +TEST(Recipe, StateCannotHaveRateThenResidual) { + ConstitutiveModel m("Both"); + auto c = m.add_parameter("c", 2.0); + auto eps = m.add_tensor_input("strain", 3, 2, roles::Strain); + auto z = m.add_scalar_state_variable( + "z", cas::make_expression(0.0)); + m.add_scalar_evolution_equation(z, c * z.current); + try { + m.add_scalar_residual_equation(z, z.current - c * trace(eps)); + FAIL() << "expected throw: state already has a rate"; + } catch (std::exception const &e) { + EXPECT_NE(std::string(e.what()).find("at most one rate or residual"), + std::string::npos) + << e.what(); + } +} + +// ...and residual first, then rate. +TEST(Recipe, StateCannotHaveResidualThenRate) { + ConstitutiveModel m("Both2"); + auto c = m.add_parameter("c", 2.0); + auto eps = m.add_tensor_input("strain", 3, 2, roles::Strain); + auto z = m.add_scalar_state_variable( + "z", cas::make_expression(0.0)); + m.add_scalar_residual_equation(z, z.current - c * trace(eps)); + try { + m.add_scalar_evolution_equation(z, c * z.current); + FAIL() << "expected throw: state already has a residual"; + } catch (std::exception const &e) { + EXPECT_NE(std::string(e.what()).find("at most one rate or residual"), + std::string::npos) + << e.what(); + } +} + +// The state may appear only as a tensor COEFFICIENT inside the t2s (e.g. +// trace(z·ε)) — the state-appears guard must still find it, i.e. collect_t2s +// recurses scalar coefficients. A false rejection here would block valid models. +TEST(Recipe, ResidualStateMayAppearAsTensorCoefficient) { + ConstitutiveModel m("CoeffState"); + auto c = m.add_parameter("c", 2.0); + auto eps = m.add_tensor_input("strain", 3, 2, roles::Strain); + auto z = m.add_scalar_state_variable( + "z", cas::make_expression(0.0)); + m.add_scalar_residual_equation(z, trace(z.current * eps) - c); + EXPECT_EQ(m.residual_equations().size(), 1u); +} + +// A residual may reference a declared SCALAR input (not only tensor inputs). +TEST(Recipe, ResidualMayReferenceScalarInput) { + ConstitutiveModel m("ScalarInputResid"); + auto temp = m.add_scalar_input("temperature"); + auto eps = m.add_tensor_input("strain", 3, 2, roles::Strain); + auto z = m.add_scalar_state_variable( + "z", cas::make_expression(0.0)); + m.add_scalar_residual_equation(z, z.current - temp * trace(eps)); + EXPECT_EQ(m.residual_equations().size(), 1u); +} + +// A state variable is bound to at most ONE evolution mechanism — adding a second +// rate to the same state is rejected (the XOR guard's evolution-side branch). +TEST(Recipe, StateRejectsSecondRate) { + ConstitutiveModel m("TwoRates"); + auto c = m.add_parameter("c", 2.0); + auto z = m.add_scalar_state_variable( + "z", cas::make_expression(0.0)); + m.add_scalar_evolution_equation(z, c * z.current); + try { + m.add_scalar_evolution_equation(z, c * z.current); + FAIL() << "expected throw: state already has a rate"; + } catch (std::exception const &e) { + EXPECT_NE(std::string(e.what()).find("already has a rate"), + std::string::npos) + << e.what(); + } +} + +// The shared handle-resolution defends against cross-recipe handle use. +TEST(Recipe, ResidualRejectsForeignHandle) { + ConstitutiveModel m1("M1"); + ConstitutiveModel m2("M2"); + auto c = m2.add_parameter("c", 2.0); + auto eps = m2.add_tensor_input("strain", 3, 2, roles::Strain); + auto z1 = m1.add_scalar_state_variable( + "z", cas::make_expression(0.0)); + try { + m2.add_scalar_residual_equation(z1, z1.current - c * trace(eps)); // m1 handle + FAIL() << "expected throw on foreign handle"; + } catch (std::exception const &e) { + EXPECT_NE(std::string(e.what()).find("different ConstitutiveModel"), + std::string::npos) + << e.what(); + } +} + } // namespace numsim::codegen