[WIP] Baseline Multi-Level Newton-Raphson (MLN)

[termi-official]

This PR is ment to be the foundation to discuss how to proceed, it is not ment to be merged as is

A brief summary on this PR is best given in the docs for the new multilevel subpackage.
This PR contains:

• a new solver (MultiLevelNewton)
• a new function (MultiLevelNonlinearFunction)
• necessary gearing for the MultiLevelNonlinearFunction to be used with the regular Newton+postcondition -- this is what most people do right now in practice

MLN as a nonlinear elimination method

I am not sure if the nonlinear preconditioning system is incomplete at this point, but I could not come up with a prototype yet. Framing the MLN as a nonlinear elimination method does not seem to be a problem per-se, but something in the framework not gearing correctly.

Fable's gives the following response:

Short version: "nonlinear elimination via precondition" is a category mismatch, elimination via postcondition is exactly our Variant-A arm, and the condensed solver won the primary spot because it's the only formulation where the rest of the framework stays both unmodified and correct.

Why not precondition (left hook, G(fu, u, p)). The left preconditioner is a residual reshaping: the solver still iterates on the full u = [ū; q] and still takes Newton steps on the full-space Jacobian of the composed map. Elimination is not a residual reshaping — it's a state correction (replace the iterate's q by q*(ū)), and G has no way to move the iterate. You could try to smuggle the condensed residual into G by having it solve the local problems internally, but then (a) G receives f(u) evaluated at the unprojected state and must throw it away — double work per evaluation; (b) the Jacobian machinery would differentiate through the inner local Newton solves, which is precisely the cost and semantics MLN exists to avoid (the whole point is that you assemble S analytically from per-point correctors); and (c) the review found that precondition + analytic jac was silently inconsistent upstream — the solver differentiates the uncomposed residual — which we turned into a hard error. MLN mandates a user-assembled S, so a precondition-based MLN would sit exactly on the refused combination. There's no formulation of elimination as a left preconditioner that isn't either wrong or a worse-dressed Variant B.

Elimination via postcondition (right hook) we did build — that's the full-space arm: MultiLevelProjection is the commit step expressed as the corrector H, plus the two pieces the hook alone can't supply (the Schur operator as the Jacobian handle, and the δq-zeroing CondensedFactorization, because the descent contract only returns full-length steps). Its structural ceiling is why it's not the default: H fires only at accepted iterates, so line-search merits and trust-region ratios are evaluated at unprojected states — measured consequence: the merit is monotonically worse in α (Armijo can fail on well-posed problems) and Dogleg's Cauchy leg bypasses the condensed solver entirely. Hence its :None-globalization restriction, and no status channel for local failure (H returns only a corrected iterate).

Why the condensed Variant B is the primary solver. Putting the local solves inside the residual makes the global problem a square, ordinary NonlinearProblem over ū — so termination modes, line search, trust region, Eisenstat–Walker forcing, linear precs on S, tracing, and the recompute_jacobian reuse signal all work unmodified and consistently, because every trial point the globalization probes is a correctly eliminated state. Failure gets real channels (Inf at trials → backtracking; commit failure → ConvergenceFailure, which later maps to a dt-cut on the ODE path). Phase 0 proved the premise empirically: stock NewtonRaphson solves the condensed problem at order 1.99, so the packaged solver is ~40 lines of orchestration around one inner cache rather than a new Newton implementation. And it keeps your taxonomy constraint clean by construction: linear preconditioning of S lives on the global solver's LinearSolve algorithm, untouched by — and composable with — the nonlinear hooks, which remain free for actual nonlinear preconditioning on top of MLN.

Missing Things

(Ordered by priority, most important first)

• Proper error+stats communication channel (vector of length num_local_problems)
• Pass outer function parameters into inner residual evaluation
• Inner solver configuration
• Adding proper citations/attributions
• Matrix-free Schur form
• Implicit AD integration
• Check that the local solutions can vary in their size (e.g. first local problem has 100 unknowns, the second local problem has 110, the third one 93, ...).

[termi-official]

I have isolated this from FerriteDiffEq.jl, see https://github.com/termi-official/FerriteDiffEq.jl/blob/master/src/ferrite_newton.jl and adjacent files.