From 90f98d9f7d45e095d0d79f1aabae9f355346c85d Mon Sep 17 00:00:00 2001 From: petlenz Date: Wed, 17 Jun 2026 21:54:53 +0200 Subject: [PATCH 1/3] Phase 2a: reject implicit-residual recipes loudly (Finding A) before Mode-B emit lands --- include/numsim_codegen/recipe.h | 17 +++++++++++++++++ src/targets/numsim_material.cpp | 17 +++++++++++++++++ tests/NumSimMaterialTargetTest.cpp | 22 ++++++++++++++++++++++ tests/RecipeTest.cpp | 21 +++++++++++++++++++++ 4 files changed, 77 insertions(+) diff --git a/include/numsim_codegen/recipe.h b/include/numsim_codegen/recipe.h index d584f22..31e1051 100644 --- a/include/numsim_codegen/recipe.h +++ b/include/numsim_codegen/recipe.h @@ -634,6 +634,23 @@ class ConstitutiveModel { // change" while the pass's mutation lives only inside this call. // ConstitutiveModel is value-typed (vectors + shared_ptr handles) // so the copy is cheap relative to the codegen itself. + // Strain-coupled implicit residuals (add_scalar_residual_equation) have NO + // emit path on the self-contained (standalone / MOOSE-local) targets: the + // self-contained pipeline below has no pass that lowers a residual into a + // Newton solve, so a residual-only recipe would otherwise emit a compute + // function that silently drops the declared state — a correctness hazard + // (reject loudly, never drop). Implicit-residual emission lands only on the + // graph-coupled NumSimMaterialTarget (Mode B: material_ref + + // solve()); see the roadmap Phase D. Reject here rather than emit a partial + // function. + if (!m_residual_equations.empty()) { + throw std::runtime_error( + "ConstitutiveModel::emit_compute_function: implicit residual " + "equations (add_scalar_residual_equation) are not supported by the " + "self-contained (standalone / MOOSE) code path — they would be " + "silently dropped. Strain-coupled residual materials are emitted only " + "by NumSimMaterialTarget (graph-coupled, Mode B)."); + } ConstitutiveModel working_copy = *this; PassContext pctx{RecipeView{working_copy}, CodeGenContext{}, std::nullopt, {}}; diff --git a/src/targets/numsim_material.cpp b/src/targets/numsim_material.cpp index bdb7ae1..c1e711f 100644 --- a/src/targets/numsim_material.cpp +++ b/src/targets/numsim_material.cpp @@ -68,6 +68,23 @@ std::string tensor_cxx_type(std::size_t dim, std::size_t rank) { // whole contract — a code generator that quietly drops a declared output is a // correctness hazard. void check_scope(ConstitutiveModel const &model) { + // Strain-coupled implicit-residual recipes (add_scalar_residual_equation) are + // a distinct contract from the rate form below: the state is defined by an + // implicit R(x, ε)=0 solved by backward_euler, not a rate integrated by + // rk_integrator. The Mode-B emission path for those lands separately; until + // then reject with a message that names the real reason — NOT the rate / + // rk_integrator / vector-solver message below, which would misdiagnose a + // residual recipe (it has one state variable but zero evolution equations, so + // it would otherwise trip the "exactly one evolution equation" check and blame + // the wrong contract). + if (!model.residual_equations().empty()) { + throw std::runtime_error( + "NumSimMaterialTarget: implicit residual equations " + "(add_scalar_residual_equation, strain-coupled state) are not yet " + "emitted — the graph-coupled Mode-B residual path (material_ref<" + "backward_euler> + solve) is a follow-up. This recipe declares a " + "residual, not a rate/rk_integrator evolution equation."); + } auto const svs = model.state_variables(); auto const eqs = model.evolution_equations(); if (svs.size() != 1 || eqs.size() != 1) { diff --git a/tests/NumSimMaterialTargetTest.cpp b/tests/NumSimMaterialTargetTest.cpp index 4740e66..67f9559 100644 --- a/tests/NumSimMaterialTargetTest.cpp +++ b/tests/NumSimMaterialTargetTest.cpp @@ -12,6 +12,8 @@ #include #include #include +#include +#include #include @@ -329,5 +331,25 @@ TEST(NumSimMaterialTarget, FloatDefaultRoundTrips) { << src; } +// Finding A (holistic review 2026-06-17): a recipe defined by an implicit +// residual (add_scalar_residual_equation) has no Mode-B emission path yet. It +// must be rejected with a message that names the RESIDUAL contract — not the +// rate/rk_integrator/vector-solver message, which would misdiagnose it (a +// residual recipe has one state variable but zero evolution equations). +TEST(NumSimMaterialTarget, RejectsResidualRecipeWithAccurateMessage) { + 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", make_expression(0.0)); + m.add_scalar_residual_equation(z, z.current - c * trace(eps)); + auto const msg = emit_throw_message(m); + // Names the residual contract... + EXPECT_NE(msg.find("residual"), std::string::npos) << msg; + // ...and does NOT misdiagnose as the rate/evolution-equation contract. + EXPECT_EQ(msg.find("exactly one scalar state variable"), std::string::npos) + << msg; +} + } // namespace } // namespace numsim::codegen diff --git a/tests/RecipeTest.cpp b/tests/RecipeTest.cpp index ebd222e..503460d 100644 --- a/tests/RecipeTest.cpp +++ b/tests/RecipeTest.cpp @@ -221,6 +221,27 @@ TEST(Recipe, StateRejectsSecondRate) { } } +// Finding A (holistic review 2026-06-17): the self-contained code path +// (standalone / MOOSE) has no pass that lowers an implicit residual into a +// Newton solve, so emit_compute_function would otherwise emit a function that +// SILENTLY DROPS the declared state. Reject loudly instead — residual emission +// lives only on the graph-coupled NumSimMaterialTarget (Mode B). +TEST(Recipe, EmitComputeFunctionRejectsResidualRecipe) { + ConstitutiveModel m("ReturnMapStandalone"); + 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 { + (void)m.emit_compute_function(); + FAIL() << "expected throw: residuals unsupported on the self-contained path"; + } catch (std::exception const &e) { + EXPECT_NE(std::string(e.what()).find("residual"), std::string::npos) + << e.what(); + } +} + // The shared handle-resolution defends against cross-recipe handle use. TEST(Recipe, ResidualRejectsForeignHandle) { ConstitutiveModel m1("M1"); From f6edc2a1790e3872a2d7396791f4b82c45aeeeef Mon Sep 17 00:00:00 2001 From: petlenz Date: Wed, 17 Jun 2026 22:06:53 +0200 Subject: [PATCH 2/3] Phase 2a: emit strain-coupled implicit-residual materials (Mode B) via NumSimMaterialTarget --- src/targets/numsim_material.cpp | 462 +++++++++++++++++- tests/CMakeLists.txt | 8 +- tests/NumSimMaterialTargetTest.cpp | 79 ++- .../generate_numsim_material_check.cpp | 22 +- .../numsim_material_check_driver.cpp | 47 ++ 5 files changed, 585 insertions(+), 33 deletions(-) diff --git a/src/targets/numsim_material.cpp b/src/targets/numsim_material.cpp index c1e711f..e3e68ca 100644 --- a/src/targets/numsim_material.cpp +++ b/src/targets/numsim_material.cpp @@ -7,6 +7,7 @@ #include #include +#include #include #include @@ -34,6 +35,15 @@ constexpr char const *state_property = "state"; constexpr char const *integrator_source_param = "integrator_source"; constexpr char const *material_base_include = "numsim-materials/core/material_base.h"; +// Mode-B (strain-coupled implicit residual) contract surface. The material +// holds a material_ref and drives the Newton loop itself via +// solve(eval). Validated against numsim-materials `materials/ +// small_strain_plasticity.h` + `solvers/backward_euler.h`. +constexpr char const *solver_source_param = "solver_source"; +constexpr char const *backward_euler_include = + "numsim-materials/solvers/backward_euler.h"; +constexpr char const *material_ref_include = + "numsim-materials/core/material_ref.h"; } // namespace contract // A recipe symbol whose name equals one of the emitted fixed members would @@ -68,23 +78,6 @@ std::string tensor_cxx_type(std::size_t dim, std::size_t rank) { // whole contract — a code generator that quietly drops a declared output is a // correctness hazard. void check_scope(ConstitutiveModel const &model) { - // Strain-coupled implicit-residual recipes (add_scalar_residual_equation) are - // a distinct contract from the rate form below: the state is defined by an - // implicit R(x, ε)=0 solved by backward_euler, not a rate integrated by - // rk_integrator. The Mode-B emission path for those lands separately; until - // then reject with a message that names the real reason — NOT the rate / - // rk_integrator / vector-solver message below, which would misdiagnose a - // residual recipe (it has one state variable but zero evolution equations, so - // it would otherwise trip the "exactly one evolution equation" check and blame - // the wrong contract). - if (!model.residual_equations().empty()) { - throw std::runtime_error( - "NumSimMaterialTarget: implicit residual equations " - "(add_scalar_residual_equation, strain-coupled state) are not yet " - "emitted — the graph-coupled Mode-B residual path (material_ref<" - "backward_euler> + solve) is a follow-up. This recipe declares a " - "residual, not a rate/rk_integrator evolution equation."); - } auto const svs = model.state_variables(); auto const eqs = model.evolution_equations(); if (svs.size() != 1 || eqs.size() != 1) { @@ -142,10 +135,445 @@ void check_scope(ConstitutiveModel const &model) { } } +// ── Mode-B strain-coupled residual emission ──────────────────────────────── +// A recipe defined by an implicit residual R(x, ε)=0 (add_scalar_residual_ +// equation) is lowered to a numsim-materials material that holds a +// material_ref and drives the Newton loop itself in ONE +// compute() (bound to the stress output — the always-pulled driver): +// +// compute() { +// auto eval = [&](value_type dz) { // trial x = x_old + dz +// value_type x = m_x.old_value() + dz; +// return std::pair{ R(x, ε), dR/dx(x, ε) }; // both t2s → scalar +// }; +// value_type dz = m_solver.get().solve(eval); +// m_x.new_value() = m_x.old_value() + dz; +// m_out_stress = σ(x, ε); ... // every output, post-solve +// } +// +// This is the robust Mode B (small_strain_plasticity's pattern): one compute() +// owns solve + state-apply + every output, so there is NO Local-edge ordering +// fragility (contrast a Mode-A material exposing residual/jacobian PROPERTIES, +// where the solve-vs-apply order rides on hash order). The consistent tangent +// dσ/dε is a follow-up (PR 2b) — tangents are rejected here. +// +// ⚠ backward_euler::solve() CLAMPS its result with std::max(x, 0) — a +// plasticity convention (a plastic-multiplier increment is physically ≥0). A +// generated residual material therefore correctly models only states whose +// solved INCREMENT is non-negative. This is surfaced in the emitted header; +// lifting it needs an unclamped solver mode upstream (numsim-materials). +std::vector emit_residual_material(ConstitutiveModel const &model) { + // ── Scope validation (residual contract) ── + auto const reqs = model.residual_equations(); + auto const svs = model.state_variables(); + if (reqs.size() != 1 || svs.size() != 1) { + throw std::runtime_error( + "NumSimMaterialTarget: a strain-coupled residual material supports " + "exactly one scalar state variable defined by one residual equation. " + "Coupled multi-state return maps need the numsim-materials vector " + "solver (numsim-materials#12)."); + } + if (!model.evolution_equations().empty()) { + // Unreachable through the public API (a state carries a rate XOR a residual, + // enforced by the recipe), but a residual recipe must never also carry a + // rate — guard against a future API that could mix them. + throw std::runtime_error( + "NumSimMaterialTarget: a residual recipe must not also declare rate " + "(evolution) equations — the state is defined by the residual alone."); + } + if (!model.tangents().empty()) { + throw std::runtime_error( + "NumSimMaterialTarget: the strain-coupled consistent tangent dσ/dε is " + "not yet emitted for a residual material — it is a follow-up (PR 2b: " + "dσ/dε = ∂σ/∂ε + ∂σ/∂x·(−∂R/∂ε / ∂R/∂x)). Drop the algorithmic-tangent " + "request for now."); + } + for (auto const &in : model.inputs()) { + if (in.kind != SymbolDecl::Kind::Tensor) { + throw std::runtime_error( + "NumSimMaterialTarget: scalar input '" + in.name + + "' is not yet wired into a residual material — a follow-up (tensor " + "inputs like strain ARE supported)."); + } + } + if (model.outputs().empty()) { + throw std::runtime_error( + "NumSimMaterialTarget: a residual material needs at least one output " + "(e.g. stress) to anchor compute() — the output's pull drives the " + "Newton solve. A state-only solve has no consumer."); + } + + auto const &req = reqs[0]; + auto const &sv = svs[req.state_variable_idx]; + auto const &cur_name = model.symbols()[sv.current_symbol_idx].name; + + // Scalar parameters the residual / outputs may reference (skip the framework + // time step — a residual material is rate-independent in this increment). + std::vector params; + std::set param_names; + for (auto const &p : model.parameters()) { + if (p.is_time_step || p.kind != SymbolDecl::Kind::Scalar) continue; + params.push_back(p); + param_names.insert(p.name); + } + + std::vector tensor_inputs; + std::set tensor_input_names; + for (auto const &in : model.inputs()) { + tensor_inputs.push_back(in); + tensor_input_names.insert(in.name); + } + + // Reserved-name guard (residual path reserves `solver_source`). + auto is_reserved_residual = [](std::string const &n) { + return n == contract::solver_source_param; + }; + if (is_reserved_residual(cur_name)) { + throw std::runtime_error( + "NumSimMaterialTarget: state variable name '" + cur_name + + "' collides with an emitted member (solver_source); rename it."); + } + for (auto const &p : params) { + if (is_reserved_residual(p.name)) { + throw std::runtime_error( + "NumSimMaterialTarget: parameter name '" + p.name + + "' collides with an emitted member (solver_source); rename it."); + } + if (p.default_value.has_value() && !std::isfinite(*p.default_value)) { + throw std::runtime_error( + "NumSimMaterialTarget: parameter '" + p.name + + "' has a non-finite default (" + fmt(*p.default_value) + + "); cannot emit as a C++ literal or JSON number."); + } + } + + // Resolve the current-state scalar holder — the diff variable for ∂R/∂x. + cas::expression_holder cur_expr; + for (auto const &[name, h] : model.scalar_symbol_map()) { + if (name == cur_name) cur_expr = h; + } + if (!cur_expr.is_valid()) { + throw std::runtime_error( + "NumSimMaterialTarget: cannot resolve the current-state scalar handle " + "for '" + sv.name + "' — state/symbol vectors out of sync."); + } + + // Register scalar symbols (state→local bare name, params→m_) and tensor + // symbols (input→local bare name) into an emit context. Shared by the residual + // /jacobian render and reused (fresh contexts) per output. + auto register_all = [&](CodeGenContext &c) { + for (auto const &[name, expr] : model.scalar_symbol_map()) { + if (param_names.contains(name)) { + c.register_symbol_scalar(expr, "m_" + name); + } else { + c.register_symbol_scalar(expr, name); + } + } + for (auto const &[name, expr] : model.tensor_symbol_map()) { + c.register_symbol_tensor(expr, name); + } + }; + + // Residual R and jacobian ∂R/∂x rendered with one shared CSE block (both live + // inside the eval lambda, so the state local and any tensor-input local they + // reference are in scope there). + CodeGenContext rc; + CodeEmitPipeline rp(rc); + register_all(rc); + rc.reset(); + auto const res_rhs = rp.t2s().apply(req.residual); + auto const jac_rhs = rp.t2s().apply(cas::diff(req.residual, cur_expr)); + auto const eval_decls = rc.render_statements(" "); // inside the lambda + + // Outputs (stress + any others), each in a fresh context so the CSE temp + // counter restarts. All are computed AFTER the solve, inside compute(). + struct EmittedOutput { + std::string name, decls, rhs; + bool is_tensor = false; + std::size_t dim = 0, rank = 0; + }; + std::vector outputs; + for (auto const &o : model.outputs()) { + if (is_reserved_residual(o.name) || o.name == cur_name || + param_names.contains(o.name) || tensor_input_names.contains(o.name)) { + throw std::runtime_error( + "NumSimMaterialTarget: output name '" + o.name + + "' collides with an emitted member, the state, a parameter, or a " + "tensor input; rename it."); + } + CodeGenContext oc; + CodeEmitPipeline op(oc); + register_all(oc); + if (o.kind == OutputDecl::Kind::Scalar) { + auto const &oexpr = + std::get>(o.expr); + LeafCollector olc; + olc.collect_scalar(oexpr); + for (auto const &n : olc.scalar_names()) { + if (n != cur_name && !param_names.contains(n)) { + throw std::runtime_error( + "NumSimMaterialTarget: scalar output '" + o.name + + "' references '" + n + "', which is neither the state '" + + cur_name + "' nor a parameter."); + } + } + oc.reset(); + auto const orhs = op.scalar().apply(oexpr); + outputs.push_back({o.name, oc.render_statements(" "), orhs, false, 0, 0}); + } else { + auto const &oexpr = + std::get>(o.expr); + LeafCollector olc; + olc.collect_tensor(oexpr); + for (auto const &n : olc.scalar_names()) { + if (n != cur_name && !param_names.contains(n)) { + throw std::runtime_error( + "NumSimMaterialTarget: tensor output '" + o.name + + "' references scalar '" + n + "', which is neither the state '" + + cur_name + "' nor a parameter."); + } + } + for (auto const &n : olc.tensor_names()) { + if (!tensor_input_names.contains(n)) { + throw std::runtime_error( + "NumSimMaterialTarget: tensor output '" + o.name + + "' references tensor '" + n + + "', which is not a declared tensor input."); + } + } + oc.reset(); + auto const orhs = op.tensor().apply(oexpr); + outputs.push_back( + {o.name, oc.render_statements(" "), orhs, true, o.dim, o.rank}); + } + } + + // Emitted-member uniqueness guard (same hazard as the rate path: synthesized + // member names can collide with recipe symbols). + { + std::set member_bases; + auto claim = [&member_bases](std::string const &base) { + if (!member_bases.insert(base).second) { + throw std::runtime_error( + "NumSimMaterialTarget: emitted member 'm_" + base + + "' would be duplicated — a recipe symbol collides with a synthesized " + "member name; rename the offending state/parameter/input/output."); + } + }; + claim(contract::solver_source_param); + claim(cur_name); + for (auto const &p : params) claim(p.name); + for (auto const &o : outputs) claim("out_" + o.name); + for (auto const &ti : tensor_inputs) { + claim(ti.name); + claim(ti.name + "_source"); + } + } + + auto const &cls = model.name(); + auto const &driver = outputs.front().name; // compute() is bound to this output + + // ── Material header ── + std::ostringstream h; + h << "// Auto-generated by numsim-codegen (NumSimMaterialTarget). Do not " + "edit.\n"; + h << "// Strain-coupled residual material for recipe \"" << cls + << "\": state " << cur_name << " solves R(" << cur_name + << ", inputs) = 0\n"; + h << "// (implicit, Mode B). Holds a material_ref and drives " + "the\n"; + h << "// Newton loop in compute() (bound to the '" << driver + << "' output); the\n"; + h << "// converged state then feeds every output.\n"; + h << "//\n"; + h << "// CONTRACT: pull the '" << driver + << "' output to drive the solve. NOTE backward_euler::solve()\n"; + h << "// clamps its increment to >= 0 (a plasticity convention), so this " + "material\n"; + h << "// models only states whose solved increment is non-negative.\n"; + h << "#pragma once\n"; + h << "#include \"" << contract::material_base_include << "\"\n"; + h << "#include \"" << contract::backward_euler_include << "\"\n"; + h << "#include \"" << contract::material_ref_include << "\"\n"; + h << "#include \n"; + h << "#include \n\n"; + h << "namespace numsim::materials::generated {\n\n"; + h << "template \n"; + h << "class " << cls << " final\n"; + h << " : public numsim::materials::material_base<" << cls + << ", Traits> {\n"; + h << "public:\n"; + h << " using base = numsim::materials::material_base<" << cls + << ", Traits>;\n"; + h << " using value_type = typename base::value_type;\n"; + h << " using input_parameter_controller =\n" + " typename base::input_parameter_controller;\n"; + h << " using solver_type = numsim::materials::backward_euler;\n\n"; + + // Constructor. + h << " template \n"; + h << " explicit " << cls << "(Args&&... args)\n"; + h << " : base(std::forward(args)...),\n"; + for (std::size_t i = 0; i < outputs.size(); ++i) { + auto const &o = outputs[i]; + std::string const ty = + o.is_tensor ? tensor_cxx_type(o.dim, o.rank) : "value_type"; + h << " m_out_" << o.name << "(base::template add_output<" << ty + << ">(\n"; + // The first output carries &compute (the solve driver); others are set as a + // side effect of compute() and have no callback of their own. + if (i == 0) { + h << " \"" << o.name << "\", &" << cls << "::compute)),\n"; + } else { + h << " \"" << o.name << "\")),\n"; + } + } + h << " m_" << cur_name + << "(base::template add_history_output(\"" << cur_name + << "\")),\n"; + for (auto const &p : params) { + h << " m_" << p.name + << "(base::template get_parameter(\"" << p.name << "\")),\n"; + } + h << " m_solver(base::template add_material_ref(\n"; + h << " base::template get_parameter(\"" + << contract::solver_source_param << "\")))"; + if (tensor_inputs.empty()) { + h << " {}\n\n"; + } else { + h << ",\n"; + for (std::size_t i = 0; i < tensor_inputs.size(); ++i) { + auto const &ti = tensor_inputs[i]; + h << " m_" << ti.name << "(base::template add_input<" + << tensor_cxx_type(ti.dim, ti.rank) << ">(\n"; + h << " base::template get_parameter(\"" << ti.name + << "_source\"),\n"; + h << " \"" << ti.name + << "\", numsim::materials::EdgeKind::Global))" + << (i + 1 < tensor_inputs.size() ? ",\n" : " {}\n\n"); + } + } + + // parameters() schema + h << " static input_parameter_controller parameters() {\n"; + h << " input_parameter_controller para{base::parameters()};\n"; + for (auto const &p : params) { + if (p.default_value.has_value()) { + h << " para.template insert(\"" << p.name + << "\").template add(value_type{" + << fmt(*p.default_value) << "});\n"; + } else { + h << " para.template insert(\"" << p.name + << "\").template add();\n"; + } + } + h << " para.template insert(\"" + << contract::solver_source_param << "\")\n"; + h << " .template add();\n"; + for (auto const &ti : tensor_inputs) { + h << " para.template insert(\"" << ti.name << "_source\")\n"; + h << " .template add();\n"; + } + h << " return para;\n"; + h << " }\n\n"; + + // compute() + h << " // Solves R(" << cur_name << ", inputs)=0 for the increment via\n"; + h << " // backward_euler::solve(), applies the state, then computes every " + "output.\n"; + h << " void compute() {\n"; + for (auto const &ti : tensor_inputs) { + h << " [[maybe_unused]] const auto& " << ti.name << " = m_" << ti.name + << ".get();\n"; + } + h << " auto eval = [&](value_type d" << cur_name + << ") -> std::pair {\n"; + h << " const value_type " << cur_name << " = m_" << cur_name + << ".old_value() + d" << cur_name << ";\n"; + if (!eval_decls.empty()) h << eval_decls; + h << " const value_type residual = " << res_rhs << ";\n"; + h << " const value_type jacobian = " << jac_rhs << ";\n"; + h << " return {residual, jacobian};\n"; + h << " };\n"; + h << " const value_type d" << cur_name + << " = m_solver.get().solve(eval);\n"; + h << " m_" << cur_name << ".new_value() = m_" << cur_name + << ".old_value() + d" << cur_name << ";\n"; + h << " [[maybe_unused]] const value_type " << cur_name << " = m_" + << cur_name << ".new_value();\n"; + for (auto const &o : outputs) { + if (!o.decls.empty()) { + // The output decls were rendered at 4-space indent; they sit directly in + // compute()'s body, which is also 4-space — emit as-is. + h << o.decls; + } + h << " m_out_" << o.name << " = " << o.rhs << ";\n"; + } + h << " }\n\n"; + + // members + h << "private:\n"; + for (auto const &o : outputs) { + std::string const ty = + o.is_tensor ? tensor_cxx_type(o.dim, o.rank) : "value_type"; + h << " " << ty << "& m_out_" << o.name << ";\n"; + } + h << " numsim_core::history_property& m_" << cur_name << ";\n"; + for (auto const &p : params) { + h << " [[maybe_unused]] const value_type& m_" << p.name << ";\n"; + } + h << " numsim::materials::material_ref& m_solver;\n"; + for (auto const &ti : tensor_inputs) { + h << " const numsim::materials::input_property<" + << tensor_cxx_type(ti.dim, ti.rank) << ",\n"; + h << " numsim::materials::property_traits>& m_" << ti.name << ";\n"; + } + h << "};\n\n"; + h << "} // namespace numsim::materials::generated\n"; + + // ── JSON config scaffold ── + std::ostringstream j; + j << "{\n"; + j << " \"_comment\": \"SCAFFOLD generated by numsim-codegen. The " + "backward_euler solver is created WITHOUT a 'function' (caller-driven " + "mode — the material drives the Newton loop). The '*_source' values are " + "PLACEHOLDER producer names; a material publishing a property of the " + "matching name/type must exist in the graph. The 'type' names must be " + "registered in the material factory.\",\n"; + j << " \"materials\": [\n"; + j << " { \"name\": \"solver\", \"type\": \"backward_euler\",\n"; + j << " \"tolerance\": 1e-12, \"max_iter\": 50 },\n"; + j << " { \"name\": \"" << cls << "\", \"type\": \"" << cls << "\",\n"; + j << " \"" << contract::solver_source_param << "\": \"solver\""; + for (auto const &ti : tensor_inputs) { + j << ",\n \"" << ti.name << "_source\": \"" << ti.name << "_producer\""; + } + for (auto const &p : params) { + if (p.default_value.has_value()) { + j << ",\n \"" << p.name << "\": " << fmt(*p.default_value); + } + } + j << " }\n"; + j << " ]\n"; + j << "}\n"; + + return { + EmittedFile{cls + ".h", h.str(), "include/materials", + EmittedFile::Kind::Header}, + EmittedFile{cls + ".config.json", j.str(), "", EmittedFile::Kind::Other}, + }; +} + } // namespace auto NumSimMaterialTarget::emit(ConstitutiveModel const &model) const -> std::vector { + // Strain-coupled implicit-residual recipes take the Mode-B path; the rate + // (rk_integrator) path below handles the rest. Routing here (before + // check_scope) keeps each path's scope validation self-contained. + if (!model.residual_equations().empty()) { + return emit_residual_material(model); + } check_scope(model); auto const &sv = model.state_variables()[0]; diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index 6b79252..d8c5225 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -194,14 +194,16 @@ if(_ncg_materials_inc AND _ncg_core_inc) ${CMAKE_CURRENT_BINARY_DIR}/generated/NonlinearDecay.h) set(GENERATED_VISCO_HEADER ${CMAKE_CURRENT_BINARY_DIR}/generated/Viscoelastic.h) + set(GENERATED_RETURNMAP_HEADER + ${CMAKE_CURRENT_BINARY_DIR}/generated/ReturnMap.h) add_custom_command( OUTPUT ${GENERATED_LINEAR_HEADER} ${GENERATED_NONLINEAR_HEADER} - ${GENERATED_VISCO_HEADER} + ${GENERATED_VISCO_HEADER} ${GENERATED_RETURNMAP_HEADER} COMMAND ${CMAKE_COMMAND} -E make_directory ${CMAKE_CURRENT_BINARY_DIR}/generated COMMAND $ ${GENERATED_LINEAR_HEADER} ${GENERATED_NONLINEAR_HEADER} - ${GENERATED_VISCO_HEADER} + ${GENERATED_VISCO_HEADER} ${GENERATED_RETURNMAP_HEADER} DEPENDS generate_numsim_material_check COMMENT "Generating numsim-materials rate materials via NumSimMaterialTarget" VERBATIM) @@ -209,7 +211,7 @@ if(_ncg_materials_inc AND _ncg_core_inc) add_executable(numsim_material_check_driver generated/numsim_material_check_driver.cpp ${GENERATED_LINEAR_HEADER} ${GENERATED_NONLINEAR_HEADER} - ${GENERATED_VISCO_HEADER}) + ${GENERATED_VISCO_HEADER} ${GENERATED_RETURNMAP_HEADER}) # numsim-materials, numsim-core, Eigen and tmech are header-only and included # as SYSTEM so their warnings stay silent under the first-party -Werror gate. # tmech (via cas's CPM) is needed for the tensor-stress material/test. diff --git a/tests/NumSimMaterialTargetTest.cpp b/tests/NumSimMaterialTargetTest.cpp index 67f9559..4555ef0 100644 --- a/tests/NumSimMaterialTargetTest.cpp +++ b/tests/NumSimMaterialTargetTest.cpp @@ -331,24 +331,81 @@ TEST(NumSimMaterialTarget, FloatDefaultRoundTrips) { << src; } -// Finding A (holistic review 2026-06-17): a recipe defined by an implicit -// residual (add_scalar_residual_equation) has no Mode-B emission path yet. It -// must be rejected with a message that names the RESIDUAL contract — not the -// rate/rk_integrator/vector-solver message, which would misdiagnose it (a -// residual recipe has one state variable but zero evolution equations). -TEST(NumSimMaterialTarget, RejectsResidualRecipeWithAccurateMessage) { +// ─── Phase 2a: Mode-B strain-coupled residual emission ─────────────────────── + +// A return-map recipe R(z, ε) = z − c·tr(ε), σ = z·ε. The state z is solved +// implicitly by backward_euler; the material drives the Newton loop itself. +auto build_return_map() -> ConstitutiveModel { 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", make_expression(0.0)); m.add_scalar_residual_equation(z, z.current - c * trace(eps)); + m.add_output("stress", z.current * eps); + return m; +} + +// The emitted residual material conforms to the Mode-B (backward_euler caller- +// driven) contract: a material_ref, a solve(eval) call, the +// residual and its jacobian inside the eval lambda, an owned history state, and +// the stress output bound to compute(). +TEST(NumSimMaterialTarget, EmitsModeBResidualMaterial) { + auto const h = header_of(NumSimMaterialTarget{}.emit(build_return_map())); + // Mode-B structural surface. + EXPECT_NE(h.find("using solver_type = numsim::materials::backward_euler"), + std::string::npos) << h; + EXPECT_NE(h.find("add_material_ref"), std::string::npos) << h; + EXPECT_NE(h.find("m_solver.get().solve(eval)"), std::string::npos) << h; + EXPECT_NE(h.find("add_history_output(\"z\")"), std::string::npos) + << h; + // The stress output drives the solve (carries &compute). + EXPECT_NE(h.find("\"stress\", &ReturnMap::compute"), std::string::npos) << h; + // The eval lambda returns {residual, jacobian}; jacobian ∂R/∂z = 1. + EXPECT_NE(h.find("return {residual, jacobian};"), std::string::npos) << h; + EXPECT_NE(h.find("const value_type jacobian = 1"), std::string::npos) << h; + // State applied as old + increment. + EXPECT_NE(h.find("m_z.new_value() = m_z.old_value() + dz"), std::string::npos) + << h; + // The solver_source param is required. + EXPECT_NE(h.find("insert(\"solver_source\")"), std::string::npos) + << h; +} + +// The strain-coupled consistent tangent is a follow-up (PR 2b): an algorithmic +// tangent on a residual material must be rejected with a message naming PR 2b, +// not silently dropped. +TEST(NumSimMaterialTarget, RejectsTangentOnResidualMaterial) { + auto m = build_return_map(); + m.add_algorithmic_tangent("dstress_dstrain", "stress", "strain"); auto const msg = emit_throw_message(m); - // Names the residual contract... - EXPECT_NE(msg.find("residual"), std::string::npos) << msg; - // ...and does NOT misdiagnose as the rate/evolution-equation contract. - EXPECT_EQ(msg.find("exactly one scalar state variable"), std::string::npos) - << msg; + EXPECT_NE(msg.find("consistent tangent"), std::string::npos) << msg; +} + +// A residual material needs at least one output to anchor compute() (the output +// pull drives the solve) — reject loudly rather than emit an un-driven solve. +TEST(NumSimMaterialTarget, RejectsResidualWithoutOutput) { + ConstitutiveModel m("NoOutput"); + 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", make_expression(0.0)); + m.add_scalar_residual_equation(z, z.current - c * trace(eps)); + EXPECT_NE(emit_throw_message(m).find("at least one output"), + std::string::npos); +} + +// Scalar inputs into a residual material are a follow-up — rejected for now. +TEST(NumSimMaterialTarget, RejectsScalarInputOnResidualMaterial) { + ConstitutiveModel m("ScalarIn"); + auto c = m.add_parameter("c", 2.0); + auto eps = m.add_tensor_input("strain", 3, 2, roles::Strain); + auto temp = m.add_scalar_input("temperature"); + auto z = + m.add_scalar_state_variable("z", make_expression(0.0)); + m.add_scalar_residual_equation(z, z.current - c * trace(eps) - temp); + m.add_output("stress", z.current * eps); + EXPECT_NE(emit_throw_message(m).find("scalar input"), std::string::npos); } } // namespace diff --git a/tests/generated/generate_numsim_material_check.cpp b/tests/generated/generate_numsim_material_check.cpp index 3825432..d91b858 100644 --- a/tests/generated/generate_numsim_material_check.cpp +++ b/tests/generated/generate_numsim_material_check.cpp @@ -20,6 +20,8 @@ #include #include #include +#include +#include #include #include @@ -51,9 +53,9 @@ bool write_header(ConstitutiveModel const& model, char const* path) { } // namespace int main(int argc, char** argv) { - if (argc < 4) { + if (argc < 5) { std::cerr << "usage: generate_numsim_material_check " - " \n"; + " \n"; return 2; } @@ -95,8 +97,24 @@ int main(int argc, char** argv) { viscoelastic.add_algorithmic_tangent("dstress_dstrain", "stress", "strain"); } + // Phase 2a: strain-coupled IMPLICIT residual material (Mode B). The state z is + // defined by R(z, ε) = z − c·tr(ε) = 0, solved by backward_euler (caller- + // driven); the stress is σ = z·ε. Unlike the rate materials above, z is solved + // implicitly inside the material's compute() — no rk_integrator. Exercises the + // material_ref + solve(eval) emission end-to-end. + ConstitutiveModel returnmap("ReturnMap"); + { + auto c = returnmap.add_parameter("c", 2.0); + auto eps = returnmap.add_tensor_input("strain", 3, 2, roles::Strain); + auto z = returnmap.add_scalar_state_variable( + "z", make_expression(0.0)); + returnmap.add_scalar_residual_equation(z, z.current - c * trace(eps)); + returnmap.add_output("stress", z.current * eps); + } + if (!write_header(linear, argv[1])) return 1; if (!write_header(nonlinear, argv[2])) return 1; if (!write_header(viscoelastic, argv[3])) return 1; + if (!write_header(returnmap, argv[4])) return 1; return 0; } diff --git a/tests/generated/numsim_material_check_driver.cpp b/tests/generated/numsim_material_check_driver.cpp index b37e435..d53f381 100644 --- a/tests/generated/numsim_material_check_driver.cpp +++ b/tests/generated/numsim_material_check_driver.cpp @@ -33,8 +33,10 @@ #if defined(__GNUC__) && !defined(__clang__) #define NCG_TENSOR_E2E 1 #include +#include "numsim-materials/solvers/backward_euler.h" #include "numsim-materials/materials/tensor_component_stepper.h" #include "Viscoelastic.h" // generated: σ = α·ε (tensor stress from scalar state) +#include "ReturnMap.h" // generated: implicit residual R(z,ε)=z−c·tr(ε), σ=z·ε #endif #include @@ -51,6 +53,7 @@ using Linear = numsim::materials::generated::LinearHardening; using Nonlinear = numsim::materials::generated::NonlinearDecay; #ifdef NCG_TENSOR_E2E using Visco = numsim::materials::generated::Viscoelastic; +using ReturnMap = numsim::materials::generated::ReturnMap; using tensor2 = tmech::tensor; #endif @@ -273,6 +276,50 @@ TEST(NumSimMaterialEndToEnd, TensorStressFromScalarStateAndStrain) { EXPECT_NEAR(C(0, 0, 1, 1), T{0}, 1e-12); // off-block zero } +// ── Strain-coupled implicit residual (Mode B) ──────────────────────────────── +// Proof the residual emit (material_ref + solve(eval)) compiles +// and runs against the REAL backward_euler in caller-driven mode: drive a strain +// producer + the generated ReturnMap; the material solves R(z,ε)=z−c·tr(ε)=0 +// internally, so z = c·tr(ε) and σ = z·ε. +TEST(NumSimMaterialEndToEnd, ResidualReturnMapSolvesAgainstBackwardEuler) { + ctx_type ctx; + param_type p; + + // strain producer: accumulates `increment` into component (0,0) per update. + p.insert("name", "stepper"); + p.insert("increment", T{0.02}); + p.insert>("indices", {0, 0}); + ctx.create>(p); + + // caller-driven backward_euler (no "function" → the material drives solve()). + p.clear(); + p.insert("name", "solver"); + p.insert("tolerance", T{1e-12}); + p.insert("max_iter", 50); + ctx.create>(p); + + p.clear(); + const T c = T{2.0}; + p.insert("name", "ReturnMap"); + p.insert("c", c); + p.insert("solver_source", "solver"); + p.insert("strain_source", "stepper"); + ctx.create(p); + + ctx.finalize(); + ctx.update(); + ctx.commit(); + + const auto eps = read_tensor(ctx, "stepper", "strain"); + const T tr = eps(0, 0) + eps(1, 1) + eps(2, 2); + const T z = ctx.get("ReturnMap", "z"); // solved scalar state + const auto sig = read_tensor(ctx, "ReturnMap", "stress"); + + EXPECT_NEAR(z, c * tr, 1e-10); // R=0 ⇒ z = c·tr(ε) + EXPECT_NEAR(sig(0, 0), z * eps(0, 0), 1e-10); // σ = z·ε + EXPECT_NEAR(sig(0, 1), T{0}, 1e-12); // off-diagonal strain is 0 +} + #endif // NCG_TENSOR_E2E } // namespace From df9fb8e347031da588f703bb0c8c1c463bb8d758 Mon Sep 17 00:00:00 2001 From: petlenz Date: Thu, 18 Jun 2026 08:27:30 +0200 Subject: [PATCH 3/3] Phase 2a review: fix multi-output CSE collision, residual leaf/reserved-name guards, jacobian-validating e2e --- src/targets/numsim_material.cpp | 103 +++++++++++++----- tests/CMakeLists.txt | 7 +- tests/NumSimMaterialTargetTest.cpp | 85 ++++++++++++++- .../generate_numsim_material_check.cpp | 24 +++- .../numsim_material_check_driver.cpp | 50 +++++++++ 5 files changed, 239 insertions(+), 30 deletions(-) diff --git a/src/targets/numsim_material.cpp b/src/targets/numsim_material.cpp index e3e68ca..64a5843 100644 --- a/src/targets/numsim_material.cpp +++ b/src/targets/numsim_material.cpp @@ -224,21 +224,29 @@ std::vector emit_residual_material(ConstitutiveModel const &model) tensor_input_names.insert(in.name); } - // Reserved-name guard (residual path reserves `solver_source`). + // Reserved-name guard. A residual material emits the fixed member `m_solver` + // (the material_ref), the schema key `solver_source`, and the compute()-local + // identifiers `eval`, `residual`, `jacobian`. A recipe symbol that renders to + // one of these bare names would produce a duplicate member / redeclared local + // (uncompilable). Applied to the state, parameters AND tensor inputs — all + // three render to bare identifiers in compute() (params are m_-prefixed in + // expressions, but `solver_source`/`solver` would still collide with the + // member/schema, so params are checked too). auto is_reserved_residual = [](std::string const &n) { - return n == contract::solver_source_param; + return n == contract::solver_source_param || n == "solver" || n == "eval" || + n == "residual" || n == "jacobian"; }; - if (is_reserved_residual(cur_name)) { - throw std::runtime_error( - "NumSimMaterialTarget: state variable name '" + cur_name + - "' collides with an emitted member (solver_source); rename it."); - } - for (auto const &p : params) { - if (is_reserved_residual(p.name)) { + auto reject_reserved = [&](std::string const &n, char const *what) { + if (is_reserved_residual(n)) { throw std::runtime_error( - "NumSimMaterialTarget: parameter name '" + p.name + - "' collides with an emitted member (solver_source); rename it."); + std::string("NumSimMaterialTarget: ") + what + " '" + n + + "' collides with an emitted member / local " + "(solver/solver_source/eval/residual/jacobian); rename it."); } + }; + reject_reserved(cur_name, "state variable name"); + for (auto const &p : params) { + reject_reserved(p.name, "parameter name"); if (p.default_value.has_value() && !std::isfinite(*p.default_value)) { throw std::runtime_error( "NumSimMaterialTarget: parameter '" + p.name + @@ -246,6 +254,18 @@ std::vector emit_residual_material(ConstitutiveModel const &model) "); cannot emit as a C++ literal or JSON number."); } } + // The Newton increment is a compute()-local named `d`; a tensor input + // (also a bare compute() local) named the same would redeclare it. + std::string const increment_local = "d" + cur_name; + for (auto const &in : tensor_inputs) { + reject_reserved(in.name, "tensor input name"); + if (in.name == increment_local) { + throw std::runtime_error( + "NumSimMaterialTarget: tensor input '" + in.name + + "' collides with the synthesized Newton-increment local '" + + increment_local + "'; rename the input or the state."); + } + } // Resolve the current-state scalar holder — the diff variable for ∂R/∂x. cas::expression_holder cur_expr; @@ -274,6 +294,35 @@ std::vector emit_residual_material(ConstitutiveModel const &model) } }; + // Residual-leaf guard: every scalar leaf of R must be the current state or a + // parameter, and every tensor leaf must be a declared tensor input. The recipe + // checks all leaves are *declared* symbols, but NOT that they map to a bound + // emit name — e.g. the previous-step state `_old` is a declared symbol + // yet has no local in compute() (the increment form x = x_old + dz is the + // emitter's own concern). Without this guard such a residual would emit an + // unbound identifier inside the eval lambda. Mirrors the rate-leaf guard. + { + LeafCollector rlc; + rlc.collect_t2s(req.residual); + for (auto const &n : rlc.scalar_names()) { + if (n != cur_name && !param_names.contains(n)) { + throw std::runtime_error( + "NumSimMaterialTarget: the residual references scalar '" + n + + "', which is neither the state '" + cur_name + + "' nor a parameter. A residual R(" + cur_name + + ", inputs) cannot depend on the previous-step state or a scalar " + "input on this path."); + } + } + for (auto const &n : rlc.tensor_names()) { + if (!tensor_input_names.contains(n)) { + throw std::runtime_error( + "NumSimMaterialTarget: the residual references tensor '" + n + + "', which is not a declared tensor input."); + } + } + } + // Residual R and jacobian ∂R/∂x rendered with one shared CSE block (both live // inside the eval lambda, so the state local and any tensor-input local they // reference are in scope there). @@ -360,7 +409,7 @@ std::vector emit_residual_material(ConstitutiveModel const &model) "member name; rename the offending state/parameter/input/output."); } }; - claim(contract::solver_source_param); + claim("solver"); // the fixed material_ref member m_solver claim(cur_name); for (auto const &p : params) claim(p.name); for (auto const &o : outputs) claim("out_" + o.name); @@ -420,13 +469,13 @@ std::vector emit_residual_material(ConstitutiveModel const &model) o.is_tensor ? tensor_cxx_type(o.dim, o.rank) : "value_type"; h << " m_out_" << o.name << "(base::template add_output<" << ty << ">(\n"; - // The first output carries &compute (the solve driver); others are set as a - // side effect of compute() and have no callback of their own. - if (i == 0) { - h << " \"" << o.name << "\", &" << cls << "::compute)),\n"; - } else { - h << " \"" << o.name << "\")),\n"; - } + // EVERY output carries &compute: pulling ANY output drives the one compute() + // that solves the state and sets all outputs. (compute() is idempotent + // within a graph evaluation — it re-solves to the same root — so a consumer + // pulling several outputs is correct, if mildly redundant.) Binding only the + // first output would leave the others stale when pulled in isolation. + h << " \"" << o.name << "\", &" << cls << "::compute)),\n"; + (void)i; } h << " m_" << cur_name << "(base::template add_history_output(\"" << cur_name @@ -502,12 +551,16 @@ std::vector emit_residual_material(ConstitutiveModel const &model) h << " [[maybe_unused]] const value_type " << cur_name << " = m_" << cur_name << ".new_value();\n"; for (auto const &o : outputs) { - if (!o.decls.empty()) { - // The output decls were rendered at 4-space indent; they sit directly in - // compute()'s body, which is also 4-space — emit as-is. - h << o.decls; - } - h << " m_out_" << o.name << " = " << o.rhs << ";\n"; + // Each output is rendered in its OWN CodeGenContext, so its CSE temporaries + // restart at t0 — and ALL outputs share this single compute() body. Without + // a nested scope, a second output with CSE temps would redeclare `t0` (a + // compile error the single-output test recipe never hit). Brace-scope every + // output so each block's temps are private. (A decl-free output gets a + // harmless empty-ish block.) + h << " {\n"; + if (!o.decls.empty()) h << o.decls; // rendered at 4-space; nested is fine + h << " m_out_" << o.name << " = " << o.rhs << ";\n"; + h << " }\n"; } h << " }\n\n"; diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index d8c5225..cfb8028 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -196,14 +196,18 @@ if(_ncg_materials_inc AND _ncg_core_inc) ${CMAKE_CURRENT_BINARY_DIR}/generated/Viscoelastic.h) set(GENERATED_RETURNMAP_HEADER ${CMAKE_CURRENT_BINARY_DIR}/generated/ReturnMap.h) + set(GENERATED_RETURNMAP_CUBIC_HEADER + ${CMAKE_CURRENT_BINARY_DIR}/generated/ReturnMapCubic.h) add_custom_command( OUTPUT ${GENERATED_LINEAR_HEADER} ${GENERATED_NONLINEAR_HEADER} ${GENERATED_VISCO_HEADER} ${GENERATED_RETURNMAP_HEADER} + ${GENERATED_RETURNMAP_CUBIC_HEADER} COMMAND ${CMAKE_COMMAND} -E make_directory ${CMAKE_CURRENT_BINARY_DIR}/generated COMMAND $ ${GENERATED_LINEAR_HEADER} ${GENERATED_NONLINEAR_HEADER} ${GENERATED_VISCO_HEADER} ${GENERATED_RETURNMAP_HEADER} + ${GENERATED_RETURNMAP_CUBIC_HEADER} DEPENDS generate_numsim_material_check COMMENT "Generating numsim-materials rate materials via NumSimMaterialTarget" VERBATIM) @@ -211,7 +215,8 @@ if(_ncg_materials_inc AND _ncg_core_inc) add_executable(numsim_material_check_driver generated/numsim_material_check_driver.cpp ${GENERATED_LINEAR_HEADER} ${GENERATED_NONLINEAR_HEADER} - ${GENERATED_VISCO_HEADER} ${GENERATED_RETURNMAP_HEADER}) + ${GENERATED_VISCO_HEADER} ${GENERATED_RETURNMAP_HEADER} + ${GENERATED_RETURNMAP_CUBIC_HEADER}) # numsim-materials, numsim-core, Eigen and tmech are header-only and included # as SYSTEM so their warnings stay silent under the first-party -Werror gate. # tmech (via cas's CPM) is needed for the tensor-stress material/test. diff --git a/tests/NumSimMaterialTargetTest.cpp b/tests/NumSimMaterialTargetTest.cpp index 4555ef0..85f0b85 100644 --- a/tests/NumSimMaterialTargetTest.cpp +++ b/tests/NumSimMaterialTargetTest.cpp @@ -361,9 +361,17 @@ TEST(NumSimMaterialTarget, EmitsModeBResidualMaterial) { << h; // The stress output drives the solve (carries &compute). EXPECT_NE(h.find("\"stress\", &ReturnMap::compute"), std::string::npos) << h; - // The eval lambda returns {residual, jacobian}; jacobian ∂R/∂z = 1. + // The eval lambda returns {residual, jacobian}. EXPECT_NE(h.find("return {residual, jacobian};"), std::string::npos) << h; - EXPECT_NE(h.find("const value_type jacobian = 1"), std::string::npos) << h; + // Pin the RESIDUAL RHS, not just the jacobian — the load-bearing output. A + // dropped coupling term or wrong sign must fail here, at the always-run unit + // layer (the tensor e2e is gcc-only). R = z − c·tr(ε): the rendered scalar + // must reference the state z and the strain trace. + EXPECT_NE(h.find("const value_type residual = "), std::string::npos) << h; + EXPECT_NE(h.find("tmech::trace(strain)"), std::string::npos) << h; + // ∂R/∂z = 1 exactly, rendered "1.0" (terminating ';' so "= 1.5" etc can't + // false-match). + EXPECT_NE(h.find("const value_type jacobian = 1.0;"), std::string::npos) << h; // State applied as old + increment. EXPECT_NE(h.find("m_z.new_value() = m_z.old_value() + dz"), std::string::npos) << h; @@ -372,6 +380,79 @@ TEST(NumSimMaterialTarget, EmitsModeBResidualMaterial) { << h; } +// Review (cpp-pro): two outputs share the single compute() body. Each is +// rendered in its own CSE context (temps restart at t0), so without a nested +// scope the second output redeclares `t0` — an uncompilable header the single- +// output recipe never exercised. Each output must be brace-scoped, and EVERY +// output must drive the solve (carry &compute), not just the first. +TEST(NumSimMaterialTarget, MultiOutputResidualScopesCseAndDrivesAll) { + ConstitutiveModel m("MultiOut"); + 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", make_expression(0.0)); + m.add_scalar_residual_equation(z, z.current - c * trace(eps)); + m.add_output("stress", sin(z.current) * eps); // CSE temp + m.add_output("stress2", cos(z.current) * eps); // CSE temp + auto const h = header_of(NumSimMaterialTarget{}.emit(m)); + // Both outputs drive the solve. + EXPECT_NE(h.find("\"stress\", &MultiOut::compute"), std::string::npos) << h; + EXPECT_NE(h.find("\"stress2\", &MultiOut::compute"), std::string::npos) << h; + // The two CSE blocks are brace-scoped (each output writes inside a `{ }`), + // so `t0` is private per block — count the opening braces of output blocks. + std::size_t braces = 0; + for (std::size_t pos = 0; + (pos = h.find("\n {\n", pos)) != std::string::npos; ++pos) + ++braces; + EXPECT_GE(braces, 2u) + << "each output must be brace-scoped to avoid CSE temp collision:\n" + << h; +} + +// Review (cpp-pro): a residual must not reference the previous-step state +// `_old` — it is a declared symbol (so the recipe accepts it) but has no +// local in the emitted compute(), which would emit an unbound identifier. Guard +// at emit time. +TEST(NumSimMaterialTarget, RejectsResidualReferencingOldState) { + ConstitutiveModel m("UsesOld"); + 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", make_expression(0.0)); + // R = z − z_old − c·tr(ε): references the previous-step state. + m.add_scalar_residual_equation(z, z.current - z.previous - c * trace(eps)); + m.add_output("stress", z.current * eps); + EXPECT_NE(emit_throw_message(m).find("previous-step state"), + std::string::npos); +} + +// Review (cpp-pro / code-reviewer): the emitter synthesizes the fixed member +// `m_solver` and the compute() locals `eval`/`residual`/`jacobian`. A recipe +// symbol named like one of these must be rejected with a rename message, not +// surface as a downstream redefinition. +TEST(NumSimMaterialTarget, RejectsSolverNameCollisionOnResidualMaterial) { + ConstitutiveModel m("SolverClash"); + auto solver = m.add_parameter("solver", 2.0); // collides with m_solver + auto eps = m.add_tensor_input("strain", 3, 2, roles::Strain); + auto z = + m.add_scalar_state_variable("z", make_expression(0.0)); + m.add_scalar_residual_equation(z, z.current - solver * trace(eps)); + m.add_output("stress", z.current * eps); + EXPECT_NE(emit_throw_message(m).find("rename"), std::string::npos); +} + +// ...and a state named like a compute() local (`residual`) is rejected too. +TEST(NumSimMaterialTarget, RejectsResidualLocalNameCollision) { + ConstitutiveModel m("LocalClash"); + auto c = m.add_parameter("c", 2.0); + auto eps = m.add_tensor_input("strain", 3, 2, roles::Strain); + auto residual = m.add_scalar_state_variable( + "residual", make_expression(0.0)); // collides + m.add_scalar_residual_equation(residual, residual.current - c * trace(eps)); + m.add_output("stress", residual.current * eps); + EXPECT_NE(emit_throw_message(m).find("rename"), std::string::npos); +} + // The strain-coupled consistent tangent is a follow-up (PR 2b): an algorithmic // tangent on a residual material must be rejected with a message naming PR 2b, // not silently dropped. diff --git a/tests/generated/generate_numsim_material_check.cpp b/tests/generated/generate_numsim_material_check.cpp index d91b858..dffbdae 100644 --- a/tests/generated/generate_numsim_material_check.cpp +++ b/tests/generated/generate_numsim_material_check.cpp @@ -53,9 +53,10 @@ bool write_header(ConstitutiveModel const& model, char const* path) { } // namespace int main(int argc, char** argv) { - if (argc < 5) { + if (argc < 6) { std::cerr << "usage: generate_numsim_material_check " - " \n"; + " " + "\n"; return 2; } @@ -112,9 +113,28 @@ int main(int argc, char** argv) { returnmap.add_output("stress", z.current * eps); } + // NONLINEAR residual to exercise the t2s-wrt-scalar jacobian (∂R/∂z, cas#285). + // R(z, ε) = z + z³ − c·tr(ε), so ∂R/∂z = 1 + 3z² is NON-constant. The linear + // ReturnMap above has ∂R/∂z ≡ 1, so a wrong jacobian would still converge to + // the right root — it cannot validate the emitted derivative. Here, with a + // tight Newton budget and c·tr(ε)=1 (root z≈0.682, where ∂R/∂z≈2.4 ≫ 1), an + // incorrect jacobian (e.g. a constant) oscillates and FAILS to converge — so + // the e2e value check actually pins the emitted derivative. + ConstitutiveModel returnmap_cubic("ReturnMapCubic"); + { + auto c = returnmap_cubic.add_parameter("c", 2.0); + auto eps = returnmap_cubic.add_tensor_input("strain", 3, 2, roles::Strain); + auto z = returnmap_cubic.add_scalar_state_variable( + "z", make_expression(0.0)); + returnmap_cubic.add_scalar_residual_equation( + z, z.current + z.current * z.current * z.current - c * trace(eps)); + returnmap_cubic.add_output("stress", z.current * eps); + } + if (!write_header(linear, argv[1])) return 1; if (!write_header(nonlinear, argv[2])) return 1; if (!write_header(viscoelastic, argv[3])) return 1; if (!write_header(returnmap, argv[4])) return 1; + if (!write_header(returnmap_cubic, argv[5])) return 1; return 0; } diff --git a/tests/generated/numsim_material_check_driver.cpp b/tests/generated/numsim_material_check_driver.cpp index d53f381..5293919 100644 --- a/tests/generated/numsim_material_check_driver.cpp +++ b/tests/generated/numsim_material_check_driver.cpp @@ -37,6 +37,7 @@ #include "numsim-materials/materials/tensor_component_stepper.h" #include "Viscoelastic.h" // generated: σ = α·ε (tensor stress from scalar state) #include "ReturnMap.h" // generated: implicit residual R(z,ε)=z−c·tr(ε), σ=z·ε +#include "ReturnMapCubic.h" // generated: NONLINEAR residual R=z+z³−c·tr(ε) #endif #include @@ -54,6 +55,7 @@ using Nonlinear = numsim::materials::generated::NonlinearDecay; #ifdef NCG_TENSOR_E2E using Visco = numsim::materials::generated::Viscoelastic; using ReturnMap = numsim::materials::generated::ReturnMap; +using ReturnMapCubic = numsim::materials::generated::ReturnMapCubic; using tensor2 = tmech::tensor; #endif @@ -320,6 +322,54 @@ TEST(NumSimMaterialEndToEnd, ResidualReturnMapSolvesAgainstBackwardEuler) { EXPECT_NEAR(sig(0, 1), T{0}, 1e-12); // off-diagonal strain is 0 } +// NONLINEAR residual R(z,ε) = z + z³ − c·tr(ε), so ∂R/∂z = 1 + 3z². This is the +// test that VALIDATES THE EMITTED JACOBIAN: with a tight Newton budget and a +// root z≈0.682 (where ∂R/∂z≈2.4 ≫ 1), a wrong derivative oscillates and the +// converged residual is NOT ~0 — unlike the linear case where any jacobian +// reaches the root. We drive c·tr(ε)=1 (strain increment 0.5 on one component). +TEST(NumSimMaterialEndToEnd, NonlinearResidualValidatesEmittedJacobian) { + ctx_type ctx; + param_type p; + + p.insert("name", "stepper"); + p.insert("increment", T{0.5}); // tr(ε) = 0.5 after one update + p.insert>("indices", {0, 0}); + ctx.create>(p); + + // Tight iteration budget: the correct jacobian (1+3z²) converges in ~6 Newton + // steps; a wrong constant jacobian would not reach tol here. + p.clear(); + p.insert("name", "solver"); + p.insert("tolerance", T{1e-13}); + p.insert("max_iter", 12); + ctx.create>(p); + + p.clear(); + const T c = T{2.0}; + p.insert("name", "ReturnMapCubic"); + p.insert("c", c); + p.insert("solver_source", "solver"); + p.insert("strain_source", "stepper"); + ctx.create(p); + + ctx.finalize(); + ctx.update(); + ctx.commit(); + + const auto eps = read_tensor(ctx, "stepper", "strain"); + const T tr = eps(0, 0) + eps(1, 1) + eps(2, 2); + const T z = ctx.get("ReturnMapCubic", "z"); + const auto sig = read_tensor(ctx, "ReturnMapCubic", "stress"); + + // The converged state must satisfy R(z,ε)=0: z + z³ = c·tr(ε). A wrong + // emitted ∂R/∂z fails to converge within the budget, so this residual ≠ 0. + const T residual = z + z * z * z - c * tr; + EXPECT_NEAR(residual, T{0}, 1e-9) << "z=" << z << " (∂R/∂z must be 1+3z²)"; + // Sanity: the root of z+z³=1 is ≈ 0.6823278. + EXPECT_NEAR(z, T{0.6823278038280193}, 1e-6); + EXPECT_NEAR(sig(0, 0), z * eps(0, 0), 1e-10); // σ = z·ε +} + #endif // NCG_TENSOR_E2E } // namespace