MIRK4 Swirling Flow III fails on default nonlinear solver with ForwardDiff Dual conversion

[ChrisRackauckas-Claude]

Reproduction

Clean checkout of SciML/BoundaryValueDiffEq.jl at current origin/master:

git clone https://github.com/SciML/BoundaryValueDiffEq.jl.git BoundaryValueDiffEq-master-failure
cd BoundaryValueDiffEq-master-failure
git fetch origin
git reset --hard origin/master

Observed HEAD:

36e9daa5f64a7e2a45ee0256559e739f027d6f8d
36e9daa5 Merge pull request #534 from ChrisRackauckas-Claude/fix-mirk-constrained-interpolation-533

Run the minimal version of lib/BoundaryValueDiffEqMIRK/test/Core/mirk_basic_tests.jl testset Swirling Flow III:

timeout 3600 /home/crackauc/.juliaup/bin/julia +1.10 --project=lib/BoundaryValueDiffEqMIRK -e 'using Pkg; Pkg.instantiate(); using BoundaryValueDiffEqMIRK; eps=0.01; function swirling_flow!(du,u,p,t); eps=p; du[1]=u[2]; du[2]=(u[1]*u[4]-u[3]*u[2])/eps; du[3]=u[4]; du[4]=u[5]; du[5]=u[6]; du[6]=(-u[3]*u[6]-u[1]*u[2])/eps; return; end; function swirling_flow_bc!(res,u,p,t); res[1]=u(0.0)[1]+1.0; res[2]=u(0.0)[3]; res[3]=u(0.0)[4]; res[4]=u(1.0)[1]-1.0; res[5]=u(1.0)[3]; res[6]=u(1.0)[4]; return; end; prob=BVProblem(swirling_flow!, swirling_flow_bc!, zeros(6), (0.0,1.0), eps); sol=solve(prob, MIRK4(); dt=0.01, abstol=1.0e-4); @show sol.retcode'

Result

This fails on clean current master with:

ERROR: MethodError: no method matching Float64(::ForwardDiff.Dual...)

The stack reaches BoundaryValueDiffEqMIRK.__perform_mirk_iteration at lib/BoundaryValueDiffEqMIRK/src/mirk.jl:334, through NonlinearSolveBase.__solve. The displayed type is enormous because the ForwardDiff.Dual tag captures the MIRK cache.

The failing resolved test environment included:

BoundaryValueDiffEqMIRK v1.17.1
ForwardDiff v1.4.1
DifferentiationInterface v0.7.18
NonlinearSolveBase v2.33.0
NonlinearSolveFirstOrder v2.1.3
SciMLBase v3.31.0
Julia 1.10.11

Additional diagnostics

The same problem succeeds if the internal nonlinear solver is explicitly set to NewtonRaphson() through the package namespace:

sol = solve(prob, MIRK4(; nlsolve = BoundaryValueDiffEqMIRK.NewtonRaphson()); dt = 0.01, abstol = 1.0e-4)
# sol.retcode = SciMLBase.ReturnCode.Success

I tried the narrow source change of making the MIRK/MIRK6I struct default match the docstring (nlsolve::N = NewtonRaphson() instead of nothing). That fixes Swirling Flow III, but it is not a valid fix: running the MIRK Core test group then produced a regression in MIRK Basic Tests, specifically Test interpolant evaluation with big defect at lib/BoundaryValueDiffEqMIRK/test/Core/mirk_basic_tests.jl:449 (SciMLBase.successful_retcode(sol) failed), with summary before interruption:

MIRK Basic Tests | 157 passed, 1 failed, 1 errored
Swirling Flow III | 1 passed
Test interpolant evaluation with big defect | 2 passed, 1 failed

Partial history check

The current failure reproduces on older recent commits too:

• 2977825c (MIRK: cache the flat nlprob u0 buffer instead of reallocating each iteration): same Float64(::ForwardDiff.Dual...) failure.
• 05541438 (Fix RAT related global error control): same Float64(::ForwardDiff.Dual...) failure.

Older comparable points were not clean bisect endpoints under today’s resolver:

• 20c40047 predates the current callable boundary-condition/subpackage setup, so the current reproducer is not comparable there.
• e0c6dafb has the current test shape but fails earlier with a PreparationMismatchError under the current dependency resolver.
• e36034da failed to load under the current resolver with UndefVarError: __concrete_nonlinearsolve_algorithm not defined.

So I did not complete a definitive bisect to a single introducing commit. The actionable narrowing is that the default nlsolve = nothing path selects the BVP-compatible nonlinear polyalgorithm; the reported case fails in its default path, while explicit NewtonRaphson() succeeds for this case but is too broad as a default replacement because it regresses another MIRK Basic test.