Array indexing error via mechanism.ro2

[kilicomu]

Hi!

The mech_converter.py script writes its error messages into the files that it is operating on, e.g. the else clause in the following snippet:

https://github.com/AtChem/AtChem2/blob/37503aefb02bb54fa3b04899c68d15afa8b1c8c0/build/mech_converter.py#L450-L466

This results in lines such as the following appearing in the mechanism.ro2 file:

0 ! error RO2 not in mechanism: CHOCOCH2O2 DIETETO2

When mechanism.ro2 is read in through readRO2species from inputFunctions.f90, this 0 value is stored in the ro2Numbers array:

https://github.com/AtChem/AtChem2/blob/37503aefb02bb54fa3b04899c68d15afa8b1c8c0/src/inputFunctions.f90#L1245-L1276

Which propagates through to the ro2sum function from outputFunctions.f90:

https://github.com/AtChem/AtChem2/blob/37503aefb02bb54fa3b04899c68d15afa8b1c8c0/src/outputFunctions.f90#L25-L42

This can result in an attempt to reference element 0 of the array y, which raises an array indexing error at runtime and crashes the model uninformatively.

Is the error message going into the mechanism.ro2 file by design to cause the model to crash on erroneous mech_converter.py output? If so, I think a more informative error message where the crash would happen would help, e.g.

if ( min(ro2Numbers) < lbound(y, 1) ) error stop "Bad RO2 numbers - check mechanism.ro2 for errors!"

If it's not intentional, it's probably worth raising an exception in the mech_converter.py script at the time when the error is encountered so as to prevent the 0 ever getting into the ro2Numbers array and forcing the user to correct the issue before running the model.

I guess the other possibility is that it isn't a deal breaking error, in which case I think the mech_generator.py script should make a clear warning to the user about the missing RO2 species and the 0 RO2 number shouldn't get into the ro2Numbers array.

What do you think? Happy to implement resulting changes, just wanted to make sure I understand any reasoning behind the error messaging before opening a PR.

Killian

[rs028]

Hi Killian, as far as I can tell/remember (@spco can correct if I am wrong) there is no particular reason why it should behave this way, especially if it can lead to errors. My guess is that it is legacy code.

The most important thing here are: 1) that the calculation of RO2 is correct 2) that there is some way to tell the user that there is a potential error (either some species is identified as RO2 but it is not, or viceversa)

If there is a better way to do this than the current implementation, yes please, will be happy to hear about it!

[spco]

Hi @kilicomu . Thanks for your debugging on this (and I hope it didn't waste too much of your time tracking it down). You're absolutely right that it should do something better :)

Catching and treating it in mech_converter.py seems the most sensible option. My memory of the details is very hazy, but I think it should be sufficient to skip adding of the line to mechanism.ro2 - readRO2species counts ro2 species from that file, so that should be robust. The only thing perhaps you could check is that we don't ingest a count of ro2 species from the mechanism file anywhere, just in case that would need modifying. I don't think that's the case though.

As to whether to raise this as a warning or error in mech_converter.py, I don't have too strong an opinion. I would likely go with a warning, because we might want to reuse a mechanism file, cut out half the reactions, and not have to fiddle around modifying the ro2 list in case any of them have been lost from the mechanism. But I'm more than happy to be overruled on that point by the users!

I'm all for defensive programming, so an error in the Fortran code as you propose to guard against/explain a zero-indexing error would also be a good addition.

If you're able to raise that PR @kilicomu that would be awesome - I don't think I will have any time to work on this atm but happy to guide it through a PR.