Hallo,
I am struggling to get the simple derivatives works. I looked at many examples of ForwardDiff but still failed to solve my problem.
I am solving a Jacobian of a Triagonal Matrix by using numerical derivatives instead of deriving the derivatives manually which is a pain (worked but lack of flexibility). This is to solve the 1D Hydrology model of the Richards’ equation by computing the derivatives numerically.
The example below fails for Ψ₂[iT,iZ] = ψ . My understanding is that there is a type problem. ForwardDiff uses ψ = Dual Type and Ψ₂[ = Array{Float64}
I Guess that an answer to my problem would be how to simply transform a Dual Type to a Float?
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# FUNCTION : ∂R∂Ψ_NUMERICAL
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
function ∂RESIDUAL∂Ψ_NUMERICAL(discret, hydro, iT, iZ, K_Aver, N_iZ, Q, Sorptivity, ΔHpond, ΔPr, ΔSink, ΔT, θ, Ψ₂)
function FUNCTION(ψ, discret, hydro, iT, iZ, K_Aver, N_iZ, Q, Sorptivity, ΔHpond, ΔPr, ΔSink, ΔT, θ, Ψ₂)
Ψ₂[iT,iZ] = ψ
∂R∂Ψ = RESIDUAL_DERIVATIVE(discret, hydro, iT, iZ, K_Aver, N_iZ, Q, Sorptivity, ΔHpond, ΔPr, ΔSink, ΔT, θ, Ψ₂)
return ∂R∂Ψ
end # function: Function
ψ =Vector{Float64}(undef, 1)
∂R∂Ψ_Numerical(ψ::Vector) = FUNCTION(abs(ψ[1]), discret, hydro, iT, iZ, K_Aver, N_iZ, Q, Sorptivity, ΔHpond, ΔPr, ΔSink, ΔT, θ, Ψ₂)
ψ[1] = Ψ₂[iT,iZ]
Func_∂R∂Ψ_Numerical = ψ -> ForwardDiff.gradient(∂R∂Ψ_Numerical, ψ)
∂R∂Ψ2 = Func_∂R∂Ψ_Numerical(ψ)[1]
return ∂R∂Ψ2
end # function: ∂RESIDUAL∂Ψ_NUMERICAL
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# FUNCTION : RESIDUAL
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
function RESIDUAL(discret, hydro, iT::Int64, iZ::Int64, Q, ΔSink, ΔT, θ, Ψ)::Float64
return Residual = discret.ΔZ[iZ] * ((θ[iT,iZ] - θ[iT-1,iZ]) - hydro.So[iZ] * (Ψ[iT,iZ] - Ψ[iT-1,iZ]) * (θ[iT,iZ] / hydro.θs[iZ])) - ΔT[iT] * (Q[iT,iZ] - Q[iT,iZ+1]) + ΔSink[iT,iZ]
end # function: RESIDUAL
The error message:
ERROR: LoadError: MethodError: no method matching Float64(::ForwardDiff.Dual{ForwardDiff.Tag{Main.residual.var"#∂R∂Ψ_Numerical#3"{Main.discretization.DISCRET,Main.hydroStruct.KOSUGI,Int64,Int64,Array{Float64,1},Int64,Array{Float64,2},Float64,Array{Float64,1},Array{Float64,1},Array{Float64,2},Array{Float64,1},Array{Float64,2},Array{Float64,2},Main.residual.var"function#2"},Float64},Float64,1})
Closest candidates are:
Float64(::Real, ::RoundingMode) where T<:AbstractFloat at rounding.jl:200
Float64(::T) where T<:Number at boot.jl:716
Float64(::Irrational{:sqrt2}) at irrationals.jl:189
Many thanks for any help you may provide,
Joseph