Numerical issues with Wilhoit-toMultiNASA conversion, apparently associated with use of indefinite analytic integrals

[gmagoon]

I've come across a tricky case that seems to cause numerical problems in the Wilhoit-to-MultiNASA conversion:

>>> w2 = thermo.Wilhoit(4.0,75.0,-54.497423456539401,303.4090028535806,-727.01542657068092,659.08598716992981,1234.5,789.0,B=3.0)
>>> ans = thermo.Wilhoit_integral_T0(w2,6.0) - thermo.Wilhoit_integral_T0(w2,.298)
>>> print ans
-27.3531739832
>>> print (thermo.Wilhoit_integral_TM1(w2,1.5) -thermo.Wilhoit_integral_TM1(w2,.298))
69.7868424873
>>> print (thermo.Wilhoit_integral2_TM1(w2,1.5) -thermo.Wilhoit_integral2_TM1(w2,.298))
3263.51989762
>>> print thermo.TintOpt_objFun_W(0.65, w2, .298, 1.5, 3)
-8.5055398813e-07
>>> print thermo.TintOpt_objFun_W(0.60, w2, .298, 1.5, 3)
-1.57702925208e-06
>>> print thermo.TintOpt_objFun_W(0.6037, w2, .298, 1.5, 3)
2.38989377976e-06
>>> print thermo.TintOpt_objFun_W(0.615, w2, .298, 1.5, 3)
8.49388015922e-07
>>> print thermo.TintOpt_objFun(0.60, w2, .298, 1.5, True, 3)
ERROR:root:Greg thought he fixed the numerical problem, but apparently it is still an issue; please e-mail him with the following results:
ERROR:root:0.6
ERROR:root:Wilhoit(cp0=4, cpInf=75, a0=-54.497423456539401, a1=303.4090028535806, a2=-727.01542657068092, a3=659.08598716992981, H0=1234.5, S0=789, B=3)
ERROR:root:0.298
ERROR:root:1.5
ERROR:root:True
ERROR:root:-1.57702925208e-06
0
>>> print thermo.Wilhoit_integral2_TM1(w2,1.5)
-574582.356995
>>> print thermo.Wilhoit_integral2_TM1(w2,.298)
-577845.876893

The above shows:
*Negative values of the objective function with fairly significant magnitude; these are inaccurate as the objective function is an integral of squared error
*Jittery, non-smooth behavior of the objective function
*Significant magnitude of the indefinite integrals, which are used to compute the definite integrals (by difference)

My investigation of this issue suggests is that it ultimately arises from numerical issues associated with the use of intermediate indefinite analytic integrals with large magnitude in the calculation. It may be possible to implement definite analytic integrals that avoid differencing of large magnitude intermediate values.

Using an older version of RMG-Py that has numerical integral capabilities, I have confirmed that use of definite numerical integrals smooths the objective function and keeps it positive.

Using the older version with the numerical option (definite numerical integrals in both objective function and NASA coefficient calculation):

>>> print thermo.Cp_TintOpt_objFun(0.60, w2, .298, 1.5, 1, 3)
8.57655311464e-005
>>> print thermo.Cp_TintOpt_objFun(0.6037, w2, .298, 1.5, 1, 3)
8.38349646131e-005
>>> print thermo.Cp_TintOpt_objFun(0.615, w2, .298, 1.5, 1, 3)
7.93030143242e-005
>>> print thermo.Cp_TintOpt_objFun(0.630, w2, .298, 1.5, 1, 3)
7.61815561034e-005
>>> print thermo.Cp_TintOpt_objFun(0.640, w2, .298, 1.5, 1, 3)
7.57683296015e-005
>>> print thermo.Cp_TintOpt_objFun(0.650, w2, .298, 1.5, 1, 3)
7.65980407849e-005
>>> print thermo.Cp_TintOpt_objFun(1.0, w2, .298, 1.5, 1, 3)
0.00101810388993
>>> print thermo.TintOpt_objFun(1.0, w2, .298, 1.5, 1, 3)
0.000430559256529

Using a modified version of the aforementioned which uses definite numerical integrals for just the objective function (indefinite analytic integrals used in matrix for NASA coefficient calculation):

>>> from rmg import thermo
>>> w2 = thermo.ThermoWilhoitData(4.0,75.0,-54.497423456539401,303.4090028535806,-727.01542657068092,659.08598716992981,1234.5,789.0,B=3.0)
>>> print (w2.integral_T0(6.0) - w2.integral_T0(.298))
-27.3531739832
>>> print thermo.TintOpt_objFun(0.60, w2, .298, 1.5, 1, 3)
7.56744975661e-006
>>> print thermo.TintOpt_objFun(1.0, w2, .298, 1.5, 1, 3)
0.00043583716888
>>> print thermo.TintOpt_objFun(0.61, w2, .298, 1.5, 1, 3)
7.57669477025e-006
>>> print thermo.TintOpt_objFun(0.62, w2, .298, 1.5, 1, 3)
7.7480199252e-006
>>> print thermo.TintOpt_objFun(0.63, w2, .298, 1.5, 1, 3)
8.12499638414e-006
>>> print thermo.TintOpt_objFun(0.64, w2, .298, 1.5, 1, 3)
8.66185655468e-006
>>> print thermo.TintOpt_objFun(0.59, w2, .298, 1.5, 1, 3)
7.68816880736e-006
>>> print thermo.TintOpt_objFun(0.58, w2, .298, 1.5, 1, 3)
8.1157850218e-006
>>> print thermo.TintOpt_objFun(0.603, w2, .298, 1.5, 1, 3)
7.51797415433e-006
>>> print thermo.TintOpt_objFun(0.607, w2, .298, 1.5, 1, 3)
7.57045654609e-006

My suspicion is that the numerical (definite) integrals are less accurate, but also less noisy/jittery, than the result computed via differencing of indefinite analytic integrals.

A problem with similar symptoms was encountered in issue #20, but this was a cython issue; this issue does not appear to be related to cython, and instead seems to be a more fundamental numerical difficulty.

[gmagoon]

Some further evidence that this is a numerical issue associated with the way in which the indefinite integrals are combined:
In the old version of RMG-Py on my laptop, I rearranged the objective function from:

    #qM1=wilhoit.integral_TM1(tint)
    #q0=wilhoit.integral_T0(tint)
    #q1=wilhoit.integral_T1(tint)
    #q2=wilhoit.integral_T2(tint)
    #q3=wilhoit.integral_T3(tint)
    #result = (wilhoit.integral2_TM1(tmax) - wilhoit.integral2_TM1(tmin) +
    #            nasa_low.integral2_TM1(tint)-nasa_low.integral2_TM1(tmin) + nasa_high.integral2_TM1(tmax) - nasa_high.integral2_TM1(tint)
    #            - 2* (b6*(wilhoit.integral_TM1(tmax)-qM1)+b1*(qM1 - wilhoit.integral_TM1(tmin))
    #            +b7*(wilhoit.integral_T0(tmax)-q0)+b2*(q0 - wilhoit.integral_T0(tmin))
    #            +b8*(wilhoit.integral_T1(tmax)-q1)+b3*(q1 - wilhoit.integral_T1(tmin))
    #            +b9*(wilhoit.integral_T2(tmax)-q2)+b4*(q2 - wilhoit.integral_T2(tmin))
    #            +b10*(wilhoit.integral_T3(tmax)-q3)+b5*(q3 - wilhoit.integral_T3(tmin))))

to

    result = (wilhoit.integral2_TM1(tmax) - wilhoit.integral2_TM1(tmin) +
                 nasa_low.integral2_TM1(tint)-nasa_low.integral2_TM1(tmin) + nasa_high.integral2_TM1(tmax) - nasa_high.integral2_TM1(tint)
                 - 2* (b6*wilhoit.integral_TM1(tmax)-b1*wilhoit.integral_TM1(tmin)+(b1-b6)*wilhoit.integral_TM1(tint)
                 +b7*wilhoit.integral_T0(tmax)-b2*wilhoit.integral_T0(tmin) +(b2-b7)*wilhoit.integral_T0(tint)
                 +b8*wilhoit.integral_T1(tmax)-b3*wilhoit.integral_T1(tmin) +(b3-b8)*wilhoit.integral_T1(tint)
                 +b9*wilhoit.integral_T2(tmax)-b4*wilhoit.integral_T2(tmin) +(b4-b9)*wilhoit.integral_T2(tint)
                 +b10*wilhoit.integral_T3(tmax)-b5*wilhoit.integral_T3(tmin)+(b5-b10)*wilhoit.integral_T3(tint)))

Results below show that the results have changed but are still jittery and still go negative:

>>> from rmg import thermo
>>> w2 = thermo.ThermoWilhoitData(4.0,75.0,-54.497423456539401,303.4090028535806,-727.01542657068092,659.08598716992981,1234.5,789.0,B=3.0)
>>> print (w2.integral_T0(6.0) - w2.integral_T0(.298))
-27.3531739832
>>> (w2.integral_T0(6.0) - w2.integral_T0(.298))
-27.353173983239685
>>> print thermo.TintOpt_objFun(0.60, w2, .298, 1.5, 1, 3)
1.8864593585e-006
>>> print thermo.TintOpt_objFun(0.61, w2, .298, 1.5, 1, 3)
8.65593392518e-007
>>> print thermo.TintOpt_objFun(0.62, w2, .298, 1.5, 1, 3)
6.189036867e-006
>>> print thermo.TintOpt_objFun(0.63, w2, .298, 1.5, 1, 3)
3.51085418515e-006
>>> print thermo.TintOpt_objFun(0.64, w2, .298, 1.5, 1, 3)
Error: :root:Greg thought he fixed the numerical problem, but apparently it is still an issue; please e-mail him with the following results:
Error: :root:0.64
Error: :root:ThermoWilhoitData(4,75,-54.5,303.4,-727,659.1,1235,789,'',B=3)
Error: :root:0.298
Error: :root:1.5
Error: :root:1
Error: :root:-4.48781247542e-006
0
>>> print thermo.TintOpt_objFun(0.59, w2, .298, 1.5, 1, 3)
2.69566317002e-006
>>> print thermo.TintOpt_objFun(0.603, w2, .298, 1.5, 1, 3)
6.85598661221e-006
>>> print thermo.TintOpt_objFun(0.607, w2, .298, 1.5, 1, 3)
Error: :root:Greg thought he fixed the numerical problem, but apparently it is still an issue; please e-mail him with the following results:
Error: :root:0.607
Error: :root:ThermoWilhoitData(4,75,-54.5,303.4,-727,659.1,1235,789,'',B=3)
Error: :root:0.298
Error: :root:1.5
Error: :root:1
Error: :root:-5.31995374331e-006
0
>>> print thermo.TintOpt_objFun(1.0, w2, .298, 1.5, 1, 3)
0.000430798881098
>>> 

[gmagoon]

Taking a closer look at the form of the integrals, it is clear that there is much room for improvement in of the grouping of the terms (which could also improve speed) even with the indefinite form, and using the definite form (grouping differences between the T components) should give further gains. However, it looks like it will be fairly time consuming, and possibly error-prone, however (and there's no guarantee it will work), so I'm trying to decide whether it is worth my time to try to address this.

Another avenue for improvement would be to divide the Wilhoit form by CpInf-Cp0 and then rescale the result by CpInf-Cp0 at the end...this would simplify things, leaving only one additive constant in the Wilhoit form to be integrated. The only reason I'm not planning to do this is that it would probably cause problems for cases where CpInf=Cp0 (monoatomic cases).

[gmagoon]

I think I've partially solved the problem, at least for this test case by using definite integrals (which took a while to manually tweak) for the Wilhoit and NASA forms. The result is still jittery, though it seems to be significantly less so, as the result doesn't go negative. My tests suggested that it was the Wilhoit changes rather than the NASA integral conversions that had the major effect.

From my old test version:

>>> from rmg import thermo
>>> w2 = thermo.ThermoWilhoitData(4.0,75.0,-54.497423456539401,303.4090028535806,-727.01542657068092,659.08598716992981,1234.5,789.0,B=3.0)
>>> print thermo.TintOpt_objFun(0.60, w2, .298, 1.5, 1, 3)
8.59392912389e-006
>>> print thermo.TintOpt_objFun(0.59, w2, .298, 1.5, 1, 3)
1.00868965092e-005
>>> print thermo.TintOpt_objFun(0.61, w2, .298, 1.5, 1, 3)
9.79530705081e-006
>>> print thermo.TintOpt_objFun(0.62, w2, .298, 1.5, 1, 3)
8.33310332382e-006
>>> print thermo.TintOpt_objFun(0.63, w2, .298, 1.5, 1, 3)
1.22686979012e-005
>>> print thermo.TintOpt_objFun(0.64, w2, .298, 1.5, 1, 3)
1.03765642052e-005
>>> print thermo.TintOpt_objFun(0.635, w2, .298, 1.5, 1, 3)
1.23651734611e-005
>>> for i in range(0,100):
    print thermo.TintOpt_objFun(0.55+i/1000., w2, .298, 1.5, 1, 3)

    
1.34421707116e-005
1.41038044603e-005
1.37015322252e-005
1.41987129609e-005
1.15972397907e-005
1.42528979268e-005
1.40551401273e-005
1.29553900479e-005
1.19883588923e-005
1.38190098369e-005
1.38874647746e-005
1.253110986e-005
8.98815051187e-006
1.20821696328e-005
1.28604960992e-005
1.31793585751e-005
1.20991844597e-005
1.21386410683e-005
1.20838321891e-005
1.19067162814e-005
1.06019297164e-005
1.12856150736e-005
1.23276258819e-005
1.0343734175e-005
1.18273310363e-005
7.88652505435e-006
1.06557554318e-005
1.14954764285e-005
8.84354994923e-006
8.9260047389e-006
1.12938914754e-005
1.14846980068e-005
1.02277108454e-005
7.1708982432e-006
1.20056065498e-005
1.07011455839e-005
7.05861475581e-006
7.96390122559e-006
1.08129443106e-005
1.09930351755e-005
1.02569028968e-005
8.14236136648e-006
1.25498600028e-005
1.05191002149e-005
9.18800651561e-006
7.94650350144e-006
1.17976205729e-005
1.04315849967e-005
1.03236498035e-005
8.20821242087e-006
8.80482184584e-006
1.04079445009e-005
8.35850187286e-006
7.4587205745e-006
7.96025869931e-006
1.1753495528e-005
1.11125600597e-005
6.77075240674e-006
7.43517557567e-006
1.05153949335e-005
8.11813879409e-006
9.75761031441e-006
7.6846781667e-006
1.10041064545e-005
1.15163629744e-005
7.72559815232e-006
8.35637183627e-006
1.03646934804e-005
8.84644123289e-006
1.02995072666e-005
8.55287908053e-006
1.20656395666e-005
1.16160363177e-005
8.92175648914e-006
8.50278502185e-006
8.40010125103e-006
1.05883291326e-005
1.15581678983e-005
8.00561883807e-006
7.82392726251e-006
1.22686979012e-005
9.31448721531e-006
8.61588159751e-006
8.66274694999e-006
1.12772413559e-005
1.23651734611e-005
8.77748789208e-006
7.7212071119e-006
1.08172935143e-005
9.9903827504e-006
1.03765642052e-005
9.34592662816e-006
1.21759076137e-005
1.14395434139e-005
9.87882322079e-006
9.46178624872e-006
1.35125619636e-005
9.95414120553e-006
1.06116413008e-005
1.08883887151e-005
>>> print thermo.TintOpt_objFun(1.0, w2, .298, 1.5, 1, 3)
0.000438772648522
>>> 

I'll push a fix shortly once I update the cython .pxd file.

[gmagoon]

Another potential avenue for improvements (which I'm not going to try at the moment due to time constraints and diminishing returns): instead of B/(T+B) for Z, use 1/(T+B); then after differencing with the two temperatures, remultiply by B (to the appropriate power).