Use multithreaded Jacobian function generated with ModelingToolKit.jl for application in DifferentialEquations.jl

General Usage
1

Hi,

I would like to complete the example PDE in the SciML News (also in the ModelingToolKit.jl documentation), section “ModelingToolkit Automatic Parallelism”, on automatically computed multithreaded Jacobian functions via ModelingToolkit.build_function with the objective to speed-up ODE simulations in DifferentialEquations.jl. I follow the steps and obtain the function multithreadedjac. I then go on with

#=
previous code, see provided links
...
multithreadedjac = eval(ModelingToolkit.build_function(vec(jac),u,multithread=true)[2])
=#
MyA = zeros(N,N)
AMx = zeros(N,N)
DA  = zeros(N,N)
# Define the initial condition as normal arrays
u0 = rand(N,N,3)
#solve 
using DifferentialEquations
prob_jac = ODEProblem(ODEFunction(f!,jac = multithreadedjac),u0,(0.0,10.0))
solve(prob_jac)

leading to a
ERROR: MethodError: no method matching (::var"#31#32")(::Array{Float64,2}, ::Array{Float64,3}, ::DiffEqBase.NullParameters, ::Float64)...

Any assistance on how to successfully translate back the multithreaded Jacobian for applications as DiffEq ODE is much appreciated.

I am using Julia v1.4, DifferentialEquations.jl v6.14.0 and ModelingToolKit.jl v3.10.2

2

Hey,
Sorry it took a bit to respond. I was redoing the documentation, and in the new version the way to do this is much more explicit:

https://mtk.sciml.ai/dev/tutorials/auto_parallel/

Cheers!

2 Likes
3

P.S. I found a sync bug when writing that so https://github.com/SciML/ModelingToolkit.jl/pull/488 is needed for the tutorial to work if you solve multiple times in rapid succession… I’ll get this patched in a few hours.

4

Thanks @ChrisRackauckas for the quick update of the documentation. It makes the application of efficient multithreading (or parallel computing) straightforward with DifferentialEquations.jl!
Regarding the actual computation of the jacobian using the function ModelingToolkit.jacobian(), is there any chance to release a version which directly considers ODE sparsity rather than applying sparse() to the dense jacobian generated (as it has been done with other ModelingToolkit.jl functions using “sparse” as a keyword argument)? There should be benefits for larger systems where jacobian calculations can take a while (my case).

5

@luis_munoz yeah… @shashi was working on that today :slight_smile:

https://github.com/JuliaSymbolics/SymbolicUtils.jl/issues/61#issuecomment-654826458

2 Likes
6

Great! I’ll keep an eye on the progress then :+1: