Duplicate reactant/product species not accumulated correctly in generic RHS

Hi,

I think there may be a generic issue in the new network RHS/Jacobian machinery when the same species appears both as a reactant and a product in a reaction.

I encountered this while independently reconstructing and validating the current aprox21 network.

One of the effective aprox21 reactions is equivalent to

He3 + He4 + H1 -> 2 He4

so He4 appears on both sides of the reaction. Its net stoichiometric coefficient in the forward direction should therefore be

-1 + 2 = +1.

However, in networks/rhs.H, the contribution of a rate to a given species is selected using a mutually exclusive if / else if chain, schematically:

if constexpr (data.species_A == species) {
    ...
}
else if constexpr (data.species_B == species) {
    ...
}
else if constexpr (data.species_C == species) {
    ...
}
else if constexpr (data.species_D == species) {
    ...
}
...

This seems to mean that if the same species occurs in more than one slot, only the first matching role contributes. In particular, a species that appears both among the reactants and products cannot have both stoichiometric contributions accumulated.

I found this independently through an analytic Jacobian comparison. I generated an aprox21 implementation that explicitly sums all reactant and product stoichiometric contributions and compared its species Jacobian against Microphysics with screening disabled.

Out of the 21 x 21 = 441 species Jacobian entries, 440 agree to about 1e-8 relative or better. The only significant discrepancy is

J(He4, He3)

associated with the effective reaction above. The Microphysics value differs by approximately a factor of two from the value obtained by directly summing the reaction stoichiometry.

The rest of the aprox21 instantaneous species RHS has also been checked against Microphysics at several states and agrees essentially to roundoff, so this appears to be localized to the duplicate-species handling rather than a general difference in the reconstructed network.

A possible generic fix would be to accumulate all matching roles independently rather than selecting exactly one role, conceptually something like

int nu = 0;

if (data.species_A == species) nu -= data.number_A;
if (data.species_B == species) nu -= data.number_B;
if (data.species_C == species) nu -= data.number_C;

if (data.species_D == species) nu += data.number_D;
if (data.species_E == species) nu += data.number_E;
if (data.species_F == species) nu += data.number_F;

with the forward and reverse contributions constructed consistently from the resulting stoichiometric coefficients.

This would also make the generic implementation robust to other reactions where a species appears in multiple reaction slots, not only this aprox21 case.

I tested against current Microphysics main around commit

3ac0bda4f00da3eed0446beadfc1cd5c51b99a40

I would be happy to provide the comparison program / numerical output if useful.

Thanks.

Start in networks/rhs.H and trace how each reaction species slot contributes to the RHS and Jacobian. Use the aprox21 He3 + He4 + H1 -> 2 He4 case and the reported J(He4, He3) discrepancy as the regression scenario. Done means all matching reactant and product roles are accumulated and the RHS and species Jacobian agree with Microphysics for this case.