Silent failure in MultiPhaseEquil ("gibbs" solver)

[GoogleCodeExporter]

In the following example, the state of the mixture is not modified, and no 
exception is generated. (Python interface, r3293)

import cantera as ct
data = """ideal_gas(name="gas", elements="C H",
                    species="nasa_gas: C2H5 C4H10,isobutane")"""
gas = ct.Solution(source=data)
gas.TPX = 300, 101325, [0.5, 0.5]
gas.equilibrate('TP','gibbs')
print gas.X # gives [0.5, 0.5], but should give approximately [0.0, 1.0]

For this initial condition, the element_potential solver fails by throwing an 
exception, and the VCS successfully finds the equilibrium composition.

Original issue reported on code.google.com by yarmond on 8 Nov 2014 at 12:58

[speth]

This bug is presumably the same as the one exercised by the (disabled) test OverconstrainedEquil.DISABLED_MultiphaseEquil in test/equil_gas.cpp, introduced in c2344377a.

[ischoegl]

Bug is present as of cantera 2.5.0a3

[speth]

No one has done any significant work on Cantera's equilibrium solvers in several years, so I certainly wouldn't expect this to have changed.

[ischoegl]

Thanks for the context: I just noticed that there are a lot of old equilibrium issues (currently 12 issues that are identified as bug, many date back to GoogleCode) that are still around. I.e. almost half of the known unresolved bugs are associated with Equilibrium.

Unfortunately, I am not familiar enough with the equilibrium solvers to make it efficient to tackle those myself. On the other hand, are there plans to address this for 2.5?

[bryanwweber]

No, there are no plans with respect to the equilibrium solvers. Feel free to dig in and learn something new 😊

[speth]

At least among the active Cantera developers, there is no one with the time, expertise, and interest to really dig into the equilibrium solvers and fix these issues. As such, no, there is currently no plan for when these issues might be addressed.

[ischoegl]

Alright. I thought I'd ask. In any case, I believe tackling the 1D code base has much higher pay-off: it's not that I'm opposed to learning something new, but it's not aligned enough to my core research.

[speth]

Updated example for the YAML era (no change to the results):

import cantera as ct
phasedef = """
phases:
- name: gas
  thermo: ideal-gas
  species: [{nasa_gas.yaml/species: [C2H5, 'C4H10,isobutane']}]
"""
gas = ct.Solution(yaml=phasedef)
gas.TPX = 300, 101325, [0.5, 0.5]
gas.equilibrate('TP','gibbs')
print(gas.X)  # gives [0.5, 0.5], but should give approximately [0.0, 1.0]

[speth]

Here's another example of what I would assume is the same bug, this time with a multiphase mixture:

import cantera as ct
solid_def = """
phases:
- name: ice
  thermo: fixed-stoichiometry
  species: [{nasa_condensed.yaml/species: [H2O(s)]}]
"""
gas_def = """
phases:
- name: vapor
  thermo: ideal-gas
  species: [{nasa_gas.yaml/species: [H2O]}]
"""

gas = ct.Solution(yaml=gas_def)
solid = ct.Solution(yaml=solid_def)
mix = ct.Mixture([gas, solid])
mix.T = 250
mix.P = 10*ct.one_atm
mix.species_moles = 'H2O:0.2, H2O(s):1.4'
mix.equilibrate('TP', solver='gibbs')
print(mix.species_moles)

which outputs the same initial composition:

[0.2 1.4]