Design (north-star): automated constitutive materials — derive the return-map residual from physics + a compile-time reduction pass #108

Status: north-star design, deferred. Not an actionable near-term task — it
sits on top of dependencies that don't exist yet (see Sequencing), and it risks
the same over-engineering the parameterization framework hit (#93, descoped at
N=2). Filed to record the target architecture and, importantly, the caveats a
critical review surfaced — so that when it is built, it's built correctly and
for the right reasons.

Vision

Let the human write the constitutive equations (elastic energy ψ, yield
function f, hardening law, flow rule) and have the code generator derive the
return-map residual, Jacobian, and consistent tangent and emit code that is
as fast as a hand-optimized material. Today the human writes the (usually
hand-reduced) residual directly via add_scalar_residual_equation; that's the
lowered IR, one layer too low.

This is the AceGen paradigm — differentiate the potentials — but with a real
symbolic CAS (numsim-cas) and generic differentiation to any order.

Architecture: a material compiler

constitutive equations (ψ, f, hardening, flow)     ← human writes physics
   → front-end
general return-map IR  (coupled residual system)   ← program derives (cas::diff + a fixed
   → ReturnMapReductionPass                            backward-Euler return-map template)
lowered residual  (scalar if reducible, else coupled)
   → LocalNewtonLowering → AlgorithmicTangent → CodeEmit
fast generated code

The SymbolValidation → TensorSpaceConsistency → LocalNewtonLowering → AlgorithmicTangent → CodeEmit pass framework already exists; the reduction is a
strength-reduction pass in the middle — a compile-time transformation, not
human work and not runtime work.

σ = ∂ψ/∂εᵉ, flow direction n = ∂f/∂σ, the Jacobian, and the consistent
tangent are all cas::diff of the human's potentials — verified reachable
(diff(‖dev σ‖, σ) → dev(σ)/‖dev σ‖, i.e. the flow direction, is computed today).

The reduction pass (the interesting middle)

For radial-return-reducible models the coupled tensor system collapses to a
scalar solve. The pass gates the reduction on four checkable preconditions:

Pass → emit the scalar radial-return recipe; fail any → fall back to the coupled
lowering (always correct, slower).

Caveats on the pass (from the critical review — do NOT hand-wave these)

• Reduction reduces the solve, not the state. Real path-dependent J2 still
  carries εᵖ (a rank-2 tensor) + α as history, updated εᵖ ← εᵖ + Δλ·n.
  The scalar collapse removes εᵖ from the unknowns, not from the persistent
  state. So a reduced-but-real J2 needs tensor history support
  (add_history_output<tensor> — the reason numsim-core#16 mattered) + the εᵖ
  update emission. It does NOT reuse the current scalar-history backend
  untouched (our present scalar golden dodges this only by assuming εᵖ_old=0,
  i.e. single-increment).
• Precondition soundness, not just conservatism. A structural match after
  canonicalize() can produce false positives (tensor-expression
  canonicalization has no confluent normal form), which would trigger a wrong
  reduction — silent-wrong material. Soundness of P1/P2 must be argued, and the
  reduction must be verified against an independent reference (below), not trusted.
• "Structure-keyed template library" is expert labor. Each reducible family
  (J2 → Drucker-Prager → cap → isotropic damage) is a new precondition set + a
  hand-authored efficient template. The pass automates dispatch and the
  derivative bookkeeping, not the per-family algorithm.

The elastic/plastic switch (KKT) — a real gap, not a footnote

The pass emits the plastic-branch residual. A real material branches:
f(σ_trial) ≤ 0 → elastic; else plastic. This is currently faked by
backward_euler's max(x,0) clamp, which is known-fragile (it caused the
Voce non-convergence near the boundary). A proper design needs: the branch, the
discontinuous consistent tangent (continuum elastic vs. algorithmic
elastoplastic), and the per-QP branch cost. This is arguably harder than the
reduction and blocks even the hand-written path from being genuinely
path-dependent.

Corrected sequencing (this is the main output of the review)

The reduction pass sits on top of things that don't exist. Build in this order:

1. Constitutive front-end — declare (ψ, f, hardening, flow); derive the
   residual system. This is the automation win, independent of any reduction,
   and removes the residual-by-hand problem for the general case.
2. Coupled return-map path (epic Epic: generated mixed-type return-map materials with consistent tangent (graph-coupled) #102, numsim-materials#11 Mandel + MOOSE backend: add runtime parity test against mock RankTwoTensor #12 vector
   solver) — unlocks the models you can't hand-reduce, and is the
   correctness reference the reduction pass will need.
3. Tensor history + the KKT switch — true gaps today (see above), needed even
   for a hand-written path-dependent J2.
4. ReturnMapReductionPass — last, evidence-gated. Only worth building once
   there are enough reducible families to amortize the template library and the
   coupled solver exists as its reference.

Correctness / acceptance (corrected)

• Self-FD is necessary but NOT sufficient. numerical_diff_checker verifies
  tangent↔stress consistency, not that σ(ε) is the correct return-mapped
  stress. A wrong reduction with a self-consistent-but-wrong σ(ε) passes self-FD.
• Acceptance for the reduction pass: the reduced material must FD-agree with the
  coupled (un-reduced) material at sampled states including the elastic/plastic
  boundary, apex/vertex, and (for spectral models) coincident-eigenvalue states —
  the measure-zero regimes a random sweep misses are exactly where reductions break.
• Where a coupled reference doesn't yet exist, cross-check against a trusted
  hand-written material (e.g. numsim-materials' own J2) — per-family, not general.

Performance framing (honest)

The scalar form eliminates the local linear solve — decisive, since generated
constitutive code is the innermost FE loop. But quote it honestly:

• baseline is a condensed general implementation (nobody solves a naive 8×8 for
  J2), not the strawman;
• the solve is a fraction of the eval (predictor, dev/norm, reconstruction are
  common), so the whole-eval speedup is ~2–10×, with the solve itself 1–2
  orders. Real, but state it at eval level.

Why record it now

The idea — reduction as a compiler pass, physics-in / fast-code-out — is the
right north-star and the same generic-differentiation thesis behind #105. The
plan needed the corrections above before it could be trusted. Filing it as
aspirational with the caveats baked in, so the near-term work (front-end, #102,
tensor history, KKT) is done first and the reduction is built last, correctly.

Related

• Epic: generated mixed-type return-map materials with consistent tangent (graph-coupled) #102 (coupled vector/mixed-type return-map — the fallback + reference)
• Verification spine: pull numerical-correctness verification (FD tangent + golden model + gates) to the front #90 (verification spine / FD net — the safety mechanism, with its limits)
• NumSimMaterialTarget: next increments (inputs, output + tangent-block properties, tensor state) + PropertySpec emission seam #92 (NumSimMaterialTarget increments — tensor history lands here)
• Formalize Target backend factory (now) + parameterization framework (deferred north star, post-verification) #93 (parameterization descope — the over-engineering precedent this must respect)
• Lower spectral nodes with a shared decomposition; MOOSE spectral-material demo #105 (spectral / isotropic tensor functions — the co-axiality reasoning that
  broadens what P2 can prove)
• numsim-materials#11 (Mandel), MOOSE backend: add runtime parity test against mock RankTwoTensor #12 (vector Newton solver)