Results sensitive on inital composition and setting parameters

[Zeracesharon]

See the testing code, and try to evaluate different initial mole fraction on N2, you will find various digree of mismatch on the results from cantera and reactorch. This phenomenon happends for reduced mechanism for IC8H18 as well as NC7H16.

Runing Solution_test with the same initial parameters, with tolerance delta=1e-4， No obvious difference is observed, which may means the kf (forward_rate_constants) kc (equilibrium_constants) kr (reverse_rate_constants) are correct.

The reasons of such mismatch are still unknown.

Testing Code:

"""
Solve a constant pressure ignition problem where the governing equations are
implemented in Python.

This demonstrates an approach for solving problems where Cantera's reactor
network model cannot be configured to describe the system in question. Here,
Cantera is used for evaluating thermodynamic properties and kinetic rates while
an external ODE solver is used to integrate the resulting equations. In this
case, the SciPy wrapper for VODE is used, which uses the same variable-order BDF
methods as the Sundials CVODES solver used by Cantera.
"""

# TODO: the reactorch class seems to be very slow here, will figure out later

import cantera as ct
import numpy as np
import reactorch as rt
from scipy.integrate import solve_ivp
import torch
from torch.autograd.functional import jacobian as jacobian
from time import perf_counter




class ReactorOde(object):
    def __init__(self, gas):
        # Parameters of the ODE system and auxiliary data are stored in the
        # ReactorOde object.
        self.gas = gas
        self.P = gas.P

    def __call__(self, t, y):
        """the ODE function, y' = f(t,y) """

        # State vector is [T, Y_1, Y_2, ... Y_K]
        self.gas.set_unnormalized_mass_fractions(y[1:])
        self.gas.TP = y[0], self.P
        rho = self.gas.density

        wdot = self.gas.net_production_rates
        dTdt = - (np.dot(self.gas.partial_molar_enthalpies, wdot) /
                  (rho * self.gas.cp))
        dYdt = wdot * self.gas.molecular_weights / rho
        #print('cantera',np.hstack((dTdt, dYdt)))

        return np.hstack((dTdt, dYdt))


class ReactorOdeRT(object):
    def __init__(self, sol):
        # Parameters of the ODE system and auxiliary data are stored in the
        # ReactorOde object.
        self.sol = sol
        self.gas = sol.gas

    def __call__(self, t, y):
        """the ODE function, y' = f(t,y) """

        TPY = torch.Tensor(y).T

        with torch.no_grad():

            self.sol.set_states(TPY)

            TYdot = self.sol.TYdot_func()
        #print('reactorch',TYdot.T)
        
        return TYdot.T.detach().cpu().numpy()

    def TYdot_jac(self, TPY):

        TPY = torch.Tensor(TPY).unsqueeze(0)

        sol.set_states(TPY)

        return sol.TYdot_func().squeeze(0)

    def jac(self, t, y):

        TPY = torch.Tensor(y)

        TPY.requires_grad = True

        jac_ = jacobian(self.TYdot_jac, TPY, create_graph=False)

        return jac_

t0_start = perf_counter()
mech_yaml = '../../data/chem.yaml'

sol = rt.Solution(mech_yaml=mech_yaml, device=None,vectorize=True)

gas = ct.Solution(mech_yaml)


# Initial condition
P = ct.one_atm * 1
T = 1800

composition = 'IC8H18:0.5,O2:12.5,N2:34.0'

gas.TPX = T, P, composition
y0 = np.hstack((gas.T, gas.Y))

# Set up objects representing the ODE and the solver
ode = ReactorOde(gas)

# Integrate the equations using Cantera
t_end = 1e-3
states = ct.SolutionArray(gas, 1, extra={'t': [0.0]})
dt = 1e-5
t = 0
ode_success = True
y = y0
t1 = perf_counter()
print('before integration','time spent {:.1e} [s]'.format(t1 - t0_start))
while ode_success and t < t_end:
    odesol = solve_ivp(ode,
                       t_span=(t, t + dt),
                       y0=y,
                       method='BDF',
                       vectorized=False, jac=None)
    # print('nfev',odesol.nfev,'njev',odesol.njev,'nlu',odesol.nlu)
    t = odesol.t[-1]
    y = odesol.y[:, -1]
    ode_successful = odesol.success

    gas.TPY = odesol.y[0, -1], P, odesol.y[1:, -1]
    states.append(gas.state, t=t)


sol.gas.TPX = T, P, composition
sol.set_pressure(sol.gas.P)
ode_rt = ReactorOdeRT(sol=sol)

t_stop1 = perf_counter()
print('finish cantera integration')
print('time spent {:.1e} [s]'.format(t_stop1 - t1))


# Integrate the equations using ReacTorch
states_rt = ct.SolutionArray(sol.gas, 1, extra={'t': [0.0]})
t = 0
ode_success = True
y = y0
dt = 1e-5

# Diable AD for jacobian seems more effient for this case.
print('reacotrch')
while ode_success and t < t_end:
    
    odesol = solve_ivp(ode_rt,
                       t_span=(t, t + dt),
                       y0=y,
                       method='BDF',
                       vectorized=True, jac=None)

    t = odesol.t[-1]
    y = odesol.y[:, -1]
    ode_successful = odesol.success

    #print('t {} T {}'.format(t, y[0]))

    # print('nfev',odesol.nfev,'njev',odesol.njev,'nlu',odesol.nlu)
    sol.gas.TPY = odesol.y[0, -1], P, odesol.y[1:, -1]
    states_rt.append(sol.gas.state, t=t)
    
t_stop = perf_counter()
print('time spent {:.1e} [s]'.format(t_stop - t_stop1))
#Plot the results
try:
    import matplotlib.pyplot as plt
    L1 = plt.plot(states.t, states.T, ls='--',
                  color='r', label='T Cantera', lw=1)
    L1_rt = plt.plot(states_rt.t, states_rt.T, ls='-',
                      color='r', label='T ReacTorch', lw=1)
    plt.xlabel('Time (s)')
    plt.ylabel('Temperature (K)')

    plt.twinx()
    L2 = plt.plot(states.t, states('OH').Y, ls='--', label='OH Cantera', lw=1)
    L2_rt = plt.plot(states_rt.t, states_rt('OH').Y,
                      ls='-',
                      label='OH ReacTorch',
                      lw=1)
    plt.ylabel('Mass Fraction')

    plt.legend(L1+L2+L1_rt+L2_rt, [line.get_label()
                                    for line in L1+L2+L1_rt+L2_rt], loc='best')

    plt.savefig('cantera_reactorch_validation.png', dpi=300)
    plt.show()
except ImportError:
    print('Matplotlib not found. Unable to plot results.')

[Zeracesharon]

Run the above testing code, you will get an results like:

[Zeracesharon]

When the vectorize of the reactorch is off, see the belown code
sol = rt.Solution(mech_yaml=mech_yaml, device=None,vectorize=False) #before vectorize=True.
The mismatch is reduced, see the fig belown

[jiweiqi]

It looks like there is something wrong in the code, the temperature difference is too large, especially the equilibrium temperature is different. Perhaps there is some issue with the thermodynamic data?

In the worse case, there is something wrong in ReacTorch.

Which mechanism is used in the above code?

[Zeracesharon]

composition = 'IC8H18:0.5,O2:12.5,N2:34.0'
mechanism is the reduced IC8H18,which could be found here(the chem.xml file inside)
https://github.com/jiweiqi/CollectionOfMechanisms/tree/master/i-Octane_C8H18/llnl_version3/ReducedByTianfengLu/sk143

[Zeracesharon]

Bug dected on Solution_test----with tolerance delta=1e-4， kf (forward_rate_constants) kc ( equilibrium_constants ) kr (reverse_rate_constants) are correct. But net rates of progress & net rates of production are considered to be wrong.
Testing code on Solution_test:

# from multiprocessing import Pool
from time import perf_counter

import cantera as ct
import torch

import reactorch as rt
import numpy as np

cpu = torch.device('cpu')

cuda = torch.device('cuda:0')

device = cpu
###################修改输入文件
mech_yaml = '../data/IC8H18_reduced.yaml'
composition = 'IC8H18:0.5,O2:12.5,N2:34.0'


sol = rt.Solution(mech_yaml=mech_yaml, device=device)

#print(sol.species_names,sol.nasa_low[4,:],sol.nasa_high[4,:])


gas = sol.gas
gas.TPX = 1800, 5 * ct.one_atm, composition


# r = ct.IdealGasReactor(gas)
r=ct.IdealGasConstPressureReactor(gas)

sim = ct.ReactorNet([r])

# time = 0.0
# t_end=10
t_end = 1e-6
idt = 0
states = ct.SolutionArray(gas, extra=['t'])


T0 = gas.T

# print('%10s %10s %10s %14s' % ('t [s]', 'T [K]', 'P [atm]', 'u [J/kg]'))

while sim.time < t_end:

    time=sim.step()
  

    states.append(r.thermo.state, t=time)

    # print('zerace',time,sim.time)
    # if r.thermo.T > T0 + 600 and idt < 1e-10:
    #     idt = sim.time
    idt+=1
    # if idt > 1e-10 and sim.time > 4 * idt:
    #     break
    # # print('zerace',time)
# print('%10.3e %10.3f %10.3f %14.6e' % (sim.time,
#                                         r.T,
#                                         r.thermo.P / ct.one_atm,
#                                         r.thermo.u))

# print('idt = {:.2e} [s] number of points {}'.format(idt, states.t.shape[0]))
print('idt = {:.2e} [s] number of points {}'.format(idt, states.t.shape))
TP = torch.stack((torch.Tensor(states.T), torch.Tensor(states.P)), dim=-1)
Y = torch.Tensor(states.Y)
# print(Y,type(Y),Y.size(),TP,TP.size())
TPY = torch.cat([TP, Y], dim=-1).to(device)

t0_start = perf_counter()

sol.set_states(TPY)

t1_stop = perf_counter()


print('sol set_states time spent {:.1e} [s]'.format(t1_stop - t0_start))


reaction_equation = gas.reaction_equations()
species_name=gas.species_names

net_ct= states.net_rates_of_progress
netrates_ct=states.net_production_rates
kf = states.forward_rate_constants
kc = states.equilibrium_constants
kr = states.reverse_rate_constants
# print('zerace',np.shape(net_ct),np.shape(netrates_ct))

net_rt= sol.qdot.detach().cpu().numpy()
netrates_rt=sol.net_production_rates.detach().cpu().numpy()
kf_rt = sol.forward_rate_constants.detach().cpu().numpy()
kc_rt = sol.equilibrium_constants.detach().cpu().numpy()
kr_rt = sol.reverse_rate_constants.detach().cpu().numpy()
# print('zerace',gas.n_reactions,gas.n_species)
# print('zerace',np.shape(kf),np.shape(kc),np.shape(kr),np.shape(net_rt),np.shape(netrates_rt))


def check_rates(i):

    eps = 1e-300
    delta = 1e-4

    ratio = (kf[:, i] + eps) / (kf_rt[:, i] + eps)

    if ratio.min() < 1 - delta or ratio.max() > 1 + delta:
        print("forward constants {} {} {:.4e} {:.4e}".format(
            i, reaction_equation[i], ratio.min(), ratio.max()))

    ratio = (kc[:, i] + eps) / (kc_rt[:, i] + eps)

    if ratio.min() < 1 - delta or ratio.max() > 1 + delta:
        print("equilibrium constants {} {} {:.4e} {:.4e}".format(
            i, reaction_equation[i], ratio.min(), ratio.max()))

    ratio = (kr[:, i] + eps) / (kr_rt[:, i] + eps)

    if ratio.min() < 1 - delta or ratio.max() > 1 + delta:
        print("reverse constants {} {} {:.4e} {:.4e}".format(
            i, reaction_equation[i], ratio.min(), ratio.max()))
        
    ratio = (net_ct[:, i] + eps) / (net_rt[:, i] + eps)  

    if ratio.min() < 1 - delta or ratio.max() > 1 + delta:
        print("net rates of progress {} {} {:.4e} {:.4e}".format(
            i, reaction_equation[i], ratio.min(), ratio.max()))
    return i

def check_rates_production(j):
    eps = 1e-300
    delta = 1e-4
    ratio = (netrates_ct[:, j] + eps) / (netrates_rt[:, j] + eps)  

    if ratio.min() < 1 - delta or ratio.max() > 1 + delta:
        print("net rates of production {} {} {:.4e} {:.4e}".format(
            j, species_name[j], ratio.min(), ratio.max()))
    return j

for i in range(gas.n_reactions):
    check_rates(i)

for j in range(gas.n_species):
    check_rates_production(j)

# if __name__ == '__main__':
#     with Pool(4) as p:
#         print(p.map(check_rates, range(gas.n_reactions)))

t1_stop = perf_counter()
print('sol check_rates time spent {:.1e} [s]'.format(t1_stop - t0_start))

with the results:

[Zeracesharon]

All dection errors on net_rate progress are all with the same feature: the ratio_min would be something like -exp(+278), ratio_max are all normal, Pointing the bug in wdot function. More details need to be explored further.

[jiweiqi]

Nice catch. I guess you can figure out the difference among those reactions shows wrong rates and the rest.

It looks like the error happens when the reaction is irreversible.

[Zeracesharon]

Testing the code_exploring further:
According to the results, molar concentration is accurate with the tolerance 1e-15, the reactant order and product_stoich_coeffs is exactly correct. But the forward rates and reverse rates show big difference.
Another funding is that as the time goes by more reactions go out of control, when testing time is 1e-10, no mismatch is found, as time goes by 1e-8: 0 reactions 1e-6: 15 reactions detected, 1e-5: 186 reactions detected 1e-3: 849 reactions detected. -------------suggesting something time-dependent goes wrong when it advances.

#from multiprocessing import Pool
from time import perf_counter

import cantera as ct
import torch

import reactorch as rt
import numpy as np

cpu = torch.device('cpu')

cuda = torch.device('cuda:0')

device = cpu
###################修改输入文件
mech_yaml = '../data/IC8H18_reduced.yaml'
composition = 'IC8H18:0.5,O2:12.5,N2:34.0'


sol = rt.Solution(mech_yaml=mech_yaml, device=device)#here vectorize seems not change the results

#print(sol.species_names,sol.nasa_low[4,:],sol.nasa_high[4,:])


gas = sol.gas
gas.TPX = 1800, 5 * ct.one_atm, composition


# r = ct.IdealGasReactor(gas)
r=ct.IdealGasConstPressureReactor(gas)

sim = ct.ReactorNet([r])

# time = 0.0
# t_end=10
t_end = 1e-5
idt = 0
states = ct.SolutionArray(gas, extra=['t'])


T0 = gas.T

# print('%10s %10s %10s %14s' % ('t [s]', 'T [K]', 'P [atm]', 'u [J/kg]'))

while sim.time < t_end:

    time=sim.step()
  

    states.append(r.thermo.state, t=time)

    # print('zerace',time,sim.time)
    # if r.thermo.T > T0 + 600 and idt < 1e-10:
    #     idt = sim.time
    idt+=1
    # if idt > 1e-10 and sim.time > 4 * idt:
    #     break
    # # print('zerace',time)
# print('%10.3e %10.3f %10.3f %14.6e' % (sim.time,
#                                         r.T,
#                                         r.thermo.P / ct.one_atm,
#                                         r.thermo.u))

# print('idt = {:.2e} [s] number of points {}'.format(idt, states.t.shape[0]))
print('idt = {:.2e} [s] number of points {}'.format(idt, states.t.shape))
TP = torch.stack((torch.Tensor(states.T), torch.Tensor(states.P)), dim=-1)
Y = torch.Tensor(states.Y)
# print(Y,type(Y),Y.size(),TP,TP.size())
TPY = torch.cat([TP, Y], dim=-1).to(device)

t0_start = perf_counter()

sol.set_states(TPY)

t1_stop = perf_counter()


print('sol set_states time spent {:.1e} [s]'.format(t1_stop - t0_start))


reaction_equation = gas.reaction_equations()
species_name=gas.species_names
frates_ct=states.forward_rates_of_progress
rrates_ct=states.reverse_rates_of_progress
net_ct= states.net_rates_of_progress
netrates_ct=states.net_production_rates
kf = states.forward_rate_constants
kc = states.equilibrium_constants
kr = states.reverse_rate_constants
# print('zerace',np.shape(net_ct),np.shape(netrates_ct))

net_rt= sol.qdot.detach().cpu().numpy()
frates_rt=sol.forward_rates_of_progress.detach().cpu().numpy()
rrates_rt=sol.reverse_rates_of_progress.detach().cpu().numpy()
netrates_rt=sol.net_production_rates.detach().cpu().numpy()
kf_rt = sol.forward_rate_constants.detach().cpu().numpy()
kc_rt = sol.equilibrium_constants.detach().cpu().numpy()
kr_rt = sol.reverse_rate_constants.detach().cpu().numpy()
# print(sol.C.size(),sol.reactant_orders.size())
# print('zerace',gas.n_reactions,gas.n_species)
# print('zerace',np.shape(kf),np.shape(kc),np.shape(kr),np.shape(net_rt),np.shape(netrates_rt))
zerace1=sol.reactant_orders.detach().cpu().numpy()-states.reactant_stoich_coeffs()
zerace2=sol.product_stoich_coeffs.detach().cpu().numpy()-states.product_stoich_coeffs()
print(np.shape(zerace1),np.nonzero(zerace1),np.nonzero(zerace2))
zerace3=sol.C.detach().cpu().numpy()-states.concentrations
print(np.shape(sol.C),np.shape(states.concentrations))
bu=np.where(abs(zerace3)>1e-10)
print(bu)
def check_rates(i):

    eps = 1e-300
    delta = 1e-4
    global counter_abnormal_reactions

    ratio = (kf[:, i] + eps) / (kf_rt[:, i] + eps)

    if ratio.min() < 1 - delta or ratio.max() > 1 + delta:
        print("forward constants {} {} {:.4e} {:.4e}".format(
            i, reaction_equation[i], ratio.min(), ratio.max()))

    ratio = (kc[:, i] + eps) / (kc_rt[:, i] + eps)

    if ratio.min() < 1 - delta or ratio.max() > 1 + delta:
        print("equilibrium constants {} {} {:.4e} {:.4e}".format(
            i, reaction_equation[i], ratio.min(), ratio.max()))

    ratio = (kr[:, i] + eps) / (kr_rt[:, i] + eps)

    if ratio.min() < 1 - delta or ratio.max() > 1 + delta:
        print("reverse constants {} {} {:.4e} {:.4e}".format(
            i, reaction_equation[i], ratio.min(), ratio.max()))
        
    ratio = (net_ct[:, i] + eps) / (net_rt[:, i] + eps)  

    if ratio.min() < 1 - delta or ratio.max() > 1 + delta:
        counter_abnormal_reactions+=1
    # if ratio.min() < 1 - delta and ratio.max() < 1:#for checking ratio_max
        # print("net rates of progress {} {:.4e} {:.4e}".format(
        #     i, ratio.min(), ratio.max()))
        
    ratio = (frates_ct[:, i] + eps) / (frates_rt[:, i] + eps)  

    if ratio.min() < 1 - delta or ratio.max() > 1 + delta:
       # pass

        print("forward rates of progress {} {:.4e} {:.4e}".format(
            i, ratio.min(), ratio.max()))
        
    ratio = (rrates_ct[:, i] + eps) / (rrates_rt[:, i] + eps)  

    if ratio.min() < 1 - delta or ratio.max() > 1 + delta:
       # pass
    # if ratio.min() < 1 - delta and ratio.max() < 1:#for checking ratio_max
        print("reverse rates of progress {} {:.4e} {:.4e}".format(
            i, ratio.min(), ratio.max()))
        
        
       
    return i

def check_rates_production(j):
    eps = 1e-300
    delta = 1e-4
    ratio = (netrates_ct[:, j] + eps) / (netrates_rt[:, j] + eps)  

    if ratio.min() < 1 - delta or ratio.max() > 1 + delta:
        pass
        # print("net rates of production {} {} {:.4e} {:.4e}".format(
        #     j, species_name[j], ratio.min(), ratio.max()))
    return j
counter_abnormal_reactions=0
for i in range(gas.n_reactions):
    check_rates(i)

for j in range(gas.n_species):
    check_rates_production(j)

# if __name__ == '__main__':
#     with Pool(4) as p:
#         print(p.map(check_rates, range(gas.n_reactions)))

t1_stop = perf_counter()
print('sol check_rates time spent {:.1e} [s]'.format(t1_stop - t0_start))
print('zerace',counter_abnormal_reactions)

Results shown belown

[Zeracesharon]

Nice catch. I guess you can figure out the difference among those reactions shows wrong rates and the rest.

It looks like the error happens when the reaction is irreversible.

Both reversible and irreversible reactions are detected. At the first begining, only something irreversbile. as the time goes by(maybe because the molar concentration and the temperature is already different from the last time states, the miss rise up to involve in more reactions)

[jiweiqi]

Can you check what the reverse rate constants for the irreversible reaction are? For both cantera and ReacTorch

[Zeracesharon]

Can you check what the reverse rate constants for the irreversible reaction are? For both cantera and ReacTorch

I have check it before in the previous coding and no abnormals are found. The constants including equilibrium &reverse constants as well as forward rate constants, they are safe. I have check the reactant orders as well as the product stoi_coefficiency, they are safe as well. After some endevors been made, I figure out some bizzar things. It is really strange but I have tried many times, it did work like what I was guessing.

[Zeracesharon]

Bizzar things occur. when the initial settings with the mechanism IC8H18_reduced.yaml are:
composition = 'IC8H18:0.5,O2:12.5,N2:34.0'
TP=1800, 5 * ct.one_atm
t_end = 1e-5
Thgrough all the camparisons, negtive of mass fraction in cantera calculations are found. See the code lying on Column 114-118 :

C_rt=sol.C.detach().cpu().numpy()
C_ct=states.concentrations
Y_rt=sol.Y.detach().cpu().numpy()
Y_ct=states.Y
print("negtive happens",np.where(Y_ct<0))

Then many negtive values of Y are detected. While for reactorch, values of Y are limited to be >=0. Relating code in solution-set_states(TPY)

       if TPY.shape[1] == self.n_species + 2:
            self.P = TPY[:, 1:2] 
            self.Y= torch.clamp(TPY[:, 2:], min=0, max=None)
        if TPY.shape[1] == self.n_species + 1:      
            self.Y= torch.clamp(TPY[:, 1:], min=0.0, max=None)

Thus, huge differences are detected for the calculation of species concentration and certainly on forward rate progress. (The normalization of Y has been commented:
#self.Y = (self.Y.T / self.Y.sum(dim=1)).T

Such differences between cantera and reactorch disappear when the code modified as:

        if TPY.shape[1] == self.n_species + 2:
            self.P = TPY[:, 1:2]
            self.Y= TPY[:, 2:]       

        if TPY.shape[1] == self.n_species + 1:
            self.P = torch.ones_like(self.T) * self.P_ref
             self.Y= TPY[:, 1:]

The testing code on Solution_test:

from time import perf_counter

import cantera as ct
import torch

import reactorch as rt
import numpy as np

cpu = torch.device('cpu')

cuda = torch.device('cuda:0')

device = cpu
###################
mech_yaml = '../data/IC8H18_reduced.yaml'
# mech_yaml = '../data/gri30.yaml'
composition = 'IC8H18:0.5,O2:12.5,N2:34.0'
# composition='CH4:0.5,O2:11,N2:40.0'


sol = rt.Solution(mech_yaml=mech_yaml, device=device)#here vectorize seems not change the results

#print(sol.species_names,sol.nasa_low[4,:],sol.nasa_high[4,:])


gas = sol.gas
gas.TPX = 1800, 5 * ct.one_atm, composition


# r = ct.IdealGasReactor(gas)
r=ct.IdealGasConstPressureReactor(gas)

sim = ct.ReactorNet([r])

# time = 0.0
# t_end=10
t_end = 1e-3
idt = 0
states = ct.SolutionArray(gas, extra=['t'])


T0 = gas.T

# print('%10s %10s %10s %14s' % ('t [s]', 'T [K]', 'P [atm]', 'u [J/kg]'))

while sim.time < t_end:

    time=sim.step()
  

    states.append(r.thermo.state, t=time)
           
    # print('zerace',time,sim.time)
    # if r.thermo.T > T0 + 600 and idt < 1e-10:
    #     idt = sim.time
    idt+=1
    # if idt > 1e-10 and sim.time > 4 * idt:
    #     break
    # # print('zerace',time)
# print('%10.3e %10.3f %10.3f %14.6e' % (sim.time,
#                                         r.T,
#                                         r.thermo.P / ct.one_atm,
#                                         r.thermo.u))

# print('idt = {:.2e} [s] number of points {}'.format(idt, states.t.shape[0]))
print('idt = {:.2e} [s] number of points {}'.format(idt, states.t.shape))
# TP = torch.stack((torch.Tensor(states.T), torch.Tensor(states.P)), dim=-1)
T=torch.Tensor(states.T)
P=torch.Tensor(states.P)
TP = torch.stack((T,P), dim=-1)

Y = torch.Tensor(states.Y)
# print(Y,type(Y),Y.size(),TP,TP.size())
TPY = torch.cat([TP, Y], dim=-1).to(device)
# print('zerace',"T",T.size(),"P",P.size(),"TP",TP.size(),'Y',Y.size(),'TPY',TPY.size())
# print('zerace',TPY.shape[1])
t0_start = perf_counter()

sol.set_states(TPY)

t1_stop = perf_counter()


print('sol set_states time spent {:.1e} [s]'.format(t1_stop - t0_start))


reaction_equation = gas.reaction_equations()
species_name=gas.species_names
frates_ct=states.forward_rates_of_progress
rrates_ct=states.reverse_rates_of_progress
net_ct= states.net_rates_of_progress
netrates_ct=states.net_production_rates
kf = states.forward_rate_constants
kc = states.equilibrium_constants
kr = states.reverse_rate_constants
# print('zerace',np.shape(net_ct),np.shape(netrates_ct))

net_rt= sol.qdot.detach().cpu().numpy()
frates_rt=sol.forward_rates_of_progress.detach().cpu().numpy()
rrates_rt=sol.reverse_rates_of_progress.detach().cpu().numpy()
netrates_rt=sol.net_production_rates.detach().cpu().numpy()
kf_rt = sol.forward_rate_constants.detach().cpu().numpy()
kc_rt = sol.equilibrium_constants.detach().cpu().numpy()
kr_rt = sol.reverse_rate_constants.detach().cpu().numpy()
# print(sol.Y.size(),sol.reactant_orders.size())
# print('zerace',gas.n_reactions,gas.n_species)
# print('zerace',np.shape(kf),np.shape(kc),np.shape(kr),np.shape(net_rt),np.shape(netrates_rt))
# zerace1=sol.reactant_orders.detach().cpu().numpy()
# zerace2=states.reactant_stoich_coeffs()
# zerace3=sol.product_stoich_coeffs.detach().cpu().numpy()
# zerace4=states.product_stoich_coeffs()

C_rt=sol.C.detach().cpu().numpy()
C_ct=states.concentrations
Y_rt=sol.Y.detach().cpu().numpy()
Y_ct=states.Y
print("negtive happens",np.where(Y_ct<0))
density_rt=sol.density_mass.detach().cpu().numpy()
density_ct=states.density_mass
# print(np.shape(sol.C),np.shape(states.concentrations))
# bu=np.where(abs(zerace3)>1e-10)
# print(bu)


def check_rates(i):

    eps = 1e-300
    delta = 1e-4
    global counter_abnormal_reactions

    ratio = (kf[:, i] + eps) / (kf_rt[:, i] + eps)

    if ratio.min() < 1 - delta or ratio.max() > 1 + delta:
        print("forward constants {} {} {:.4e} {:.4e}".format(
            i, reaction_equation[i], ratio.min(), ratio.max()))

    ratio = (kc[:, i] + eps) / (kc_rt[:, i] + eps)

    if ratio.min() < 1 - delta or ratio.max() > 1 + delta:
        print("equilibrium constants {} {} {:.4e} {:.4e}".format(
            i, reaction_equation[i], ratio.min(), ratio.max()))

    ratio = (kr[:, i] + eps) / (kr_rt[:, i] + eps)

    if ratio.min() < 1 - delta or ratio.max() > 1 + delta:
        print("reverse constants {} {} {:.4e} {:.4e}".format(
            i, reaction_equation[i], ratio.min(), ratio.max()))
        
    ratio = (net_ct[:, i] + eps) / (net_rt[:, i] + eps)  

    if ratio.min() < 1 - delta or ratio.max() > 1 + delta:
        counter_abnormal_reactions+=1
  
        print("net rates of progress {} {:.4e} {:.4e}".format(
            i, ratio.min(), ratio.max()))
        
    ratio = (frates_ct[:, i] + eps) / (frates_rt[:, i] + eps)  

    if ratio.min() < 1 - delta or ratio.max() > 1 + delta:
        # pass

        print("forward rates of progress {}{} {:.4e} {:.4e}".format(
            i, reaction_equation[i],ratio.min(), ratio.max()))
        
    ratio = (rrates_ct[:, i] + eps) / (rrates_rt[:, i] + eps)  

    if ratio.min() < 1 - delta or ratio.max() > 1 + delta:
        # pass
   
        print("reverse rates of progress{} {} {:.4e} {:.4e}".format(
            i, reaction_equation[i],ratio.min(), ratio.max()))
        
        
       
    return i

def check_rates_production(j):
    eps = 1e-300
    delta = 1e-4
    ratio = (netrates_ct[:, j] + eps) / (netrates_rt[:, j] + eps)  

    if ratio.min() < 1 - delta or ratio.max() > 1 + delta:
        # pass
        print("net rates of production {} {} {:.4e} {:.4e}".format(
            j, species_name[j], ratio.min(), ratio.max()))
    
    ratio = (C_ct[:, j] + eps) / (C_rt[:, j] + eps) 
    if ratio.min() < 1 - delta or ratio.max() > 1 + delta:
        # pass
        print("concentration {} {} {:.4e} {:.4e}".format(
            j, species_name[j], ratio.min(), ratio.max()))
        
    ratio = (Y_ct[:, j] + eps) / (Y_rt[:, j] + eps) 
    
    if ratio.min() < 1 - delta or ratio.max() > 1 + delta:
        # pass
        print("mass fraction {} {} {:.4e} {:.4e}".format(
            j, species_name[j], ratio.min(), ratio.max()))
        
    return j


    
    
    
counter_abnormal_reactions=0
for i in range(gas.n_reactions):
    check_rates(i)

for j in range(gas.n_species):
    check_rates_production(j)   
    
eps = 1e-300
delta = 1e-4
ratio = (density_ct[:] + eps) / (density_rt[:] + eps) 
if ratio.min() < 1 - delta or ratio.max() > 1 + delta:
    # pass
    print("mass density {:.4e} {:.4e}".format(
             ratio.min(), ratio.max()))


t1_stop = perf_counter()
print('sol check_rates time spent {:.1e} [s]'.format(t1_stop - t0_start))
print('abnormal reactions detected',counter_abnormal_reactions)

Before modified results are:

After the modification:

[jiweiqi]

For irreversible reaction, the rate constants should be zero. What's the value you observed?

[jiweiqi]

I think Cantera may also clips the concentration to be within zero and one.

[Zeracesharon]

Something needs to be explored further:
The formular used in reactorch for the calculation is

But the testing shows that the density mass is different from what cantera has been calculated.

[Zeracesharon]

I think Cantera may also clips the concentration to be within zero and one.

The mass fraction but not the concentration. Give a look at the code, you will find at this setting,

Y_ct=states.Y
print("negtive happens",np.where(Y_ct<0))

You will see many negtive numbers have been detected

[jiweiqi]

I guess allowing negative mass fraction has a lot of other consequences, for example, computing density

[Zeracesharon]

I guess allowing negative mass fraction has a lot of other consequences, for example, computing density

Yes, so when I modified the code in reactorch, the mismatch of reactorch and cantera in terms of concentration and mass fraction as well as abnormal reactions all disappear.

[Zeracesharon]

For irreversible reaction, the rate constants should be zero. What's the value you observed?

I only check the difference between cantera and reactorch. No special attentions on irreversible constants are paid yet. But I would try.

[Zeracesharon]

In terms of the mass density mismatch, I am thinking and trying to evaluate the concentration and mass density according to the following formula, to see if there is any difference.

[jiweiqi]

I agree with you the problem is probably mainly related to the tiny negative mass fractions

[Zeracesharon]

For irreversible reaction, the rate constants should be zero. What's the value you observed?

I only check the difference between cantera and reactorch. No special attentions on irreversible constants are paid yet. But I would try.

Testing code on Solution_test:
For both those without clips or those with clips in Y, no nonzero reverse constants for irrevisible has been detected which means the reverse constants are safe for both cases.

reverse_reaction=sol.is_reversible
for i in range(gas.n_reactions):
    check_rates(i)
    if reverse_reaction[i]==0:
        # krr=kr_rt[:,i]
        # krc=kr[:,i]
        krr=np.nonzero(kr_rt[:,i])
        krc=np.nonzero(kr[:,i])
        if np.nonzero(krr)[0].size !=0:
            print('nonzero irreverse constants detected in reactorch')
        if np.nonzero(krc)[0].size !=0:
            print('nonzero irreverse constants detected in reactorch')

[Zeracesharon]

I agree with you the problem is probably mainly related to the tiny negative mass fractions

I would continue to commit my time figuring out the mismatch in mass density for Solution_test since they shows obvious (nonnegligible) difference for cantera and reactorch. Would this be Okay? If not, what are your suggestions.

[jiweiqi]

Sure. But I think you had better keep the normalization of the mass fraction, eventually. Although you can relax it for debugging code.

[Zeracesharon]

Testings on density mass and mean malecular weights are passed. No further modifications are needed. Density mass calculation is safe. Besides, many other calculations are evluated in the code. They all perform well in comparision with the counterparts in cantera including Tdot Ydot.
Testing variables are :
forward_rates_of_progress
reverse_rates_of_progress
net_rates_of_progress
net_production_rates
forward_rate_constants
equilibrium_constants
reverse_rate_constants
product_stoich_coeffs()
reactant_orders
concentrations
density_mass
mean_molecular_weight
TYdot
Testing code

# from multiprocessing import Pool
from time import perf_counter

import cantera as ct
import torch

import reactorch as rt
import numpy as np

cpu = torch.device('cpu')

cuda = torch.device('cuda:0')

device = cpu
###################修改输入文件
mech_yaml = '../data/IC8H18_reduced.yaml'
composition = 'IC8H18:0.5,O2:12.5,N2:34.0'
# composition='CH4:0.5,O2:11,N2:40.0'
# composition = 'IC8H18:0.08,O2:1.0,N2:3.76'


sol = rt.Solution(mech_yaml=mech_yaml, device=device)#here vectorize seems not change the results

#print(sol.species_names,sol.nasa_low[4,:],sol.nasa_high[4,:])


gas = sol.gas
gas.TPX = 1800, 5 * ct.one_atm, composition


# r = ct.IdealGasReactor(gas)
r=ct.IdealGasConstPressureReactor(gas)

sim = ct.ReactorNet([r])

# time = 0.0
# t_end=10
t_end = 1e-3
idt = 0
states = ct.SolutionArray(gas, extra=['t'])


T0 = gas.T

# print('%10s %10s %10s %14s' % ('t [s]', 'T [K]', 'P [atm]', 'u [J/kg]'))

while sim.time < t_end:

    time=sim.step()
  

    states.append(r.thermo.state, t=time)   
           

    # if r.thermo.T > T0 + 600 and idt < 1e-10:
    #     idt = sim.time
    idt+=1
    # if idt > 1e-10 and sim.time > 4 * idt:
    #     break
    # # print('zerace',time)
# print('%10.3e %10.3f %10.3f %14.6e' % (sim.time,
#                                         r.T,
#                                         r.thermo.P / ct.one_atm,
#                                         r.thermo.u))

# print('idt = {:.2e} [s] number of points {}'.format(idt, states.t.shape[0]))
print('idt = {:.2e} [s] number of points {}'.format(idt, states.t.shape))
# TP = torch.stack((torch.Tensor(states.T), torch.Tensor(states.P)), dim=-1)
T=torch.Tensor(states.T)
P=torch.Tensor(states.P)
TP = torch.stack((T,P), dim=-1)

Y = torch.Tensor(states.Y)
# print(Y,type(Y),Y.size(),TP,TP.size())
TPY = torch.cat([TP, Y], dim=-1).to(device)
# print('zerace',"T",T.size(),"P",P.size(),"TP",TP.size(),'Y',Y.size(),'TPY',TPY.size())
# print('zerace',TPY.shape[1])
t0_start = perf_counter()

sol.set_states(TPY)

t1_stop = perf_counter()


print('sol set_states time spent {:.1e} [s]'.format(t1_stop - t0_start))


reaction_equation = gas.reaction_equations()
species_name=gas.species_names
frates_ct=states.forward_rates_of_progress
rrates_ct=states.reverse_rates_of_progress
net_ct= states.net_rates_of_progress
netrates_ct=states.net_production_rates
kf = states.forward_rate_constants
kc = states.equilibrium_constants
kr = states.reverse_rate_constants
# print('zerace',np.shape(net_ct),np.shape(netrates_ct))

net_rt= sol.qdot.detach().cpu().numpy()
frates_rt=sol.forward_rates_of_progress.detach().cpu().numpy()
rrates_rt=sol.reverse_rates_of_progress.detach().cpu().numpy()
netrates_rt=sol.net_production_rates.detach().cpu().numpy()
kf_rt = sol.forward_rate_constants.detach().cpu().numpy()
kc_rt = sol.equilibrium_constants.detach().cpu().numpy()
kr_rt = sol.reverse_rate_constants.detach().cpu().numpy()
# print(sol.Y.size(),sol.reactant_orders.size())
# print('zerace',gas.n_reactions,gas.n_species)
# print('zerace',np.shape(kf),np.shape(kc),np.shape(kr),np.shape(net_rt),np.shape(netrates_rt))
# zerace1=sol.reactant_orders.detach().cpu().numpy()
# zerace2=states.reactant_stoich_coeffs()
# zerace3=sol.product_stoich_coeffs.detach().cpu().numpy()
# zerace4=states.product_stoich_coeffs()

C_rt=sol.C.detach().cpu().numpy()
C_ct=states.concentrations
Y_rt=sol.Y.detach().cpu().numpy()
Y_ct=states.Y
# print("negtive happens",np.where(Y_ct<0))
# print("negtive happens",np.where(C_ct<0))
density_rt=sol.density_mass.detach().cpu().numpy()
density_ct=states.density_mass
reverse_reaction=sol.is_reversible
# print('size match',np.shape(density_ct),np.shape(density_rt))

Ydot_rt=sol.Ydot.detach().cpu().numpy()
Tdot_rt=sol.Tdot.detach().cpu().numpy()
mean_W_rt=sol.mean_molecular_weight.detach().cpu().numpy()
mean_W_ct=states.mean_molecular_weight
# print('size match',np.shape(mean_W_ct),np.shape(mean_W_rt))
#################################cantera TYdot calculation
rho = states.DP[:][0]
wdot = netrates_ct
temp_entro=np.dot(states.partial_molar_enthalpies, wdot.T)
Tdot_ct = - ( np.diag(temp_entro)/
                  (rho * states.cp_mass))
temp_ydot=wdot * states.molecular_weights
Ydot_ct=np.zeros_like(temp_ydot)
for i in range(idt):

    Ydot_ct[i,:] = temp_ydot[i,:] / rho[i]





def check_rates(i):

    eps = 1e-300
    delta = 1e-4
    global counter_abnormal_reactions

    ratio = (kf[:, i] + eps) / (kf_rt[:, i] + eps)

    if ratio.min() < 1 - delta or ratio.max() > 1 + delta:
        # pass
        print("forward constants {} {} {:.4e} {:.4e}".format(
            i, reaction_equation[i], ratio.min(), ratio.max()))

    ratio = (kc[:, i] + eps) / (kc_rt[:, i] + eps)

    if ratio.min() < 1 - delta or ratio.max() > 1 + delta:
        # pass
        print("equilibrium constants {} {} {:.4e} {:.4e}".format(
            i, reaction_equation[i], ratio.min(), ratio.max()))

    ratio = (kr[:, i] + eps) / (kr_rt[:, i] + eps)

    if ratio.min() < 1 - delta or ratio.max() > 1 + delta:
        # pass
        print("reverse constants {} {} {:.4e} {:.4e}".format(
            i, reaction_equation[i], ratio.min(), ratio.max()))
        
    ratio = (net_ct[:, i] + eps) / (net_rt[:, i] + eps)  

    if ratio.min() < 1 - delta or ratio.max() > 1 + delta:
        counter_abnormal_reactions+=1
  
        print("net rates of progress {} {:.4e} {:.4e}".format(
            i, ratio.min(), ratio.max()))
        
    ratio = (frates_ct[:, i] + eps) / (frates_rt[:, i] + eps)  

    if ratio.min() < 1 - delta or ratio.max() > 1 + delta:
        # pass

        print("forward rates of progress {}{} {:.4e} {:.4e}".format(
            i, reaction_equation[i],ratio.min(), ratio.max()))
        
    ratio = (rrates_ct[:, i] + eps) / (rrates_rt[:, i] + eps)  

    if ratio.min() < 1 - delta or ratio.max() > 1 + delta:
        # pass
   
        print("reverse rates of progress{} {} {:.4e} {:.4e}".format(
            i, reaction_equation[i],ratio.min(), ratio.max()))           
    return i

def check_rates_production(j):
    eps = 1e-300
    delta = 1e-5
    ratio = (netrates_ct[:, j] + eps) / (netrates_rt[:, j] + eps)  

    if ratio.min() < 1 - delta or ratio.max() > 1 + delta:
        # pass
        print("net rates of production {} {} {:.4e} {:.4e}".format(
            j, species_name[j], ratio.min(), ratio.max()))
    
    ratio = (C_ct[:, j] + eps) / (C_rt[:, j] + eps) 
    if ratio.min() < 1 - delta or ratio.max() > 1 + delta:
        # pass
        print("concentration {} {} {:.4e} {:.4e}".format(
            j, species_name[j], ratio.min(), ratio.max()))
        
    ratio = (Y_ct[:, j] + eps) / (Y_rt[:, j] + eps) 
    
    if ratio.min() < 1 - delta or ratio.max() > 1 + delta:
        # pass
        print("mass fraction {} {} {:.4e} {:.4e}".format(
            j, species_name[j], ratio.min(), ratio.max()))
        
    ratio = (Ydot_ct[:, j] + eps) / (Ydot_rt[:, j] + eps) 
    
    if ratio.min() < 1 - delta or ratio.max() > 1 + delta:
        # pass
        print("Ydot {} {} {:.4e} {:.4e}".format(
            j, species_name[j], ratio.min(), ratio.max()))
        
    return j


    
    
    
counter_abnormal_reactions=0
for i in range(gas.n_reactions):
    check_rates(i)
    if reverse_reaction[i]==0:
        # krr=kr_rt[:,i]
        # krc=kr[:,i]
        krr=np.nonzero(kr_rt[:,i])
        krc=np.nonzero(kr[:,i])
        if np.nonzero(krr)[0].size !=0:
            print('nonzero irreverse constants detected in reactorch')
        if np.nonzero(krc)[0].size !=0:
            print('nonzero irreverse constants detected in reactorch')
        # print('nonzero irreverse constants_reactorch',np.nonzero(krr)[0].size)
        # print('nonzero irreverse constants_cantera',np.nonzero(krc)[0].size)

for j in range(gas.n_species):
    check_rates_production(j)   
    
eps = 1e-300
delta = 1e-5
ratio_density = (density_ct + eps) / (density_rt[:,0] + eps) 
if ratio_density.min() < 1 - delta or ratio_density.max() > 1 + delta:
    # pass
    print("mass density {:.4e} {:.4e}".format(
             ratio_density.min(), ratio_density.max()))
# print('min index value',np.argmin(ratio),ratio[(np.argmin(ratio))])
# print('max index value',np.argmax(ratio),ratio[(np.argmax(ratio))]) 
ratio = (Tdot_ct + eps) / (Tdot_rt + eps) 
if ratio.min() < 1 - delta or ratio.max() > 1 + delta:
    # pass
    print("Tdot {:.4e} {:.4e}".format(
             ratio.min(), ratio.max()))
    
ratio_mean = (mean_W_ct + eps) / (mean_W_rt[:,0] + eps) 
if ratio_mean.min() < 1 - delta or ratio_mean.max() > 1 + delta:
    # pass
    print("mean W {:.4e} {:.4e}".format(
             ratio_mean.min(), ratio_mean.max()))


t1_stop = perf_counter()
print('sol check_rates time spent {:.1e} [s]'.format(t1_stop - t0_start))
print('abnormal reactions detected',counter_abnormal_reactions)

Results show that all testing variables pass the accuracy checking.

[Zeracesharon]

While new issue arises:
When the code in set_TPY has been modified as:

        if TPY.shape[1] == self.n_species + 2:
            self.P = TPY[:, 1:2]          
            self.Y = TPY[:, 2:]
        if TPY.shape[1] == self.n_species + 1:
            self.P = torch.ones_like(self.T) * self.P_ref
            self.Y = TPY[:, 1:]

Runing the code in exaples/constant_pressure_autoignition: bug comes again, reactorch cannot calculate but cantera succeeded. Detailed information is still exploring.
Testing code in auto_ignition:

"""
Solve a constant pressure ignition problem where the governing equations are
implemented in Python.

This demonstrates an approach for solving problems where Cantera's reactor
network model cannot be configured to describe the system in question. Here,
Cantera is used for evaluating thermodynamic properties and kinetic rates while
an external ODE solver is used to integrate the resulting equations. In this
case, the SciPy wrapper for VODE is used, which uses the same variable-order BDF
methods as the Sundials CVODES solver used by Cantera.
"""

# TODO: the reactorch class seems to be very slow here, will figure out later

import cantera as ct
import numpy as np
import reactorch as rt
from scipy.integrate import solve_ivp
import torch
from torch.autograd.functional import jacobian as jacobian
from time import perf_counter




class ReactorOde(object):
    def __init__(self, gas):
        # Parameters of the ODE system and auxiliary data are stored in the
        # ReactorOde object.
        self.gas = gas
        self.P = gas.P

    def __call__(self, t, y):
        """the ODE function, y' = f(t,y) """

        # State vector is [T, Y_1, Y_2, ... Y_K]
        self.gas.set_unnormalized_mass_fractions(y[1:])
        self.gas.TP = y[0], self.P
        rho = self.gas.density

        wdot = self.gas.net_production_rates
        dTdt = - (np.dot(self.gas.partial_molar_enthalpies, wdot) /
                  (rho * self.gas.cp))
        dYdt = wdot * self.gas.molecular_weights / rho
        #print('cantera',np.hstack((dTdt, dYdt)))

        return np.hstack((dTdt, dYdt))


class ReactorOdeRT(object):
    def __init__(self, sol):
        # Parameters of the ODE system and auxiliary data are stored in the
        # ReactorOde object.
        self.sol = sol
        self.gas = sol.gas

    def __call__(self, t, y):
        """the ODE function, y' = f(t,y) """

        TPY = torch.Tensor(y).T

        with torch.no_grad():

            self.sol.set_states(TPY)

            TYdot = self.sol.TYdot_func()
        #print('reactorch',TYdot.T)
        
        return TYdot.T.detach().cpu().numpy()

    def TYdot_jac(self, TPY):

        TPY = torch.Tensor(TPY).unsqueeze(0)

        sol.set_states(TPY)

        return sol.TYdot_func().squeeze(0)

    def jac(self, t, y):

        TPY = torch.Tensor(y)

        TPY.requires_grad = True

        jac_ = jacobian(self.TYdot_jac, TPY, create_graph=False)

        return jac_

t0_start = perf_counter()
mech_yaml = '../../data/chem.yaml'

sol = rt.Solution(mech_yaml=mech_yaml, device=None,vectorize=False)

gas = ct.Solution(mech_yaml)


# Initial condition
P = ct.one_atm * 5
T = 1800

composition = 'IC8H18:0.5,O2:12.5,N2:34.0'

gas.TPX = T, P, composition
y0 = np.hstack((gas.T, gas.Y))

# Set up objects representing the ODE and the solver
ode = ReactorOde(gas)

# Integrate the equations using Cantera
t_end = 1e-3
states = ct.SolutionArray(gas, 1, extra={'t': [0.0]})
dt = 1e-5
t = 0
ode_success = True
y = y0
t1 = perf_counter()
print('before integration','time spent {:.1e} [s]'.format(t1 - t0_start))
while ode_success and t < t_end:
    odesol = solve_ivp(ode,
                       t_span=(t, t + dt),
                       y0=y,
                       method='BDF',
                       vectorized=False, jac=None)
    # print('nfev',odesol.nfev,'njev',odesol.njev,'nlu',odesol.nlu)
    t = odesol.t[-1]
    y = odesol.y[:, -1]
    ode_successful = odesol.success

    gas.TPY = odesol.y[0, -1], P, odesol.y[1:, -1]
    states.append(gas.state, t=t)


sol.gas.TPX = T, P, composition
sol.set_pressure(sol.gas.P)
ode_rt = ReactorOdeRT(sol=sol)

t_stop1 = perf_counter()
print('finish cantera integration')
print('time spent {:.1e} [s]'.format(t_stop1 - t1))


# Integrate the equations using ReacTorch
states_rt = ct.SolutionArray(sol.gas, 1, extra={'t': [0.0]})
t = 0
ode_success = True
y = y0
dt = 1e-5

# Diable AD for jacobian seems more effient for this case.
print('reacotrch')
while ode_success and t < t_end:
    
    odesol = solve_ivp(ode_rt,
                       t_span=(t, t + dt),
                       y0=y,
                       method='BDF',
                       vectorized=True, jac=None)

    t = odesol.t[-1]
    y = odesol.y[:, -1]
    ode_successful = odesol.success

    #print('t {} T {}'.format(t, y[0]))

    # print('nfev',odesol.nfev,'njev',odesol.njev,'nlu',odesol.nlu)
    sol.gas.TPY = odesol.y[0, -1], P, odesol.y[1:, -1]
    states_rt.append(sol.gas.state, t=t)
    
t_stop = perf_counter()
print('time spent {:.1e} [s]'.format(t_stop - t_stop1))
#Plot the results
try:
    import matplotlib.pyplot as plt
    L1 = plt.plot(states.t, states.T, ls='--',
                  color='r', label='T Cantera', lw=1)
    L1_rt = plt.plot(states_rt.t, states_rt.T, ls='-',
                      color='r', label='T ReacTorch', lw=1)
    plt.xlabel('Time (s)')
    plt.ylabel('Temperature (K)')

    plt.twinx()
    L2 = plt.plot(states.t, states('OH').Y, ls='--', label='OH Cantera', lw=1)
    L2_rt = plt.plot(states_rt.t, states_rt('OH').Y,
                      ls='-',
                      label='OH ReacTorch',
                      lw=1)
    plt.ylabel('Mass Fraction')

    plt.legend(L1+L2+L1_rt+L2_rt, [line.get_label()
                                    for line in L1+L2+L1_rt+L2_rt], loc='best')

    plt.savefig('cantera_reactorch_validation.png', dpi=300)
    plt.show()
except ImportError:
    print('Matplotlib not found. Unable to plot results.')

[Zeracesharon]

Abnormal things detected again, the reasons for bug are figured out. But the way of clear is still under investigation. It is because of negtive C, when exp log runs in wdot function, nan and inf occurs. But python would let it happen without warning. When it is calculated through solve_ivp, the bug shows out.

[Zeracesharon]

Testing cantera clippings and normalization, results show:
the function set_unnormalized_mass_fractions : no clipping & no normalization for Y, only normalization for X
Setting by TYX: clipping & both normalization for Y and X
Stting by TPY: clipping, no normalization for Y, normalization for X

Testing code:

import torch 
import numpy as np
import cantera as ct

mech_yaml = '../../data/IC8H18_reduced.yaml'
gas = ct.Solution(mech_yaml)

# Initial condition
P = ct.one_atm * 1
T = 1800
composition = 'IC8H18:0.5,O2:12.5,N2:34.0'

gas.TPX = T, P, composition
print('cantera setting by TPX')
if (gas.Y<0).any():
    print("no clipping for Y",np.where(gas.Y<0))
if (gas.concentrations<0).any():
    print("no clipping for concentration",np.where(gas.concentrations<0))
if (np.sum(gas.Y)!=1):
    print('unnormalization happen for Y')
if (np.sum(gas.X)!=1):
    print('unnormalization happen for X') 
    
    
print('cantera setting by unnormalized way')
Y=np.zeros(gas.n_species)
Y[0]=-0.4
Y[5]=-0.3
Y[7]=0.1
Y[10]=1.4
gas.set_unnormalized_mass_fractions(Y)
if (gas.Y<0).any():
    print("no clipping for Y",np.where(gas.Y<0))
if (gas.concentrations<0).any():
    print("no clipping for concentration",np.where(gas.concentrations<0))
if (gas.forward_rates_of_progress<0).any():
    print('forward rates of progress calculation without clipping',np.where(gas.forward_rates_of_progress<0))
if (np.sum(gas.Y)!=1):
    print('mass fraction unnormalization happen')
if (np.sum(gas.X)!=1):
    print('unnormalization happen for X') 

print('cantera setting by TPY')
gas.Y = Y
if (gas.Y<0).any():
    print("no clipping for Y",np.where(gas.Y<0))
if (gas.concentrations<0).any():
    print("no clipping for concentration",np.where(gas.concentrations<0))
if (np.sum(gas.Y)!=1):
    print('unnormalizations for Y')
if (np.sum(gas.X)!=1):
    print('unnormalization happen for X') 

Results:

[Zeracesharon]

Personal thought about this investigation:
the normalizaiton in rt : self.Y = (self.Y.T / self.Y.sum(dim=1)).T needs to be ignored. Here is the reason. Because when calculating jac numerically, it performs tiny mass fraction perturbation for a certain species k, then Yk_p=Yk+eps. Then the total mass fraction sum is not equal to 1. If you enforce it to be added to one, the remaining Y would change accordingly. Then the jac can't be calculated correctly since we only want to perturb a certain Yk in order to get its deriavative.
For exp log calculation: The reason why it fails to be correct in rt is because of when the Yk=0, Ck(the concentration of x) =0. if you perturbs Yk, Ck changes accordingly. When the solve_ivp tries to evaluate the jac, large tiny negative Ck would occur in order to obtain its deriative. At this time exp log would be large NAN & INF matrix because log(C) C<=0 evenrthough you implement C+eps, it still cannot ensure C+eps>0 since C would be something -1e-30. IF you enfore C.where(C<=0,eps,C), this would solve the problem of NAN & INF, but it is still problematic. Since when C =0, it perturbs to be C_p= -1e-30, the calculation would not change at all.But sometimes as the time goes, certain species mass fraction change from 0-tiny values. So this is also problematic. I am wondering, it may induced by its calculation method, solve_ivp use finite differences for calculating jac function.--------------details need to be explored and discussed further

[jiweiqi]

Thanks a lot for looking into this. It is a very interesting analysis. I think an easy way to workaround is

• adding an option to solution_update to choose whether doing normalization and clipping.
• add an option to wdot to choose whether evaluating wdot via log-exp or via iteration.

Then it should be identical to Cantera. Allowing tiny negative Y and X should be alright in the code.