possible bug in default third-body efficiency

[michael-a-hansen]

I think there's a bug in how pyJac handles three-body and falloff reactions when the default third-body efficiency is not 1. In verifying a Jacobian code I'm working on, I've tested against pyJac on 3-body, Lindemann, and Troe reactions, with the N-th species as a reactant or not, and with the N-th species as a special third-body or not. In this test matrix there are persistent diffs between pyJac and cantera/my-code (which do not diff) in the source term (RHS) calculation only when the default third-body efficiency is not zero. Due to this, I don't think it's an N-th species thing.

The most serious impact would be the RHS. Here's code to show pyJac vs Cantera. Attached are two toy mechanisms with only one reaction (reversible, three-body, N-th is not a reactant, N-th is not a special TB). Swap out the mech variable in the code to see the difference.

I couldn't upload an XML to github, so they're in txt format, but just change the extension and it should work fine. The code expects each pyJac-cython binary file to be put in a folder named like the mechanism.

rev_3body_noN_noN_nonUnity.txt rev_3body_noN_noN_unity.txt

from numpy import hstack, max, abs, ones, zeros
from importlib import util
from cantera import Solution, one_atm
from os.path import join, abspath

mech = 'rev_3body_noN_noN_unity'  # pyjac will match cantera
#mech = 'rev_3body_noN_noN_nonUnity'  # pyjac will not match cantera

xml =  mech + '.xml'
dest = mech + '/'  # put the pyjac-generated binary file in a folder named like the mechanism

magic_name = 'pyjacob.cpython-36m-darwin.so'

pyjac_abs_path = join(abspath(dest), magic_name)
pyjacob_module_spec = util.spec_from_file_location('pyjacob', pyjac_abs_path)
try:
  pyjac_module = util.module_from_spec(pyjacob_module_spec)
  pyjacob_module_spec.loader.exec_module(pyjac_module)
except ImportError:
  raise ImportError('You have not built a pyJac wrapper for the specified mechanism!')

pyjac_rhs = pyjac_module.py_dydt
pyjac_jac = pyjac_module.py_eval_jacobian

gas = Solution(xml)

ns = gas.n_species
p = one_atm
T = 1000.
gas.TPX = T, p, ones(ns)

state = hstack((T, gas.Y))

rhsPJ = zeros(ns)
jacPJ = zeros(ns * ns)
pyjac_rhs(0., p, state, rhsPJ)
pyjac_jac(0., p, state, jacPJ)

rhsCN = gas.net_production_rates * gas.molecular_weights / gas.density

print('RHS of the species equations:')
print('pyjac:', rhsPJ[1:])
print('cantera:', rhsCN[:-1])

print('pyjac vs cantera:', max(abs(rhsCN[:-1] - rhsPJ[1:])))

[skyreflectedinmirrors]

@michael-a-hansen thanks for the report. Looking at the two mechanisms, I think I see the issue.

On line 1031 (and 991 for that matter) we check for a default efficiency of 0, and if so treat the first non-zero efficiency as a single species third-body efficiency (i.e., with an efficiency of 1, and all other species with and efficiency of 0). Of course, we have a species "C" with efficiency = 2

I think the most robust way to handle this case would be to treat this third-body, would be as a "mixture" and then go set the third body efficiencies of the rest of the species in the mechanism to zero -- this would also fix the case (which we obviously haven't encountered yet) that the default efficiency was not 0 or 1, or if you set the default efficiency to 0 and then had 2+ species with specified efficiencies.

I'm not sure this specification is "correct" strictly speaking -- in chemkin wouldn't you be forced to go in and manually specify the third-body efficiency? None the less, if Cantera supports this sort of specification we probably should as well.

When I get a chance today, I'll test this out (by setting the "C" efficiency to 1) to see if that fixes the problem

[skyreflectedinmirrors]

@michael-a-hansen I think the latest commits (in pull request currently, will merge shortly) should fix this, my updated test script:

from numpy import hstack, max, abs, ones, zeros
from importlib import util
from cantera import Solution, one_atm, gas_constant
from os.path import join, abspath
from sys import argv

if argv[1] == 'unity':
  print('unity')
  mech = 'rev_3body_noN_noN_unity'  # pyjac will match cantera
else:
  print('non_unity')
  mech = 'rev_3body_noN_noN_nonUnity'  # pyjac will not match cantera

xml =  mech + '.xml'
dest = mech + '/'  # put the pyjac-generated binary file in a folder named like the mechanism

magic_name = 'pyjacob.cpython-36m-x86_64-linux-gnu.so'

pyjac_abs_path = join(abspath(dest), magic_name)
pyjacob_module_spec = util.spec_from_file_location('pyjacob', pyjac_abs_path)
try:
  pyjac_module = util.module_from_spec(pyjacob_module_spec)
  pyjacob_module_spec.loader.exec_module(pyjac_module)
except ImportError:
  raise ImportError('You have not built a pyJac wrapper for the specified mechanism!')

pyjac_rhs = pyjac_module.py_dydt
pyjac_jac = pyjac_module.py_eval_jacobian
pyjac_pmod = pyjac_module.py_get_rxn_pres_mod

gas = Solution(xml)

ns = gas.n_species
p = one_atm
T = 1000.
gas.TPX = T, p, ones(ns)

state = hstack((T, gas.Y))

rhsPJ = zeros(ns)
jacPJ = zeros(ns * ns)
pyjac_rhs(0., p, state, rhsPJ)
pyjac_jac(0., p, state, jacPJ)

rhsCN = gas.net_production_rates * gas.molecular_weights / gas.density

print('RHS of the species equations:')
print('pyjac:', rhsPJ[1:])
print('cantera:', rhsCN[:-1])

print('pyjac vs cantera:', max(abs(rhsCN[:-1] - rhsPJ[1:])))

print('pres mod')
pmodPJ = zeros(1)
pmodAN = zeros(1)
pyjac_pmod(T, p, gas.concentrations, pmodPJ)
rxn = gas.reactions()[0]
pmodAN[0] = rxn.default_efficiency * p / (gas_constant * T) + sum([(eff - rxn.default_efficiency) * gas.concentrations[gas.species_index(sp)] for sp, eff in rxn.efficiencies.items()])
print('pyjac:', pmodPJ[0])
print('analytical:', pmodAN[0])
print('pyjac vs analytical:', max(abs(pmodPJ[0] - pmodAN)))

and output:

non_unity
RHS of the species equations:
pyjac: [ 0.32158177 -0.16079088  0.          5.10458496 -5.26537584]
cantera: [ 0.32158177 -0.16079088  0.          5.10458496 -5.26537584]
pyjac vs cantera: 7.9936057773e-14
pres mod
pyjac: 0.00406219904472
analytical: 0.00406219904472
pyjac vs analytical: 3.46944695195e-18

Can you give it a shot when you get a chance?