ChemEquil solver fails or generates SIGFPE

[GoogleCodeExporter]

Certain initial conditions cause unexpected failures in the ChemEquil 
equilibrium solver. For example, with the attached mechanism:

    import Cantera
    gas = Cantera.IdealGasMix('chemC2H4_100.xml')

    def f(i, Yf):
        gas.set(P=100002, T=299.999, Y='C2H4:%f, N2:0.594407, O2:0.179732' % Yf)
        gas.equilibrate('HP', i)

f(0, 0.225861) # Fails after exceeding the maximum number of iterations
f(1, 0.225861) # Succeeds
f(2, 0.225861) # Succeeds
f(0, 0.226861) # Succeeds

In addition, if trapping of floating point exceptions is enabled, e.g. by 
placing a call to:

    #include <fenv.h>
    feenableexcept(FE_OVERFLOW | FE_UNDERFLOW);

Somewhere before the call to equilibrate (I put it in the constructor of Phase, 
just to test), an exception will be caught in ChemEquil::equilibrate (approx. 
line 876 of ChemEquil.cpp), when trying to compute:

    jac.leftMult(DATA_PTR(res_trial), DATA_PTR(grad));

Where 'jac' looks like:

-3.4455e-08,  5.0889e-09,  1.7723e-08, 4.3784e+147, -3.8413e-06
 1.7243e-07, -3.0886e-07, -3.8914e-07, 6.0754e+154,  3.5798e-05
 2.0702e-07, -2.2973e-07, -5.4029e-07, 6.6276e+154,  4.3018e-05
 0.0000e+00,  5.1571e-08,  5.4633e-08, 1.0121e+148, -3.0169e-06
 1.6591e-06, -1.1083e-06, -3.0991e-06, 3.4299e+155,  1.1981e+01

The overflow is caused by the entries in the 4th column.

I think this has something to do with the fact that N2 is the only 
nitrogen-containing species defined in the mechanism, and N is the 4th element 
defined. This doesn't explain why the error only occurs for certain mixture 
compositions, though.

Original issue reported on code.google.com by yarmond on 29 Mar 2013 at 5:30

Attachments:

• chemC2H4_100.xml

[GoogleCodeExporter]

Note: this issue was originally reported in the Cantera User's Group:

https://groups.google.com/forum/?fromgroups=#!topic/cantera-users/7QUmPpNmRfA

Original comment by yarmond on 29 Mar 2013 at 5:31

[GoogleCodeExporter]

The changes in r2550 and r2551 eliminate the floating point exceptions. 
However, there are still cases where the solver unexpectedly fails to converge. 
The following script (for the 'new' Python module) which tests random (but 
repeatable) mixture compositions has a failure rate of 5.8%:

import cantera
from random import random as rand
import random
random.seed(0)

gas = cantera.Solution('gri30.xml')
#gas = cantera.Solution('chemC2H4_100.xml')

fail_count = 0
def g():
    a,b,c = rand(), 2*rand(), rand()
    gas.TPX = 299.999, 100002, 'C2H4:%f, N2:%f, O2:%f' % (a,b,c)
    h0,p0 = gas.HP
    print a,b,c
    try:
        gas.equilibrate('HP', 0)
    except:
        print 'failed at', gas.TPX
        global fail_count
        fail_count += 1
    h1,p1 = gas.HP
    print '        ',h1-h0, p1-p0

for i in range(500):
    print '======', i, '======'
    g()
print 'Finished with %d failures.' % fail_count

Original comment by yarmond on 30 Sep 2013 at 2:04

• Changed state: InProgress