Derivatives in multicomponent multiparameter mixtures returning NaN

[fefiedler]

teqp returns NaNs in multicomponent (more than 2 components) multiparameter models, if the molefraction of one of the components is set to 1.

import teqp
import numpy as np

model = teqp.build_multifluid_model(['hydrogen','argon','methane'], teqp.get_datapath())

temp = 300
dens = 5930.22080526693
derivs = [0,1,2,3]
molefrac = np.array([1.0,0.0,0.0])
for NT in derivs:
    for ND in derivs:
        if NT+ND < 4: print(temp,dens,NT,ND,model.get_Arxy(NT,ND,temp,dens,molefrac))

I have also tested a four component mixture, and additionally, everytime the first two components are set to zero, a NaN is returned.

import teqp
import numpy as np

model = teqp.build_multifluid_model(['hydrogen','argon','methane','nitrogen'], teqp.get_datapath())

temp = 300
dens = 5930.22080526693
derivs = [0,1,2,3]
molefrac = np.array([0.0,0.0,0.8,0.2])
for NT in derivs:
    for ND in derivs:
        if NT+ND < 4 :print(temp,dens,NT,ND,model.get_Arxy(NT,ND,temp,dens,molefrac))

[ianhbell]

I found the problem, it is

teqp/include/teqp/models/multifluid_reducing.hpp

Line 215 in 9729a32

sum2 = sum2 + 2.0 * z[i] * z[j] * (z[i] + z[j]) / (pow2(beta(i, j)) * z[i] + z[j]) * Yij(i, j);

when z[i] = 0 and z[j] = 0, so the term pow2(beta(i, j)) * z[i] + z[j] is also zero and you divide by zero. Fixing it looks trivial (just don’t add the term in the summation when z[i]=0 and z[j]=0) because the numerator is also zero, but I have dealt with similar problems and it is not quite as simple as it looks. The trivial solution is:

typename MoleFractions::value_type sum2 = 0.0;
for (auto i = 0U; i < N - 1; ++i) {
    for (auto j = i + 1; j < N; ++j) {
        auto den = (pow2(beta(i, j)) * z[i] + z[j]);
        if (getbaseval(den) != 0){
            sum2 = sum2 + 2.0 * z[i] * z[j] * (z[i] + z[j]) / den * Yij(i, j);
        }
    }
}

[ianhbell]

The trivial solution is also fine for first derivatives with respect to composition, but does NOT work for second derivatives. This represents a fundamental problem with the GERG reducing function.