Phase D: emit consistent tangent dσ/dε as a rank-4 tensor property

src/targets/numsim_material.cpp

@@ -6,6 +6,7 @@
 #include <numsim_codegen/recipe.h>
 
 #include <numsim_cas/core/diff.h>
+#include <numsim_cas/tensor/tensor_diff.h>
 
 #include <cmath>
 #include <format>
@@ -89,13 +90,28 @@ void check_scope(ConstitutiveModel const &model) {
   // properties with their own update callbacks. Only tensor INTERNAL STATE is
   // still blocked (Mandel + vector solver) — a tensor stress from a scalar state
   // needs neither, so it is supported here.
-  if (!model.tangents().empty()) {
-    // Tensor stress outputs ARE emitted now, so a tangent's `of_output` resolves
-    // — but the consistent tangent dσ/dε itself is Phase D (it emits the
-    // tangent-block properties the solver assembles), so reject it here.
-    throw std::runtime_error(
-        "NumSimMaterialTarget: algorithmic tangents are out of scope for the "
-        "rate material (Phase D emits tangent-block properties).");
+  // Algorithmic/consistent tangents dσ/dε ARE emitted (Phase D) as rank-4 tensor
+  // properties via cas::diff. In THIS target's scope the rate cannot depend on
+  // strain (rate-leaf guard), so the integrated state is independent of the
+  // current strain (dx/dε ≡ 0) and dσ/dε = ∂σ/∂ε exactly — no J⁻¹ solve / no
+  // rank-4 inverse, just diff(stress, strain). Validate each spec resolves.
+  for (auto const &t : model.tangents()) {
+    bool of_ok = false;
+    for (auto const &o : model.outputs())
+      if (o.name == t.of_output && o.kind == OutputDecl::Kind::Tensor) of_ok = true;
+    if (!of_ok) {
+      throw std::runtime_error(
+          "NumSimMaterialTarget: tangent '" + t.name + "' differentiates '" +
+          t.of_output + "', which is not a declared tensor output (stress).");
+    }
+    bool wrt_ok = false;
+    for (auto const &in : model.inputs())
+      if (in.name == t.wrt_input && in.kind == SymbolDecl::Kind::Tensor) wrt_ok = true;
+    if (!wrt_ok) {
+      throw std::runtime_error(
+          "NumSimMaterialTarget: tangent '" + t.name + "' differentiates w.r.t. '" +
+          t.wrt_input + "', which is not a declared tensor input (strain).");
+    }
   }
   // Tensor inputs (e.g. strain) ARE wired (Global-edge input_property); scalar
   // inputs are a separate small follow-up — reject those loudly for now.
@@ -303,6 +319,44 @@ auto NumSimMaterialTarget::emit(ConstitutiveModel const &model) const
     }
   }
 
+  // Consistent tangents dσ/dε: a rank-4 tensor property = cas::diff(stress, strain).
+  // Folded into `outputs` as rank-4 tensor outputs, so the ctor/callback/member
+  // emission below handles them uniformly. The strain's roles::Strain symmetric
+  // space makes diff return the minor-symmetric rank-4 identity (P_sym), so the
+  // tangent is minor-symmetric as a stress-strain tangent must be.
+  for (auto const &t : model.tangents()) {
+    cas::expression_holder<cas::tensor_expression> sigma;
+    std::size_t tdim = 0;
+    for (auto const &o : model.outputs()) {
+      if (o.name == t.of_output && o.kind == OutputDecl::Kind::Tensor) {
+        sigma = std::get<cas::expression_holder<cas::tensor_expression>>(o.expr);
+        tdim = o.dim;
+      }
+    }
+    cas::expression_holder<cas::tensor_expression> eps;
+    for (auto const &[name, h] : model.tensor_symbol_map()) {
+      if (name == t.wrt_input) eps = h;
+    }
+    // check_scope already verified both resolve; guard belt-and-braces so a
+    // future check_scope/emit drift surfaces as a clear error, not cas::diff UB.
+    if (!sigma.is_valid() || !eps.is_valid()) {
+      throw std::runtime_error(
+          "NumSimMaterialTarget: tangent '" + t.name +
+          "' could not resolve its stress output / strain input handle.");
+    }
+    auto const tangent_expr = cas::diff(sigma, eps); // rank-4 ∂σ/∂ε
+
+    CodeGenContext tc;
+    CodeEmitPipeline tp(tc);
+    register_scalars(tc);
+    for (auto const &[name, expr] : model.tensor_symbol_map()) {
+      tc.register_symbol_tensor(expr, name);
+    }
+    tc.reset();
+    auto const trhs = tp.tensor().apply(tangent_expr);
+    outputs.push_back({t.name, tc.render_statements("    "), trhs, true, tdim, 4});
+  }
+
   bool has_tensor = !tensor_inputs.empty();
   for (auto const &o : outputs)
     if (o.is_tensor) has_tensor = true;

tests/NumSimMaterialTargetTest.cpp

@@ -253,18 +253,40 @@ TEST(NumSimMaterialTarget, RejectsReservedStateName) {
             std::string::npos);
 }
 
-// The stress TENSOR output is now emitted, so a tangent-bearing recipe is no
-// longer caught by an output guard — the consistent tangent dσ/dε itself is
-// Phase D, so the tangents() guard is now the load-bearing rejection.
-TEST(NumSimMaterialTarget, RejectsTangentBearingRecipe) {
+// Phase D: the consistent tangent dσ/dε is emitted as a rank-4 tensor property
+// via cas::diff(stress, strain). With a strain-INDEPENDENT rate (dα/dε=0) this is
+// the exact consistent tangent. For σ=2μ·ε the tangent is 2μ·P_sym, so the
+// emitted rank-4 RHS is minor-symmetric (otimesu + otimesl).
+TEST(NumSimMaterialTarget, EmitsConsistentTangent) {
   ConstitutiveModel m("WithTangent");
   auto mu = m.add_parameter("mu", 0.5);
   auto a = m.add_scalar_state_variable("a", make_expression<scalar_constant>(0.0));
   m.add_scalar_evolution_equation(a, mu * a.current);
-  auto eps = m.add_tensor_input("eps", 3, 2, roles::Strain);
+  auto eps = m.add_tensor_input("strain", 3, 2, roles::Strain);
   m.add_output("stress", 2 * mu * eps, roles::Stress);
-  m.add_algorithmic_tangent("dstress_deps", "stress", "eps");
-  EXPECT_NE(emit_throw_message(m).find("algorithmic tangents"),
+  m.add_algorithmic_tangent("dstress_dstrain", "stress", "strain");
+  auto const h = header_of(NumSimMaterialTarget{}.emit(m));
+  // rank-4 tensor property + its own callback
+  EXPECT_NE(h.find("add_output<tmech::tensor<value_type, 3, 4>>(\n"
+                   "            \"dstress_dstrain\", "
+                   "&WithTangent::update_dstress_dstrain"),
+            std::string::npos)
+      << h;
+  EXPECT_NE(h.find("void update_dstress_dstrain() {"), std::string::npos) << h;
+  // minor-symmetric: both otimesu AND otimesl present in the tangent RHS
+  EXPECT_NE(h.find("tmech::otimesu"), std::string::npos) << h;
+  EXPECT_NE(h.find("tmech::otimesl"), std::string::npos) << h;
+}
+
+// A tangent whose `of_output` is not a declared tensor output is rejected.
+TEST(NumSimMaterialTarget, RejectsTangentOfMissingStress) {
+  ConstitutiveModel m("BadTangent");
+  auto mu = m.add_parameter("mu", 0.5);
+  auto a = m.add_scalar_state_variable("a", make_expression<scalar_constant>(0.0));
+  m.add_scalar_evolution_equation(a, mu * a.current);
+  m.add_tensor_input("strain", 3, 2, roles::Strain);
+  m.add_algorithmic_tangent("dstress_dstrain", "stress", "strain"); // no "stress"
+  EXPECT_NE(emit_throw_message(m).find("not a declared tensor output"),
             std::string::npos);
 }
 

tests/generated/generate_numsim_material_check.cpp

@@ -90,6 +90,9 @@ int main(int argc, char** argv) {
     viscoelastic.add_scalar_evolution_equation(alpha, K * alpha.current);
     auto eps = viscoelastic.add_tensor_input("strain", 3, 2, roles::Strain);
     viscoelastic.add_output("stress", alpha.current * eps);
+    // Phase D: consistent tangent dσ/dε. Since the rate (K·α) is strain-
+    // independent, dα/dε=0, so dσ/dε = ∂σ/∂ε = α·P_sym (minor-symmetric).
+    viscoelastic.add_algorithmic_tangent("dstress_dstrain", "stress", "strain");
   }
 
   if (!write_header(linear, argv[1])) return 1;

tests/generated/numsim_material_check_driver.cpp

@@ -260,6 +260,17 @@ TEST(NumSimMaterialEndToEnd, TensorStressFromScalarStateAndStrain) {
   const auto sig = read_tensor(ctx, "Viscoelastic", "stress");
   EXPECT_NEAR(sig(0, 0), alpha * eps(0, 0), 1e-12); // σ = α·ε on the driven comp
   EXPECT_NEAR(sig(0, 1), T{0}, 1e-12);              // off-diagonal strain is 0
+
+  // Phase D: consistent tangent dσ/dε = α·P_sym (rank-4, minor-symmetric).
+  using tensor4 = tmech::tensor<T, 3, 4>;
+  auto* cp = dynamic_cast<
+      numsim_core::property<tensor4, numsim::materials::property_traits>*>(
+      ctx.find_property("Viscoelastic", "dstress_dstrain"));
+  ASSERT_NE(cp, nullptr);
+  const auto C = cp->get();
+  EXPECT_NEAR(C(0, 0, 0, 0), alpha, 1e-12);       // P_sym(0,0,0,0)=1
+  EXPECT_NEAR(C(0, 1, 0, 1), alpha * 0.5, 1e-12); // minor-symmetric ½
+  EXPECT_NEAR(C(0, 0, 1, 1), T{0}, 1e-12);        // off-block zero
 }
 
 #endif // NCG_TENSOR_E2E