Maintainer triage (June 29, 2026)
- Status: Verified bug against the current
master source in this fork.
- Severity: Critical.
- Area labels: area:drying.
- Tackle batch: P1 - critical blockers and severe unit errors.
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: critical | Category: units | Location: PharmaPy/Drying_Model.py lines 302-305
What's wrong
dsat_dt = -sum_dry / dens_liq / self.porosity, where sum_dry = dry_rate.sum(axis=1). dry_rate comes from get_drying_rate as k_ya_V(Δy) with k_y in mol/s/m^2 and a_V in m^2/m^3, so it is a MOLAR volumetric rate [mol/s/m^3]. dens_liq (self.rho_liq, line 250) is a MASS density [kg/m^3]. Dividing a molar rate by a mass density yields [mol/(kg·s)], not the [1/s] required for a saturation (volume-fraction) time derivative.
Why it's incorrect
Saturation S = liquid volume / pore volume. Its rate of change from evaporation is -(Σ_i N_i·MW_i)/(rho_liq·porosity), where N_i [mol/s/m^3] is the molar evaporation rate and MW_i [kg/mol] converts it to mass, then /rho_liq [kg/m^3] gives volume rate, then /porosity converts cake basis to pore basis. The code omits the MW_i factor entirely, so a molar quantity is divided by a mass-based density. The result is dimensionally inconsistent and numerically wrong by a factor of MW (order 0.02-0.1 kg/mol for typical solvents), badly mis-scaling the entire drying transient.
Suggested regression test
Multiply each component molar rate by its molecular weight before summing: sum_dry_mass = (dry_rate[:, idx_volatiles] * mw_volatiles).sum(axis=1); then dsat_dt = -sum_dry_mass/dens_liq/porosity. Confirm get_drying_rate truly returns mol/s/m^3 (k_y units) and check the units of self.rho_liq.
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: critical | Category: units | Location:
PharmaPy/Drying_Model.pylines 302-305What's wrong
dsat_dt = -sum_dry / dens_liq / self.porosity, where sum_dry = dry_rate.sum(axis=1). dry_rate comes from get_drying_rate as k_ya_V(Δy) with k_y in mol/s/m^2 and a_V in m^2/m^3, so it is a MOLAR volumetric rate [mol/s/m^3]. dens_liq (self.rho_liq, line 250) is a MASS density [kg/m^3]. Dividing a molar rate by a mass density yields [mol/(kg·s)], not the [1/s] required for a saturation (volume-fraction) time derivative.
Why it's incorrect
Saturation S = liquid volume / pore volume. Its rate of change from evaporation is -(Σ_i N_i·MW_i)/(rho_liq·porosity), where N_i [mol/s/m^3] is the molar evaporation rate and MW_i [kg/mol] converts it to mass, then /rho_liq [kg/m^3] gives volume rate, then /porosity converts cake basis to pore basis. The code omits the MW_i factor entirely, so a molar quantity is divided by a mass-based density. The result is dimensionally inconsistent and numerically wrong by a factor of MW (order 0.02-0.1 kg/mol for typical solvents), badly mis-scaling the entire drying transient.
Suggested regression test
Multiply each component molar rate by its molecular weight before summing: sum_dry_mass = (dry_rate[:, idx_volatiles] * mw_volatiles).sum(axis=1); then dsat_dt = -sum_dry_mass/dens_liq/porosity. Confirm get_drying_rate truly returns mol/s/m^3 (k_y units) and check the units of self.rho_liq.
Verification: 2/2 adversarial reviewers confirmed (refuted: 0).