InterfaceKinetics::solvePseudoSteadyStateProblem does not enforce coverage constraints

[skyreflectedinmirrors]

Cantera version

Dev

Operating System

Ubuntu

Python/MATLAB version

3.6

Expected Behavior

Solution given by InterfaceKinetics::solvePsudeoSteadyStateProblem should at least satisfy the constraint sum(coverages) = 1. Ideally, matching the output of InterfaceKinetics::advanceCoverages.

Actual Behavior

It appears to find a multiple of the steady-state solution that does not satisfy the constraint.

Steps to reproduce

1. Download the SiF4_NH3_mec.cti model from the 1D PFR notebook

2. Save the following as test.cpp

#include "cantera/zerodim.h"
#include "cantera/IdealGasMix.h"
#include "cantera/IncompressibleSolid.h"
#include "cantera/Interface.h"

using namespace Cantera;

void run(bool advance_coverages=false, int ifunc=-1)
{
    IdealGasMix* gas = new IdealGasMix("SiF4_NH3_mec.cti", "gas");
    gas->setState_TPX(1713, 2 * 101325 / 760.0, "NH3:6, SiF4:1");
    IncompressibleSolid* NBulk = new IncompressibleSolid("SiF4_NH3_mec.cti", "NBulk");
    NBulk->setState_TP(1713, 2 * 101325 / 760.0);
    IncompressibleSolid* SiBulk = new IncompressibleSolid("SiF4_NH3_mec.cti", "SiBulk");
    SiBulk->setState_TP(1713, 2 * 101325 / 760.0);
    std::vector<ThermoPhase*> phases({gas, NBulk, SiBulk});
    Interface* interface = importInterface("SiF4_NH3_mec.cti", "SI3N4", phases);
    interface->setState_TP(1713, 2 * 101325 / 760.0);

    std::cout << (advance_coverages ? "Advance" : "Pseudo") << std::endl;
    std::cout << "ifunc = " << ifunc << std::endl;
    if (advance_coverages)
    {
        interface->advanceCoverages(100);
    }
    else
    {
        try
        {
            interface->solvePseudoSteadyStateProblem(ifunc);
        }
        catch(Cantera::CanteraError &e)
        {
            std::cerr << e.what() << std::endl;
        }
    }

    std::vector<double> cov(interface->nSpecies());
    interface->getCoverages(&cov[0]);
    std::cout << "cov = {";
    for (size_t i = 0; i < interface->nSpecies(); ++i)
    {
        if (i)
            std::cout << ", ";
        std::cout << cov[i];
    }
    std::cout << "}" << std::endl;
}

int main()
{
    run(true);
    for (int i = -1; i <= 4; ++i)
    {
        run(false, i);
    }

}

3. Compile, e.g., g++ -std=c++11 test.cpp -lcantera -I$CONDA_PREFIX/include -L$CONDA_PREFIX/lib -lpthread

4. Obtain the following output:

./a.out 
Advance
ifunc = -1
cov = {0.0625701, 0.915542, 0.000314168, 0.0208512, 0.00024098, 0.00048196}
Pseudo
ifunc = -1
cov = {0.125201, 1.83197, 0.000628639, 0.0417225, 0.00024098, 0.00048196}
Pseudo                                                                                                                                                                  
ifunc = 0                                                                                                                                                               
                                                                                                                                                                        
***********************************************************************                                                                                                 
CanteraError thrown by solve(DenseMatrix& A, double* b):                                                                                                                
Matrix appears to be rank-deficient: non-zero pivots = 5; columns = 6                                                                                                   
***********************************************************************                                                                                                 
                                                                                                                                                                        
cov = {1, 0, 0, 0, 0, 4e-10}                                                                                                                                            
Pseudo
ifunc = 1
cov = {0.125201, 1.83197, 0.000628639, 0.0417225, 0.00024098, 0.00048196}
Pseudo
ifunc = 2

***********************************************************************
CanteraError thrown by solve(DenseMatrix& A, double* b):
Matrix appears to be rank-deficient: non-zero pivots = 5; columns = 6
***********************************************************************

cov = {1, 0, 0, 0, 0, 4e-10}
Pseudo
ifunc = 3

***********************************************************************
CanteraError thrown by solve(DenseMatrix& A, double* b):
Matrix appears to be rank-deficient: non-zero pivots = 5; columns = 6
***********************************************************************

cov = {1, 0, 0, 0, 0, 4e-10}
Pseudo
ifunc = 4
cov = {0.125201, 1.83197, 0.000628639, 0.0417225, 0.00024098, 0.00048196}

It appears that all the various solution methods either fail, or find the same (incorrect) answer at least.