SciML / OrdinaryDiffEq.jl

High performance ordinary differential equation (ODE) and differential-algebraic equation (DAE) solvers, including neural ordinary differential equations (neural ODEs) and scientific machine learning (SciML)
https://diffeq.sciml.ai/latest/
Other
563 stars 211 forks source link

Long Compile Times for StaticArrays Arguments to `solve` #2535

Open cadojo opened 1 week ago

cadojo commented 1 week ago

Describe the bug 🐞

I think there may be a regression in compilation times solve(::ODEProblem) with StaticArrays types. I have a minimum reproducible example below. First a Vector type is used, and that compiles quickly. Next an MVector type is used, and the solve call takes approximately one hour to complete; I haven't confirmed this locally (the call hangs, and I wait a few minutes before killing the process), but in CI the call does finish, but it takes around one hour.

I don't have any specific package versions for when this worked, but I was able to run solve with MVector types in the last year. Here is one example from 5 months ago; the full job completed in 13 minutes.

Expected behavior

The solve call with MVector should compile and evaluate within an order of magnitude of Vector.

Minimal Reproducible Example 👇

import Pkg
Pkg.activate(tempdir())
Pkg.add("StaticArrays")
Pkg.add("OrdinaryDiffEq")

using StaticArrays
using OrdinaryDiffEq

μ = 0.012150584395829193
u = [
    0.8222791798525636
    0.0
    0.0
    0.0
    0.13799313228400034
    0.0
]

function r2bp(du, u, p, t)
    x, y, z, ẋ, ẏ, ż = u
    μ, = p

    r = sqrt(x^2 + y^2 + z^2)

    du[1] = ẋ
    du[2] = ẏ
    du[3] = ż
    du[4] = μ * x / r^3
    du[5] = μ * y / r^3
    du[6] = μ * z / r^3

    return nothing
end

@info "Solving Problem with Vector Input"
@warn "This should compile in a couple of seconds."
problem = ODEProblem(r2bp, u, (0.0, 3.0), (μ,))
@time solution = solve(problem, Vern9(), reltol = 1e-12, abstol = 1e-12)

@info "Solving Problem with MVector State"
@error "This will take very long to compile!"
problem = ODEProblem(r2bp, MVector{6}(u), (0.0, 3.0), (μ,))
@time solution = solve(problem, Vern9(), reltol = 1e-12, abstol = 1e-12)

Environment (please complete the following information):

Status `~/.julia/environments/v1.11/Project.toml`
  [1dea7af3] OrdinaryDiffEq v6.90.1
  [90137ffa] StaticArrays v1.9.8
Status `~/.julia/environments/v1.11/Manifest.toml`
  [47edcb42] ADTypes v1.10.0
  [7d9f7c33] Accessors v0.1.38
  [79e6a3ab] Adapt v4.1.1
  [ec485272] ArnoldiMethod v0.4.0
  [4fba245c] ArrayInterface v7.17.0
  [4c555306] ArrayLayouts v1.10.4
  [62783981] BitTwiddlingConvenienceFunctions v0.1.6
  [70df07ce] BracketingNonlinearSolve v1.1.0
  [2a0fbf3d] CPUSummary v0.2.6
  [7057c7e9] Cassette v0.3.14
  [d360d2e6] ChainRulesCore v1.25.0
  [fb6a15b2] CloseOpenIntervals v0.1.13
  [38540f10] CommonSolve v0.2.4
  [bbf7d656] CommonSubexpressions v0.3.1
  [f70d9fcc] CommonWorldInvalidations v1.0.0
  [34da2185] Compat v4.16.0
  [a33af91c] CompositionsBase v0.1.2
  [2569d6c7] ConcreteStructs v0.2.3
  [187b0558] ConstructionBase v1.5.8
  [adafc99b] CpuId v0.3.1
  [9a962f9c] DataAPI v1.16.0
  [864edb3b] DataStructures v0.18.20
  [e2d170a0] DataValueInterfaces v1.0.0
  [2b5f629d] DiffEqBase v6.160.0
  [163ba53b] DiffResults v1.1.0
  [b552c78f] DiffRules v1.15.1
  [a0c0ee7d] DifferentiationInterface v0.6.22
  [ffbed154] DocStringExtensions v0.9.3
  [4e289a0a] EnumX v1.0.4
  [f151be2c] EnzymeCore v0.8.6
  [d4d017d3] ExponentialUtilities v1.26.1
  [e2ba6199] ExprTools v0.1.10
⌅ [6b7a57c9] Expronicon v0.8.5
  [7034ab61] FastBroadcast v0.3.5
  [9aa1b823] FastClosures v0.3.2
  [29a986be] FastLapackInterface v2.0.4
  [a4df4552] FastPower v1.1.1
  [1a297f60] FillArrays v1.13.0
  [6a86dc24] FiniteDiff v2.26.0
  [f6369f11] ForwardDiff v0.10.38
  [f62d2435] FunctionProperties v0.1.2
  [069b7b12] FunctionWrappers v1.1.3
  [77dc65aa] FunctionWrappersWrappers v0.1.3
⌅ [46192b85] GPUArraysCore v0.1.6
  [c145ed77] GenericSchur v0.5.4
  [86223c79] Graphs v1.12.0
  [3e5b6fbb] HostCPUFeatures v0.1.17
  [615f187c] IfElse v0.1.1
  [d25df0c9] Inflate v0.1.5
  [3587e190] InverseFunctions v0.1.17
  [92d709cd] IrrationalConstants v0.2.2
  [82899510] IteratorInterfaceExtensions v1.0.0
  [692b3bcd] JLLWrappers v1.6.1
  [ef3ab10e] KLU v0.6.0
  [ba0b0d4f] Krylov v0.9.8
  [10f19ff3] LayoutPointers v0.1.17
  [5078a376] LazyArrays v2.2.2
  [87fe0de2] LineSearch v0.1.4
  [d3d80556] LineSearches v7.3.0
  [7ed4a6bd] LinearSolve v2.37.0
  [2ab3a3ac] LogExpFunctions v0.3.28
  [bdcacae8] LoopVectorization v0.12.171
  [d8e11817] MLStyle v0.4.17
  [1914dd2f] MacroTools v0.5.13
  [d125e4d3] ManualMemory v0.1.8
  [bb5d69b7] MaybeInplace v0.1.4
  [46d2c3a1] MuladdMacro v0.2.4
  [d41bc354] NLSolversBase v7.8.3
  [77ba4419] NaNMath v1.0.2
  [8913a72c] NonlinearSolve v4.1.0
  [be0214bd] NonlinearSolveBase v1.3.3
  [5959db7a] NonlinearSolveFirstOrder v1.0.0
  [9a2c21bd] NonlinearSolveQuasiNewton v1.0.0
  [26075421] NonlinearSolveSpectralMethods v1.0.0
  [6fe1bfb0] OffsetArrays v1.14.1
  [bac558e1] OrderedCollections v1.6.3
  [1dea7af3] OrdinaryDiffEq v6.90.1
  [89bda076] OrdinaryDiffEqAdamsBashforthMoulton v1.1.0
  [6ad6398a] OrdinaryDiffEqBDF v1.1.2
  [bbf590c4] OrdinaryDiffEqCore v1.11.0
  [50262376] OrdinaryDiffEqDefault v1.1.0
  [4302a76b] OrdinaryDiffEqDifferentiation v1.2.0
  [9286f039] OrdinaryDiffEqExplicitRK v1.1.0
  [e0540318] OrdinaryDiffEqExponentialRK v1.1.0
  [becaefa8] OrdinaryDiffEqExtrapolation v1.2.1
  [5960d6e9] OrdinaryDiffEqFIRK v1.4.0
  [101fe9f7] OrdinaryDiffEqFeagin v1.1.0
  [d3585ca7] OrdinaryDiffEqFunctionMap v1.1.1
  [d28bc4f8] OrdinaryDiffEqHighOrderRK v1.1.0
  [9f002381] OrdinaryDiffEqIMEXMultistep v1.1.0
  [521117fe] OrdinaryDiffEqLinear v1.1.0
  [1344f307] OrdinaryDiffEqLowOrderRK v1.2.0
  [b0944070] OrdinaryDiffEqLowStorageRK v1.2.1
  [127b3ac7] OrdinaryDiffEqNonlinearSolve v1.2.4
  [c9986a66] OrdinaryDiffEqNordsieck v1.1.0
  [5dd0a6cf] OrdinaryDiffEqPDIRK v1.1.0
  [5b33eab2] OrdinaryDiffEqPRK v1.1.0
  [04162be5] OrdinaryDiffEqQPRK v1.1.0
  [af6ede74] OrdinaryDiffEqRKN v1.1.0
  [43230ef6] OrdinaryDiffEqRosenbrock v1.3.1
  [2d112036] OrdinaryDiffEqSDIRK v1.1.0
  [669c94d9] OrdinaryDiffEqSSPRK v1.2.0
  [e3e12d00] OrdinaryDiffEqStabilizedIRK v1.1.0
  [358294b1] OrdinaryDiffEqStabilizedRK v1.1.0
  [fa646aed] OrdinaryDiffEqSymplecticRK v1.1.0
  [b1df2697] OrdinaryDiffEqTsit5 v1.1.0
  [79d7bb75] OrdinaryDiffEqVerner v1.1.1
  [65ce6f38] PackageExtensionCompat v1.0.2
  [d96e819e] Parameters v0.12.3
  [f517fe37] Polyester v0.7.16
  [1d0040c9] PolyesterWeave v0.2.2
  [d236fae5] PreallocationTools v0.4.24
  [aea7be01] PrecompileTools v1.2.1
  [21216c6a] Preferences v1.4.3
  [3cdcf5f2] RecipesBase v1.3.4
  [731186ca] RecursiveArrayTools v3.27.3
  [f2c3362d] RecursiveFactorization v0.2.23
  [189a3867] Reexport v1.2.2
  [ae029012] Requires v1.3.0
  [7e49a35a] RuntimeGeneratedFunctions v0.5.13
  [94e857df] SIMDTypes v0.1.0
  [476501e8] SLEEFPirates v0.6.43
  [0bca4576] SciMLBase v2.61.0
  [19f34311] SciMLJacobianOperators v0.1.1
  [c0aeaf25] SciMLOperators v0.3.12
  [53ae85a6] SciMLStructures v1.5.0
  [efcf1570] Setfield v1.1.1
  [727e6d20] SimpleNonlinearSolve v2.0.0
  [699a6c99] SimpleTraits v0.9.4
  [ce78b400] SimpleUnPack v1.1.0
  [47a9eef4] SparseDiffTools v2.23.0
  [0a514795] SparseMatrixColorings v0.4.10
  [e56a9233] Sparspak v0.3.9
  [276daf66] SpecialFunctions v2.4.0
  [aedffcd0] Static v1.1.1
  [0d7ed370] StaticArrayInterface v1.8.0
  [90137ffa] StaticArrays v1.9.8
  [1e83bf80] StaticArraysCore v1.4.3
  [10745b16] Statistics v1.11.1
  [7792a7ef] StrideArraysCore v0.5.7
  [2efcf032] SymbolicIndexingInterface v0.3.35
  [3783bdb8] TableTraits v1.0.1
  [bd369af6] Tables v1.12.0
  [8290d209] ThreadingUtilities v0.5.2
  [a759f4b9] TimerOutputs v0.5.25
  [d5829a12] TriangularSolve v0.2.1
  [781d530d] TruncatedStacktraces v1.4.0
  [3a884ed6] UnPack v1.0.2
  [3d5dd08c] VectorizationBase v0.21.71
  [19fa3120] VertexSafeGraphs v0.2.0
  [1d5cc7b8] IntelOpenMP_jll v2024.2.1+0
  [856f044c] MKL_jll v2024.2.0+0
  [efe28fd5] OpenSpecFun_jll v0.5.5+0
  [1317d2d5] oneTBB_jll v2021.12.0+0
  [0dad84c5] ArgTools v1.1.2
  [56f22d72] Artifacts v1.11.0
  [2a0f44e3] Base64 v1.11.0
  [ade2ca70] Dates v1.11.0
  [8ba89e20] Distributed v1.11.0
  [f43a241f] Downloads v1.6.0
  [7b1f6079] FileWatching v1.11.0
  [9fa8497b] Future v1.11.0
  [b77e0a4c] InteractiveUtils v1.11.0
  [4af54fe1] LazyArtifacts v1.11.0
  [b27032c2] LibCURL v0.6.4
  [76f85450] LibGit2 v1.11.0
  [8f399da3] Libdl v1.11.0
  [37e2e46d] LinearAlgebra v1.11.0
  [56ddb016] Logging v1.11.0
  [d6f4376e] Markdown v1.11.0
  [a63ad114] Mmap v1.11.0
  [ca575930] NetworkOptions v1.2.0
  [44cfe95a] Pkg v1.11.0
  [de0858da] Printf v1.11.0
  [9a3f8284] Random v1.11.0
  [ea8e919c] SHA v0.7.0
  [9e88b42a] Serialization v1.11.0
  [1a1011a3] SharedArrays v1.11.0
  [6462fe0b] Sockets v1.11.0
  [2f01184e] SparseArrays v1.11.0
  [fa267f1f] TOML v1.0.3
  [a4e569a6] Tar v1.10.0
  [8dfed614] Test v1.11.0
  [cf7118a7] UUIDs v1.11.0
  [4ec0a83e] Unicode v1.11.0
  [e66e0078] CompilerSupportLibraries_jll v1.1.1+0
  [deac9b47] LibCURL_jll v8.6.0+0
  [e37daf67] LibGit2_jll v1.7.2+0
  [29816b5a] LibSSH2_jll v1.11.0+1
  [c8ffd9c3] MbedTLS_jll v2.28.6+0
  [14a3606d] MozillaCACerts_jll v2023.12.12
  [4536629a] OpenBLAS_jll v0.3.27+1
  [05823500] OpenLibm_jll v0.8.1+2
  [bea87d4a] SuiteSparse_jll v7.7.0+0
  [83775a58] Zlib_jll v1.2.13+1
  [8e850b90] libblastrampoline_jll v5.11.0+0
  [8e850ede] nghttp2_jll v1.59.0+0
  [3f19e933] p7zip_jll v17.4.0+2
Info Packages marked with ⌅ have new versions available but compatibility constraints restrict them from upgrading. To see why use `status --outdated -m`
Julia Version 1.11.1
Commit 8f5b7ca12ad (2024-10-16 10:53 UTC)
Build Info:
  Official https://julialang.org/ release
Platform Info:
  OS: macOS (arm64-apple-darwin22.4.0)
  CPU: 8 × Apple M1
  WORD_SIZE: 64
  LLVM: libLLVM-16.0.6 (ORCJIT, apple-m1)
Threads: 4 default, 0 interactive, 2 GC (on 4 virtual cores)
Environment:
  JULIA_NUM_THREADS = 4
  JULIA_SHELL = /bin/bash
ChrisRackauckas commented 1 week ago

I think there may be a regression in compilation times for very stiff ODE functions with StaticArrays types

But the MWE is non-stiff?

ChrisRackauckas commented 1 week ago

And it took less than 3 seconds?

2.841622 seconds (15.83 M allocations: 969.635 MiB, 5.62% gc time, 99.99% compilation time)

cadojo commented 1 week ago

And it took less than 3 seconds?

When the ODEProblem has a state vector type of Vector (the default in the example), it compiles quickly (a couple of seconds). When you uncomment the MVector line, though, it takes tens of minutes to compile (around one hour in CI). I've replaced the MRE with a modified version that shows both input types. Could you try running that again?

But the MWE is non-stiff?

Ah, sorry about that. I misunderstood the definition of stiff. The reason I tried to characterize the dynamics at all (e.g. stiff) is I think I have seen some equations which compile quickly with MVector state vector types, and so I thought there was an "ODE sensitivity" check somewhere in the stack that triggered extra compilation for some problems and not others. I thought that would be information that might help flag the issue, but I think I'm mistaken.

ChrisRackauckas commented 1 week ago

When the ODEProblem has a state vector type of Vector (the default in the example), it compiles quickly (a couple of seconds). When you uncomment the MVector line, though, it takes tens of minutes to compile (around one hour in CI). I've replaced the MRE with a modified version that shows both input types. Could you try running that again?

I did that. It compiled just fine with MVector.

cadojo commented 1 week ago

I did that. It compiled just fine with MVector.

Hm, not for me. I added some Pkg.activate lines at the top just to make sure. It never finishes running for me, and I am assuming it's stuck on compilation because I've solved these problems with the same initial conditions before.

The linked issue was reported by another user too. I solved the issue by changing all code examples from MVector to Vector.

Sorry this is an odd case. Maybe it's an architecture issue? I'm on an M1 MacBook Air. The output from me running the MRE is below; it just hangs on the solve call with the MVector input.

  Activating project at `/var/folders/vc/rjqggbq92mj_strb625fy2g40000gn/T`
    Updating registry at `~/.julia/registries/General.toml`
   Resolving package versions...
  No Changes to `/private/var/folders/vc/rjqggbq92mj_strb625fy2g40000gn/T/Project.toml`
  No Changes to `/private/var/folders/vc/rjqggbq92mj_strb625fy2g40000gn/T/Manifest.toml`
   Resolving package versions...
  No Changes to `/private/var/folders/vc/rjqggbq92mj_strb625fy2g40000gn/T/Project.toml`
  No Changes to `/private/var/folders/vc/rjqggbq92mj_strb625fy2g40000gn/T/Manifest.toml`
[ Info: Solving Problem with Vector Input
┌ Warning: This should compile in a couple of seconds.
└ @ Main ~/Projects/Astrodynamics/GeneralAstrodynamics.jl/check.jl:37
  2.234717 seconds (21.28 M allocations: 1.192 GiB, 9.82% gc time, 99.99% compilation time)
[ Info: Solving Problem with MVector State
┌ Error: This will take very long to compile!
└ @ Main ~/Projects/Astrodynamics/GeneralAstrodynamics.jl/check.jl:42
cadojo commented 1 week ago

I think I found the issue for me. Changing the algorithm from Vern9 to Tsit5 solves the issue for me. I have a new MRE below which tries out the two algorithms. For me, the final solve call hangs.

Modified MRE ```julia import Pkg Pkg.activate(tempdir()) Pkg.add("StaticArrays") Pkg.add("OrdinaryDiffEq") using StaticArrays using OrdinaryDiffEq μ = 0.012150584395829193 u = [ 0.8222791798525636 0.0 0.0 0.0 0.13799313228400034 0.0 ] function r2bp(du, u, p, t) x, y, z, ẋ, ẏ, ż = u μ, = p r = sqrt(x^2 + y^2 + z^2) du[1] = ẋ du[2] = ẏ du[3] = ż du[4] = μ * x / r^3 du[5] = μ * y / r^3 du[6] = μ * z / r^3 return nothing end @info "Solving Problem with Vector Input using Tsit5" @warn "This should compile in a couple of seconds." problem = ODEProblem(r2bp, u, (0.0, 3.0), (μ,)) @time solution = solve(problem, Tsit5(), reltol = 1e-12, abstol = 1e-12) @info "Solving Problem with MVector Input using Tsit5" @warn "This should also compile in a couple of seconds." problem = ODEProblem(r2bp, MVector{6}(u), (0.0, 3.0), (μ,)) @time solution = solve(problem, Tsit5(), reltol = 1e-12, abstol = 1e-12) @info "Solving Problem with Vector Input using Vern9" @warn "This should compile in a couple of seconds." problem = ODEProblem(r2bp, u, (0.0, 3.0), (μ,)) @time solution = solve(problem, Vern9(), reltol = 1e-12, abstol = 1e-12) @info "Solving Problem with MVector Input using Vern9" @error "This takes a long time to compile for me; upwards of one hour." problem = ODEProblem(r2bp, MVector{6}(u), (0.0, 3.0), (μ,)) @time solution = solve(problem, Vern9(), reltol = 1e-12, abstol = 1e-12) ```
ChrisRackauckas commented 1 week ago

Still only took:

2.199512 seconds (12.69 M allocations: 758.995 MiB, 8.69% gc time, 99.99% compilation time)

How did you install Julia, using juliaup? What version?

cadojo commented 1 week ago

Strange. Yeah, I installed it through juliaup. Version 1.11.1.

topolarity commented 1 week ago

I can reproduce this on my x86_64-linux machine.

It looks like an issue with a subtyping query in the broadcast exploding: https://github.com/JuliaLang/julia/issues/56606

topolarity commented 5 days ago

Good news! Testing with https://github.com/JuliaLang/julia/pull/56640 seems to resolve the problem 👍 (huge thanks to @N5N3):

[ Info: Solving Problem with Vector Input
┌ Warning: This should compile in a couple of seconds.
└ @ Main ~/repos/julia/mwe.jl:36
  4.455485 seconds (20.61 M allocations: 1.135 GiB, 16.08% gc time, 100.00% compilation time)
[ Info: Solving Problem with MVector State
┌ Error: This will take very long to compile!
└ @ Main ~/repos/julia/mwe.jl:41
  4.371111 seconds (18.77 M allocations: 902.214 MiB, 15.37% gc time, 100.00% compilation time)