BoundaryValueDiffEqFIRK: Limit function range explored by jac_alg = BVPJacobianAlgorithm(AutoFiniteDiff())?

[caMi11er]

After updating to latest Julia 1.11.4 and BoundaryValueDiffEq v5.16.0, for the first time I see a glimmer of hope for BoundaryValueDiffEqFIRK, [despite remaining unavailability of solution interpolation of first derivative (in bc! function) , and rejection of RadauIIa5(autodiff=false or ... =AutoFiniteDiff().] (hint, hint!)
It turns out that the "jac_alg" fix that EricQQY [provided] (#304 (comment)), for the lastest Rodas4P solver also enables solver RadauIIa5 to run now:

prob = BVProblem{true}(BVPFunction(shtructure!, bcond2!; bcresid_prototype = zeros(5) ), initguess, rspan, (eos, rspan, Pr0, pf) ) 

sol = solve(prob, RadauIIa5(jac_alg = BVPJacobianAlgorithm(AutoFiniteDiff())); dt = 0.5, dtmax=0.5, callback=cb, abstol=1e-4, reltol=1e-4, verbose=false) 

But it still isn't useful yet, because of continuous warning messages from IntervalNonlinearSolve, which is called by my derivative function. These warnings arise from the limited range in function space in which it is known or defined. The problem seems to be that the function finite_difference_jacobian! is exploring values outside the available range:

  [1] error(s::String)
    @ Base ./error.jl:35
  [2] #15
    @ ./REPL[67]:3 [inlined]
  [3] handle_message(::TransformerLogger{TransformerLogger{TransformerLogger{…}, var"#13#14"}, var"#15#16"}, ::LogLevel, ::Vararg{Any}; kwargs::@Kwargs{})
    @ LoggingExtras ~/.julia/packages/LoggingExtras/cFgEq/src/CompositionalLoggers/transformer.jl:22
  [4] handle_message(::TransformerLogger{TransformerLogger{…}, var"#15#16"}, ::LogLevel, ::String, ::Module, ::Symbol, ::Symbol, ::String, ::Int64)
    @ LoggingExtras ~/.julia/packages/LoggingExtras/cFgEq/src/CompositionalLoggers/transformer.jl:20
  [5] #invokelatest#2
    @ ./essentials.jl:1055 [inlined]
  [6] invokelatest
    @ ./essentials.jl:1052 [inlined]
  [7] macro expansion
    @ ./logging/logging.jl:353 [inlined]
  [8] solve(::IntervalNonlinearProblem{…}, ::ITP{…}; maxiters::Int64, abstol::Nothing, verbose::Bool, kwargs::@Kwargs{})
    @ BracketingNonlinearSolve ~/.julia/packages/BracketingNonlinearSolve/1Lfiz/src/itp.jl:86
  [9] solve
    @ ~/.julia/packages/BracketingNonlinearSolve/1Lfiz/src/itp.jl:59 [inlined]
 [10] nb
    @ ~/Applications/Julia/RealisticSeismology-main/src/sim.jl:153 [inlined]
 [11] ε(sim::Sim{…}, p::Float64)
    @ Main ~/Applications/Julia/RealisticSeismology-main/scripts/preprofile.jl:22
 [12] shtructure!(dy_dr::Vector{Float64}, y::Vector{Float64}, param::Tuple{Sim{…}, Tuple{…}, Float64, Float64}, r::Float64)
    @ Main ~/Applications/Julia/RealisticSeismology-main/src/RealMod/src/shtarMfreeFIRK.jl:59
 [13] BVPFunction
    @ ~/.julia/packages/SciMLBase/LGLmH/src/scimlfunctions.jl:2524 [inlined]
 [14] Φ!(residual::Vector{…}, fᵢ_cache::BoundaryValueDiffEqCore.FakeDiffCache{…}, k_discrete::Vector{…}, f!::BVPFunction{…}, TU::BoundaryValueDiffEqFIRK.FIRKTableau{…}, y::Vector{…}, u::Vector{…}, p::Tuple{…}, mesh::Vector{…}, mesh_dt::Vector{…}, stage::Int64)
    @ BoundaryValueDiffEqFIRK ~/.julia/packages/BoundaryValueDiffEqFIRK/xca6M/src/collocation.jl:33
 [15] Φ!
    @ ~/.julia/packages/BoundaryValueDiffEqFIRK/xca6M/src/collocation.jl:2 [inlined]
 [16] __firk_loss_collocation!(resid::Vector{…}, u::Vector{…}, p::Tuple{…}, y::Vector{…}, mesh::Vector{…}, residual::Vector{…}, cache::BoundaryValueDiffEqFIRK.FIRKCacheExpand{…})
    @ BoundaryValueDiffEqFIRK ~/.julia/packages/BoundaryValueDiffEqFIRK/xca6M/src/firk.jl:759
 [17] #117
    @ ~/.julia/packages/BoundaryValueDiffEqFIRK/xca6M/src/firk.jl:411 [inlined]
 [18] FixTail
    @ ~/.julia/packages/DifferentiationInterface/7eD1K/src/utils/context.jl:169 [inlined]
 [19] finite_difference_jacobian!(J::Matrix{…}, f::DifferentiationInterface.FixTail{…}, x::Vector{…}, cache::FiniteDiff.JacobianCache{…}, f_in::Nothing; relstep::Float64, absstep::Float64, colorvec::UnitRange{…}, sparsity::Nothing, dir::Bool)
    @ FiniteDiff ~/.julia/packages/FiniteDiff/EBPBu/src/jacobians.jl:444
 [20] finite_difference_jacobian! (repeats 2 times)
    @ ~/.julia/packages/FiniteDiff/EBPBu/src/jacobians.jl:406 [inlined]
 [21] jacobian
    @ ~/.julia/packages/DifferentiationInterface/7eD1K/ext/DifferentiationInterfaceFiniteDiffExt/twoarg.jl:375 [inlined]
 [22] __construct_nlproblem(cache::BoundaryValueDiffEqFIRK.FIRKCacheExpand{…}, y::Vector{…}, loss_bc::BoundaryValueDiffEqFIRK.var"#115#121"{…}, loss_collocation::BoundaryValueDiffEqFIRK.var"#117#123"{…}, loss::BoundaryValueDiffEqFIRK.var"#119#125"{…}, ::SciMLBase.StandardBVProblem)
    @ BoundaryValueDiffEqFIRK ~/.julia/packages/BoundaryValueDiffEqFIRK/xca6M/src/firk.jl:483
 [23] __construct_nlproblem
    @ ~/.julia/packages/BoundaryValueDiffEqFIRK/xca6M/src/firk.jl:426 [inlined]
 [24] __perform_firk_iteration(cache::BoundaryValueDiffEqFIRK.FIRKCacheExpand{…}, abstol::Float64, adaptive::Bool; nlsolve_kwargs::@NamedTuple{}, kwargs::@Kwargs{…})
    @ BoundaryValueDiffEqFIRK ~/.julia/packages/BoundaryValueDiffEqFIRK/xca6M/src/firk.jl:349
 [25] __perform_firk_iteration
    @ ~/.julia/packages/BoundaryValueDiffEqFIRK/xca6M/src/firk.jl:347 [inlined]
 [26] solve!(cache::BoundaryValueDiffEqFIRK.FIRKCacheExpand{…})
    @ BoundaryValueDiffEqFIRK ~/.julia/packages/BoundaryValueDiffEqFIRK/xca6M/src/firk.jl:302
 [27] __solve(::BVProblem{…}, ::RadauIIa5{…}; kwargs::@Kwargs{…})
    @ BoundaryValueDiffEqCore ~/.julia/packages/BoundaryValueDiffEqCore/Q7n7W/src/BoundaryValueDiffEqCore.jl:37
 [28] __solve
    @ ~/.julia/packages/BoundaryValueDiffEqCore/Q7n7W/src/BoundaryValueDiffEqCore.jl:34 [inlined]
 [29] #solve_call#35
    @ ~/.julia/packages/DiffEqBase/zYZst/src/solve.jl:635 [inlined]
 [30] solve_call
    @ ~/.julia/packages/DiffEqBase/zYZst/src/solve.jl:592 [inlined]
 [31] #solve_up#44
    @ ~/.julia/packages/DiffEqBase/zYZst/src/solve.jl:1142 [inlined]
 [32] solve_up
    @ ~/.julia/packages/DiffEqBase/zYZst/src/solve.jl:1120 [inlined]
 [33] solve(prob::BVProblem{…}, args::RadauIIa5{…}; sensealg::Nothing, u0::Nothing, p::Nothing, wrap::Val{…}, kwargs::@Kwargs{…})
    @ DiffEqBase ~/.julia/packages/DiffEqBase/zYZst/src/solve.jl:1057
 [34] Shtar(eos::Sim{…}, r0i::Float64, Pr0::Float64, pf::Float64)
    @ Main ~/Applications/Julia/RealisticSeismology-main/src/RealMod/src/shtarMfreeFIRK.jl:256
 [35] top-level scope
    @ REPL[82]:1
Some type information was truncated. Use `show(err)` to see complete types.

(I managed to convert those warnings into Tracebacking errors via a trick using TransformerLogger in packages LoggingExtras, Logging.)

So my question is:

Is there some (convenient) way to limit the function range of exploration by finite_difference_jacobian ?
Placing hard limits or unphysical extrapolations on my derivative function would seem to risk corrupting the jacobian.

thanks

[caMi11er]

One further observation — this variant is also accepted, but produces a different TraceBack

sol = solve(prob, RadauIIa5(jac_alg =  BVPJacobianAlgorithm(nothing)); dt = 0.5, dtmax=0.5, callback=cb, abstol=1e-4, reltol=1e-4, verbose=false) 

  [1] error(s::String)
    @ Base ./error.jl:35
  [2] #15
    @ ./REPL[67]:3 [inlined]
  [3] handle_message(::TransformerLogger{TransformerLogger{TransformerLogger{…}, var"#13#14"}, var"#15#16"}, ::LogLevel, ::Vararg{Any}; kwargs::@Kwargs{})
    @ LoggingExtras ~/.julia/packages/LoggingExtras/cFgEq/src/CompositionalLoggers/transformer.jl:22
  [4] handle_message(::TransformerLogger{TransformerLogger{…}, var"#15#16"}, ::LogLevel, ::String, ::Module, ::Symbol, ::Symbol, ::String, ::Int64)
    @ LoggingExtras ~/.julia/packages/LoggingExtras/cFgEq/src/CompositionalLoggers/transformer.jl:20
  [5] #invokelatest#2
    @ ./essentials.jl:1055 [inlined]
  [6] invokelatest
    @ ./essentials.jl:1052 [inlined]
  [7] macro expansion
    @ ./logging/logging.jl:353 [inlined]
  [8] solve(::IntervalNonlinearProblem{…}, ::ITP{…}; maxiters::Int64, abstol::Nothing, verbose::Bool, kwargs::@Kwargs{})
    @ BracketingNonlinearSolve ~/.julia/packages/BracketingNonlinearSolve/1Lfiz/src/itp.jl:86
  [9] solve
    @ ~/.julia/packages/BracketingNonlinearSolve/1Lfiz/src/itp.jl:59 [inlined]
 [10] nb
    @ ~/Applications/Julia/RealisticSeismology-main/src/sim.jl:153 [inlined]
 [11] ε(sim::Sim{…}, p::Float64)
    @ Main ~/Applications/Julia/RealisticSeismology-main/scripts/preprofile.jl:22
 [12] shtructure!(dy_dr::Base.ReinterpretArray{…}, y::Base.ReinterpretArray{…}, param::Tuple{…}, r::Float64)
    @ Main ~/Applications/Julia/RealisticSeismology-main/src/RealMod/src/shtarMfreeFIRK.jl:59
 [13] BVPFunction
    @ ~/.julia/packages/SciMLBase/LGLmH/src/scimlfunctions.jl:2524 [inlined]
 [14] Φ!(residual::Vector{…}, fᵢ_cache::PreallocationTools.DiffCache{…}, k_discrete::Vector{…}, f!::BVPFunction{…}, TU::BoundaryValueDiffEqFIRK.FIRKTableau{…}, y::Vector{…}, u::Vector{…}, p::Tuple{…}, mesh::Vector{…}, mesh_dt::Vector{…}, stage::Int64)
    @ BoundaryValueDiffEqFIRK ~/.julia/packages/BoundaryValueDiffEqFIRK/xca6M/src/collocation.jl:33
 [15] Φ!
    @ ~/.julia/packages/BoundaryValueDiffEqFIRK/xca6M/src/collocation.jl:2 [inlined]
 [16] __firk_loss_collocation!
    @ ~/.julia/packages/BoundaryValueDiffEqFIRK/xca6M/src/firk.jl:759 [inlined]
 [17] #117
    @ ~/.julia/packages/BoundaryValueDiffEqFIRK/xca6M/src/firk.jl:411 [inlined]
 [18] compute_ydual_twoarg(f!::BoundaryValueDiffEqFIRK.var"#117#123"{…}, y::Vector{…}, prep::DifferentiationInterfaceForwardDiffExt.ForwardDiffTwoArgPushforwardPrep{…}, x::Vector{…}, tx::NTuple{…}, contexts::DifferentiationInterface.Constant{…})
    @ DifferentiationInterfaceForwardDiffExt ~/.julia/packages/DifferentiationInterface/7eD1K/ext/DifferentiationInterfaceForwardDiffExt/twoarg.jl:59
 [19] pushforward!
    @ ~/.julia/packages/DifferentiationInterface/7eD1K/ext/DifferentiationInterfaceForwardDiffExt/twoarg.jl:122 [inlined]
 [20] _sparse_jacobian_aux!(f_or_f!y::Tuple{…}, jac::SparseArrays.SparseMatrixCSC{…}, prep::DifferentiationInterfaceSparseMatrixColoringsExt.SMCPushforwardSparseJacobianPrep{…}, backend::AutoSparse{…}, x::Vector{…}, contexts::DifferentiationInterface.Constant{…})
    @ DifferentiationInterfaceSparseMatrixColoringsExt ~/.julia/packages/DifferentiationInterface/7eD1K/ext/DifferentiationInterfaceSparseMatrixColoringsExt/jacobian.jl:308
 [21] jacobian!
    @ ~/.julia/packages/DifferentiationInterface/7eD1K/ext/DifferentiationInterfaceSparseMatrixColoringsExt/jacobian.jl:237 [inlined]
 [22] jacobian(f!::BoundaryValueDiffEqFIRK.var"#117#123"{…}, y::Vector{…}, prep::DifferentiationInterfaceSparseMatrixColoringsExt.SMCPushforwardSparseJacobianPrep{…}, backend::AutoSparse{…}, x::Vector{…}, contexts::DifferentiationInterface.Constant{…})
    @ DifferentiationInterfaceSparseMatrixColoringsExt ~/.julia/packages/DifferentiationInterface/7eD1K/ext/DifferentiationInterfaceSparseMatrixColoringsExt/jacobian.jl:250
 [23] __construct_nlproblem(cache::BoundaryValueDiffEqFIRK.FIRKCacheExpand{…}, y::Vector{…}, loss_bc::BoundaryValueDiffEqFIRK.var"#115#121"{…}, loss_collocation::BoundaryValueDiffEqFIRK.var"#117#123"{…}, loss::BoundaryValueDiffEqFIRK.var"#119#125"{…}, ::SciMLBase.StandardBVProblem)
    @ BoundaryValueDiffEqFIRK ~/.julia/packages/BoundaryValueDiffEqFIRK/xca6M/src/firk.jl:483
 [24] __construct_nlproblem
    @ ~/.julia/packages/BoundaryValueDiffEqFIRK/xca6M/src/firk.jl:426 [inlined]
 [25] __perform_firk_iteration(cache::BoundaryValueDiffEqFIRK.FIRKCacheExpand{…}, abstol::Float64, adaptive::Bool; nlsolve_kwargs::@NamedTuple{}, kwargs::@Kwargs{…})
    @ BoundaryValueDiffEqFIRK ~/.julia/packages/BoundaryValueDiffEqFIRK/xca6M/src/firk.jl:349
 [26] __perform_firk_iteration
    @ ~/.julia/packages/BoundaryValueDiffEqFIRK/xca6M/src/firk.jl:347 [inlined]
 [27] solve!(cache::BoundaryValueDiffEqFIRK.FIRKCacheExpand{…})
    @ BoundaryValueDiffEqFIRK ~/.julia/packages/BoundaryValueDiffEqFIRK/xca6M/src/firk.jl:302
 [28] __solve(::BVProblem{…}, ::RadauIIa5{…}; kwargs::@Kwargs{…})
    @ BoundaryValueDiffEqCore ~/.julia/packages/BoundaryValueDiffEqCore/Q7n7W/src/BoundaryValueDiffEqCore.jl:37
 [29] __solve
    @ ~/.julia/packages/BoundaryValueDiffEqCore/Q7n7W/src/BoundaryValueDiffEqCore.jl:34 [inlined]
 [30] #solve_call#35
    @ ~/.julia/packages/DiffEqBase/zYZst/src/solve.jl:635 [inlined]
 [31] solve_call
    @ ~/.julia/packages/DiffEqBase/zYZst/src/solve.jl:592 [inlined]
 [32] #solve_up#44
    @ ~/.julia/packages/DiffEqBase/zYZst/src/solve.jl:1142 [inlined]
 [33] solve_up
    @ ~/.julia/packages/DiffEqBase/zYZst/src/solve.jl:1120 [inlined]
 [34] solve(prob::BVProblem{…}, args::RadauIIa5{…}; sensealg::Nothing, u0::Nothing, p::Nothing, wrap::Val{…}, kwargs::@Kwargs{…})
    @ DiffEqBase ~/.julia/packages/DiffEqBase/zYZst/src/solve.jl:1057
 [35] Shtar(eos::Sim{…}, r0i::Float64, Pr0::Float64, pf::Float64)
    @ Main ~/Applications/Julia/RealisticSeismology-main/src/RealMod/src/shtarMfreeFIRK.jl:256
 [36] top-level scope
    @ REPL[93]:1
Some type information was truncated. Use `show(err)` to see complete types.

[caMi11er]

More information about the range of target function values explored by the solver.
I installed a push onto a stack inside my derivative function.

The range explored is at least 10,000 times larger than the values provided (and observably retrieved) in my initialguess function covering the whole span in t, including wrong-sign. So the solver is wasting most of its time strolling abroad.

I seem to be missing some fundamental constraint on the "universe of discourse".
Did I miss some fundamental advice in the documents about how the solver is supposed to learn about the relevant target-function domains?

If provided in the boundarycondition function, would the solver respect them at all stages in the process?

Should I do this by specifying such as…

residual[1] = (sol[1,:] > maxvalue 

[caMi11er]

The plot of this comedy has thickened further.

I tried a form of the boundary condition function called by RadauIIa5 that might better constrain the range of target function values explored by the solver…

    124 function bcond2!(resid, sol, par, r)     
    125     eos, rspan, Peak, pf = par  
    126     rin, rout = rspan 
    127     δ = 0.1    # Proximity test
    128     Pmax = 1e-4    # km^(-2)
    129 
    130 
    131 
    132     global typeofr
    133     typeofr = typeof(r) 
    134 
    135     resid[1] = (sol(r)[2] < 0) + (sol(r)[2] > Pmax)   # Universe of discourse! 
    136     if ( (rin - r) < δ)
    137         resid[2] = sol(r)[2]  #  p = 0 at span ends
    138         resid[3] = 10*(sol(r)[1] - sol(r+0.1)[1])    #  rmβ'(inner)=0
    139     elseif ( (rout - r) < δ)
    140         resid[4] = sol(r)[2]  #  p = 0 at span ends
    141         resid[5] = 10*(sol(r)[1] - sol(r-0.1)[1])    #  rmβ'(outer)=0
    142     end
    143 end 

but this faIled because, I now discover, the t argument (called r here) is provided as a vector.

julia> typeofr
Vector{Float64} (alias for Array{Float64, 1})

This is unlike the scalar t argument of the target-derivative function, and the initiaL guess function.
I raised such questions about this t argument some weeks ago, but I don't recall seeing an answer.
Is this because the answer is different for the various solvers? ! Did I miss this fact in the documentation?

Anyway, please clarify the implications.

Since the solver probably cannot accept a 2-dimensional matrix of residuals containing one t-span vector for each of the residuals coded above, I wondered if I am expected to create a residual vector with the same length as the provided t vector, and combine all constraints that apply at each t value into one corresponding residual vector element.

On the other hand, I should specify the length of the residual vector in the problem definition; e.g., bcresid_prototype = zeros(5).

So my only option seems to be to provide a few residual values, each element of which applies to the entire t-vector.

But this seems like a weak constraint, and certaInly not what is meant by the concept of "boundary condition".

So now I totally confused. I wonder about the concept of the authors of RadauIIa5.

Surely I can't be the first person to confront this enigma! Really?

[ChrisRackauckas]

Please share a runnable example.

[caMi11er]

Yes, a request for a running example is usually fine.

In this case, browsing the FIRK code helped me to better understand the documentation for the standard multi-point BCProblem, so I made progress to the next question.

I now understand that the solver creates a NonlinearSolve minimization (not root-finding?) problem with a loss function calculated presumably as the L2 norm of the resid vector from bc! . (If true, might be helpful to add this to the docs.)

So I provided a resid vector in bc! that combined specific end-span conditions with hard discontinuous penalties over the entire span in an effort to keep the derivative calculations within the defined calculable domain.
(I suppose such discontinuities may seriously degrade the jacobian benefits? Some such limits closely approach the domain of interest in this problem.)
Outside these limits, the derivatives are unphysical. Will I need to fake some extrapolations?

function bcond2!(resid, sol, par, r)     
    eos, rspan, Peak, pf = par  
    rin, rout = rspan 

    Pmax = 1e-4    # sol[2] span of validity (0.0, Pmax) 

    global x   # x is last element of initial guess at t = span[2]

#  Target function domain limits:
    resid =  1.0*( (sol[2] .< 0.0) .+ (sol[2] .> Pmax) )  # While declaring resid vector as Float64

# Add in further (positive definite) penalty contributions at span-ends:
    resid[1] += 10.0*abs(sol[2,begin]) + 10.0*(sol[1,begin] > 0.0)  
    resid[end] += 10.0*(abs(sol[2,end]) + abs(sol[1,end] - x[1]) ) 
end 

Is this the general idea?
I worry that the several-dozen hard domain limits may dilute the influence of the end-span boundary conditions.
Hence the 10x weights.

It still doesn't work.
Values of the target functions are explored by the solver that are now more than twice both the initial guess functions and the hard discontinuous limits defined in the resid vector.

Hmmm…

Thanks for reading this.