Fix type instability in solve pathway (#454)

[ChrisRackauckas-Claude]

Summary

Addresses #454 by adding function barriers at key points in the solve chain where type instability propagates.

Changes

BoundaryValueDiffEqCore:

• Add Val-dispatched function barrier in __internal_nlsolve_problem for the nlls=Nothing case, preventing Union{NonlinearProblem, NonlinearLeastSquaresProblem} from propagating

MIRK/FIRK/MIRKN/Ascher/Shooting solvers:

• Extract internal problem solve into __solve_internal_problem barrier functions, ensuring the algorithm selection and solve happens in a type-stable context

Tests:

• Add BVProblem Constructor test section verifying BVProblem{true}(...), BVProblem{false}(...), TwoPointBVProblem{true/false}(...) are @inferred-clean
• Add __init type stability test

Type stability status

Component	Status	Notes
BVProblem{true}(f!, bc!, u0, tspan)	Type-stable	Explicit iip avoids runtime isinplace check
BVProblem(f!, bc!, u0, tspan)	Type-unstable	Same as ODEProblem — isinplace is runtime. Needs upstream SciMLBase fix
solve(prob, MIRK5(; jac_alg=...); nlls=Val(false))	Type-stable	With explicit jac_alg and nlls
solve(prob, MIRK5(); dt=...)	Type-unstable	Default AutoSparse jac_alg has complex inference

Upstream work needed

The convenience constructor BVProblem(f, bc, u0, tspan) type instability lives in SciMLBase (same pattern as ODEProblem). A companion PR to SciMLBase with function barriers for the BVP convenience constructors would complete the fix. The BoundaryValueDiffEq changes in this PR are independent and beneficial regardless.

Test plan

• Existing type stability tests pass (@inferred solve(...) with explicit jac_alg and nlls=Val(false))
• New constructor type stability tests pass
• Issue Possible type instability in BVP.BVProblem() and BVP.solve() #454 reproducer works correctly
• CI

🤖 Generated with Claude Code

[ChrisRackauckas-Claude]

Implementation notes

What this PR fixes:
The solve pathway had type instabilities from two sources:

1. __internal_nlsolve_problem returning Union{NonlinearProblem, NonlinearLeastSquaresProblem} when nlls=Nothing (the default)
2. The algorithm selection and solve cascading from that union through the entire solve chain

Function barriers at these points ensure each branch compiles to type-stable code. With nlls=Val(false) and an explicit jac_alg, the full solve call is @inferred-clean.

Upstream SciMLBase changes (not included):
I also have function barriers ready for SciMLBase's BVProblem, BVPFunction, SecondOrderBVProblem, and DynamicalBVPFunction convenience constructors. These add Val{iip}() dispatch so BVPFunction{iip}(f, bc) construction happens in a type-stable context. Should I open a companion PR to SciMLBase?

Workaround for users (no code changes needed):
Users who need type stability today can use:

prob = BVProblem{true}(f!, bc!, u0, tspan; nlls=Val(false))
sol = solve(prob, MIRK5(; jac_alg=BVPJacobianAlgorithm(AutoForwardDiff(; chunksize=N))); dt=h)

[github-actions]

Benchmark Results

Click to check benchmark results

	master	a893ce8...	master / a893ce8...
Simple Pendulum/IIP/BoundaryValueDiffEqMIRK.MIRK2()	0.591 ± 0.017 s	0.612 ± 0.025 s	0.966 ± 0.048
Simple Pendulum/IIP/BoundaryValueDiffEqMIRK.MIRK3()	12.4 ± 0.56 ms	12.8 ± 0.72 ms	0.971 ± 0.07
Simple Pendulum/IIP/BoundaryValueDiffEqMIRK.MIRK4()	2.39 ± 0.15 ms	2.42 ± 0.16 ms	0.988 ± 0.09
Simple Pendulum/IIP/BoundaryValueDiffEqMIRK.MIRK5()	2.89 ± 0.27 ms	3 ± 0.5 ms	0.962 ± 0.18
Simple Pendulum/IIP/BoundaryValueDiffEqMIRK.MIRK6()	1.13 ± 0.24 ms	1.16 ± 0.26 ms	0.978 ± 0.3
Simple Pendulum/IIP/MultipleShooting(10, Tsit5; grid_coarsening = false)	2.49 ± 0.74 ms	2.61 ± 0.75 ms	0.954 ± 0.39
Simple Pendulum/IIP/MultipleShooting(10, Tsit5; grid_coarsening = true)	4.89 ± 1.3 ms	5.33 ± 1.3 ms	0.917 ± 0.33
Simple Pendulum/IIP/MultipleShooting(100, Tsit5; grid_coarsening = false)	0.0674 ± 0.0078 s	0.0718 ± 0.012 s	0.939 ± 0.19
Simple Pendulum/IIP/MultipleShooting(100, Tsit5; grid_coarsening = true)	0.0999 ± 0.017 s	0.111 ± 0.026 s	0.901 ± 0.26
Simple Pendulum/IIP/Shooting(Tsit5())	0.301 ± 0.094 ms	0.308 ± 0.097 ms	0.979 ± 0.44
Simple Pendulum/OOP/BoundaryValueDiffEqMIRK.MIRK2()	0.723 ± 0.005 s	0.745 ± 0.013 s	0.971 ± 0.018
Simple Pendulum/OOP/BoundaryValueDiffEqMIRK.MIRK3()	14.9 ± 4.2 ms	15.3 ± 6.2 ms	0.969 ± 0.48
Simple Pendulum/OOP/BoundaryValueDiffEqMIRK.MIRK4()	2.77 ± 0.15 ms	2.78 ± 0.13 ms	0.997 ± 0.07
Simple Pendulum/OOP/BoundaryValueDiffEqMIRK.MIRK5()	3.45 ± 0.2 ms	3.54 ± 0.25 ms	0.974 ± 0.089
Simple Pendulum/OOP/BoundaryValueDiffEqMIRK.MIRK6()	1.33 ± 0.11 ms	1.35 ± 0.15 ms	0.988 ± 0.13
Simple Pendulum/OOP/MultipleShooting(10, Tsit5; grid_coarsening = false)	3.45 ± 3.3 ms	3.84 ± 3.4 ms	0.898 ± 1.2
Simple Pendulum/OOP/MultipleShooting(10, Tsit5; grid_coarsening = true)	6.9 ± 6 ms	7.55 ± 7 ms	0.915 ± 1.2
Simple Pendulum/OOP/MultipleShooting(100, Tsit5; grid_coarsening = false)	0.0819 ± 0.0016 s	0.0861 ± 0.0031 s	0.951 ± 0.039
Simple Pendulum/OOP/MultipleShooting(100, Tsit5; grid_coarsening = true)	0.125 ± 0.0027 s	0.135 ± 0.0056 s	0.925 ± 0.043
Simple Pendulum/OOP/Shooting(Tsit5())	0.599 ± 0.047 ms	0.612 ± 0.048 ms	0.979 ± 0.11
time_to_load	6.69 ± 0.068 s	6.94 ± 0.098 s	0.963 ± 0.017

### Benchmark Plots A plot of the benchmark results has been uploaded as an artifact to the workflow run for this PR. Go to "Actions"->"Benchmark a pull request"->[the most recent run]->"Artifacts" (at the bottom).

[ChrisRackauckas-Claude]

Updated approach — fixing the true instability

After further investigation, the function barrier approach was a band-aid. The true root cause is a one-line bug in SciMLBase.

Root cause

In SciMLBase/src/problems/bvp_problems.jl:189:

_nlls = Nothing # Cannot reliably infer

This is hit for StandardBVProblem without bcresid_prototype. But the solver always creates a residual matching u0 size (via zero(u) in __get_bcresid_prototype), so nlls is definitionally false in this case.

This Nothing type parameter propagates into __internal_nlsolve_problem which returns Union{NonlinearProblem, NonlinearLeastSquaresProblem}, poisoning the entire solve chain with type instability.

Fix

SciMLBase one-liner:

- _nlls = Nothing # Cannot reliably infer
+ _nlls = false   # Without bcresid_prototype, residual matches u0 size

With this fix, @inferred solve(prob, MIRK5(; jac_alg=...); dt=0.2) passes without any nlls=Val(false) workaround.

I've updated this PR to just update the tests (removing the now-unnecessary nlls=Val(false) for StandardBVProblem). The function barriers have been reverted. Will open a companion SciMLBase PR.