# Solving ODE inside a double for loop

**URL:** <https://discourse.julialang.org/t/solving-ode-inside-a-double-for-loop/101822>\
**Category:** Performance\
**Created:** [July 20, 2023, 8:29am UTC](https://discourse.julialang.org/t/solving-ode-inside-a-double-for-loop/101822 "2023-07-20T08:29:16Z")\
**Posts on this page:** 10\
**Page:** 1

<div class="post-metadata">

**Author:** ![David\_F](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/david_f/32/42706_2.png) [@David\_F](https://discourse.julialang.org/u/David_F)\
**Post date:** [July 20, 2023, 8:29am UTC](https://discourse.julialang.org/t/solving-ode-inside-a-double-for-loop/101822/1 "2023-07-20T08:29:16Z")

</div>

Hello,

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:

```julia
    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)
            else
               k = Int(ii)-1
            end
                
            ## 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])
        end
              
        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)
            break
        end

```

The ODE function is defined as:

```julia
    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        
    end

```

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [July 20, 2023, 6:33pm UTC](https://discourse.julialang.org/t/solving-ode-inside-a-double-for-loop/101822/2 "2023-07-20T18:33:50Z")

</div>

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

> **[Integrator Interface · DifferentialEquations.jl](https://docs.sciml.ai/DiffEqDocs/stable/basics/integrator/)**
>
> Documentation for DifferentialEquations.jl.

> **[Event Handling and Callback Functions · DifferentialEquations.jl](https://docs.sciml.ai/DiffEqDocs/stable/features/callback_functions/)**
>
> Documentation for DifferentialEquations.jl.

> [@David\_F](#):
>
> ```julia
> T = M\RHS       
> push!(Tsurf, T[1])
> 
> ```

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

> **[Differential Algebraic Equations · DifferentialEquations.jl](https://docs.sciml.ai/DiffEqDocs/stable/tutorials/dae_example/)**
>
> Documentation for DifferentialEquations.jl.

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

---

<div class="post-metadata">

**Author:** ![David\_F](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/david_f/32/42706_2.png) [@David\_F](https://discourse.julialang.org/u/David_F)\
**Post date:** [July 21, 2023, 5:56am UTC](https://discourse.julialang.org/t/solving-ode-inside-a-double-for-loop/101822/3 "2023-07-21T05:56:19Z")

</div>

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](https://encyclopediaofmath.org/wiki/Crank-Nicolson_method) 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.

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [July 21, 2023, 12:31pm UTC](https://discourse.julialang.org/t/solving-ode-inside-a-double-for-loop/101822/4 "2023-07-21T12:31:44Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![David\_F](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/david_f/32/42706_2.png) [@David\_F](https://discourse.julialang.org/u/David_F)\
**Post date:** [July 21, 2023, 6:58pm UTC](https://discourse.julialang.org/t/solving-ode-inside-a-double-for-loop/101822/5 "2023-07-21T18:58:22Z")

</div>

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.

```julia
#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
end

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:

> [IDAS ERROR] IDACalcIC  
> Newton/Linesearch algorithm failed to converge.  
> retcode: InitialFailure

---

<div class="post-metadata">

**Author:** ![contradict](https://avatars.discourse-cdn.com/v4/letter/c/ac91a4/32.png) [@contradict](https://discourse.julialang.org/u/contradict)\
**Post date:** [July 21, 2023, 11:21pm UTC](https://discourse.julialang.org/t/solving-ode-inside-a-double-for-loop/101822/6 "2023-07-21T23:21:49Z")

</div>

I think the `kinetics` function should always compute a residual:

```julia
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

    nothing
end

```

![chemical](https://global.discourse-cdn.com/julialang/original/3X/2/c/2c5227a2bc827f9d18820ee5c82072ef35fccf82.png)

---

<div class="post-metadata">

**Author:** ![David\_F](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/david_f/32/42706_2.png) [@David\_F](https://discourse.julialang.org/u/David_F)\
**Post date:** [July 22, 2023, 12:53pm UTC](https://discourse.julialang.org/t/solving-ode-inside-a-double-for-loop/101822/7 "2023-07-22T12:53:47Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![Oscar\_Smith](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/oscar_smith/32/25343_2.png) [@Oscar\_Smith](https://discourse.julialang.org/u/Oscar_Smith)\
**Post date:** [July 22, 2023, 1:00pm UTC](https://discourse.julialang.org/t/solving-ode-inside-a-double-for-loop/101822/8 "2023-07-22T13:00:16Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![David\_F](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/david_f/32/42706_2.png) [@David\_F](https://discourse.julialang.org/u/David_F)\
**Post date:** [July 22, 2023, 1:06pm UTC](https://discourse.julialang.org/t/solving-ode-inside-a-double-for-loop/101822/9 "2023-07-22T13:06:28Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![contradict](https://avatars.discourse-cdn.com/v4/letter/c/ac91a4/32.png) [@contradict](https://discourse.julialang.org/u/contradict)\
**Post date:** [July 24, 2023, 1:54pm UTC](https://discourse.julialang.org/t/solving-ode-inside-a-double-for-loop/101822/10 "2023-07-24T13:54:56Z")

</div>

Re-creating the problem from scratch can be expensive, can you use [`remake`](https://docs.sciml.ai/DiffEqDocs/stable/basics/problem/#Modification-of-problem-types)? Have you [profiled](https://docs.julialang.org/en/v1/manual/profile/) to see exactly which part of the code is slow?