diff --git a/CMakeLists.txt b/CMakeLists.txt index 8a3324c..d7b5556 100644 --- a/CMakeLists.txt +++ b/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 diff --git a/include/numsim_codegen/passes/internal/algorithmic_tangent.h b/include/numsim_codegen/passes/internal/algorithmic_tangent.h index 25f7469..652ec9c 100644 --- a/include/numsim_codegen/passes/internal/algorithmic_tangent.h +++ b/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 const &expr, cas::expression_holder const &arg) -> cas::expression_holder { -#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() return cas::diff(expr, arg); #else (void)expr; diff --git a/src/targets/numsim_material.cpp b/src/targets/numsim_material.cpp index 64a5843..5f74fff 100644 --- a/src/targets/numsim_material.cpp +++ b/src/targets/numsim_material.cpp @@ -3,11 +3,18 @@ #include #include #include +#include #include #include +#include +#include #include +#include +#include #include +#include +#include #include #include @@ -181,13 +188,6 @@ std::vector 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 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>(src->expr); + cas::expression_holder 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(-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 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 \n"; + h << "#include \n"; h << "#include \n\n"; h << "namespace numsim::materials::generated {\n\n"; h << "template \n"; @@ -546,6 +654,16 @@ std::vector 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_" diff --git a/tests/AlgorithmicTangentTest.cpp b/tests/AlgorithmicTangentTest.cpp index db1cb80..9d58754 100644 --- a/tests/AlgorithmicTangentTest.cpp +++ b/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("eps", 3, std::size_t{2}); auto x = make_expression("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 diff --git a/tests/NumSimMaterialTargetTest.cpp b/tests/NumSimMaterialTargetTest.cpp index 85f0b85..6b03e33 100644 --- a/tests/NumSimMaterialTargetTest.cpp +++ b/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>"), + std::string::npos) + << h; + EXPECT_NE(h.find("\"dstress_dstrain\", &ReturnMap::compute"), + std::string::npos) + << h; + EXPECT_NE(h.find("tmech::tensor& 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<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()"), + std::string::npos) + << h; + EXPECT_NE(h.find("tmech::otimesl(tmech::eye()"), + 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("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 diff --git a/tests/generated/generate_numsim_material_check.cpp b/tests/generated/generate_numsim_material_check.cpp index dffbdae..0ab7d0c 100644 --- a/tests/generated/generate_numsim_material_check.cpp +++ b/tests/generated/generate_numsim_material_check.cpp @@ -111,6 +111,9 @@ int main(int argc, char** argv) { "z", make_expression(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; diff --git a/tests/generated/numsim_material_check_driver.cpp b/tests/generated/numsim_material_check_driver.cpp index 5293919..f0ea941 100644 --- a/tests/generated/numsim_material_check_driver.cpp +++ b/tests/generated/numsim_material_check_driver.cpp @@ -322,6 +322,61 @@ TEST(NumSimMaterialEndToEnd, ResidualReturnMapSolvesAgainstBackwardEuler) { EXPECT_NEAR(sig(0, 1), T{0}, 1e-12); // off-diagonal strain is 0 } +// Phase 2b: the strain-coupled CONSISTENT TANGENT of the generated ReturnMap, +// verified numerically through the real backward_euler solve. For σ = z·ε with +// z solving R = z − c·tr(ε) = 0: +// dσ/dε = ∂σ/∂ε + ∂σ/∂z ⊗ dz/dε = z·I⁴ˢ + ε ⊗ (c·I). +// The second term is the implicit coupling a naive explicit ∂σ/∂ε drops. The +// load-bearing assertion is the off-block C_{0011} = c·ε₀₀ ≠ 0 — exactly the +// component the rate-path (strain-only) Viscoelastic tangent above has as ZERO. +TEST(NumSimMaterialEndToEnd, ResidualReturnMapConsistentTangentHasCouplingTerm) { + ctx_type ctx; + param_type p; + + p.insert("name", "stepper"); + p.insert("increment", T{0.02}); + p.insert>("indices", {0, 0}); + ctx.create>(p); + + 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 e00 = eps(0, 0); + const T z = ctx.get("ReturnMap", "z"); // = c·tr(ε) = c·e00 + + using tensor4 = tmech::tensor; + auto* cp = dynamic_cast< + numsim_core::property*>( + ctx.find_property("ReturnMap", "dstress_dstrain")); + ASSERT_NE(cp, nullptr); + const auto C = cp->get(); + + // Coupling term (the whole point): C_{0011} = z·I⁴ˢ_{0011} + ε₀₀·(c·I)_{11} + // = 0 + c·e00. Naive ∂σ/∂ε ⇒ 0. + EXPECT_NEAR(C(0, 0, 1, 1), c * e00, 1e-10); + EXPECT_GT(std::abs(C(0, 0, 1, 1)), 1e-6) << "coupling term must be nonzero"; + // Explicit base term (correction vanishes here): C_{0101} = z·I⁴ˢ_{0101} = z·½. + EXPECT_NEAR(C(0, 1, 0, 1), z * T{0.5}, 1e-10); + // Diagonal: C_{0000} = z·I⁴ˢ_{0000} + c·e00 = z + c·e00. + EXPECT_NEAR(C(0, 0, 0, 0), z + c * e00, 1e-10); +} + // 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 @@ -368,6 +423,24 @@ TEST(NumSimMaterialEndToEnd, NonlinearResidualValidatesEmittedJacobian) { // 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·ε + + // Phase 2b: the tangent on the NONLINEAR residual VALIDATES THE DIVISION by + // ∂R/∂z. Here ∂R/∂z = 1+3z² ≠ 1, so dz/dε = c·I/(1+3z²) and the coupling term + // C_{0011} = ε₀₀·c/(1+3z²). If the emitter dropped the /∂R/∂z factor (which the + // LINEAR ReturnMap's ∂R/∂z≡1 cannot detect), C_{0011} would be ε₀₀·c instead — + // off by the ~2.4× jacobian here. This is the load-bearing division check. + using tensor4 = tmech::tensor; + auto* cp = dynamic_cast< + numsim_core::property*>( + ctx.find_property("ReturnMapCubic", "dstress_dstrain")); + ASSERT_NE(cp, nullptr); + const auto C = cp->get(); + const T dRdz = T{1} + T{3} * z * z; // = ∂R/∂z at the converged z + EXPECT_GT(dRdz, T{2}) << "must be ≫1 so the division is observable"; + EXPECT_NEAR(C(0, 0, 1, 1), c * eps(0, 0) / dRdz, 1e-9); + // Contrast: without the division it would be c·ε₀₀ — assert we are NOT that. + EXPECT_GT(std::abs(C(0, 0, 1, 1) - c * eps(0, 0)), 1e-3) + << "tangent must divide by ∂R/∂z, not omit it"; } #endif // NCG_TENSOR_E2E