Bad performance of DifferentialEquations.jl

New to Julia
#1

I am trying to solve a ODE with ModelingToolkit.jl, unfortunately I have been waiting for the result for more than 11 hours.
I would to use Multithread, or some CUDA package in order to increase this awful performance, but I am pretty new to this (not only to Julia).

using ModelingToolkit, DifferentialEquations;
using Plots;
using DataFrames, CSV;

@parameters t ξ;
@variables φ(t) ψ(t);

N = 2;
n = 4;
V₀ = 0.002;
B = φ^N;
V = V₀ * φ^n;

D = Differential(t);

dB = expand_derivatives(D(B)/D(φ));
dV = expand_derivatives(D(V)/D(φ));
ddB = expand_derivatives(D(dB)/D(φ));
ddV = expand_derivatives(D(dV)/D(φ));

H1 = (3 * ξ * ψ * dB - sqrt(9 * ξ^2 * φ^2 * dB^2 + (1 - 6 * ξ * B)*(ψ^2 + 2 * V)) )/(1 - 6 * ξ * B)

dψ = -3 * H1 * ψ +((3 * ξ * dB *(1 - 3 * ξ * ddB))/(1 - 6 * ξ * B + 9 *ξ ^2 * dB^2)) * ψ^2 - (dV + 6 * ξ * (2 * V * dB - dV * B))/(1 - 6 * ξ * B + 9 * ξ^2 * dB^2)

eqs1 = [D(φ) ~ ψ,
D(ψ) ~ dψ]
sys1 = ODESystem(eqs1)

#green
u01 = [φ => 0.000,
ψ => 0.010]
p1 = [ξ => -1.0]
tspan1 = (0.0,150.0)
prob1 = ODEProblem(sys1,u01,tspan1,p1)
sol1 = solve(prob1,AutoTsit5(Rosenbrock23()))
df1 = DataFrame(sol1)
CSV.write(“greenxi-1fit0,010.csv”,df1)

plot(sol1,
vars = (t,φ),
title = “φ vs t”,
xaxis = “Time (t)”,
yaxis = “φ(t)”,
color = green,
legend = false)

#2

There is a few reasons why your code is running slow. Mainly that all your variables are in the global scope and so the compiler can not generate efficient code for it. The first thing I’d do is wrap your entire code into a function. That should help a lot. Then consider reading the Performance Tips.

That being said, I am not sure why that piece of code is taking more than 11 hours to run. What part of the code is getting hung?

Lastly, please provide julia version versioninfo() and package versions by julia> ] st

1 Like
#3

versioninfo():

Commit 6aaedec (2021-04-23 05:59 UTC)
Platform Info:
OS: Linux (x86_64-redhat-linux)
CPU: Intel(R) Core™ i7-10750H CPU @ 2.60GHz
WORD_SIZE: 64
LIBM: libopenlibm
LLVM: libLLVM-11.0.1 (ORCJIT, skylake)

] status:

 Status `~/.julia/environments/v1.6/Project.toml`

[336ed68f] CSV v0.8.4
[a93c6f00] DataFrames v1.1.1
[0c46a032] DifferentialEquations v6.16.0
[28b8d3ca] GR v0.57.4
[7073ff75] IJulia v1.23.2
[4138dd39] JLD v0.12.3
[961ee093] ModelingToolkit v5.16.0
[e7bfaba1] NumericalIntegration v0.3.3
[47be7bcc] ORCA v0.5.0
[8314cec4] PGFPlotsX v1.2.10
[f0f68f2c] PlotlyJS v0.14.1
[91a5bcdd] Plots v1.13.2
[d330b81b] PyPlot v2.9.0
[24249f21] SymPy v1.0.44

This is the part which is getting hung:

u01 = [φ => 0.000,
ψ => 0.010]
p1 = [ξ => -1.0]
tspan1 = (0.0,150.0)
prob1 = ODEProblem(sys1,u01,tspan1,p1)
sol1 = solve(prob1,AutoTsit5(Rosenbrock23()))
df1 = DataFrame(sol1)
CSV.write(“greenxi-1fit0,010.csv”,df1)

Thank you so much @affans, I will check the performance tips

#4

I do not think it is because of that. The ModelingToolkit converts system of equations into a function by itself and calls with given initial state u01 and p1. @leonandonayre you can try a small time step(if it is oscillating fast then use a fraction of the period of oscillation) for tspan to see if it really shows something or if it stuck in somewhere in compilation. If it still takes much time then interrupt and show us some stacktrace from the interruption.

1 Like
#5

This example is broken running H1 = (3ξψdB - sqrt(9ξ^2 * φ^2 * dB^2 + (1 - 6ξB)(ψ^2 + 2V)) )/(1 - 6ξB) gives
ERROR: UndefVarError: ξψdB not defined.

#6

There are some * missing when I copied the code from the notebook

Screenshot_2021-05-10 Toporensky 1404 6226 - Jupyter Notebook
Screenshot_2021-05-10 Toporensky 1404 6226 - Jupyter Notebook1148×970 111 KB

#7

Can you please update the original post to be actually runable? Otherwise it is very difficult to debug this.

2 Likes
#8

Sorry. I think now it will work

#9

The problem appears to be that the model completely blows up between x=14 and x=15. If you use Rodas5() as your solver, you will see that it finishes in about 5 seconds, but stops after around x=14 because the solution blows up to infinity.

2 Likes
#10

Just run the model with solvers from 3 languages: Julia, C++, and Fortran.

using OrdinaryDiffEq, Sundials, ODEInterfaceDiffEq
sol1 = solve(prob1,Rodas5(),abstol=1e-12,reltol=1e-12)
sol2 = solve(prob1,CVODE_BDF(),abstol=1e-12,reltol=1e-12)
sol3 = solve(prob1,radau(),abstol=1e-12,reltol=1e-12)

@show sol1.t[end], sol2.t[end], sol3.t[end]
@show sol1[end], sol2[end], sol3[end]

(sol1.t[end], sol2.t[end], sol3.t[end]) = (14.57381167838683, 14.57381163596623, 14.573811679098052)
(sol1[end], sol2[end], sol3[end]) = ([213178.70363115682, 4.8383636378845696e17], [1.868111482362633e36, 2.930079193719137e117], [349223.9235893567, 2.3765960233233695e18])

You see that every robust solver from each language all blows up to infinity at almost precisely the same point in time. That is very strong evidence that the true solution of the model you had written down goes to infinity at that time.

It seems like the automatic stiffness detection falls into a loop (which we should fix), but at least your issue is that the differential equation you wrote down is divergent and you should probably double check for sign issues.

12 Likes
#11

Should I report the issue in github?

#12

Yes