bugfix and ignition example

example/ignition/ignition.jl

@@ -0,0 +1,70 @@
+using Arrhenius
+using LinearAlgebra
+using DifferentialEquations
+using ForwardDiff
+using DiffEqSensitivity
+using Sundials
+using Plots
+using DelimitedFiles
+using Profile
+
+cd("./example/ignition")
+
+cantera_data = readdlm("ignition.dat")
+ct_ts = cantera_data[:, 1]
+ct_T = cantera_data[:, 2]
+ct_Y = cantera_data[:, 3:end]
+
+gas = CreateSolution("../../mechanism/gri30.yaml")
+ns = gas.n_species
+
+T0 = 980.0
+P = one_atm * 15.0
+
+X0 = zeros(ns);
+X0[species_index(gas, "CH4")] = 1.0 / 2.0
+X0[species_index(gas, "O2")] = 1.0
+X0[species_index(gas, "N2")] = 3.76
+X0 = X0 ./ sum(X0);
+Y0 = X2Y(gas, X0, dot(X0, gas.MW));
+
+# Y0 = ones(ns) ./ns
+u0 = vcat(Y0, T0)
+@inbounds function dudt!(du, u, p, t)
+    T = u[end]
+    Y = @view(u[1:ns])
+    mean_MW = 1.0 / dot(Y, 1 ./ gas.MW)
+    ρ_mass = P / R / T * mean_MW
+    X = Y2X(gas, Y, mean_MW)
+    C = Y2C(gas, Y, ρ_mass)
+    cp_mole, cp_mass = get_cp(gas, T, X, mean_MW)
+    h_mole = get_H(gas, T, Y, X)
+    S0 = get_S(gas, T, P, X)
+    # qdot = wdot_func(gas.reaction, T, C, S0, h_mole; get_qdot=true)
+    wdot = wdot_func(gas.reaction, T, C, S0, h_mole)
+    Ydot = wdot / ρ_mass .* gas.MW
+    Tdot = -dot(h_mole, wdot) / ρ_mass / cp_mass
+    du .= vcat(Ydot, Tdot)
+end
+
+tspan = [0.0, 0.1];
+prob = ODEProblem(dudt!, u0, tspan);
+@time sol = solve(prob, CVODE_BDF(), reltol = 1e-6, abstol = 1e-9)
+
+tmin = 1.e-2
+plt = plot(
+    sol.t .+ tmin,
+    sol[species_index(gas, "H"), :],
+    lw = 2,
+    label = "Arrhenius.jl",
+);
+plot!(plt, ct_ts .+ tmin, ct_Y[:, species_index(gas, "H")], label = "Cantera")
+ylabel!(plt, "Mass Fraction of H")
+xlabel!(plt, "Time [s]")
+pltT = plot(sol.t .+ tmin, sol[end, :], lw = 2, label = "Arrhenius.jl");
+plot!(pltT, ct_ts .+ tmin, ct_T, label = "Cantera")
+ylabel!(pltT, "Temperature [K]")
+xlabel!(pltT, "Time [s]")
+title!(plt, "GRI30 CH4 Ignition @980K/15atm")
+pltsum = plot(plt, pltT, xscale = :identity, legend = :left, framestyle = :box)
+png(pltsum, "CH4_ignition.png")

example/ignition/ignition.py

@@ -0,0 +1,46 @@
+import sys
+import os
+import numpy as np
+import matplotlib.pyplot as plt
+from skopt.space import Space
+from skopt.sampler import Lhs
+
+import cantera as ct
+
+np.random.seed(1)
+
+gas = ct.Solution('../../mechanism/gri30.yaml')
+gas.TPX = 980, 15*ct.one_atm, 'CH4:1.0,O2:2.0,N2:7.52'
+
+r = ct.IdealGasConstPressureReactor(gas, energy='on')
+sim = ct.ReactorNet([r])
+time = 0.0
+states = ct.SolutionArray(gas, extra=['t'])
+
+print('%10s %10s %10s %14s' % ('t [s]', 'T [K]', 'P [Pa]', 'u [J/kg]'))
+while sim.time < 1.e-1:
+    states.append(r.thermo.state, t=sim.time)
+    sim.step()
+
+nodedata = np.vstack((states.t, states.T, states.Y.T)).T
+np.savetxt('ignition.dat', nodedata)
+
+plt.clf()
+plt.subplot(2, 2, 1)
+plt.plot(states.t, states.T, 'o')
+plt.xlabel('Time (ms)')
+plt.ylabel('Temperature (K)')
+plt.subplot(2, 2, 2)
+plt.plot(states.t, states.Y[:, gas.species_index('H2')])
+plt.xlabel('Time (ms)')
+plt.ylabel('CH4 Mole Fraction')
+plt.subplot(2, 2, 3)
+plt.plot(states.t, states.Y[:, gas.species_index('O2')])
+plt.xlabel('Time (ms)')
+plt.ylabel('O2 Mole Fraction')
+plt.subplot(2, 2, 4)
+plt.plot(states.t, states.Y[:, gas.species_index('H2O')])
+plt.xlabel('Time (ms)')
+plt.ylabel('CO2 Mole Fraction')
+plt.tight_layout()
+plt.show()

src/Arrhenius.jl

@@ -19,7 +19,7 @@ Following codes are used during development phase only.
 # using Profile
 # using ForwardDiff
 # using Zygote
-
+#
 # gas = Arrhenius.CreateSolution("./mechanism/gri30.yaml")
 # ns = gas.n_species
 # @show gas.thermo.nasa_high[Arrhenius.species_index(gas, "O2"), :]
@@ -30,17 +30,6 @@ Following codes are used during development phase only.
 # P = Arrhenius.one_atm
 # wdot = Arrhenius.set_states(gas, T0, P, Y0)
 # @show wdot
-
-# function profile_test(n)
-#     for i = 1:n
-#         Arrhenius.set_states(gas, T0, P, Y0)
-#     end
-# end
-# @time profile_test(1000)
-#
-# # Profile.clear
-# # @profile profile_test(1000)
-# # Juno.profiler(; C = false)
 #
 # u0 = vcat(Y0, T0);
 # R = Arrhenius.R
@@ -58,12 +47,23 @@ Following codes are used during development phase only.
 #     wdot = Arrhenius.wdot_func(gas.reaction, T, C, S0, h_mole)
 #     return wdot[end]
 # end
+# @time grad = ForwardDiff.gradient(f, u0)
+
+# function profile_test(n)
+#     for i = 1:n
+#         Arrhenius.set_states(gas, T0, P, Y0)
+#     end
+# end
+# @time profile_test(1000)
+#
+# # Profile.clear
+# # @profile profile_test(1000)
+# # Juno.profiler(; C = false)
+
 # u0[end] = 1200 + rand()
 # @time f(u0)
 # @time f(u0)
 #
-# @time grad = ForwardDiff.gradient(f, u0)
-# @time grad = ForwardDiff.gradient(f, u0)
 
 # Profile.clear
 # @profile grad = ForwardDiff.gradient(x -> f(x), u0)

src/DataStructure.jl

@@ -4,13 +4,8 @@ struct Reaction
     reactant_orders::Array{Float64, 2}
     is_reversible::Array{Bool, 1}
     Arrhenius_coeffs::Array{Float64, 2}
-    Arrhenius_A0::Array{Float64, 1}
-    Arrhenius_b0::Array{Float64, 1}
-    Arrhenius_Ea0::Array{Float64, 1}
-    Troe_A::Array{Float64, 1}
-    Troe_T1::Array{Float64, 1}
-    Troe_T2::Array{Float64, 1}
-    Troe_T3::Array{Float64, 1}
+    Arrhenius_0::Array{Float64, 2}
+    Troe_::Array{Float64, 2}
     index_three_body::Array{Int64, 1}
     index_falloff::Array{Int64, 1}
     index_falloff_troe::Array{Int64, 1}
@@ -25,6 +20,8 @@ end
 struct Thermo
     nasa_low::Array{Float64, 2}
     nasa_high::Array{Float64, 2}
+    Trange::Array{Float64, 2}
+    isTcommon::Bool
 end
 
 struct Solution

src/Kinetics.jl

@@ -10,29 +10,24 @@ function wdot_func(reaction, T, C, S0, h_mole; get_qdot=false)
         @inbounds _kf[i] *= dot(@view(reaction.efficiencies_coeffs[:, i]), C)
     end
 
-    j = 1
-    for i in reaction.index_falloff
-        @inbounds A0 = reaction.Arrhenius_A0[j]
-        @inbounds b0 = reaction.Arrhenius_b0[j]
-        @inbounds Ea0 = reaction.Arrhenius_Ea0[j]
+    for (j, i) in enumerate(reaction.index_falloff)
+        @inbounds A0, b0, Ea0 = reaction.Arrhenius_0[j, :]
         @inbounds k0 = A0 * exp(b0 * log(T) - Ea0 * 4184.0 / R / T)
         @inbounds Pr =
             k0 * dot(@view(reaction.efficiencies_coeffs[:, i]), C) / _kf[i]
         lPr = log10(Pr)
         _kf[i] *= (Pr / (1 + Pr))
 
-        if (i in reaction.index_falloff_troe) & (reaction.Troe_A[j] > 1.e-12)
-            # TODO: (reaction.Troe_A[j] > 1.e-12) is for compatity, will be removed
+        if (reaction.Troe_[j, 1] > 1.e-12)
             @inbounds F_cent =
-                (1 - reaction.Troe_A[j]) * exp(-T / reaction.Troe_T3[j]) +
-                reaction.Troe_A[j] * exp(-T / reaction.Troe_T1[j]) +
-                exp(-reaction.Troe_T2[j] / T)
+                (1 - reaction.Troe_[j, 1]) * exp(-T / reaction.Troe_[j, 4]) +
+                reaction.Troe_[j, 1] * exp(-T / reaction.Troe_[j, 2]) +
+                exp(-reaction.Troe_[j, 3] / T)
             lF_cent = log10(F_cent)
             _C = -0.4 - 0.67 * lF_cent
             N = 0.75 - 1.27 * lF_cent
             @inbounds f1 = (lPr + _C) / (N - 0.14 * (lPr + _C))
             @inbounds _kf[i] *= exp(log(10.0) * lF_cent / (1 + f1^2))
-            j = j + 1
         end
     end
 

src/Solution.jl

@@ -13,13 +13,16 @@ function CreateSolution(mech)
 
     nasa_low = zeros(n_species, 7)
     nasa_high = zeros(n_species, 7)
+    Trange = zeros(n_species, 3)
 
-    _species_names = [yaml["species"][i]["name"] for i in 1:length(yaml["species"])]
+    _species_names =
+        [yaml["species"][i]["name"] for i = 1:length(yaml["species"])]
 
     for (i, species) in enumerate(species_names)
         spec = yaml["species"][findfirst(x -> x == species, _species_names)]
         nasa_low[i, :] = spec["thermo"]["data"][1]
         nasa_high[i, :] = spec["thermo"]["data"][2]
+        Trange[i, :] .= spec["thermo"]["temperature-ranges"]
 
         for j = 1:n_elements
             if haskey(spec["composition"], elements[j])
@@ -28,7 +31,9 @@ function CreateSolution(mech)
         end
     end
 
-    thermo = Thermo(nasa_low, nasa_high)
+    isTcommon = (maximum(Trange[:, 2]) - minimum(Trange[:, 2])) < 0.01
+
+    thermo = Thermo(nasa_low, nasa_high, Trange, isTcommon)
 
     npz = npzread("$mech.npz")
     MW = npz["molecular_weights"]
@@ -42,23 +47,31 @@ function CreateSolution(mech)
         Arrhenius_A0 = npz["Arrhenius_A0"]
         Arrhenius_b0 = npz["Arrhenius_b0"]
         Arrhenius_Ea0 = npz["Arrhenius_Ea0"]
+    else
+        Arrhenius_A0 = []
+        Arrhenius_b0 = []
+        Arrhenius_Ea0 = []
+    end
+
+    if haskey(npz, "Troe_A")
         Troe_A = npz["Troe_A"]
         Troe_T1 = npz["Troe_T1"]
         Troe_T2 = npz["Troe_T2"]
         Troe_T3 = npz["Troe_T3"]
     else
-        Arrhenius_A0 = []
-        Arrhenius_b0 = []
-        Arrhenius_Ea0 = []
         Troe_A = []
         Troe_T1 = []
         Troe_T2 = []
         Troe_T3 = []
     end
+    Arrhenius_0 = hcat(Arrhenius_A0, Arrhenius_b0, Arrhenius_Ea0)
+    Troe_ = hcat(Troe_A, Troe_T1, Troe_T2, Troe_T3)
 
     index_three_body = []
     index_falloff = []
     index_falloff_troe = []
+    _ind_troe = []
+    j = 0
     for i = 1:n_reactions
         reaction = yaml["reactions"][i]
         if haskey(reaction, "type")
@@ -67,13 +80,21 @@ function CreateSolution(mech)
             end
             if yaml["reactions"][i]["type"] == "falloff"
                 push!(index_falloff, i)
+                j = j + 1
                 if haskey(yaml["reactions"][i], "Troe")
                     push!(index_falloff_troe, i)
+                    push!(_ind_troe, j)
                 end
             end
         end
     end
 
+    if length(index_falloff) > length(index_falloff_troe)
+        _Troe = zeros(length(index_falloff), 4) .+ 1.e-16
+        _Troe[_ind_troe, :] .= Troe_
+        Troe_ = _Troe
+    end
+
     i_reactant = []
     i_product = []
     for i = 1:n_reactions
@@ -82,21 +103,16 @@ function CreateSolution(mech)
     end
 
     vk = product_stoich_coeffs - reactant_stoich_coeffs
-    vk_sum = sum(vk, dims=1)[1, :]
+    vk_sum = sum(vk, dims = 1)[1, :]
 
     reaction = Reaction(
         product_stoich_coeffs,
         reactant_stoich_coeffs,
         reactant_orders,
         is_reversible,
         Arrhenius_coeffs,
-        Arrhenius_A0,
-        Arrhenius_b0,
-        Arrhenius_Ea0,
-        Troe_A,
-        Troe_T1,
-        Troe_T2,
-        Troe_T3,
+        Arrhenius_0,
+        Troe_,
         index_three_body,
         index_falloff,
         index_falloff_troe,

src/Thermo.jl

@@ -1,9 +1,16 @@
 "get specific of heat capacity"
 function get_cp(gas, T, X, mean_MW)
+    cp_T = [1.0, T, T^2, T^3, T^4]
     if T <= 1000.0
-        cp = @view(gas.thermo.nasa_low[:, 1:5]) * [1.0, T, T^2, T^3, T^4]
+        cp = @view(gas.thermo.nasa_low[:, 1:5]) * cp_T
     else
-        cp = @view(gas.thermo.nasa_high[:, 1:5]) * [1.0, T, T^2, T^3, T^4]
+        cp = @view(gas.thermo.nasa_high[:, 1:5]) * cp_T
+    end
+    # TODO: not sure if inplace operation will be an issue for AD
+    if !gas.thermo.isTcommon
+        ind_correction = @. (T > 1000.0) & (T < gas.thermo.Trange[:, 2])
+        cp[ind_correction] .=
+            @view(gas.thermo.nasa_low[ind_correction, 1:5]) * cp_T
     end
     cp_mole = dot(cp, X) * R
     cp_mass = cp_mole / mean_MW
@@ -19,6 +26,11 @@ function get_H(gas, T, Y, X)
     else
         h_mole = @view(gas.thermo.nasa_high[:, 1:6]) * H_T * R * T
     end
+    if !gas.thermo.isTcommon
+        ind_correction = @. (T > 1000.0) & (T < gas.thermo.Trange[:, 2])
+        h_mole[ind_correction] .=
+            @view(gas.thermo.nasa_low[ind_correction, 1:6]) * H_T * R * T
+    end
     # H_mole = dot(h_mole, X)
     return h_mole
 end
@@ -37,6 +49,13 @@ function get_S(gas, T, P, X)
     else
         S0 = @view(gas.thermo.nasa_high[:, [1, 2, 3, 4, 5, 7]]) * S_T * R
     end
+    if !gas.thermo.isTcommon
+        ind_correction = @. (T > 1000.0) & (T < gas.thermo.Trange[:, 2])
+        S0[ind_correction] .=
+            @view(gas.thermo.nasa_low[ind_correction, [1, 2, 3, 4, 5, 7]]) *
+            S_T *
+            R
+    end
     # _X = @. S0 - R * log(clamp(X, 1.e-30, Inf))
     # s_mole = _X .- R * (P / one_atm)
     # S_mole = dot(s_mole, X)