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

• Added labels: ****
• Removed labels: ****

[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
• Added labels: ****
• Removed labels: ****

[speth]

The frequency of failures seems to increase as both the initial temperature and pressure are increased. Setting the initial conditions in the above example to 1200 K and 5 bar increases the failure rate to 111 out of 500 samples (22%).

[ischoegl]

File does not open - may no longer be reproducible.

[bryanwweber]

The input file is linked on the relevant Users' Group thread: https://groups.google.com/g/cantera-users/c/7QUmPpNmRfA/m/CJJTRzbMg-wJ

[speth]

The solver does not fail at the single point that was originally demonstrated from that input file in any case. The behavior for the wider range of test cases using the script I provided and the GRI 3.0 mechanism is still the same. Here is an updated version for modern Python and Cantera versions:

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

gas = ct.Solution('gri30.yaml')

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

    if abs(h1 - h0) > 1e-3 or abs(p1 - p0) > 1e-4:
        print('        ', h1-h0, p1-p0)

for i in range(500):
    g()
print(f'Finished with {fail_count} failures.')

As before, this fails on 29/500 cases as written and with T and P increased to 1200 K and 5 bar, the number of failures increases to 111/500.