search

Cantera / cantera

Chemical kinetics, thermodynamics, and transport tool suite
https://cantera.org
Other
615 stars 349 forks source link

Bug in handling of surface reactions imported from YAML #902

Closed ischoegl closed 4 years ago

ischoegl commented 4 years ago

Problem description

The switch from methane_pox_on_pt.cti to methane_pox_on_pt.yaml input apparently broke the Python example surf_pfr.py (see original issue report #877).

Edit: The example catalytic_combustion.py is affected by a closely related issue (see #873).

Steps to reproduce

A run of the Python example surf_pfr.py currently fails (see #877). Switching from YAML back to CTI input restores the old behavior. It can be easily verified that re-creating YAML using cti2yaml does not fix the issue. The same issue also arises if the original file is converted from XML (via ctml_writer and ctml2yaml).

Behavior

Time integration reduces CO levels to negative values until the density itself becomes negative and an error is raised (see #877), which appears to be either due to a faulty cti2yaml / ctml2yaml conversion of reaction parameters, or a faulty initialization of a cantera object from YAML. The context points at an issue with surface reactions.

System information

Additional context

The behavior became apparent after 1866e35 which switched example input files to YAML.

bryanwweber commented 4 years ago

Spent parts of three days debugging this... man, this was tricky. The intermediate cause appears to be because the thermo phases are stored in a different order inside the kinetics object for YAML vs. CTI files. This means that the order of species in the internal arrays is different, which results in a different rates of production. This cause can be verified by changing the order of the ThermoPhase arrays in the _SolutionBase.__cinit__() method:

        if isinstance(self, Kinetics):
            for phase in adjacent:
                # adjacent bulk phases for a surface phase
                v.push_back(phase.thermo)
            v.push_back(self.thermo)  # Move this after the adjacent phases to fix the example

With that change, both surf_pfr.py and catalytic_combustor.py work again.

I don't believe this is the root cause, or the appropriate fix, because the underlying calculation of the rates of production is intended to be independent of the order that phases are added. More digging is required to find where this assumption is violated in the current code.

Debugging details follow, just in case anyone finds it interesting, and also because I'm so frustrated at how long this took that I want to write it all out!

First, I verified for surf_pfr.py that the mole fractions of several species were becoming negative. After that, I printed the rates of production of the species and compared between the YAML file and CTI file (it was very fortunate that the example still worked with the CTI file). That verified that the rates of production were indeed different between the mechanisms. I verified by hand that all the rate constant parameters had been converted properly from CTI to YAML and ruled out that as the cause. The next step was to use LLDB (I'm on a mac, gdb would work equivalently on Linux) to step into the C++ code for computing the rates. At that point, I noticed that the creation rates stored in the internal arrays had all the same numbers, they were just rearranged. When I looked at the documentation for the creation rates, I saw that it lists all the species of all the phases, and I guessed that the ThermoPhases stored within each of the surface kinetics objects were in different orders. I verified this by accessing the `Interface._phase_indices` property in Python, which showed that the `Surf` and `IdealGas` phases indeed had different internal indices. Actually, I found that the `_phase_indices` property was different between the two input files by writing a loop in Python to assert that each property of the `Interface` instance was equal. The last debugging step so far was to check the change I mentioned earlier, and put the `Surf` phase second, to match how the CTI file is initialized. Interestingly, even though the CTI initializer in Python has the same code that the YAML file does, the C++ code re-runs the loop, since the XML definition must list the associated phases, see: https://github.com/Cantera/cantera/blob/cf1c0816e7d535a1fc385063aebb8b8e93a85233/src/kinetics/importKinetics.cpp#L155-L160
ischoegl commented 4 years ago

:+1: whoa ... that indeed sounds like a tricky one to track down! It’s really surprising that this hadn’t come up thus far, and it was curious to find that switching things back to CTI appeared to ‘fix’ things ... anyhow, kudos.