Skip to content

Phase D PR1: recipe API for implicit residual equations (strain-coupled scalar state) #99

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
petlenz merged 3 commits into from
Jun 16, 2026
Merged

Phase D PR1: recipe API for implicit residual equations (strain-coupled scalar state) #99

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
9ae7ba4
Phase D PR1: recipe API for implicit residual equations (strain-coupl…
Jun 16, 2026
284c42d
Phase D PR1 review: lock the second-rate rejection + cover scalar-inp…
Jun 16, 2026
e1b7121
Phase D PR1 round-2 review: accurate mutual-exclusion message; docume…
Jun 16, 2026
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
246 changes: 178 additions & 68 deletions include/numsim_codegen/recipe.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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<cas::tensor_to_scalar_expression> 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 /
Expand Down Expand Up @@ -481,74 +498,11 @@ class ConstitutiveModel {
ScalarStateVariableHandle const &state_var,
cas::expression_holder<cas::scalar_expression> 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<void const *>(state_var.model_token),
static_cast<void const *>(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<cas::scalar>()` 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<cas::scalar const *>(
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
Expand All @@ -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<cas::tensor_to_scalar_expression> 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,
Expand Down Expand Up @@ -769,6 +753,11 @@ class ConstitutiveModel {
return m_evolution_equations;
}

[[nodiscard]] auto residual_equations() const noexcept
-> std::span<ResidualEquation const> {
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
Expand Down Expand Up @@ -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<void const *>(state_var.model_token),
static_cast<void const *>(this)));
}
auto const *typed =
dynamic_cast<cas::scalar const *>(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<cas::tensor_to_scalar_expression> const &residual,
std::string_view sv_name) const {
LeafCollector lc;
lc.collect_t2s(residual);
std::vector<std::string> 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
Expand Down Expand Up @@ -1043,6 +1152,7 @@ class ConstitutiveModel {
std::vector<OutputDecl> m_outputs;
std::vector<StateVariable> m_state_variables; // Phase 2.1
std::vector<EvolutionEquation> m_evolution_equations; // Phase 2.2
std::vector<ResidualEquation> 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<TangentSpec> m_tangents; // Phase 3b-1 (issue #35)
Expand Down
Loading
Loading