Phase 2b: emit the strain-coupled consistent tangent dσ/dε for Mode-B residual materials

CMakeLists.txt

@@ -170,6 +170,13 @@ target_include_directories(numsim_codegen_headers
 target_link_libraries(numsim_codegen_headers INTERFACE NumSim_CAS)
 target_compile_features(numsim_codegen_headers INTERFACE cxx_std_23)
 
+# NOTE: the ∂(tensor)/∂(scalar) seam (detail::diff_tensor_wrt_scalar, cas#275)
+# self-detects the overload via __has_include of the CAS visitor header — see
+# include/numsim_codegen/passes/internal/algorithmic_tangent.h. No compile
+# definition is set here: an unconditional macro would make that inline
+# function's body differ between TUs that do and don't see it (an ODR hazard),
+# and would also lie if the pin ever drifted to a CAS without the overload.
+
 add_library(numsim_codegen STATIC
     src/targets/standalone_cxx.cpp
     src/targets/moose_material.cpp

include/numsim_codegen/passes/internal/algorithmic_tangent.h

@@ -22,9 +22,9 @@
 // coefficient product-rule term). Tracked upstream as numsim-cas#275.
 //
 // `diff_tensor_wrt_scalar` is the single seam through which that derivative is
-// taken. Until #275 lands it throws with a precise diagnostic; when the
-// overload ships, define NUMSIM_CODEGEN_HAVE_DIFF_TENSOR_WRT_SCALAR (or wire it
-// to a CAS feature macro) and the body collapses to a one-line `cas::diff`.
+// taken. It self-detects the cas#275 overload from the pin (see the
+// __has_include below) and collapses to a one-line `cas::diff`; absent the
+// overload it throws with a precise diagnostic.
 //
 // NOTE: this term only fires for a strain-coupled (t2s) residual. With the
 // current scalar-residual Newton machinery dx/dε ≡ 0, so AlgorithmicTangentPass
@@ -44,7 +44,15 @@ namespace numsim::codegen::detail {
     cas::expression_holder<cas::tensor_expression> const &expr,
     cas::expression_holder<cas::scalar_expression> const &arg)
     -> cas::expression_holder<cas::tensor_expression> {
-#ifdef NUMSIM_CODEGEN_HAVE_DIFF_TENSOR_WRT_SCALAR
+// Capability detection. The cas#275 overload — diff(tensor, scalar) — ships as
+// the `tensor_differentiation_wrt_scalar` visitor header; detect its PRESENCE
+// directly rather than relying on an externally-defined compile macro. This
+// keeps `diff_tensor_wrt_scalar` a SINGLE inline definition across every TU
+// (the body no longer depends on whether a TU links `numsim_codegen_headers`),
+// which would otherwise be an ODR hazard. `NUMSIM_CODEGEN_HAVE_DIFF_TENSOR_WRT_SCALAR`
+// is kept only as an explicit override/escape hatch.
+#if defined(NUMSIM_CODEGEN_HAVE_DIFF_TENSOR_WRT_SCALAR) ||                       \
+    __has_include(<numsim_cas/tensor/visitors/tensor_differentiation_wrt_scalar.h>)
   return cas::diff(expr, arg);
 #else
   (void)expr;

src/targets/numsim_material.cpp

@@ -3,11 +3,18 @@
 #include <numsim_codegen/code_emit/code_emit_pipeline.h>
 #include <numsim_codegen/code_emit/codegen_context.h>
 #include <numsim_codegen/code_emit/leaf_collector.h>
+#include <numsim_codegen/passes/internal/algorithmic_tangent.h>
 #include <numsim_codegen/recipe.h>
 
 #include <numsim_cas/core/diff.h>
+#include <numsim_cas/scalar/scalar_constant.h>
+#include <numsim_cas/scalar/scalar_operators.h>
 #include <numsim_cas/tensor/tensor_diff.h>
+#include <numsim_cas/tensor/tensor_functions.h>
+#include <numsim_cas/tensor/tensor_operators.h>
 #include <numsim_cas/tensor_to_scalar/tensor_to_scalar_diff.h>
+#include <numsim_cas/tensor_to_scalar/tensor_to_scalar_functions.h>
+#include <numsim_cas/tensor_to_scalar/tensor_to_scalar_operators.h>
 
 #include <cmath>
 #include <format>
@@ -181,13 +188,6 @@ std::vector<EmittedFile> emit_residual_material(ConstitutiveModel const &model)
         "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(
@@ -397,6 +397,105 @@ std::vector<EmittedFile> emit_residual_material(ConstitutiveModel const &model)
     }
   }
 
+  // ── Consistent tangent(s) dσ/dε (strain-coupled, Mode B) ──
+  // For a requested tangent of a stress output σ w.r.t. a strain input ε:
+  //   dσ/dε = ∂σ/∂ε + ∂σ/∂z ⊗ dz/dε ,   dz/dε = −∂R/∂ε / ∂R/∂z
+  // where z is the scalar Newton state. ∂R/∂z is the scalar jacobian (a scalar —
+  // NO matrix inverse, since the state is scalar), so the whole tangent is a
+  // rank-4 tensor EXPRESSION in (z, ε, params). It is emitted as an extra
+  // post-solve output, evaluated at the converged z exactly like the stress
+  // (and, like every output, bound to &compute). ∂σ/∂z is taken through the
+  // detail::diff_tensor_wrt_scalar seam (cas#275).
+  for (auto const &t : model.tangents()) {
+    // of_output must be an emitted TENSOR output (the stress).
+    OutputDecl const *src = nullptr;
+    for (auto const &o : model.outputs()) {
+      if (o.name == t.of_output) {
+        src = &o;
+        break;
+      }
+    }
+    if (src == nullptr) {
+      throw std::runtime_error(
+          "NumSimMaterialTarget: consistent tangent '" + t.name +
+          "' references output '" + t.of_output +
+          "', which is not a declared output.");
+    }
+    if (src->kind != OutputDecl::Kind::Tensor) {
+      throw std::runtime_error(
+          "NumSimMaterialTarget: consistent tangent '" + t.name +
+          "' differentiates output '" + t.of_output +
+          "', which is scalar — dσ/dε requires a tensor (stress) output.");
+    }
+    // wrt_input must be a declared tensor input (the strain), and share σ's dim.
+    SymbolDecl const *arg = nullptr;
+    for (auto const &ti : tensor_inputs) {
+      if (ti.name == t.wrt_input) {
+        arg = &ti;
+        break;
+      }
+    }
+    if (arg == nullptr) {
+      throw std::runtime_error(
+          "NumSimMaterialTarget: consistent tangent '" + t.name +
+          "' differentiates w.r.t. '" + t.wrt_input +
+          "', which is not a declared tensor input.");
+    }
+    if (arg->dim != src->dim) {
+      throw std::runtime_error(
+          "NumSimMaterialTarget: consistent tangent '" + t.name + "': output '" +
+          t.of_output + "' (dim " + std::to_string(src->dim) + ") and input '" +
+          t.wrt_input + "' (dim " + std::to_string(arg->dim) +
+          ") have different tensor dimensions.");
+    }
+    // The tangent output name shares the same collision surface as any output.
+    if (is_reserved_residual(t.name) || t.name == cur_name ||
+        param_names.contains(t.name) || tensor_input_names.contains(t.name)) {
+      throw std::runtime_error(
+          "NumSimMaterialTarget: consistent-tangent name '" + t.name +
+          "' collides with an emitted member, the state, a parameter, or a "
+          "tensor input; rename it.");
+    }
+    for (auto const &o : outputs) {
+      if (o.name == t.name) {
+        throw std::runtime_error(
+            "NumSimMaterialTarget: consistent-tangent name '" + t.name +
+            "' collides with output '" + o.name + "'; rename it.");
+      }
+    }
+
+    // Resolve σ and ε expression handles.
+    auto const &sigma =
+        std::get<cas::expression_holder<cas::tensor_expression>>(src->expr);
+    cas::expression_holder<cas::tensor_expression> eps;
+    for (auto const &[name, h] : model.tensor_symbol_map()) {
+      if (name == t.wrt_input) eps = h;
+    }
+    if (!eps.is_valid()) {
+      throw std::runtime_error(
+          "NumSimMaterialTarget: cannot resolve the tensor-input handle for '" +
+          t.wrt_input + "' — input/symbol maps out of sync.");
+    }
+
+    // dσ/dε = ∂σ/∂ε + otimes(∂σ/∂z, (−1/∂R/∂z)·∂R/∂ε).
+    auto const dsig_deps = cas::diff(sigma, eps);                     // rank-4
+    auto const dsig_dz = detail::diff_tensor_wrt_scalar(sigma, cur_expr); // rank-2
+    auto const dR_deps = cas::diff(req.residual, eps);               // rank-2
+    auto const dR_dz = cas::diff(req.residual, cur_expr);            // scalar
+    auto const neg1 = cas::make_expression<cas::scalar_constant>(-1.0);
+    auto const dz_deps = (neg1 / dR_dz) * dR_deps;                   // rank-2
+    auto const tangent = dsig_deps + cas::otimes(dsig_dz, dz_deps);  // rank-4
+
+    CodeGenContext tc;
+    CodeEmitPipeline tp(tc);
+    register_all(tc);
+    tc.reset();
+    auto const trhs = tp.tensor().apply(tangent);
+    // dim from the operands (checked equal); rank = rank(σ) + rank(ε).
+    outputs.push_back({t.name, tc.render_statements("    "), trhs, true,
+                       src->dim, src->rank + arg->rank});
+  }
+
   // Emitted-member uniqueness guard (same hazard as the rate path: synthesized
   // member names can collide with recipe symbols).
   {
@@ -439,12 +538,21 @@ std::vector<EmittedFile> emit_residual_material(ConstitutiveModel const &model)
     << "' 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 << "// models only states whose solved increment is non-negative. This "
+       "caveat\n";
+  h << "// applies to the CONSISTENT TANGENT dσ/dε too (if emitted): it is\n";
+  h << "// evaluated at the solved state via the implicit-function theorem, so "
+       "on a\n";
+  h << "// clamped (elastic) increment where the raw root is negative the "
+       "tangent is\n";
+  h << "// outside its domain of validity — the correct value there is ∂σ/∂ε "
+       "alone.\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 <tmech/tmech.h>\n";
+  h << "#include <stdexcept>\n";
   h << "#include <utility>\n\n";
   h << "namespace numsim::materials::generated {\n\n";
   h << "template <typename Traits>\n";
@@ -546,6 +654,16 @@ std::vector<EmittedFile> emit_residual_material(ConstitutiveModel const &model)
   h << "    };\n";
   h << "    const value_type d" << cur_name
     << " = m_solver.get().solve(eval);\n";
+  // Reject loudly if the Newton solve did not converge: every output (stress AND
+  // the consistent tangent) is evaluated at this state, so a non-converged root
+  // would silently feed a wrong stress / wrong stiffness into the graph. (Note:
+  // this does NOT catch backward_euler's std::max(x,0) clamp returning a
+  // non-root on the elastic branch — see the header caveat; that needs an
+  // upstream was_clamped() accessor, numsim-materials follow-up.)
+  h << "    if (!m_solver.get().converged())\n";
+  h << "      throw std::runtime_error(\"" << cls
+    << "::compute: backward_euler did not converge solving R(" << cur_name
+    << ", inputs)=0; stress/tangent would be evaluated at a non-root state.\");\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_"

tests/AlgorithmicTangentTest.cpp

@@ -393,19 +393,45 @@ TEST(AlgorithmicTangent, PlainInputTangentStaysNonSymmetrized) {
   EXPECT_EQ(src.find("tmech::otimesl"), std::string::npos) << src;
 }
 
-// Pins the scaffold's flip-on point: the ∂σ/∂x seam throws a precise
-// numsim-cas#275 diagnostic until the upstream diff(tensor, scalar) lands.
-TEST(AlgorithmicTangent, DiffTensorWrtScalarSeamThrowsUntilCas275) {
+// Flip-on point (PR 2b): cas#275 — diff(tensor_expression, scalar_expression) —
+// is in the pin, and the seam is enabled via
+// NUMSIM_CODEGEN_HAVE_DIFF_TENSOR_WRT_SCALAR. The seam now RETURNS the real
+// ∂(tensor)/∂(scalar), including the scalar-coefficient product-rule term the
+// stub could not compute: ∂(x·ε)/∂x = ε.
+TEST(AlgorithmicTangent, DiffTensorWrtScalarSeamComputesProductRuleTerm) {
   using namespace numsim::cas;
   auto eps = make_expression<tensor>("eps", 3, std::size_t{2});
   auto x = make_expression<scalar>("x");
-  try {
-    (void)detail::diff_tensor_wrt_scalar(eps, x);
-    FAIL() << "expected the diff(tensor,scalar) seam to throw";
-  } catch (std::runtime_error const &e) {
-    EXPECT_NE(std::string(e.what()).find("numsim-cas#275"), std::string::npos)
-        << e.what();
-  }
+
+  auto d = detail::diff_tensor_wrt_scalar(x * eps, x); // ∂(x·ε)/∂x = ε
+  ASSERT_TRUE(d.is_valid());
+
+  CodeGenContext ctx;
+  CodeEmitPipeline p(ctx);
+  ctx.register_symbol_tensor(eps, "eps");
+  ctx.register_symbol_scalar(x, "x");
+  ctx.reset();
+  // Product rule: 1·ε + x·0 = ε. The stub threw here; the flip must render ε.
+  EXPECT_EQ(p.tensor().apply(d), "eps");
+
+  // Non-trivial coefficient: ∂((x·x)·ε)/∂x = 2x·ε. This exercises the
+  // scalar-coefficient product-rule term that a degenerate ∂(x·ε)/∂x (factor 1)
+  // does not — a diff that returned ε instead of 2x·ε would pass the check above
+  // but fail here.
+  auto d2 = detail::diff_tensor_wrt_scalar((x * x) * eps, x);
+  ASSERT_TRUE(d2.is_valid());
+  CodeGenContext ctx2;
+  CodeEmitPipeline p2(ctx2);
+  ctx2.register_symbol_tensor(eps, "eps");
+  ctx2.register_symbol_scalar(x, "x");
+  ctx2.reset();
+  auto const r2 = p2.tensor().apply(d2);
+  // 2x·ε needs CSE temps, so apply() returns a temp ref and the real work is in
+  // the rendered statements. The whole program is the statements + final expr.
+  auto const prog2 = ctx2.render_statements() + r2;
+  EXPECT_NE(prog2, "eps") << prog2;                        // NOT the trivial term
+  EXPECT_NE(prog2.find("x"), std::string::npos) << prog2;  // carries the 2x factor
+  EXPECT_NE(prog2.find("eps"), std::string::npos) << prog2;
 }
 
 // Round-2 review (test-quality MAJOR-5): the pass running ALONGSIDE local Newton

tests/NumSimMaterialTargetTest.cpp

@@ -456,11 +456,80 @@ TEST(NumSimMaterialTarget, RejectsResidualLocalNameCollision) {
 // 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) {
+// Phase 2b: the strain-coupled consistent tangent dσ/dε for a Mode-B residual
+// material. σ = z·ε, R = z − c·tr(ε) ⇒
+//   dσ/dε = ∂σ/∂ε + ∂σ/∂z ⊗ dz/dε = z·I⁴ˢ + ε ⊗ (c·I),   dz/dε = −∂R/∂ε/∂R/∂z = c·I.
+// The correction term ε⊗(c·I) is what a naive ∂σ/∂ε alone would drop.
+TEST(NumSimMaterialTarget, EmitsConsistentTangentForResidualMaterial) {
   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;
+  auto const h = header_of(NumSimMaterialTarget{}.emit(m));
+
+  // Emitted as a rank-4 output, bound to compute() like every output.
+  EXPECT_NE(h.find("add_output<tmech::tensor<value_type, 3, 4>>"),
+            std::string::npos)
+      << h;
+  EXPECT_NE(h.find("\"dstress_dstrain\", &ReturnMap::compute"),
+            std::string::npos)
+      << h;
+  EXPECT_NE(h.find("tmech::tensor<value_type, 3, 4>& m_out_dstress_dstrain;"),
+            std::string::npos)
+      << h;
+
+  // The implicit correction ∂σ/∂z ⊗ dz/dε: ∂σ/∂z = ε (the strain), assembled as
+  // an outer product. Pinning this is the whole point — a dropped coupling term
+  // would still leave a plausible-but-wrong ∂σ/∂ε.
+  EXPECT_NE(h.find("tmech::outer_product<tmech::sequence<1, 2>, "
+                   "tmech::sequence<3, 4>>(strain,"),
+            std::string::npos)
+      << h;
+  // dz/dε = −∂R/∂ε/∂R/∂z carries a negation — pin the sign (a dropped/flipped
+  // sign flips the coupling term, which the compiler-independent layer must
+  // catch, not only the gcc-gated numeric e2e).
+  EXPECT_NE(h.find("-1.0 *"), std::string::npos) << h;
+  // The explicit base term ∂σ/∂ε = z·I⁴ˢ: minor-symmetric identity (BOTH otimesu
+  // AND otimesl — mirroring the rate-path tangent test) scaled by the state z.
+  EXPECT_NE(h.find("tmech::otimesu(tmech::eye<double, 3, 2>()"),
+            std::string::npos)
+      << h;
+  EXPECT_NE(h.find("tmech::otimesl(tmech::eye<double, 3, 2>()"),
+            std::string::npos)
+      << h;
+  EXPECT_NE(h.find("z * "), std::string::npos) << h; // the z·I⁴ˢ coefficient
+  EXPECT_NE(h.find("m_out_dstress_dstrain ="), std::string::npos) << h;
+
+  // The tangent is evaluated AFTER the solve (uses the converged z).
+  auto const solve_at = h.find("m_solver.get().solve(eval)");
+  auto const tangent_at = h.find("m_out_dstress_dstrain =");
+  ASSERT_NE(solve_at, std::string::npos);
+  ASSERT_NE(tangent_at, std::string::npos);
+  EXPECT_LT(solve_at, tangent_at) << "tangent must use the converged state";
+}
+
+// A tangent can only differentiate a TENSOR (stress) output.
+TEST(NumSimMaterialTarget, RejectsTangentOfScalarOutputOnResidualMaterial) {
+  auto m = build_return_map();
+  m.add_output("scalar_out", 2.0 * make_expression<scalar>("c"));
+  m.add_algorithmic_tangent("dsc_dstrain", "scalar_out", "strain");
+  EXPECT_NE(emit_throw_message(m).find("is scalar"), std::string::npos);
+}
+
+// wrt_input must be a declared tensor input.
+TEST(NumSimMaterialTarget, RejectsTangentWrtUnknownInputOnResidualMaterial) {
+  auto m = build_return_map();
+  m.add_algorithmic_tangent("dstress_dghost", "stress", "ghost");
+  EXPECT_NE(emit_throw_message(m).find("not a declared tensor input"),
+            std::string::npos);
+}
+
+// The tangent output name shares the emitted-identifier collision surface. A
+// name like "solver" is not a recipe symbol (so the request-time
+// availability check passes) but collides with the emitted m_solver member —
+// caught by the emit-time guard.
+TEST(NumSimMaterialTarget, RejectsTangentNameCollidingWithEmittedMember) {
+  auto m = build_return_map();
+  m.add_algorithmic_tangent("solver", "stress", "strain");
+  EXPECT_NE(emit_throw_message(m).find("collides"), std::string::npos);
 }
 
 // A residual material needs at least one output to anchor compute() (the output

tests/generated/generate_numsim_material_check.cpp

@@ -111,6 +111,9 @@ int main(int argc, char** argv) {
         "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);
+    // Phase 2b: the strain-coupled consistent tangent dσ/dε. The driver checks
+    // the off-block coupling term the naive ∂σ/∂ε alone misses (C_{0011}=c·ε₀₀).
+    returnmap.add_algorithmic_tangent("dstress_dstrain", "stress", "strain");
   }
 
   // NONLINEAR residual to exercise the t2s-wrt-scalar jacobian (∂R/∂z, cas#285).
@@ -129,6 +132,10 @@ int main(int argc, char** argv) {
     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);
+    // Phase 2b: a tangent on the NONLINEAR residual — here ∂R/∂z = 1+3z² ≠ 1, so
+    // dz/dε = −∂R/∂ε/∂R/∂z genuinely exercises the division by the jacobian
+    // (the linear ReturnMap has ∂R/∂z≡1, where dropping the divisor is invisible).
+    returnmap_cubic.add_algorithmic_tangent("dstress_dstrain", "stress", "strain");
   }
 
   if (!write_header(linear, argv[1])) return 1;