Phase 2a: emit strain-coupled implicit-residual materials (Mode 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<backward_euler> +
+    // 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, {}};

tests/CMakeLists.txt

@@ -194,22 +194,29 @@ 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)
+    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_VISCO_HEADER} ${GENERATED_RETURNMAP_HEADER}
+               ${GENERATED_RETURNMAP_CUBIC_HEADER}
         COMMAND ${CMAKE_COMMAND} -E make_directory
                 ${CMAKE_CURRENT_BINARY_DIR}/generated
         COMMAND $<TARGET_FILE:generate_numsim_material_check>
                 ${GENERATED_LINEAR_HEADER} ${GENERATED_NONLINEAR_HEADER}
-                ${GENERATED_VISCO_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)
 
     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}
+        ${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.

tests/NumSimMaterialTargetTest.cpp

@@ -12,6 +12,8 @@
 #include <numsim_cas/scalar/scalar_operators.h>
 #include <numsim_cas/scalar/scalar_std.h>
 #include <numsim_cas/tensor/tensor_definitions.h>
+#include <numsim_cas/tensor_to_scalar/tensor_to_scalar_functions.h>
+#include <numsim_cas/tensor_to_scalar/tensor_to_scalar_operators.h>
 
 #include <gtest/gtest.h>
 
@@ -329,5 +331,163 @@ TEST(NumSimMaterialTarget, FloatDefaultRoundTrips) {
       << src;
 }
 
+// ─── 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<scalar_constant>(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<backward_euler>, 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<Traits>"),
+            std::string::npos) << h;
+  EXPECT_NE(h.find("add_material_ref<solver_type>"), std::string::npos) << h;
+  EXPECT_NE(h.find("m_solver.get().solve(eval)"), std::string::npos) << h;
+  EXPECT_NE(h.find("add_history_output<value_type>(\"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}.
+  EXPECT_NE(h.find("return {residual, jacobian};"), 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;
+  // The solver_source param is required.
+  EXPECT_NE(h.find("insert<std::string>(\"solver_source\")"), std::string::npos)
+      << 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<scalar_constant>(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
+// `<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<scalar_constant>(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<scalar_constant>(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<scalar_constant>(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.
+TEST(NumSimMaterialTarget, RejectsTangentOnResidualMaterial) {
+  auto m = build_return_map();
+  m.add_algorithmic_tangent("dstress_dstrain", "stress", "strain");
+  auto const msg = emit_throw_message(m);
+  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<scalar_constant>(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<scalar_constant>(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
 } // namespace numsim::codegen

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<cas::scalar_constant>(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");

tests/generated/generate_numsim_material_check.cpp

@@ -20,6 +20,8 @@
 #include <numsim_cas/scalar/scalar_std.h>
 #include <numsim_cas/tensor/tensor_definitions.h>
 #include <numsim_cas/tensor/tensor_operators.h>
+#include <numsim_cas/tensor_to_scalar/tensor_to_scalar_functions.h>
+#include <numsim_cas/tensor_to_scalar/tensor_to_scalar_operators.h>
 
 #include <fstream>
 #include <iostream>
@@ -51,9 +53,10 @@ bool write_header(ConstitutiveModel const& model, char const* path) {
 } // namespace
 
 int main(int argc, char** argv) {
-  if (argc < 4) {
+  if (argc < 6) {
     std::cerr << "usage: generate_numsim_material_check <linear.h> "
-                 "<nonlinear.h> <viscoelastic.h>\n";
+                 "<nonlinear.h> <viscoelastic.h> <returnmap.h> "
+                 "<returnmapcubic.h>\n";
     return 2;
   }
 
@@ -95,8 +98,43 @@ 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<backward_euler> + 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<scalar_constant>(0.0));
+    returnmap.add_scalar_residual_equation(z, z.current - c * trace(eps));
+    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<scalar_constant>(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;
 }