I am using the IDA solver from Sundials.jl to solve a system of DAEs. It works fine, but each time step needs a lot of allocations, even though my callback function does not allocate at all.
MWE:
using KiteModels
kps4::KPS4 = KPS4(KCU(se()))
integrator = KiteModels.init_sim!(kps4; stiffness_factor=0.035, prn=false)
@timev next_step!(kps4, integrator, v_ro = 0, dt=0.05)
nothing
Output:
julia> include("mwes/mwe_02.jl")
0.003074 seconds (16.70 k allocations: 1.230 MiB)
elapsed time (ns): 3073747
gc time (ns): 0
bytes allocated: 1290000
pool allocs: 16703
non-pool GC allocs: 0
minor collections: 0
full collections: 0
This is for one high-level time-step of 50ms, which needs 50 to 1500 callbacks to find a solution.
Any idea how to reduce the allocations?
I know, I should use ModelingToolkit, but re-writing my model for ModelingToolkit would be a huge task. Is there anything else I can do to improve the speed or reduce the allocations?
You can find the function init_sim!() here: KiteModels.jl/src/KiteModels.jl at 0abe18f6de4b6ba62df145c13423b05584c25817 · ufechner7/KiteModels.jl · GitHub
and the function next_step!() here: KiteModels.jl/src/KiteModels.jl at 0abe18f6de4b6ba62df145c13423b05584c25817 · ufechner7/KiteModels.jl · GitHub
The state vector is has 61 elements of type Float64.
The residual function looks like this:
"""
residual!(res, yd, y::MVector{S, SimFloat}, s::KPS4, time) where S
N-point tether model, four points for the kite on top:
Inputs:
State vector y = pos1, pos2, ... , posn, vel1, vel2, . .., veln, length, v_reel_out
Derivative yd = posd1, posd2, ..., posdn, veld1, veld2, ..., veldn, lengthd, v_reel_outd
Output:
Residual res = res1, res2 = vel1-posd1, ..., veld1-acc1, ...,
Additional parameters:
s: Struct with work variables, type KPS4
S: The dimension of the state vector
The number of the point masses of the model N = S/6, the state of each point
is represented by two 3 element vectors.
"""
function residual!(res, yd, y::Vector{SimFloat}, s::KPS4, time)
S = length(y)
y_ = MVector{S, SimFloat}(y)
yd_ = MVector{S, SimFloat}(yd)
residual!(res, yd_, y_, s, time)
end
function residual!(res, yd, y::MVector{S, SimFloat}, s::KPS4, time) where S
...
end
. But yes, in a serious project you’d probably just fix that if you’re worried about robustness and performance.