Very high memory use with DelayDifferentialEquations

Specific Domains Modelling & Simulations
1

I have a large PDE system that I originally was solving with an explicit ODE solver via ODEProblem with DiffEq.jl. I have added a component to that system that requires the time derivative info so I changed it to a DDE and swapped over to the DDEProblem solver. The issue I am having now is that when using the DDE solver I am running out of memory, whereas my original ODE problem stayed at a fixed memory usage.

I have the following MWE below. When running on my machine the ODE solve takes about 1 GB of RAM whereas the DDE solve quickly jumps to 10GB+ of RAM. Is there some kind of issue with the DDE suite or is there a simple optimization that I need to be aware of?

*Note that my actual code follows all of the performance tips and guidelines as far as not allocating variables and performing calculations in-place as much as possible.

using DifferentialEquations, LinearAlgebra, SparseArrays
using Logging: global_logger
using TerminalLoggers: TerminalLogger
global_logger(TerminalLogger())

const N = 32
const xyd_brusselator = range(0, stop=1, length=N)
brusselator_f(x, y, t) = (((x - 0.3)^2 + (y - 0.6)^2) <= 0.1^2) * (t >= 1.1) * 5.0
limit(a, N) = a == N + 1 ? 1 : a == 0 ? N : a
function brusselator_2d_loop(du, u, p, t)
    A, B, alpha, dx = p
    alpha = alpha / dx^2
    @inbounds for I in CartesianIndices((N, N))
        i, j = Tuple(I)
        x, y = xyd_brusselator[I[1]], xyd_brusselator[I[2]]
        ip1, im1, jp1, jm1 = limit(i + 1, N), limit(i - 1, N), limit(j + 1, N),
        limit(j - 1, N)
        du[i, j, 1] = alpha * (u[im1, j, 1] + u[ip1, j, 1] + u[i, jp1, 1] + u[i, jm1, 1] -
                               4u[i, j, 1]) +
                      B + u[i, j, 1]^2 * u[i, j, 2] - (A + 1) * u[i, j, 1] +
                      brusselator_f(x, y, t)
        du[i, j, 2] = alpha * (u[im1, j, 2] + u[ip1, j, 2] + u[i, jp1, 2] + u[i, jm1, 2] -
                               4u[i, j, 2]) +
                      A * u[i, j, 1] - u[i, j, 1]^2 * u[i, j, 2]
    end
end
p = (3.4, 1.0, 10.0, step(xyd_brusselator))

function init_brusselator_2d(xyd)
    N = length(xyd)
    u = zeros(N, N, 2)
    for I in CartesianIndices((N, N))
        x = xyd[I[1]]
        y = xyd[I[2]]
        u[I, 1] = 22 * (y * (1 - y))^(3 / 2)
        u[I, 2] = 27 * (x * (1 - x))^(3 / 2)
    end
    u
end
u0 = init_brusselator_2d(xyd_brusselator)
prob_ode_brusselator_2d = ODEProblem(brusselator_2d_loop, u0, (0.0, 20.0), p)

@time solve(prob_ode_brusselator_2d, SSPRK43(), saveat=0.5, progress=true, progress_steps=1);

function brusselator_2d_loop_h(du, u, h, p, t)
    A, B, alpha, dx = p
    alpha = alpha / dx^2
    @inbounds for I in CartesianIndices((N, N))
        i, j = Tuple(I)
        x, y = xyd_brusselator[I[1]], xyd_brusselator[I[2]]
        ip1, im1, jp1, jm1 = limit(i + 1, N), limit(i - 1, N), limit(j + 1, N),
        limit(j - 1, N)
        du[i, j, 1] = alpha * (u[im1, j, 1] + u[ip1, j, 1] + u[i, jp1, 1] + u[i, jm1, 1] -
                               4u[i, j, 1]) +
                      B + u[i, j, 1]^2 * u[i, j, 2] - (A + 1) * u[i, j, 1] +
                      brusselator_f(x, y, t)
        du[i, j, 2] = alpha * (u[im1, j, 2] + u[ip1, j, 2] + u[i, jp1, 2] + u[i, jm1, 2] -
                               4u[i, j, 2]) +
                      A * u[i, j, 1] - u[i, j, 1]^2 * u[i, j, 2]
    end
end
D_dt = similar(u0)
h(D_dt, p, t, ::Type{Val{1}}) = (D_dt .= 0.0)
prod_dde = DDEProblem(brusselator_2d_loop_h, u0, h, (0.0, 20.0), p)
alg = MethodOfSteps(SSPRK43())

@time sol = solve(prod_dde, alg, saveat=0.5, progress=true, progress_steps=1);
2

I don’t think you need a DDE to represent this model. Instead, you can use a DAE (in mass matrix form)

3

No, in general DDEs require a lot more memory because you need to store history for being able to handle past reconstructions based on interpolations. We will likely add some more options in the near future based on maximum delay bounds to better drop history, but even then you should expect the DDE solver to take more memory than the ODE solver because of this. If you think about it, DDEs use memory of the system, so of course we need to have that memory stored (though we could make use of mmap’d arrays to improve this).

I presume this is just an MWE to show that the DDE solver takes more memory than the ODE solver, which of course has to be true given the above.

4

I was basing off of

I have added a component to that system that requires the time derivative info so I changed it to a DDE

Which implies the reason for the DDE is trying to get a derivative that you could get from a DAE.

5

Oh indeed, if it’s just for current time derivative then it should be a DAE, not only for memory but also because it would then be implicit and you’d get a more stable solve.