code optimization

example/pyrolysis.jl

@@ -25,8 +25,7 @@ P = one_atm
 
 set_states(gas, mgas, T0, P, Y0)
 
-@show mgas.ρ_mass
-@show mgas.wdot
+@show mgas.ρ_mass mgas.wdot
 
 u0 = vcat(Y0, T0)
 

src/Arrhenius.jl

@@ -16,26 +16,26 @@ end
 # using DifferentialEquations
 # using Sundials
 # using Profile
-
+#
 # gas = Arrhenius.CreateSolution("./mechanism/JP10skeletal.yaml")
 # mgas = Arrhenius.CreateMSolution(gas)
-
+#
 # ns = gas.n_species
-
+#
 # Y0 = zeros(ns)
 # Y0[1] = 1.0
 # T0 = 1200.0
 # P = Arrhenius.one_atm
 # Arrhenius.set_states(gas, mgas, T0, P, Y0)
 # @show mgas.ρ_mass mgas.wdot
-
+#
 # function profile_test(n)
 #     for i = 1:n
 #         Arrhenius.set_states(gas, mgas, T0, P, Y0)
 #     end
 # end
 # @time profile_test(1000)
-
+#
 # Profile.clear
 # @profile profile_test(1000)
 # Juno.profiler(; C = false)

src/Kinetics.jl

@@ -1,13 +1,13 @@
 function wdot_func(gas, mgas)
     T = mgas.T
     C = mgas.C
-    reaction = gas.reaction
 
-    mgas.kf = @. reaction.Arrhenius_coeffs[:, 1] * exp(reaction.Arrhenius_coeffs[:, 2] * log(T) -
-                                       reaction.Arrhenius_coeffs[:, 3] * (4184.0 / R / T))
+    mgas.kf = @. @view(reaction.Arrhenius_coeffs[:, 1]) *
+                                       exp(@view(reaction.Arrhenius_coeffs[:, 2]) * log(T) -
+                                       @view(reaction.Arrhenius_coeffs[:, 3]) * (4184.0 / R / T))
 
     for i in reaction.index_three_body
-        mgas.kf[i] *= dot(reaction.efficiencies_coeffs[:, i], C)
+        mgas.kf[i] *= dot(@view(reaction.efficiencies_coeffs[:, i]), C)
     end
 
     jj = 1
@@ -17,7 +17,7 @@ function wdot_func(gas, mgas)
         b0 = reaction.Arrhenius_b0[j]
         Ea0 = reaction.Arrhenius_Ea0[j]
         k0 = A0 * exp(b0 * log(T) - Ea0 * 4184.0 / R / T)
-        Pr = k0 * dot(reaction.efficiencies_coeffs[:, i], C) / kinf
+        Pr = k0 * dot(@view(reaction.efficiencies_coeffs[:, i]), C) / kinf
         lPr = log10(Pr)
 
         mgas.kf[i] *= (Pr / (1 + Pr))
@@ -40,24 +40,40 @@ function wdot_func(gas, mgas)
         mgas.kf[i] *= exp(log(10.0) * lF_cent / (1 + f1^2))
     end
 
-    vk = reaction.product_stoich_coeffs - reaction.reactant_stoich_coeffs
+    # vk = reaction.product_stoich_coeffs - reaction.reactant_stoich_coeffs
     ΔS_R = vk' * mgas.S0 / R
     ΔH_RT = vk' * mgas.h_mole / (R * T)
     Keq = exp.(ΔS_R .- ΔH_RT .+ log(one_atm / R / T) .* sum(vk, dims=1)[1, :])
     mgas.kr = @. mgas.kf / Keq * reaction.is_reversible
-    for i in 1:gas.n_reactions
-        rop_f = mgas.kf[i]
-        for j in reaction.i_reactant[i]
-            rop_f *= C[j]^reaction.reactant_orders[j, i]
-        end
 
-        rop_r = mgas.kr[i]
-        for j in reaction.i_product[i]
-            rop_r *= C[j]^reaction.product_stoich_coeffs[j, i]
+    # TODO: computing qdot seems to be the bottle neck
+    function qdot_func!(rop_f, rop_r, i, kfi, kri)
+        if length(reaction.i_reactant[i]) == 1
+            ind = reaction.i_reactant[i][1]
+            rop_f = kfi * C[ind]^ind
+        elseif length(reaction.i_reactant[i]) == 2
+            ind = reaction.i_reactant[i]
+            rop_f = kfi * C[ind[1]]^ind[1] * C[ind[2]]^ind[2]
+        end
+        if reaction.is_reversible[i]
+            if length(reaction.i_product[i]) == 1
+                ind = reaction.i_product[i][1]
+                rop_r = kri * C[ind]^ind
+            elseif length(reaction.i_product[i]) == 2
+                ind = reaction.i_product[i]
+                rop_r = kri * C[ind[1]]^ind[1] * C[ind[2]]^ind[2]
+            end
+        else
+            rop_r = 0.0
         end
+    end
 
+    rop_f = 0.0
+    rop_r = 0.0
+    for i in 1:gas.n_reactions
+        qdot_func!(rop_f, rop_r, i, mgas.kf[i], mgas.kr[i])
         mgas.qdot[i] = rop_f - rop_r
     end
 
-    mgas.wdot = (reaction.product_stoich_coeffs - reaction.reactant_orders) * mgas.qdot
+    mgas.wdot = vk * mgas.qdot
 end

src/Solution.jl

@@ -78,5 +78,7 @@ function CreateSolution(mech)
                         index_three_body, index_falloff, efficiencies_coeffs,
                         i_reactant, i_product)
 
+    global vk = product_stoich_coeffs - reactant_stoich_coeffs
+    global reaction
     gas = Solution(n_species, n_reactions, molecular_weights, species_names, thermo, reaction)
 end