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.

6

Hey, I recently stumbled on this thread as I have similar problem as jarias had: changing from ODEProblem to DDEProblem to get the current (first and second) time derivative. I looked up DAEs but could not find any info how to utilize them in this kind of situation. Could you please explain a bit further how one would do this? Thank you.

7

Can you give an MWE of what you’re trying to do?

8

Here’s an MWE of my current DDE system. Essentially I’m using the DDE part to calculate the first and second derivative of rho0vac through a few hoops.

using DifferentialEquations

function solve_equations()
    k = logrange(1, 10^4, 400)
    kav = (k[1:end - 1] + k[2:end]) / 2
    kavlength = length(kav)

    iconds = initialConditions(kavlength)
    h(out, p, t) = (out .= iconds)
    tspan = (0, 4)

    out1 = zeros(size(iconds))
    out2 = zeros(size(iconds))
    tau = 1e-3
    dddeltarho0 = zeros(kavlength)
    rho0vac = zeros(kavlength)
    rho0vacprev = zeros(kavlength)
    rho0vacprev2 = zeros(kavlength)
    omega2 = zeros(kavlength)

    prob = DDEProblem(diffEquations!, iconds, h, tspan, (kav, kavlength, out1, out2, tau,
                      dddeltarho0, rho0vac, rho0vacprev, rho0vacprev2, omega2), constant_lags=[tau])
    sol = solve(prob, MethodOfSteps(Vern7()), abstol=1e-15, reltol=1e-15)
end

function initialConditions(kavlength)
    deltarho0 = zeros(kavlength)
    ddeltarho0 = zeros(kavlength)
    arr = [15, -0.8, 1, 0.001, 0]
    initial = vcat(arr, deltarho0, ddeltarho0)
    return initial
end

function diffEquations!(du, u, h, args, t)
    kav, kavlength, out1, out2, tau, dddeltarho0, rho0vac, rho0vacprev, rho0vacprev2, omega2 = args
    h(out1, args, t - tau)
    h(out2, args, t - 2 * tau)
    
    infl, dinfl, a, sigma, dsigma = @view u[1:5]
    deltarho0 = @view u[6:6 + kavlength - 1]
    ddeltarho0 = @view u[6 + kavlength:6 + 2 * kavlength - 1]

    H = sqrt((dinfl^2 + infl^2) / 6)
    R = dinfl^2 - 2 * infl^2
    Rprev = out1[2]^2 - 2 * out1[1]^2
    Rprev2 = out2[2]^2 - 2 * out2[1]^2

    m2Eff = a^2 * ((1e-11)^2 + R * (1 / 6 - 0.2))
    m2Effprev = out1[3]^2 * ((1e-11)^2 + Rprev * (1 / 6 - 0.2))
    m2Effprev2 = out2[3]^2 * ((1e-11)^2 + Rprev2 * (1 / 6 - 0.2))
    
    for i in eachindex(rho0vac)
        omega2[i] = kav[i] * kav[i] + m2Eff
        rho0vac[i] = 1 / (2 * sqrt(omega2[i]))
        rho0vacprev[i] = 1 / (2 * sqrt(kav[i] * kav[i] + m2Effprev))
        rho0vacprev2[i] = 1 / (2 * sqrt(kav[i] * kav[i] + m2Effprev2))
    end

    for i in eachindex(rho0vac)
        drho0vac = (rho0vac[i] - rho0vacprev[i]) / tau
        ddrho0vac = (rho0vac[i] - 2 * rho0vacprev[i] + rho0vacprev2[i]) / tau^2
        deltarho2 = omega2[i] * deltarho0[i]^2 +
                    a^2 * 0.25 * (drho0vac + ddeltarho0[i])^2 / (2 * (rho0vac[i] + deltarho0[i]))
        dddeltarho0[i] = 4 * (deltarho2 - omega2[i] * deltarho0[i]) / a^2 -
                         H * (drho0vac + ddeltarho0[i]) - ddrho0vac
    end

    ddinfl = -3 * H * dinfl - infl
    da = a * H
    ddsigma = - H * dsigma - m2Eff * sigma / a^2

    du[1] = dinfl
    du[2] = ddinfl
    du[3] = da
    du[4] = dsigma
    du[5] = ddsigma
    du[6:6 + kavlength - 1] .= ddeltarho0
    du[6 + kavlength:6 + 2 * kavlength - 1] .= dddeltarho0
end

solve_equations()
9

What is the system of equations here?

10

The system of DEs is

\begin{cases} \ddot\phi=-3H\dot{\phi}-\phi\\ \dot{a}=Ha\\ \ddot\sigma=-H\dot{\sigma}-\frac{m_{Eff}^2}{a^2}\sigma\\ \ddot{\delta\rho}_0=\frac{4\left(\delta\rho_2-\omega^2\delta\rho_0\right)}{a^2}-H\left(\dot{\rho}_{0{,}vac}+\dot{\delta\rho}_0\right)-\ddot\rho_{0{,}vac} \end{cases}

where

\begin{cases} H=\sqrt{\frac{\dot{\phi }^2+\phi ^2}{6}}\\ R=\dot{\phi }^2-2\phi\\ m_{Eff}^2=a^2\left(m_R^2+R\left(\frac{1}{6}-\xi \right)\right)\\ \omega ^2=k^2+m_{Eff}^2\\ \rho_{0,vac}=\frac{1}{2\omega }\\ \delta \rho _2=\frac{1}{2\left(\rho _{0{,}vac}+\delta \rho _0\right)}\left(\omega ^2\delta \rho _0^2+\frac{1}{4}a^2\left(\dot{\rho }_{0{,}vac}+\dot{\delta \rho }_0\right)^2\right). \end{cases}

The fourth DE is evaluated with different values of k (in the above code for 400 values).

11

Instead of getting your derivatives by making this a delay equation, you can get them just by adding variables.

so specifically, instead of

\ddot\phi=-3H\dot{\phi}-\phi

you will have

\begin{cases} \dot{d\phi}=-3H\dot{\phi}-\phi\\ \dot{\phi} = d\phi \end{cases}

etc

1 Like
12

Does this even need to be a DAE? Looks like a normal 2nd order ODE to me.

1 Like
13

Thank you for your replies. I feel like I might not have articulated the problem well enough. I am using the delay functions only to calculate the derivatives of \rho_{0,vac}: \dot{\rho}_{0,vac} and \ddot{\rho}_{0,vac}. The reason I (believe I) can’t use the analytical derivatives is that I’ll use parameters for which m_{Eff}^2 is negative and |m_{Eff}^2| grows with time. Thus \ddot{\rho}_{0,vac}, for example, would have a term \frac{\Theta(\omega^2)}{\omega^5} where \omega=0 for some k s inside the region where the DEs are solved. I don’t know if such DEs with singularities are solvable. I have a way of limiting \frac{1}{2\omega} so that it doesn’t blow up but don’t think it extends to higher powers, not with reasonable accuracy at least.