Solving ODE inside a double for loop


I have two for loops, say the outer one iterates on a time variable, and the inner one iterates on a position variable.

Inside the position loop, I call ODEproblem() and solve() to get the solution for a system of ODE. The ODE function takes in a vector (T[rr]), which is calculated in the outer (time) loop, then solves the system for every value of the position in the inner loop, and outputs a vector, which is needed in the next iteration of the outer (time) loop.

Even when I set maxiters=1 and use a high reltol for the solve() function, there’s a huge slow-down of the overall code. This probably means that some allocations or other internal operations take a lot of time. Can anyone suggest a way of doing this better?

The relevant part of the code is the following:

    for ii = (1:nₜ)
        ra = sol(ii*dt)[2]                  
        rf = b_radius(sol(ii*dt)[2])[1]      
        r = range(ra, stop=rf, length=nᵣ)
        r = Vector(r)
        dr = (rf-ra)/(nᵣ - 1)                
        A = αD*dt/(dr^2)                     
        B = αD*dt/(2*dr)      
        C = (sol(ii*dt, Val{1})[2].*sol(ii*dt)[2].^2) * dt/(4*dr)
        MainD = (1+A).*ones(nᵣ-1)
        insert!(MainD, 1, 1+A)
        downD = -(A/2 .+ B./r[2:nᵣ-1] - C./r[2:nᵣ-1].^2).*ones(nᵣ-2)
        insert!(downD, nᵣ-1, -A)
        upD = -(A/2 .- B./r[2:nᵣ-1] + C./r[2:nᵣ-1].^2).*ones(nᵣ-2)
        insert!(upD, 1, -A)

        M = spdiagm(
             0 => MainD,
            -1 => downD,
             1 => upD,
        m = Int(ii)
        for rr = (1:nᵣ)    
            if Int(ii) == 1
               k = Int(ii)
               k = Int(ii)-1
            ## this part causes huge slowdown ##
            chemprob = ODEProblem(kinetics!, [Y1_vs_R_t[k, rr]; Y2_vs_R_t[k, rr]; Y3_vs_R_t[k, rr]; Y4_vs_R_t[k, rr]], ii*dt, T[rr])
            chemsol = solve(chemprob, KenCarp4(), reltol=1.0e-1, maxiters = 1)

            Y1_vs_R_t[m, rr] = chemsol.u[end][1]
            Y2_vs_R_t[m, rr] = chemsol.u[end][2]
            Y3_vs_R_t[m, rr] = chemsol.u[end][3]
            Y4_vs_R_t[m, rr] = chemsol.u[end][4]
            R1_vs_R_t[rr] = exp(lnZ₁).*exp(-Eₐ₁./(R_gas.*T[rr])).*Y1_vs_R_t[rr]
            R2_vs_R_t[rr] = exp(lnZ₂).*exp(-Eₐ₂./(R_gas.*T[rr])).*Y2_vs_R_t[rr]
            R3_vs_R_t[rr] = exp(lnZ₃).*exp(-Eₐ₃./(R_gas.*T[rr])).*Y3_vs_R_t[rr]
            Q_vs_R_t[rr] = (1.0/(ρ.*Cᵥ)).*(Q₁.*R1_vs_R_t[rr] .+ Q₂.*R2_vs_R_t[rr] .+ Q₃.*R3_vs_R_t[rr])
        Q1 = Q_vs_R_t[2:nᵣ-1]
        Q2 = Q_vs_R_t[1]
        Q3 = Q_vs_R_t[nᵣ]
        RHS = (A/2 .- B./r[2:nᵣ-1] + C./r[2:nᵣ-1].^2).*T[1:nᵣ-2] .+ (1-A).*T[2:nᵣ-1] .+ (A/2 .+ B./r[2:nᵣ-1] - C./r[2:nᵣ-1].^2).*T[3:nᵣ] .+ dt.*dissipation(ii*dt, r[2:nᵣ-1], Y_new[2:nᵣ-1]) .+ dt.*Q1
        insert!(RHS, 1,  (1-A).*T[1]  .+ (A .+ 2*B./r[1] - C./r[1].^2).*(2*dr*φa/λₛ)  .+ A.*T[2] .+ dt.*dissipation(ii*dt, r[1], Y_new[1]) .+ dt.*Q2)                                                   
        insert!(RHS, nᵣ, (1-A).*T[nᵣ] .+ (A .+ 2*B./r[nᵣ] - C./r[nᵣ].^2).*(2*dr*φb/λₛ) .+ A.*T[nᵣ-1] .+ dt.*dissipation(ii*dt, r[nᵣ], Y_new[nᵣ]) .+ dt.*Q3)                                             

        T = M\RHS       
        push!(Tsurf, T[1])

        if T[1] >= 2000.0
            time_Pcrit = Int(ii)

The ODE function is defined as:

    function kinetics!(dY, Y, p, t)
      R1 = (exp(lnZ₁).*exp(-Eₐ₁/(R_gas.*p))).*Y[1]       
      R2 = (exp(lnZ₂).*exp(-Eₐ₂/(R_gas.*p))).*Y[2]       
      R3 = (exp(lnZ₃).*exp(-Eₐ₃/(R_gas.*p))).*Y[3].^2    
      dY[1] = (-1.0/ρ).*R1
      dY[2] = (1.0/ρ).*(R1 - R2)
      dY[3] = (1.0/ρ).*(R2 - R3)
      dY[4] = (1.0/ρ).*R3        

Why are you solving one step at a time? You may want to consider using a callback, or the integrator interface.

It looks like you’re trying to define some variables algebraically? Your system would be a DAE.

Only imposing the DAE condition at dt points will give only O(dt) convergence, so you’ll want to do this more correctly.

Thanks! I will try to give some context, so hopefully the issue will be easier to help with.

Why are you solving one step at a time? You may want to consider using a callback, or the integrator interface.

Completely ignoring the inner for loop, the code solves the transient thermal conduction problem with heat source terms using the Crank-Nicolson method. This provides the temperature field (in 1D) for every instant in time from 1:nₜ.

Since one of the heat sources (Q_vs_R_t[rr]) comes from a chemical reaction (which in turn depends on the temperature calculated previously), I need to solve the kinetics for each position (which is taken care of in the inner loop), and add the resultant heat back to the heat conduction equation to update the temperature field for the next time step.

The T = M\RHS part provides the temperature field in the given time step.
push!(Tsurf, T[1]) takes just the temperature value in the leftmost position because it is of physical interest, but this line is irrelevant for the issue.

In short, the two calculations are coupled (temperature and kinetics of heat generation), which is why I need to solve the kinetics at each new time step.

I hope this explains the situation.

Yes, that’s a DAE, and Crank-Nicolson isn’t a good way to do that. Have you checked using DAEProblem with IDA?

I’m not sure how to properly set up a DAE problem for my case. I have tried the following setup just to see what I get for the kinetic part, but I get an error right at the start.

#display(HTML("<style>.container { width:100% !important; }</style>"))
using Revise
using DifferentialEquations

ρ = 1670.0                                   
Eₐ₁ = 198.74                     
Eₐ₂ = 188.28                      
Eₐ₃ = 143.09                      
lnZ₁ = 45.2                         
lnZ₂ = 40.0                           
lnZ₃ = 34.9                         
R_gas = 0.00831447             

function kinetics(out, dY, Y, p, t)
    TT = p
    R1 = (exp(lnZ₁).*exp(-Eₐ₁/(R_gas.*TT))).*Y[1]       
    R2 = (exp(lnZ₂).*exp(-Eₐ₂/(R_gas.*TT))).*Y[2]       
    R3 = (exp(lnZ₃).*exp(-Eₐ₃/(R_gas.*TT))).*Y[3].^2    

    out[1] = (-1.0/ρ).*R1
    out[2] = (1.0/ρ).*(R1 - R2)
    out[3] = (1.0/ρ).*(R2 - R3)
    out[4] = Y[1] + Y[2] + Y[3] + Y[4] - 1.0

u₀ = [1.0, 0.0, 0.0, 0.0]
du₀ = ones(length(u₀))
tspan = (0.0, 1.0e-6)

differential_vars = [true, true, true, false]
chemprob = DAEProblem(kinetics, du₀, u₀, tspan, 1000.0, differential_vars = differential_vars)

using Sundials
sol = solve(chemprob, IDA())

The above code results in the following error:

Newton/Linesearch algorithm failed to converge.
retcode: InitialFailure

I think the kinetics function should always compute a residual:

function kinetics(out, dY, Y, p, t)
    TT = p
    R1 = (exp(lnZ₁).*exp(-Eₐ₁/(R_gas.*TT))).*Y[1]
    R2 = (exp(lnZ₂).*exp(-Eₐ₂/(R_gas.*TT))).*Y[2]
    R3 = (exp(lnZ₃).*exp(-Eₐ₃/(R_gas.*TT))).*Y[3].^2

    out[1] = dY[1] - (-1.0/ρ).*R1
    out[2] = dY[2] - (1.0/ρ).*(R1 - R2)
    out[3] = dY[3] - (1.0/ρ).*(R2 - R3)
    out[4] = Y[1] + Y[2] + Y[3] + Y[4] - 1.0



1 Like

Thanks, it works now.
However, this way the calculation is much slower compared with the original method.

your code is now type understand because you have a ton of non combat globals

I don’t think that declaring variable types is really the issue here. The code runs pretty fast as long as I don’t call the ODEProblem(…) and solve(…) functions inside the loops even if I use maxiter=1, which means that it’s not the ODE solution process itself which causes the slow-down but something else.

Re-creating the problem from scratch can be expensive, can you use remake? Have you profiled to see exactly which part of the code is slow?