Error in reaction path flux calculation

The reaction path flux calculation appears to fail for reactions such as:

H + HO2 <=> 2 OH

For example, running the following code:

import cantera as ct
g = ct.Solution('h2o2.xml')
g.TPX = 900, ct.one_atm, 'H:0.1, HO2:0.1, AR:6'
diagram = ct.ReactionPathDiagram(g, 'H')
ropf = g.forward_rates_of_progress
ropr = g.reverse_rates_of_progress
for i in range(g.n_reactions):
    if (ropf[i] > 1e-8 or ropr[i] > 1e-8) and ('H' in g.reaction(i).reactants or 'H' in g.reaction(i).products):
        print('{:2d}  {:25s}  {:11.5g} {:11.5g}'.format(i, g.reaction_equation(i), ropf[i], ropr[i]))
print(diagram.get_data())
print(g['H'].net_production_rates)

Produces the following output:

 5  H + O2 + M <=> HO2 + M               0  1.3514e-06
 8  AR + H + O2 <=> AR + HO2             0  1.5244e-05
10  2 H + M <=> H2 + M             0.46067           0
14  H + HO2 <=> H2O + O             130.12           0
15  H + HO2 <=> H2 + O2             1176.1           0
16  H + HO2 <=> 2 OH                2809.2           0

H HO2 H2 H2O H2O2 
H HO2 0 -1.65958e-05
H H2 1176.99 -0
H H2O 130.122 -0
H H2O2 0 -0
HO2 H2 1176.07 -0
HO2 H2O 130.122 -0
HO2 H2O2 79.6986 -0
H2 H2O 0 -0
H2 H2O2 0 -0
H2O H2O2 0 -0

[-4116.29642908]

From this, we see:

Reactions 5+8 give the correct flux for HO2 -> H
Reactions 10+15 give the correct flux for H -> H2
Reaction 14 gives the correct flux for H -> H2O
Reaction 16 is not contributing to the apparent flux of H -> OH
Reaction 16 is not contributing to the apparent flux of HO2 -> OH

In addition, the sum of all fluxes from H (-1.6e-5 - 1176 - 130 = -1306) do not equal the net production rate of H (-4116).

In Cantera 2.1.2, this calculation is differently wrong. Instead of giving zero flux for H -> OH and HO2 ->OH, it gives twice the correct value. A slightly modified version of the above script:

import cantera as ct
g = ct.Solution('h2o2.xml')
g.TPX = 900, ct.one_atm, 'H:0.1, HO2:0.1, AR:6'
diagram = ct.ReactionPathDiagram(g, 'H')
ropf = g.forward_rates_of_progress
ropr = g.reverse_rates_of_progress
for i in range(g.n_reactions):
    if (ropf[i] > 1e-8 or ropr[i] > 1e-8) and ('H' in g.reaction(i).reactants or 'H' in g.reaction(i).products):
        print('{:2d}  {:25s}  {:11.5g} {:11.5g}'.format(i, g.reaction_equation(i), ropf[i], ropr[i]))
print(diagram.get_data())
print(g['H'].net_production_rates)

gives the output:

 7  H + O2 + AR <=> HO2 + AR             0  1.5244e-05
 9  2 H + M <=> H2 + M             0.46067           0
13  H + HO2 <=> O + H2O             130.12           0
14  H + HO2 <=> O2 + H2             1176.1           0
15  H + HO2 <=> 2 OH                2809.2           0

H HO2 H2 H2O OH H2O2 
H HO2 0 -1.52444e-05
H H2 1176.99 -0
H H2O 130.122 -0
H OH 5618.36 -0
H H2O2 0 -0
HO2 H2 1176.07 -0
HO2 H2O 130.122 -0
HO2 OH 5618.36 -0
HO2 H2O2 79.6986 -0
H2 H2O 0 -0
H2 OH 0 -0
H2 H2O2 0 -0
H2O OH 0 -0
H2O H2O2 0 -0
OH H2O2 0 -0

[-4116.29643043]

speth@f2ff344 is a partial fix for this issue, and restores the double-counting behavior seen before 37f71bd9.

speth@12f2823 is an attempt to fix the double-counting seen when computing the flux of H->OH and HO2->OH in the reaction H + HO2 <=> 2 OH. It does achieve this goal. However, this change seems to then cause the flux of HO2->H2O2 from the reaction 2 HO2 <=> H2O2 + O2 not to be correct (where I think this flux should be twice the ROP for the reaction, since two H atoms are being transferred). As an example:

import cantera as ct
g = ct.Solution('h2o2.xml')

g.TPX = 900, ct.one_atm, 'H:0.1, HO2:0.1, AR:6'
diagram = ct.ReactionPathDiagram(g, 'H')
ropf = g.forward_rates_of_progress
ropr = g.reverse_rates_of_progress
for i in range(g.n_reactions):
    if (ropf[i] > 1e-8 or ropr[i] > 1e-8):
        print('{:2d}  {:25s}  {:11.5g} {:11.5g}'.format(
            i, g.reaction_equation(i), ropf[i], ropr[i]))

print()
for i,line in enumerate(diagram.get_data().split('\n')):
    s = line.split()
    if len(s) != 4:
        continue
    print('{:>6s} -> {:6s} {: 9.2f} {: 9.2f}'.format(
        s[0], s[1], float(s[2]), float(s[3])))

gives the output:

 5  H + O2 + M <=> HO2 + M               0  1.3514e-06
 8  AR + H + O2 <=> AR + HO2             0  1.5244e-05
10  2 H + M <=> H2 + M             0.46067           0
14  H + HO2 <=> H2O + O             130.12           0
15  H + HO2 <=> H2 + O2             1176.1           0
16  H + HO2 <=> 2 OH                2809.2           0
25  2 HO2 <=> H2O2 + O2             15.426           0
26  2 HO2 <=> H2O2 + O2             24.423           0

     H -> HO2         0.00     -0.00
     H -> H2       1176.53     -0.00
     H -> H2O       130.12     -0.00
     H -> OH       2809.18     -0.00
     H -> H2O2        0.00     -0.00
   HO2 -> H2       1176.07     -0.00
   HO2 -> H2O       130.12     -0.00
   HO2 -> OH       2809.18     -0.00
   HO2 -> H2O2       39.85     -0.00
    H2 -> H2O         0.00     -0.00
    H2 -> OH          0.00     -0.00
    H2 -> H2O2        0.00     -0.00
   H2O -> OH          0.00     -0.00
   H2O -> H2O2        0.00     -0.00
    OH -> H2O2        0.00     -0.00

Cantera / cantera

Error in reaction path flux calculation #377