Maintainer triage (June 29, 2026)
- Status: Verified bug against the current
master source in this fork.
- Severity: High.
- Area labels: area:kinetics, area:reactors, area:thermo.
- Tackle batch: P4 - high-impact numerical/model correctness.
The source review confirmed that the cited code path and failure mode are still present. The original audit details are preserved below for reproduction notes and suggested fixes.
Original audit details
Part of the model-correctness audit — see #17.
Severity: high | Category: conservation | Location: PharmaPy/Kinetics.py + PharmaPy/Reactors.py + PharmaPy/ThermoModule.py lines Kinetics.py:275-281, 593-614; Reactors.py:680-694; ThermoModule.py:258-282
What's wrong
The per-reaction rate returned by get_rxn_rates(overall_rates=False) is temp_terms*f_terms (line 609). Species production rates are formed as rate . normalized_stoich.T (line 611), where normalized_stoich = stoich_matrix.T / abs(ref_stoich) (line 281) and ref_stoich is the FIRST reactant coefficient of each reaction. Thus the reaction-extent variable that multiplies into the species balances is defined on a basis where the first reactant has coefficient 1. The energy balance (Reactors.py:694) computes source_term = -(deltah_rxn * rates).sum()*vol, but deltah_rxn comes from getHeatOfRxn fed the RAW (un-normalized) stoich_matrix (Reactors.py:651,681; ThermoModule.py:276 delta_cp = integral . stoich_matrix.T) and delta_hrxn reference is per-reaction-as-written. So the heat term mixes a raw-stoichiometry enthalpy with a first-reactant-normalized extent rate.
Why it's incorrect
For a reaction like '2A --> B', normalization sets ref_stoich=2, normalized nu_A=-1, so the material balance is dCA/dt=-r and r is the rate of A consumption. The thermodynamically consistent extent rate xi satisfies dCA/dt = nu_A_raw * xi = -2*xi, i.e. xi = r/2. The heat released per unit time is deltaH_as_written * xi = deltaH * r/2. The code instead computes deltaH * r, over-counting the heat by a factor equal to |ref_stoich| (=2 here) whenever the first reactant's stoichiometric coefficient differs from 1. This violates energy conservation: the heat generated is not consistent with the moles actually reacting.
Suggested regression test
Set up an adiabatic batch reactor with reaction '2A --> B' and a nonzero delta_hrxn, and the equivalent reaction written as 'A --> 0.5 B' with the same per-mole-of-A enthalpy. The two should give identical temperature trajectories. With the current code the '2A-->B' form releases twice the heat. Alternatively, divide deltah_rxn by abs(ref_stoich) (or feed getHeatOfRxn the normalized_stoich) and confirm the trajectories match.
Verification: 2/2 adversarial reviewers confirmed (refuted: 0).
Maintainer triage (June 29, 2026)
mastersource in this fork.The source review confirmed that the cited code path and failure mode are still present. The original audit details are preserved below for reproduction notes and suggested fixes.
Original audit details
Part of the model-correctness audit — see #17.
Severity: high | Category: conservation | Location:
PharmaPy/Kinetics.py + PharmaPy/Reactors.py + PharmaPy/ThermoModule.pylines Kinetics.py:275-281, 593-614; Reactors.py:680-694; ThermoModule.py:258-282What's wrong
The per-reaction rate returned by get_rxn_rates(overall_rates=False) is temp_terms*f_terms (line 609). Species production rates are formed as rate . normalized_stoich.T (line 611), where normalized_stoich = stoich_matrix.T / abs(ref_stoich) (line 281) and ref_stoich is the FIRST reactant coefficient of each reaction. Thus the reaction-extent variable that multiplies into the species balances is defined on a basis where the first reactant has coefficient 1. The energy balance (Reactors.py:694) computes source_term = -(deltah_rxn * rates).sum()*vol, but deltah_rxn comes from getHeatOfRxn fed the RAW (un-normalized) stoich_matrix (Reactors.py:651,681; ThermoModule.py:276 delta_cp = integral . stoich_matrix.T) and delta_hrxn reference is per-reaction-as-written. So the heat term mixes a raw-stoichiometry enthalpy with a first-reactant-normalized extent rate.
Why it's incorrect
For a reaction like '2A --> B', normalization sets ref_stoich=2, normalized nu_A=-1, so the material balance is dCA/dt=-r and r is the rate of A consumption. The thermodynamically consistent extent rate xi satisfies dCA/dt = nu_A_raw * xi = -2*xi, i.e. xi = r/2. The heat released per unit time is deltaH_as_written * xi = deltaH * r/2. The code instead computes deltaH * r, over-counting the heat by a factor equal to |ref_stoich| (=2 here) whenever the first reactant's stoichiometric coefficient differs from 1. This violates energy conservation: the heat generated is not consistent with the moles actually reacting.
Suggested regression test
Set up an adiabatic batch reactor with reaction '2A --> B' and a nonzero delta_hrxn, and the equivalent reaction written as 'A --> 0.5 B' with the same per-mole-of-A enthalpy. The two should give identical temperature trajectories. With the current code the '2A-->B' form releases twice the heat. Alternatively, divide deltah_rxn by abs(ref_stoich) (or feed getHeatOfRxn the normalized_stoich) and confirm the trajectories match.
Verification: 2/2 adversarial reviewers confirmed (refuted: 0).