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:
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:
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
