NonlinearSolve uses AD Jacobian even when I specify it analytically

Hm, I recently tried specification of the Jacobian in this way - however, NonlinearSolve still tries to attempt to construct the Jacobian using AD. Has anything changed here?

That’s a different topic, so please make a separate thread with an MWE.

So here is an MWE that demonstrates unexpected behaviour (for me) at multple places:

using NonlinearSolve

function r(gamma, params_nonlinear)
    println("r") # Does show up
    return gamma^params_nonlinear.exponent
end

function dr(gamma, params_nonlinear)
    println("dr") # Does not show up
    return params_nonlinear.exponent * gamma^(params_nonlinear.exponent - 1.0)
end

h = NonlinearFunction(r; jac = dr)
gamma_0 = 1.0
params_nonlinear = (exponent = 2.0, )
sol = NonlinearSolve.solve(NonlinearProblem(h, gamma_0, params_nonlinear), 
                           NewtonRaphson(;concrete_jac = true, autodiff = false), 
                           abstol = 2 * eps(Float64))
println(sol.stats) # SciMLBase.NLStats(27, 26, 0, 26, 26) => 26 Jacobians created

You seem to be using an old version of NonlinearSolve. Passing false to autodiff will fail in the new releases. In v4 I am getting:

julia> sol = NonlinearSolve.solve(NonlinearProblem(h, gamma_0, params_nonlinear), 
                                  NewtonRaphson(;concrete_jac = true), 
                                  abstol = 2 * eps(Float64))
r
dr
r
dr
r
dr
r
dr
r
dr
r
dr
r
dr
r
dr
r
dr
r
dr
r
dr
r
dr
r
dr
r
dr
r
dr
r
dr
r
dr
r
dr
r
dr
r
dr
r
dr
r
dr
r
dr
r
dr
r
dr
r
dr
r
r
r
retcode: Success
u: 1.4901161193847656e-8
1 Like