NonlinearSolveAlg DAE stage residual: daenlf/oopdaenlf do not match the _compute_rhs DAEFunction methods

[ChrisRackauckas-Claude]

Summary

Once the nlstep_data field read of #4127 is guarded, NonlinearSolveAlg on a DAEProblem
still cannot take a step: the DAE stage-residual closures daenlf/oopdaenlf do not match the
_compute_rhs!/_compute_rhs DAEFunction methods they call, and the parameter tuple that
feeds them has three different shapes in three different places. init throws a MethodError
at the first residual evaluation, so the NonlinearSolveAlg DAE path has never run.

This is a separate defect from #4127 (which is a field lookup); filing it so the umbrella
complaint there — "NonlinearSolveAlg cannot be used with a DAEProblem at all" — is not lost
when that issue closes.

Reproduction

The Robertson DAE from #4127, on top of the nlstep_data guard:

using OrdinaryDiffEqBDF
using OrdinaryDiffEqNonlinearSolve: NonlinearSolveAlg
using NonlinearSolve: NewtonRaphson
using ADTypes

function daef!(out, du, u, p, t)
    out[1] = -0.04u[1] + 1.0e4 * u[2] * u[3] - du[1]
    out[2] = 0.04u[1] - 3.0e7 * u[2]^2 - 1.0e4 * u[2] * u[3] - du[2]
    out[3] = u[1] + u[2] + u[3] - 1.0
    return nothing
end
prob = DAEProblem(
    daef!, [-0.04, 0.04, 0.0], [1.0, 0.0, 0.0], (0.0, 1.0e5);
    differential_vars = [true, true, false]
)
init(prob, DFBDF(nlsolve = NonlinearSolveAlg(NewtonRaphson(; autodiff = AutoForwardDiff()))))

MethodError: no method matching _compute_rhs!(::Vector{Float64}, ::Vector{Float64},
::Vector{Float64}, ::Float64, ::Float64, ::Float64, ::Vector{Float64}, ::Float64,
::SciMLBase.NullParameters, ::Float64, ::DAEFunction{...}, ::Vector{Float64})
    daenlf @ utils.jl:435
    NonlinearFunction @ scimlfunctions.jl:2857
    evaluate_f @ utils.jl:196
    #__init#4 @ solve.jl:190

DImplicitEuler fails identically. Out-of-place fails the same way one method over, through
oopdaenlf → _compute_rhs.

The mismatches (master, 5cf0839)

1. The in-place closure calls a method that does not exist.
daenlf (lib/OrdinaryDiffEqNonlinearSolve/src/utils.jl:428) is a copy of odenlf with
method dropped:

function daenlf(ztmp, z, p)
    tmp, ustep, γ, α, tstep, k, invγdt, _p, dt, f = p
    return _compute_rhs!(tmp, ztmp, ustep, γ, α, tstep, k, invγdt, _p, dt, f, z)[1]
end

but the DAEFunction methods of _compute_rhs! (newton.jl:726, newton.jl:764) are

_compute_rhs!(tmp, ztmp, ustep, α, tstep, k, invγdt, p, uprev, f::TF, z) where {TF <: DAEFunction}

— eleven arguments, no γ, no dt, and a uprev that the parameter tuple never carries
(the DAE stage state is ustep = uprev + z).

2. Three different parameter shapes for the same tuple.

• build_nlsolver builds ten entries for isdae (utils.jl:805):
  (tmp, ustep, γ, α, tstep, k, invγdt, p, dt, f)
• initialize! rebuilds eleven for f isa DAEFunction (newton.jl:153):
  (tmp, ztmp, ustep, γ, α, tstep, k, invγdt, p, dt, f)
• daenlf destructures ten.

Whatever the right shape is, these three have to agree.

3. The residual never reaches NonlinearSolve's buffer.
Even with the arity fixed, the DAE _compute_rhs! writes the du estimate into its ztmp
argument and the residual into k, returning _vec(k):

@.. broadcast = false ztmp = (tmp + α * z) * invγdt
@.. ustep = uprev + z
f(k, ztmp, ustep, p, tstep)
return _vec(k), ustep

daenlf passes NonlinearSolve's residual buffer as that ztmp, so the buffer would come back
holding the du estimate while the residual sits in k (an in-place NonlinearFunction
ignores the return value). The DAE path needs a scratch buffer for the du estimate that is
distinct from the residual output — NonlinearSolveCache has no spare one today, and
nlsolver.ztmp is not free (compute_step! stores the inner solution there).

4. Out-of-place passes the wrong p plus an extra argument.
oopdaenlf (utils.jl:946) forwards the whole parameter tuple as p and keeps dt:

function oopdaenlf(z, p)
    tmp, α, tstep, invγdt, _p, dt, uprev, f = p
    return _compute_rhs(tmp, α, tstep, invγdt, p, dt, uprev, f, z)[1]   # p, and dt
end

against _compute_rhs(tmp, α, tstep, invγdt, p, uprev, f::TF, z) (newton.jl:680). This one
is a one-line fix (_p, drop dt), unlike the in-place path.

5. The convergence measure is the ODE one.
compute_step!(::NLSolver{<:NonlinearSolveAlg, true}, integrator) forms
ustep = compute_ustep!(ustep, tmp, γ, z, method) (i.e. tmp + γz, or z for
COEFFICIENT_MULTISTEP), while the DAE stage state is uprev + z. calculate_residuals! is
weighted by that ustep, so the outer convergence test would be scaled wrongly even once the
residual is correct. Compare compute_step!(::NLSolver{<:NLNewton, true}, integrator, γW),
which branches on isdae for exactly this.

Notes

• use_w_reuse is !isdae && ..., so cache.W === nothing on the DAE path and
  _update_nlsolvealg_W! is not reached — worth remembering, since it reads f.mass_matrix,
  which a DAEFunction also does not have.
• The isdae branches (daenlf/oopdaenlf, DAEResidualJacobianWrapper, the isdae guard on
  use_w_reuse) show the path is meant to be supported, so this looks like never-exercised
  plumbing rather than an intentional restriction. If it is not meant to be supported, the
  cleaner fix is an ArgumentError from build_nlsolver for isdae, the way
  HomotopyNonlinearSolveAlg already rejects DAEs.

/cc @ChrisRackauckas