使用 fortran 包装的右侧时，Julia 的 DifferentialEquations 包失败

Question

我正在尝试使用 Julia 的 DifferentialEquations 方法求解常微分方程组。我的 ODE 的右侧被包裹 Fortran 90。这是我的 Julia 代码：

using DifferentialEquations
function rhs(dNdt,N,p,t)
    ccall((:__atmos_MOD_rhs, "./EvolveAtmFort.so"), Cvoid,(Ref{Float64}, Ref{Float64},Ref{Float64}),t,N,dNdt)
end

N0 = [0.0,298.9,0.0562,22.9,0.0166,35.96,0.0,0.0,0.0,0.0]*6.022e23
tspan = [0.0,1.0e6*365.0*24.0*60.0*60.0]
prob = ODEProblem(rhs,N0,tspan)
sol = solve(prob,Rodas5());

这会产生以下长错误，这与计算右侧的 derivative/jacobian 有关。下面，我只包含 Stacktrace 中看起来很重要的一些部分。

MethodError: no method matching Float64(::ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.TimeGradientWrapper{ODEFunction{true,typeof(rhs),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Array{Float64,1},DiffEqBase.NullParameters},Float64},Float64,1})
Closest candidates are:
  Float64(::Real, !Matched::RoundingMode) where T<:AbstractFloat at rounding.jl:200
  Float64(::T) where T<:Number at boot.jl:715
  Float64(!Matched::Int8) at float.jl:60
  ...

Stacktrace:
[1] convert(::Type{Float64},...
[2] Base.RefValue{Float64}...
[3] convert(::Type{Ref{Float64}}, ...
[4] cconvert(::Type{T} where T, 
[5] rhs(::Array{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.TimeGradientWrapper{...
[6] (::ODEFunction{true,typeof(rhs),LinearAlgebra...
[7] (::DiffEqBase.TimeGradientWrapper{ODEFunction{true,typeof(rhs),LinearAlgebra...
[8] derivative!(::Array{Float64,1},...
[9] calc_tderivative!(::OrdinaryDiffEq...
[10] calc_rosenbrock_differentiation! at...
[etc...]

当我使用无 jacobian 方法时，如 Tsit5()，集成工作正常。只有需要雅可比计算的方法才会失败。我究竟做错了什么？如何调整我的 Fortran 包装程序以便我可以使用隐式方法？谢谢！

Answer 1

这是由于某些隐式求解器中的隐式自动微分。您需要将其关闭，即 Rodas5(autodiff=false).

使用 fortran 包装的右侧时，Julia 的 DifferentialEquations 包失败

Julia's DifferentialEquations package fails when using fortran-wrapped right-hand-side

julia

jit

differential-equations