GitLab

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