Maintainer triage (June 29, 2026)
- Status: Verified bug against the current
master source in this fork.
- Severity: High.
- Area labels: area:thermo.
- Tackle batch: P2 - shared thermo, phase, and basis 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
Severity: High · Part of #67 (round-2 audit batch)
phi_ij numerator uses (M_i/M_j)^0.25 instead of the required (M_j/M_i)^0.25
Location: PharmaPy/ThermoModule.py:363-365 (ThermoPhysicalManager.getViscosityMix, vapor branch)
The Wilke (1950) interaction factor is phi_ij = [1 + (mu_i/mu_j)^0.5 * (M_j/M_i)^0.25]^2 / [sqrt(8)*(1 + M_i/M_j)^0.5]. Line 363 builds mw_term = np.outer(self.mw, 1/self.mw), so mw_term[i,j] = M_i/M_j. Line 365 uses mw_term0.25 in the numerator: (1 + visc_term0.5 * mw_term**0.25)*2 / np.sqrt(8(1+mw_term)), i.e. it uses (M_i/M_j)^0.25 instead of (M_j/M_i)^0.25. The denominator (1+M_i/M_j)^0.5 is correct. The downstream np.dot(mole_frac, phi_mix.T) uses the correct textbook orientation and does NOT cancel the inversion (verified numerically: applying (1/mw_term)**0.25 reproduces the exact textbook Wilke mixture viscosity to machine precision). For i=j phi_ii=1 either way, so the error is silent for self-terms and corrupts only cross terms.
Related existing issues (referenced, not modified): #18, #28, #45, #62, #73 (other thermo issues -- enthalpy latent-heat weighting, average MW from mass fractions, VaporPhase molar-density units, vectorized Keq Henry masking, VaporPhase constructor basis) -- reference only; none touch getViscosityMix or the Wilke MW exponent.
Found by a multi-agent source audit (round 2); confirmed by re-reading the source and cross-checked against all 80 existing open issues. Each defect was verified by an independent code-truth skeptic and a novelty skeptic.
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
Severity: High · Part of #67 (round-2 audit batch)
phi_ij numerator uses (M_i/M_j)^0.25 instead of the required (M_j/M_i)^0.25
Location:
PharmaPy/ThermoModule.py:363-365 (ThermoPhysicalManager.getViscosityMix, vapor branch)The Wilke (1950) interaction factor is phi_ij = [1 + (mu_i/mu_j)^0.5 * (M_j/M_i)^0.25]^2 / [sqrt(8)*(1 + M_i/M_j)^0.5]. Line 363 builds mw_term = np.outer(self.mw, 1/self.mw), so mw_term[i,j] = M_i/M_j. Line 365 uses mw_term0.25 in the numerator: (1 + visc_term0.5 * mw_term**0.25)*2 / np.sqrt(8(1+mw_term)), i.e. it uses (M_i/M_j)^0.25 instead of (M_j/M_i)^0.25. The denominator (1+M_i/M_j)^0.5 is correct. The downstream np.dot(mole_frac, phi_mix.T) uses the correct textbook orientation and does NOT cancel the inversion (verified numerically: applying (1/mw_term)**0.25 reproduces the exact textbook Wilke mixture viscosity to machine precision). For i=j phi_ii=1 either way, so the error is silent for self-terms and corrupts only cross terms.
Related existing issues (referenced, not modified): #18, #28, #45, #62, #73 (other thermo issues -- enthalpy latent-heat weighting, average MW from mass fractions, VaporPhase molar-density units, vectorized Keq Henry masking, VaporPhase constructor basis) -- reference only; none touch getViscosityMix or the Wilke MW exponent.
Found by a multi-agent source audit (round 2); confirmed by re-reading the source and cross-checked against all 80 existing open issues. Each defect was verified by an independent code-truth skeptic and a novelty skeptic.