Mixed-rank coupled LocalNewtonSystem (scalar + vector + tensor unknowns in one solve)

[petlenz]

Summary

Generalize the local Newton solver to a single coupled system over heterogeneous, mixed-rank unknowns — scalars, vectors, and tensors solved together. This is the full LocalNewtonSystem IR (deferred from Phase 2.4); the scalar-only coupled solver in PR #83 (3b-2b) is its first slice.

Motivation

Real local constitutive return-mapping rarely solves for one kind of quantity. Typical coupled local systems mix ranks:

• J2, isotropic + kinematic hardening: Δγ (scalar) + back-stress X (tensor) + α (scalar).
• Crystal plasticity: a vector of slip rates γ̇ⱼ (+ scalars).
• Damage–plasticity: damage d (scalar) + plastic strain εᵖ (tensor) + hardening (scalar).
• Gradient / micromorphic models: scalar + vector + tensor micro-variables.

All such unknowns are mutually coupled and must be solved in one monolithic Newton iteration, not staggered (unless the author chooses to — see the auto-coupling policy in #83).

Design: flatten → block Jacobian → dense solve → scatter

A mixed-rank coupled Newton is solved by flattening all unknowns into one component vector x ∈ ℝ^M:

Unknown	Components
scalar (Δγ, α, d)	1
vector / rank-1 (slip rates, directors)	dim
symmetric rank-2 (εᵖ, sym. X)	6 (Mandel)
non-symmetric rank-2 (F, full X)	9

Each Newton iteration: assemble residual R ∈ ℝ^M; assemble the block Jacobian J ∈ ℝ^{M×M}, where J_ab = ∂R_a/∂x_b is a sub-block sized by the two unknowns' component counts (scalar-scalar 1×1, scalar-tensor 1×6, tensor-tensor 6×6, …); solve the dense M×M; scatter the update back to the typed unknowns. This is the classic computational-plasticity return-mapping residual vector.

Flattening convention — Mandel, not Voigt

• Symmetric rank-2 → Mandel-6 (√2 on off-diagonals): orthonormal, so the 6-vector / 6×6 block are isometric to the tensors, the Jacobian block is symmetric iff the operator is, and there is no stress-vs-strain factor-of-2 bookkeeping. Voigt's engineering-shear asymmetry is a notorious bug source — worse in machine-generated index math.
• Non-symmetric rank-2 → 9-component flat (Voigt/Mandel-6 is invalid — it drops the skew part).
• The representation is chosen automatically per unknown from tensor_space (codegen already tracks symmetry — see the PR Minor-symmetric consistent tangents (#35 follow-up) #81 minor-symmetry work where roles::Strain carries Symmetric).
• Internal-solve representation only — the tangent/output stays full-tensor (current behavior).

What exists vs what's new

• The LinearAlgebraEmitter boundary is already correct (PR Phase 3b-2b: coupled vector Newton (#35) #83): the dense solve is rank-agnostic — it solves an M×M of doubles. Mixed-rank handling lives above it (the IR + flatten/scatter + block-diff), not inside the solve. No change to that seam.
• 3b-2b is the scalar-only first slice: today NewtonSystem holds vector<string> unknowns (all scalar) and a scalar N×N Jacobian. The generalization replaces this with typed unknowns (rank / dim / symmetry), per-unknown component offsets, and a block-assembled Jacobian. (The NewtonSystem doc-comment already marks "scalar-unknowns only — 3b-2d needs a sibling type / block generalization".)

Dependencies

The off-diagonal Jacobian blocks need exactly the upstream CAS diff overloads already on the critical path:

• ∂(tensor residual)/∂(scalar unknown) → numsim-cas#275 diff(tensor, scalar)
• ∂(scalar/t2s residual)/∂(tensor unknown) → numsim-cas#285 diff(t2s, scalar) (+ the existing diff(t2s, tensor))
• ∂(tensor)/∂(tensor) → rank-4, exists
• building tensor-tensor blocks may use rank-4 tensor_inv → tensor_inv rank-4 emit (blocks on numsim-cas exposing the contraction-index pair) #43 / numsim-cas#276

Note: flattening to a real M×M solve sidesteps a literal rank-4 tensor inverse (you invert the assembled real matrix, not a rank-4 tensor), but rank-4 diff is still needed to build the tensor-tensor blocks.

Scope

• LocalNewtonSystem IR: an ordered list of typed unknowns (name, rank, dim, symmetry) + residual expressions; per-unknown component count + offset.
• Flatten/scatter codegen: tensor↔component mapping (Mandel-6 / 9-flat / dim-vector), chosen from tensor_space.
• Block Jacobian assembly: per (a,b) pair, diff(R_a, x_b) at the right rank, flattened into the J(offset_a.., offset_b..) block.
• LocalNewtonLoweringPass: group coupled mixed-rank equations (generalize the current scalar union-find), emit one M×M dense solve via the existing LinearAlgebraEmitter.
• Interplay with the divergence/robustness work (Local-Newton robustness: detect non-convergence + guard singular Jacobian (scalar + coupled) #85): the M×M solve has the same singular/non-convergence concerns.

Acceptance

• A recipe with a coupled {scalar, symmetric-tensor} system (e.g. {Δγ, εᵖ}) emits a single dense Newton over the flattened 1+6=7 vector and converges (numerically verified, mirroring the 3b-2b scalar compile-and-run test).
• A {scalar, vector} system (crystal-plasticity-like) flattens correctly.
• A non-symmetric tensor unknown flattens to 9 components, not 6.
• Existing scalar-only coupled recipes stay byte-identical.

Related

• PR Phase 3b-2b: coupled vector Newton (#35) #83 — 3b-2b scalar coupled vector Newton (the first slice; LinearAlgebraEmitter, auto-coupling).
• numsim-cas#275 / #285 — the mixed-rank diff overloads this needs.
• tensor_inv rank-4 emit (blocks on numsim-cas exposing the contraction-index pair) #43 / numsim-cas#276 — rank-4 tensor_inv (tensor-tensor blocks).
• Local-Newton robustness: detect non-convergence + guard singular Jacobian (scalar + coupled) #85 — local-Newton divergence/singularity robustness (applies to the M×M solve too).
• Scalar-type abstraction in generated code (enable float / fixed-point / autodiff) #84 — scalar-type abstraction (orthogonal; both touch LinearAlgebraEmitter).
• Phase 2.4 LocalNewtonSystem IR (originally deferred) / Phase 3b-2 (Phase 3b — Algorithmic tangent + Kuhn-Tucker lowering #35).