libantioch / antioch

C++ Chemical Kinetics, Thermodynaimics, and Transport Library
https://libantioch.github.io/
Other
22 stars 17 forks source link

Antioch fails to read gri30.xml #211

Closed roystgnr closed 6 years ago

roystgnr commented 8 years ago

We get an antioch_not_implemented() error in read_reaction_set_data.C, when the code can't decide what sort of reaction "falloff" is.

In the XML file, there's a <falloff type="Lindemann"/> tag, but our code seems to be assuming that we'll have specified that as type=LindemannFalloff instead.

SylvainPlessis commented 8 years ago

Antioch is much more explicit than the gri30.xml, indeed...

But I wonder why I haven't coded a gri-compatible parser. The fix is to make XMLParser<NumericType>::reaction_chemical_process() aware of this definition, and to make the XMLParser<NumericType>::Troe_*_parameter(...) compatible with the non-specific block definition.

I'm on it.

pbauman commented 8 years ago

We should have a more helpful error message too, e.g. write out what was found (if anything) and then list what is currently understood by Antioch.

roystgnr commented 8 years ago

Thanks, @SylvainPlessis ! But don't feel bad if this gets put on the back-burner - this is something I'd like to fix for other users' benefit, not something I need for myself. Changing the XML file to be Antioch-compatible was pretty trivial, and I'm actually going to be using a different reaction dataset in the end anyway.

pbauman commented 8 years ago

I should also note that, if you enjoy pain, I'm pretty sure the ChemKin GRI3 input file is parsed correctly (I'd tested this several months ago).

roystgnr commented 8 years ago

@pbauman, we get the "list what is currently understood by Antioch" part right; but the error message omits the "what was found" part, presumably because if you had verbose=true that part was already just printed.

SylvainPlessis commented 8 years ago

This is linked to issue #184, du to #172. This is where I explain that indeed, we were not GRI (Cantera) compatible for falloff. It was an answer to #163. Here's the whole history for it.

PR #212 is where I fixed it.

To sum it up, the GRI/Cantera way was not explicit/rigorous enough for me (and I guess I was lazy enough not to look for the order of the Troe parameters), so I used a more explicit description.

roystgnr commented 8 years ago

Looks like we're not quite there yet. #212 doesn't seem to handle three-body reactions?

SylvainPlessis commented 8 years ago

Now this one is unexpected...

The only supposed difference is threeBody (GRI) instead of ThreeBody. Which is supposedly taken care of at the line 163 of _src/parsing/src/read_reaction_setdata.C. There is even a comment to explicitely say it!

Should I start crying or there's something else?

roystgnr commented 8 years ago

When parsing reaction 0085, I now get

(gdb) p my_rxn->type()
$1 = Antioch::ReactionType::TROE_FALLOFF

Looks like the problem here is that there's again no explicit 'threeBody' or 'ThreeBody' keyword in the file. The reaction type="falloff", there's a falloff type="Troe", but the only clue that it should really be TROE_FALLOFF_THREE_BODY instead is that there's a "(+ M)" in the equation and there's an efficiencies section in the rateCoeff.

SylvainPlessis commented 8 years ago

Ah... Indeed I haven't even though about those falloff three-body reactions... So const std::string XMLParser<NumericType>::reaction_chemical_process() const should check for a (+ M) string in the equation and (or?) the presence of the efficiencies tag.

I suggest that the efficiencies tag presence is the test here. How are you feeling about this?

roystgnr commented 8 years ago

I'm feeling annoyed that they didn't put the metadata into the metadata.

My personal preference would be to check for \( *\+ *M *\) in the equation, use that to determine the reaction is three-body, but then continue to verify that only three-body reactions have efficiencies.

That way we're still correct in the case where a three-body equation is specified but the efficiencies tag was left out as a way of implying all efficiencies are 1.0.

I'd also move the verification out of an assert... parsing isn't in an inner loop, and these files are clearly ill-defined enough that we want to be extra careful when reading them.

SylvainPlessis commented 8 years ago

...these files are clearly ill-defined...

This one made me smile from ear to ear. Welcome into a chemicaler world :smile:

SylvainPlessis commented 8 years ago

Can I use #include <regex> on this one? That would greatly simplify the code.

roystgnr commented 8 years ago

I take it this means your vote is "yes, mandate at least C++11" for the question in #206?

Works for me, then.

SylvainPlessis commented 8 years ago

Indeed.
I'll go ahead with regex then. @pbauman, I'd like to have your approval on this one though.

SylvainPlessis commented 8 years ago

My personal preference would be to check for \( *\+ *M *\) in the equation, use that to determine the reaction is three-body, but then continue to verify that only three-body reactions have efficiencies.

I don't know about checking that there are no extra information in the description. This is typically so much more painful. As of now, there is no check of that sort (apart from the Arrhenius in place of Kooij stuff, because this is way too irritating).

That way we're still correct in the case where a three-body equation is specified but the efficiencies tag was left out as a way of implying all efficiencies are 1.0.

This is ensured by the ThreeBodyReaction and FalloffThreeBodyReaction constructors, they fill up the efficiencies array with 1., we should be safe here.

pbauman commented 7 years ago

Taken care of in #228.

roystgnr commented 7 years ago

And now I'm hitting another parsing failure, as soon as I add nitrous oxide into the mix:

Reaction "0185":
 eqn: N2O (+ M) [=] N2 + O (+ M)
 type: LindemannFalloffThreeBody
reversible: 1

   reactants: N2O:1,1, 
    N2O 1
   products: N2:1,1, O:1,1, 
    N2 1
    O 1
 rate: Arrhenius model
  High pressure limit rate constant
   A: 7.91e+10 s-1
   b: 0
   E: 56020 cal/mol
Assertion `3 == data.size()' failed.
3 = 3
data.size() = 4
SylvainPlessis commented 7 years ago

I think I found the reason. This is a Kooij rate constant labeled Arrhenius with a zero power parameter. So a very (very) badly described Arrhenius (unless you add uncertainties on the power parameter, in this case it's a classic Kooij, but I didn't consider this option as we do not describe uncertainties in the input files).

I think the problem lies between the zero test in xml_parser.C and read_reaction_set_data.C. They don't seem to agree.

If I'm correct, you should have no warning on this reaction, the expected (and uncorrect by the way, the warning should not be printed in an Arrhenius case) behavior would be having the warning of the second zero test.

The simplest way to verify this is to print the data vector juste before build_rate which is where the error is found. If the second number is zero, this is the faulty parameter that should not be there. Too many for Arrhenius because no power parameter is supposed to be there, too few for Kooij because no Tref. This could explain this schizophrenic reaction.

I can't verify it right now, cppunit v 1.13.1 does not agree with make check, and it's a little late to install a newer version, verify it works, and open an issue about it. I'll look deeper tomorrow.

Hope it helps

SylvainPlessis commented 7 years ago

Ok I've found the problem, it's not in the == 0 part, it's because of the double agent nature of this Arrhenius that is coded as a Kooij that is equivalent to an Arrhenius. For performance reason, Antioch takes it for an Arrhenius (we won't compute useless power of zero).

Now the thing is, the test for the power parameter is whether the parser finds a power parameter, which in this case will turn out to be true. The filtering of the Arrhenius badly written as a Kooij should operate at the == 0 test, same criteria as before, except if we are in a falloff, because who says the second rate constant won't have a non-null power parameter? So in this case it takes the useless power parameter, but when it comes to Tref, the test is about the kinetics model, and Arrhenius doesn't care about this. So we basically have a kinetics superposition of models, but Antioch isn't quantic...

The quick solution:

We will lose the power optimisation for the falloff in the case of both rate constants are the double agent Arrhenius, but it'll work.

Ok, I've posted the PR (#230), and I didn't post this one before... I though I did...

This is the post the PR refer to.

pbauman commented 6 years ago

I believe this is good. I'm hitting GRI3 XML much harder in #254 and it's parsing fine. So I'll go ahead and close.

Looks like this is related to #253 in @SylvainPlessis comment above.