Possible mistake in omega calculation in pyrolysisChemistryModel.C
Hi,
I was trying to use solidChemistry and pyrolysisChemistryModel in a solver, and I was getting weird results or runtime errors (undefined behavior) for a simple case where reactions took place in a solid medium! When I tracked down the runtime error, I found that the problem arises from the line 267 of pyrolysisChemistryModel.C, where we have
const scalar exp = R.lhs()[si].exponent;
where we get strange values for the "exp" variable! If you take a look at the whole block of this code:
const label Nl = R.lhs().size();
for (label s=0; s<Nl; s++)
{
label si = R.lhs()[s].index;
const scalar exp = R.lhs()[si].exponent; // This is line 267, the problematic line IMO
kf *=
pow(c1[si]/Ys0_[si][celli], exp)
*(Ys0_[si][celli]);
}
return kf;
I believe the "si" at line 267, should be replaced by "s", so the new code should be:
const scalar exp = R.lhs()[s].exponent;
I tried this and it solved the problem in my code! Because as far as I understand, "s" is the species index in the R.lhs() (left hand side of the reaction) and "si" is the species index in the species list! In case it is not clear or it is not convincing, please take a look at the following reaction file as an example:
species
(
A
B
C
);
gaseousSpecies
(
D
E
F
);
reactions
{
reaction1
{
type irreversibleArrheniusSolidReaction;
reaction "B = C";
A 1.5e14;
beta 0;
Ta 23653.7493709109;
Tcrit 500;
}
}
So in this case, to calculate the omega of the reaction, "s" will be 0 but "si" will be 1! Therefore, trying to read R.lhs()[si].exponent leads to undefined behavior, because R.lhs() has a size of 1!
I tried to explain the problem as clear as I could, so sorry for the long post, and please let me know if there is still something missing.
Best regards, Morteza