Add check whether relaxation is better than the full step

[SouthEndMusic]

Fixes #2442

Not sure whether this should be upstreamed to LineSearches.jl

[oscardssmith]

Presumably it should be upstreamed.

[SouthEndMusic]

@oscardssmith I still feel some work on this side is required as well. Most line-search algorithms from LineSearches.jl error when ϕ'(0) < 0 is not satsiified, but that shouldn't stop the ODE solver. The relaxation can be refactored so that line-search based relaxation is only performed when that condition is met. I think that fixes the problem I'm encountering. Then upstream this check should be added as wel to BackTracking where it is currently missing.

[SouthEndMusic]

Hm, on closer inspection the condition for BackTracking is somewhat more complicated, so it can not be done that generically. I suppose putting the line-search in try-catch is not desirable?

[SouthEndMusic]

Maybe the method for computing the derivative of the residual with respect to the relative step size is just inaccurate for my particular use case, which would justify my current version of this PR again.

[oscardssmith]

one other question: If you use the (still somewhat undocumented NonlinearSolveAlg) does this bug still appear?

[ChrisRackauckas]

We really should just be using NonlinearSolveAlg and pushing that forward rather than small improvements to the other.

[SouthEndMusic]

I'm trying to run the ODE solver with nlsolve = NonlinearSolveAlg(NewtonRaphson(; linesearch = BackTracking(), autodiff = AutoForwardDiff()) but I get an error because my rhs is called with eltype(du) = Float64 and eltype(u) = ForwardDiff.Dual{...}.

[SouthEndMusic]

Stacktrace:

ERROR: MethodError: no method matching water_balance!(::ComponentArrays.ComponentVector{…}, ::ComponentArrays.ComponentVector{…}, ::Ribasim.Parameters{…}, ::Float64)

Closest candidates are:
  water_balance!(::V, ::V, ::Ribasim.Parameters, ::Number) where V<:(ComponentArrays.ComponentVector)   @ Ribasim c:\Users\konin_bt\Ribasim_development\Ribasim\core\src\solve.jl:4

Stacktrace:
  [1] (::SciMLBase.ODEFunction{…})(::ComponentArrays.ComponentVector{…}, ::Vararg{…})
    @ SciMLBase C:\Users\konin_bt\.julia\packages\SciMLBase\vhP5T\src\scimlfunctions.jl:2299
  [2] _compute_rhs!(tmp::ComponentArrays.ComponentVector{…}, ztmp::ComponentArrays.ComponentVector{…}, ustep::ComponentArrays.ComponentVector{…}, γ::Float64, α::Float64, tstep::Float64, k::ComponentArrays.ComponentVector{…}, invγdt::Float64, method::OrdinaryDiffEq.MethodType, p::Ribasim.Parameters{…}, 
dt::Float64, f::SciMLBase.ODEFunction{…}, z::ComponentArrays.ComponentVector{…})
    @ OrdinaryDiffEq C:\Users\konin_bt\.julia\packages\OrdinaryDiffEq\s27pa\src\nlsolve\newton.jl:330  [3] (::OrdinaryDiffEq.var"#178#181")(ztmp::ComponentArrays.ComponentVector{…}, z::ComponentArrays.ComponentVector{…}, p::Tuple{…})
    @ OrdinaryDiffEq C:\Users\konin_bt\.julia\packages\OrdinaryDiffEq\s27pa\src\nlsolve\utils.jl:220 
  [4] NonlinearFunction
    @ C:\Users\konin_bt\.julia\packages\SciMLBase\vhP5T\src\scimlfunctions.jl:2300 [inlined]
  [5] JacobianWrapper
    @ C:\Users\konin_bt\.julia\packages\SciMLBase\vhP5T\src\function_wrappers.jl:108 [inlined]       
  [6] vector_mode_dual_eval!
    @ C:\Users\konin_bt\.julia\packages\ForwardDiff\PcZ48\src\apiutils.jl:31 [inlined]
  [7] vector_mode_jacobian!(result::Matrix{…}, f!::SciMLBase.JacobianWrapper{…}, y::ComponentArrays.ComponentVector{…}, x::ComponentArrays.ComponentVector{…}, cfg::ForwardDiff.JacobianConfig{…})        
    @ ForwardDiff C:\Users\konin_bt\.julia\packages\ForwardDiff\PcZ48\src\jacobian.jl:153
  [8] jacobian!
    @ C:\Users\konin_bt\.julia\packages\ForwardDiff\PcZ48\src\jacobian.jl:78 [inlined]
  [9] jacobian!
    @ C:\Users\konin_bt\.julia\packages\ForwardDiff\PcZ48\src\jacobian.jl:76 [inlined]
 [10] sparse_jacobian!
    @ C:\Users\konin_bt\.julia\packages\SparseDiffTools\hj4gQ\src\highlevel\forward_mode.jl:65 [inlined]
 [11] JacobianCache
    @ C:\Users\konin_bt\.julia\packages\NonlinearSolve\82B4C\src\internal\jacobian.jl:138 [inlined]  
 [12] JacobianCache
    @ C:\Users\konin_bt\.julia\packages\NonlinearSolve\82B4C\src\internal\jacobian.jl:115 [inlined]  
 [13] __step!(cache::NonlinearSolve.GeneralizedFirstOrderAlgorithmCache{…}; recompute_jacobian::Nothing, kwargs::@Kwargs{})
    @ NonlinearSolve C:\Users\konin_bt\.julia\packages\NonlinearSolve\82B4C\src\core\generalized_first_order.jl:219
 [14] __step!
    @ C:\Users\konin_bt\.julia\packages\NonlinearSolve\82B4C\src\core\generalized_first_order.jl:215 
[inlined]
 [15] #step!#218
    @ C:\Users\konin_bt\.julia\packages\NonlinearSolve\82B4C\src\core\generic.jl:50 [inlined]        
 [16] step!
    @ C:\Users\konin_bt\.julia\packages\NonlinearSolve\82B4C\src\core\generic.jl:45 [inlined]        
 [17] compute_step!(nlsolver::OrdinaryDiffEq.NLSolver{…}, integrator::OrdinaryDiffEq.ODEIntegrator{…})
    @ OrdinaryDiffEq C:\Users\konin_bt\.julia\packages\OrdinaryDiffEq\s27pa\src\nlsolve\newton.jl:106 [18] nlsolve!(nlsolver::OrdinaryDiffEq.NLSolver{…}, integrator::OrdinaryDiffEq.ODEIntegrator{…}, cache::OrdinaryDiffEq.OrdinaryDiffEqBDF.QNDFCache{…}, repeat_step::Bool)
    @ OrdinaryDiffEq C:\Users\konin_bt\.julia\packages\OrdinaryDiffEq\s27pa\src\nlsolve\nlsolve.jl:54 [19] perform_step!(integrator::OrdinaryDiffEq.ODEIntegrator{…}, cache::OrdinaryDiffEq.OrdinaryDiffEqBDF.QNDFCache{…}, repeat_step::Bool)
    @ OrdinaryDiffEq.OrdinaryDiffEqBDF C:\Users\konin_bt\.julia\packages\OrdinaryDiffEq\s27pa\lib\OrdinaryDiffEqBDF\src\bdf_perform_step.jl:912
 [20] perform_step!
    @ C:\Users\konin_bt\.julia\packages\OrdinaryDiffEq\s27pa\lib\OrdinaryDiffEqBDF\src\bdf_perform_step.jl:847 [inlined]
 [21] solve!(integrator::OrdinaryDiffEq.ODEIntegrator{…})
    @ OrdinaryDiffEq C:\Users\konin_bt\.julia\packages\OrdinaryDiffEq\s27pa\src\solve.jl:557
 [22] solve!(model::Ribasim.Model{OrdinaryDiffEq.ODEIntegrator{…}})
    @ Ribasim c:\Users\konin_bt\Ribasim_development\Ribasim\core\src\model.jl:237
 [23] macro expansion
    @ .\timing.jl:279 [inlined]
 [24] top-level scope
    @ c:\Users\konin_bt\Ribasim_development\runners\runner.jl:26

[oscardssmith]

oh good point, I need to fix that...

[SouthEndMusic]

Looking at how NonlinearSolve.jl deals with this here, that looks a lot like what I propose above. So indeed, using NonlinearSolveAlg will probably solve the problem. I'll close the issue + PR once I have been able to test it.

[SouthEndMusic]

@oscardssmith btw, currently my run fails because remake is not defined here:

OrdinaryDiffEq.jl/lib/OrdinaryDiffEqNonlinearSolve/src/newton.jl

Line 71 in 1eef9db

new_prob = remake(cache.prob, p = nlp_params, u0 = z)