Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Add check whether relaxation is better than the full step #2443

Open
wants to merge 1 commit into
base: master
Choose a base branch
from

Conversation

SouthEndMusic
Copy link
Member

@SouthEndMusic SouthEndMusic commented Aug 29, 2024

Fixes #2442

Not sure whether this should be upstreamed to LineSearches.jl

@oscardssmith
Copy link
Contributor

Presumably it should be upstreamed.

@SouthEndMusic
Copy link
Member Author

SouthEndMusic commented Aug 29, 2024

@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
Copy link
Member Author

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
Copy link
Member Author

SouthEndMusic commented Aug 29, 2024

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
Copy link
Contributor

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

@ChrisRackauckas
Copy link
Member

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

@SouthEndMusic
Copy link
Member Author

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
Copy link
Member Author

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
Copy link
Contributor

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

@SouthEndMusic
Copy link
Member Author

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
Copy link
Member Author

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

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

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

[OrdinaryDiffEqNonLinearSolve] relax = BackTracking() leads to worse residual in NLNewton
3 participants