Broken extract_submechanism.py example

[ischoegl]

System information

• Cantera version: 2.5.0a4 (current master / commit c16c4b0)
• OS: Ubuntu
• Python/MATLAB version: Python 3.6.8

Expected behavior

The example kinetics/extract_submechanism.py runs.

Actual behavior

In [1]: %run extract_submechanism.py
Species: H2, H, O, O2, OH, H2O, HO2, H2O2, C, CO, CO2, N2
---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
/src/build/python/cantera/examples/kinetics/extract_submechanism.py in <module>()
     37 
     38 # Filter reactions, keeping only those that only involve the selected species
---> 39 all_reactions = ct.Reaction.listFromFile('gri30.yaml')
     40 reactions = []
     41 

interfaces/cython/cantera/reaction.pyx in cantera._cantera.Reaction.listFromFile()

ValueError: A Kinetics object is required.

To Reproduce

See above.

[ischoegl]

This bug will also affect a range of unit tests once xml input is migrated to yaml. A fix will require discussion ...

[sabasafdar]

hey, i'm interested in working on it

[bryanwweber]

@sabasafdar This may not be a very good "first issue" to resolve. It requires some pretty deep knowledge of the new input format. Feel free to try it out if you want, and if you get stuck feel free to ask questions here.

[ischoegl]

@sabasafdar ... From what I can tell, there are several ways to address this; the new YAML philosophy is somewhat different than the old XML, so some guidance from developers is certainly needed: a fix would potentially reach deep into the C++ sub layer.

[kyleniemeyer]

I'm going to try tackling this

[bryanwweber]

FWIW, this was introduced in 1866e35 which switched the input file to YAML. The cause of the error is that constructing a Reaction from a YAML entry requires a Kinetics object, to be able to tell whether or not the reaction is a surface reaction, and to get the units for the activity coefficients. See:

cantera/src/kinetics/Reaction.cpp

Lines 1071 to 1084 in cf1c081

	if (kin.thermo(kin.reactionPhaseIndex()).nDim() < 3) {
	// See if this is an electrochemical reaction
	Reaction testReaction(0);
	parseReactionEquation(testReaction, node["equation"], kin);
	if (isElectrochemicalReaction(testReaction, kin)) {
	unique_ptr<ElectrochemicalReaction> R(new ElectrochemicalReaction());
	setupElectrochemicalReaction(*R, node, kin);
	return unique_ptr<Reaction>(move(R));
	} else {
	unique_ptr<InterfaceReaction> R(new InterfaceReaction());
	setupInterfaceReaction(*R, node, kin);
	return unique_ptr<Reaction>(move(R));
	}
	}

and

cantera/src/kinetics/Reaction.cpp

Line 308 in cf1c081

Units rxn_phase_units = kin.thermo(kin.reactionPhaseIndex()).standardConcentrationUnits();

#845 formalized the idea that Kinetics should not be able to be instantiated as their own object, at least in Python. Instead, users should pass an instance of Solution to Reaction.listFromFile(), since Solution is an instance of Kinetics via inheritance, or an instance of Interface (I think) for the same reason, depending on the use case.

The simplest fix for the error here is to replace the offending line with

all_reactions = ct.Solution('gri30.yaml').reactions()

rather than using listFromFile(), since the return value of Solution.reactions() and listFromFile() will be identical in this case, and using listFromFile() would require the same Solution object.

[ischoegl]

@bryanwweber

The cause of the error is that constructing a Reaction from a YAML entry requires a Kinetics object, to be able to tell whether or not the reaction is a surface reaction, and to get the units for the activity coefficients.

I see this in a similar way, but am not sure that creating a temporary Kinetics object (as part of a Solution object in Python) is the only way to resolve this. Presumably, Reaction objects can exist without being part of a Kinetics object (at least in Python - see old behavior), so there should be a way to address the root cause of this bug. The main issue appears that Reaction.listFromFile is not implemented correctly for YAML, i.e. it should be either fixed in the main code (not just in the example), or be eliminated from the YAML interface in favor of your proposed alternative.

this was introduced in 1866e35 which switched the input file to YAML.

... the switch of input file formats made a pre-existing bug apparent.

@kyleniemeyer ... what's your opinion on this?

[bryanwweber]

@ischoegl I will let @speth jump in as necessary, but my understanding is that the YAML format is a little more flexible in terms of the units, such that it is necessary to determine the units of the activity concentration when the reaction is created (see the definition of readArrhenius:

cantera/src/kinetics/Reaction.cpp

Line 332 in cf1c081

Arrhenius readArrhenius(const Reaction& R, const AnyValue& rate,

). I don't think the requirement for a Kinetics object is a bug, but a change in behavior that changes how this function must be used.

That said, there is a case where this function could be useful, which is if you have a Solution generated from one input file, and want to get reactions from another input file that use the same species definitions. In this case, you could use Reaction.listFromFIle and pass the second input file, plus the instance of Solution from the first input file.

Nonetheless, for this particular example, since the reactions and species are in the same file, the fix I had earlier seems like the best bet. I'll push a branch shortly, once I see if any other examples need fixes (aside from those in #876 😄 )

[ischoegl]

@bryanwweber ... fwiw, it would be nice if the Kinetics object would be hidden from the Python interface (E.g. temporarily allocate it in C++). That would fix the root cause of the bug in Reactions.ListFromFile? I.e. no examples or other user code would need fixing.

PS: the name of the method should likely be deprecated in favor of list_from_file as the old name is not pythonic ...

[bryanwweber]

@ischoegl

temporarily allocate it in C++

Which ThermoPhase definition should be used for this allocation? What if the input file has no ThermoPhase definition, only a reactions section? How do we know if it should be bulk-phase Kinetics or InterfaceKinetics? Those are rhetorical questions (so no need to answer), and I'm pointing them out just to say that changing the behavior of Reaction.listFromFile() will not really be a simple change, and will not necessarily be any better than the current case, just different. My preference is to fix the example and defer other updates to the Reaction.listFromFile() interface until after 2.5.0 is released. People have to opt in to the new behavior by updating to use YAML files, so it should not come up that often, and I think the new behavior is correct (i.e., not a bug). In case it does come up before other changes are made, we can refer people to the workaround that I mentioned (create a Solution, pass that to Reaction.listFromFile()) as needed.

[ischoegl]

@bryanwweber ... I'm not sure that we're talking about the same thing. Let's use a unit test as an example:

cantera/test/kinetics/kineticsFromYaml.cpp

Lines 10 to 18 in cf1c081

	TEST(Reaction, ElementaryFromYaml)
	{
	auto sol = newSolution("gri30.yaml");
	AnyMap rxn = AnyMap::fromYamlString(
	"{equation: N + NO <=> N2 + O,"
	" rate-constant: [-2.70000E+13 cm^3/mol/s, 0, 355 cal/mol],"
	" negative-A: true}");
	auto R = newReaction(rxn, *(sol->kinetics()));

The newReaction call implies that instantiation of Reaction objects via YAML requires Kinetics, and the latter are associated with a ThermoPhase via the member variable m_thermo. So, at the face of it, reaction definitions (see YAML input) can exist by themselves, whereas underlying objects need context. I.e. YAML files without phase definitions cannot be used to create valid Reaction objects.

I'm still on-board with your approach to load the entire solution, but would argue that doing this within Reaction.list_from_file at the C++ layer (i.e. use the C++ version of newSolution) would be much preferable as it fixes the actual root cause. As it stands, the YAML version of the listFromFile method is broken and should imho be fixed prior to 2.5.

PS: the change of behavior (from XML input) is summarized here ...

cantera/include/cantera/kinetics/Reaction.h

Lines 264 to 271 in cf1c081

	//! Create a new Reaction object for the reaction defined in `rxn_node`
	//!
	//! @deprecated The XML input format is deprecated and will be removed in
	//! Cantera 3.0.
	shared_ptr<Reaction> newReaction(const XML_Node& rxn_node);
	//! Create a new Reaction object using the specified parameters
	unique_ptr<Reaction> newReaction(const AnyMap& rxn_node, const Kinetics& kin);

i.e. the XML version didn't need Kinetics, whereas YAML does.

[ischoegl]

@bryanwweber ... on further inspection, a simple fix including the original call would be

gas = ct.Solution('gri30.yaml')
all_reactions = ct.Reaction.listFromFile('gri30.yaml', gas)

I still think that it would be an intuitive implementation if Reaction.listFromFile defaulted to the kinetics object found in the file if gas (or the Kinetics portion thereof) is None.

PS: After looking at this some more, there's probably too much guesswork in creating a temporary Kinetics object. So I've come around: fixing the example(s) is fine, but I believe it would help if the signature of Reaction.listFromFile would be clarified as

    @staticmethod
    def listFromFile(filename, _SolutionBase sol=None, section='reactions'):

and the error thrown would refer to Solution instead of Kinetics (and potentially provide more insights on the changed behavior due to the transition from CTI to YAML).

[bryanwweber]

@ischoegl Just to clarify a few things, recognizing that we're in agreement 😄 :

YAML files without phase definitions cannot be used to create valid Reaction objects.

They can be used to create Reaction objects, as long as you have a valid Kinetics object constructed somehow, e.g., from a different input file.

I believe it would help if the signature of Reaction.listFromFile would be clarified as...

I agree, except InterfaceKinetics are also a valid option, so a Solution object is not the only choice to pass there (i.e., one could pass an instance of Interface). So some clarification is probably warranted, but I really want to defer that discussion until after 2.5.0 (which is very, very close! We are burning down the last few critical issues, and this doesn't seem critical, regardless of how small the change is).

on further inspection, a simple fix including the original call would be

I don't see a reason to include the original call. It doesn't demonstrate a use case for Reaction.listFromFile that is all that useful, since Solution.reactions() gives the equivalent list in one fewer line, and without the duplicate processing. It would be good if we had an example that did show a good use case for this function, but I don't think this is it. Anyways, we can discuss this more if necessary on my forthcoming PR.

[ischoegl]

Ad 2.5.0 - I realize that it's close ... when is feature freeze / beta scheduled? There are a couple of PR's that appear to be close but haven't been updated in a while.

Myself, I don't have any plans for further PR's anytime soon (am planning to wait until the counter goes up to at least #1k)

Regarding listFromFile ... I don't like the new behavior whatsoever as it is neither intuitive nor elegant. I guess this is related to #747 and #750, where autodetection of reaction types was a clear preference. There may be ways around this to have things both ways, but none are simple.

[ischoegl]

One last comment ... InterfaceKinetics is already included in my signature change suggestion above (as part of Interface), but needs to be reflected in the error message. Fwiw, Interface and other _SolutionBase-derived objects should likely also be included in the error messages of #845.