Preallocate u and du in DifferentialEquations.jl

General Usage
1

Hi all,

we are currently writing a simulation engine which is a massive ODE with >>100k elements and therefore need to be very memory conscious.
We need to pre-allocate the u and du vectors and have DifferentialEquations.jl use them in-place. Currently there seems to be no way to supply u and du, only u0, which is not used internally.

This is how we would like to use it:

using OrdinaryDiffEq
function lorenz!(du,u,p,t)
 du[1] = 10.0(u[2]-u[1])
 du[2] = u[1]*(28.0-u[3]) - u[2]
 du[3] = u[1]*u[2] - (8/3)*u[3]
end
u = [1.0;0.0;0.0]
du = [0.0; 0.0; 0.0]
tspan = (0.0,100.0)

# ODEProblem doesnt accept u and du, only u0
prob = ODEProblem(lorenz!,u, du, tspan)

integrator = init(prob,Tsit5())
step!(integrator,1)

Is there a way to make this work?

2

That allocates all of the necessary vectors. Then step!(integrator) is non-allocating.

How else would you know what vectors to allocate? Different methods need different vectors of different things, so that’s the general form to figure out what needs to be cached for a method.

3

I am not sure whether that’s the problem but Tsit5 will probably need to allocate a lot more memory than just u and du since it is a higher Order method right?
So it might be that you have to resort to lower order methods like BS3 or perhaps even Euler?
There are also multistep methods such as VCABM although I don’t know how much they need to allocate maybe someone more knowledgeable will be able to estimate this

4

Thanks for the quick answer!

I know internally that the lorenz! function needs du and u to be a vector with three elements - those are the elements I use, regardless of the method use.

There should be a way for me to allocate those vectors?

5

But how many vectors? For Tsit5 you need 7, for Ven7 you need over 20. init is how you can allocate the cache you find in integrator.cache.

1 Like
6

yes, I see that. I could still save some allocations when pre-allocating du and u.

Maybe I just don’t fully understand what you are doing there.

7

Which allocations?

8

we could avoid the allocation of u0.

9
  • alias_u0: allows the solver to alias the initial condition array that is contained in the problem struct. Defaults to false.
1 Like