SLACKHA / pyJac

Creates C and CUDA analytical Jacobians for chemical kinetics ODE systems
http://slackha.github.io/pyJac/
MIT License
52 stars 23 forks source link

Jacobian creation with SRI Falloff Function #12

Closed TJaravel closed 7 years ago

TJaravel commented 7 years ago

An error occurs in create_jacobian.py when treating a reaction with SRI Falloff function.

The following fix may work l249 of create_jacobian.py Current:

    _elif rxn.sri:
        jline += '* pow({:.6} * '.format(reac.sri_par[0])
        # Need to check for negative parameters, and
        # skip "-" sign if so.
        if reac.sri_par[1] > 0.0:
            line += 'exp(-{:.6} / T)'.format(reac.sri_par[1])
        else:
            line += 'exp({:.6} / T)'.format(abs(reac.sri_par[1]))

        if reac.sri_par[2] > 0.0:
            line += ' + exp(-T / {:.6}), X) '.format(reac.sri_par[2])
        else:
            line += ' + exp(T / {:.6}), X) '.format(abs(reac.sri_par[2]))

        if (len(reac.sri_par) == 5 and
                reac.sri_par[3] != 1.0 and reac.sri_par[4] != 0.0):
            line += ('* {:.8e} * '.format(reac.sri_par[3]) +
                     'pow(T, {:.6}) '.format(reac.sri_par[4])_

Suggested fix:

    elif rxn.sri:
        jline += '* pow({:.6} * '.format(rxn.sri_par[0])
        # Need to check for negative parameters, and
        # skip "-" sign if so.
        if rxn.sri_par[1] > 0.0:
            jline += 'exp(-{:.6} / T)'.format(rxn.sri_par[1])
        else:
            jline += 'exp({:.6} / T)'.format(abs(rxn.sri_par[1]))

        if rxn.sri_par[2] > 0.0:
            jline += ' + exp(-T / {:.6}), X) '.format(rxn.sri_par[2])
        else:
            jline += ' + exp(T / {:.6}), X) '.format(abs(rxn.sri_par[2]))

        if (len(rxn.sri_par) == 5 and
                rxn.sri_par[3] != 1.0 and rxn.sri_par[4] != 0.0):
            jline += ('* {:.8e} * '.format(rxn.sri_par[3]) +
                     'pow(T, {:.6}) '.format(rxn.sri_par[4])
skyreflectedinmirrors commented 7 years ago

Thanks for the bug report, the current version is definitely incorrect. I will test out your fix / incorporate ASAP

kyleniemeyer commented 7 years ago

Strange that we didn't encounter this when testing the model with SRI reactions—does this only occur when one of the optional parameters is used (or not used)?

skyreflectedinmirrors commented 7 years ago

@kyleniemeyer It appears that the SRI reaction in the Sarathy i-pentanol model does not have any non-default third-body efficiencies, hence the containing if statement on line 234 is not triggered

@TJaravel would you mind sharing the model you were using when this error occurred (or really, just the SRI reaction in question)? It would help test that the fix is correct (as compared to the automatic differentiation soln.)

TJaravel commented 7 years ago

The SRI reaction (chemkin):

H+CH3(+M)=CH4(+M) 1.2000e+15 -0.400 0.00 LOW/ 6.40e+23 -1.800 0.0/ SRI/ 0.4500 797.0000 979.0000 1.0000 0.0000/ H2/ 2.00/ CO/ 2.00/ CO2/ 3.00/ H2O/ 5.00/

skyreflectedinmirrors commented 7 years ago

Ok, running the functional tester using one of the GRI PaSR outputs with this reaction inserted in place of the similar Troe reaction in GRIMech 3.0 results in a max error of 0.436%

I think this fix is good.

Thanks again for the report @TJaravel