Units-and-basis correctness strategy: canonical registry, boundary validation, and balance tests (not Pint in the core)

[bernalde]

Goal

Establish a repository-wide units-and-basis correctness strategy that prevents the largest recurring defect class in this codebase — dimensional and physical-basis errors — without paying runtime cost in the ODE/DAE hot path. This is proposed infrastructure under #3, and it ties together the two refactors already planned in #118 (StateLayout) and #12 (UnitOperation base + constants module). It is a design/architecture issue; individual bug fixes remain tracked under epic #17.

Motivation

Dimensional and basis mistakes (mass vs molar vs volume, mass-fraction vs mole-fraction, kg vs g, mol/m³ vs kmol/m³, latent-heat basis) are by far the most common and most severe verified defects in the audit. A non-exhaustive cluster:

• Mass/molar-basis mixing: Critical: [Phases] VaporPhase mass-basis enthalpy weights latent heat by mole fraction #18, Critical: [Drying] Liquid saturation balance mixes molar drying rate with mass density #20, Critical: [Drying] Gas species balance adds molar source terms to mass-fraction states #21, High: [Drying] Energy terms multiply molar drying rates by mass-basis heat quantities #48, High: [Drying] Condensed-phase energy balance doubles latent heat #24, High: [Crystallizers] MSMPR concentration balance uses solid density in the liquid transfer term #46, High: [Containers/Phases] DynamicCollector divides mass-basis enthalpy by mole-basis Cp #84, High: [Phases] Multi-temperature VaporPhase mass-basis latent heat drops component columns #101
• Fraction-type confusion (mole vs mass fraction, both dimensionless [-]): High: [Drying] Gas molecular weight is computed from mass fractions as mole fractions #28, Medium: [Distillation] Outlet streams pass mole fractions as mole_conc #58, Medium: [Phases] SolidPhase mass-fraction distribution path is unreachable and dimensionally wrong #64, High: [Crystallizers] Mass-fraction conversion is dead code in Batch/Semibatch crystal balances #47, Medium: [Drying] Volatile liquid mass-fraction column is forced to zero each RHS evaluation #42
• Missing/duplicated unit conversions: High: [Phases] SolidPhase moles use stale mass and miss kg-to-g conversion #50 (kg→g), Low: [Phases] VaporPhase.getDensity returns mol/m3 while molar density convention is kmol/m3 #45 (mol/m³ vs kmol/m³), Medium: [Drying] Gas energy balance multiplies convective enthalpy by density twice #37, High: [Containers] Mixer slurry energy balance weights liquid enthalpy by total slurry Cp #25, Medium: [Containers] Mixer energy balance uses inconsistent inlet and outlet enthalpy references #38
• Index/basis handling in conversions: High: [ThermoModule] Solvent index 0 is treated as missing in concentration conversions #57, High: [Phases] VaporPhase constructor maps mole_conc into LiquidPhase mass_conc #73

These are not one-off slips; they recur because the physical basis of a quantity is carried only by convention and re-derived by hand at every call site. The same root cause drives #118: state vectors rely on memorized positional order.

Why not simply adopt Pint everywhere (the obvious first idea): Pint cannot live in the numerical core, and it would not even catch the worst bugs here:

1. The integrator boundary strips units. Every dynamic unit op exposes a unit_model(self, time, states, ...) RHS/residual (e.g. Reactors.py:369, Drying_Model.py:283, Evaporators.py:719) that Assimulo (CVode/IDA) and SciPy call thousands of times per solve and require plain float64 ndarrays. Any Pint quantity would have to be unwrapped before the call and rewrapped after — so units can never actually flow through the hot path.
2. State is a heterogeneous-unit flat vector. balances = np.concatenate((material, energy)) plus index-slice unpacking (NameAnalysis.py, and the dryer example in State vectors rely on memorized positional order across unit operations — introduce a single source of truth (StateLayout) #118) put concentrations, temperature, and holdup in one array. Pint arrays are single-unit; they cannot represent that vector.
3. Performance. Pint-wrapped array arithmetic is 10–100× slower per operation — unacceptable inside an RHS.
4. Pint's blind spot is our worst bug. mol/s vs kg/s differ in dimension (Pint catches them), but mole-fraction vs mass-fraction are both [-] and Pint is silent — yet that confusion drives High: [Drying] Gas molecular weight is computed from mass fractions as mole fractions #28, Medium: [Distillation] Outlet streams pass mole fractions as mole_conc #58, Medium: [Phases] SolidPhase mass-fraction distribution path is unreachable and dimensionally wrong #64, High: [Crystallizers] Mass-fraction conversion is dead code in Batch/Semibatch crystal balances #47, Medium: [Drying] Volatile liquid mass-fraction column is forced to zero each RHS evaluation #42. Pint alone would give false confidence on exactly the hardest cases.

Conclusion: the fix is a canonical-basis discipline enforced at boundaries and in tests, not units-in-the-core. Pint is at most an optional backend for the boundary validator, where it is cheap.

Proposed strategy (three layers)

1. Canonical unit + basis registry (single source of truth)

Define one authoritative unit and basis per physical quantity / state variable in the shared constants module introduced by #12 (e.g. PharmaPy/constants.py or PharmaPy/units.py). Every module imports canonical symbols instead of restating conventions or hard-coding factors like gas_ct = 8.314 or ad-hoc kg→g conversions.

• Enumerate every state variable and model quantity with its canonical SI(-derived) unit and basis (mass/molar/volume/number).
• Centralize physical constants and conversion factors (source, units, basis, valid range documented per AGENTS.md).
• Replace scattered magic factors (e.g. the kg↔g handling behind High: [Phases] SolidPhase moles use stale mass and miss kg-to-g conversion #50, the density convention behind Low: [Phases] VaporPhase.getDensity returns mol/m3 while molar density convention is kmol/m3 #45) with named references to the registry.

2. Validate and normalize once, at the public boundary

Do all unit/basis validation and conversion at construction / ingestion — Phase/Stream setup, inlet-condition setters, property-data loading — then work in the canonical internal representation as plain arrays. Inner functions and RHS code assume the invariant; no interior re-validation.

• Add boundary validators at phase/stream constructors and inlet setters that check units/basis and convert to canonical form, raising specific, actionable exceptions on mismatch.
• Explicitly validate finite-choice options (activity models, modes) so an unknown value never silently selects an unintended model (relates to High: [Evaporators] Activity model is passed positionally and VLE falls back to ideal gamma #27).
• Optional Pint pilot: back one boundary validator with Pint to parse/convert user-supplied units, since ingestion is not hot. Evaluate before spreading.

3. Named state layout + balance-based units tests

• Adopt the StateLayout single-source-of-truth from State vectors rely on memorized positional order across unit operations — introduce a single source of truth (StateLayout) #118 so packing/unpacking no longer relies on positional memory (this is where most basis-at-the-boundary bugs originate).
• Add balance-based units tests in the suite from Expand test coverage across modules (unit + integration) #7: assert conservation on asymmetric fixtures (mass-in = mass-out, mole/mass consistency, energy closure). These catch basis errors — including the fraction-type confusion Pint cannot — at zero runtime cost, because they assert a known physical balance rather than a dimension.
• Give each such test module a docstring stating scope and basis assumptions (per AGENTS.md testing section).

Acceptance criteria

• A documented canonical unit+basis registry exists and is imported by at least one unit-operation family end to end.
• At least one unit operation (pilot) validates/normalizes units at its public boundary and carries no interior re-validation.
• Balance-based units tests exist for the pilot op and fail if a basis error is reintroduced (verified by mutation, per AGENTS.md).
• A short docs/ note (or AGENTS.md cross-reference) records the canonical-basis convention and the "no Pint in the core" decision with its rationale.

Suggested sequencing

1. Fix the confirmed basis defects first (P1/P2 in [Epic] Verified model and numerical correctness defects (#18-#66) #17: Critical: [Phases] VaporPhase mass-basis enthalpy weights latent heat by mole fraction #18, Critical: [Drying] Liquid saturation balance mixes molar drying rate with mass density #20, Critical: [Drying] Gas species balance adds molar source terms to mass-fraction states #21, High: [Drying] Gas molecular weight is computed from mass fractions as mole fractions #28, High: [Drying] Condensed-phase energy balance doubles latent heat #24, High: [Drying] Energy terms multiply molar drying rates by mass-basis heat quantities #48, High: [Phases] SolidPhase moles use stale mass and miss kg-to-g conversion #50, High: [ThermoModule] Solvent index 0 is treated as missing in concentration conversions #57, Medium: [Distillation] Outlet streams pass mole fractions as mole_conc #58, Medium: [Phases] SolidPhase mass-fraction distribution path is unreachable and dimensionally wrong #64, …) so the pilot is built on correct balances, not on top of live bugs.
2. Land the infrastructure this depends on: Introduce a shared UnitOperation base class and constants module #12 (constants module → registry home) and State vectors rely on memorized positional order across unit operations — introduce a single source of truth (StateLayout) #118 (StateLayout → named packing).
3. Pilot on one unit operation (recommend a reactor or the dryer, since the dryer already carries the densest basis-bug history) to prove payoff before spreading.
4. Roll out per unit-operation family, adding balance tests under Expand test coverage across modules (unit + integration) #7 as each is converted.

Relationships

• Parent roadmap: [Epic] Restore CI, raise the SE baseline, and extend modeling (bioreactors, crystallizers, scheduling) #3
• Depends on / complements: Introduce a shared UnitOperation base class and constants module #12 (constants + base class), State vectors rely on memorized positional order across unit operations — introduce a single source of truth (StateLayout) #118 (StateLayout), Expand test coverage across modules (unit + integration) #7 (unit-test coverage)
• Motivated by (prevents this defect class): epic [Epic] Verified model and numerical correctness defects (#18-#66) #17 and the basis cluster Critical: [Phases] VaporPhase mass-basis enthalpy weights latent heat by mole fraction #18, Critical: [Drying] Liquid saturation balance mixes molar drying rate with mass density #20, Critical: [Drying] Gas species balance adds molar source terms to mass-fraction states #21, High: [Drying] Condensed-phase energy balance doubles latent heat #24, High: [Drying] Gas molecular weight is computed from mass fractions as mole fractions #28, Medium: [Drying] Gas energy balance multiplies convective enthalpy by density twice #37, Medium: [Containers] Mixer energy balance uses inconsistent inlet and outlet enthalpy references #38, Medium: [Drying] Volatile liquid mass-fraction column is forced to zero each RHS evaluation #42, Low: [Phases] VaporPhase.getDensity returns mol/m3 while molar density convention is kmol/m3 #45, High: [Crystallizers] MSMPR concentration balance uses solid density in the liquid transfer term #46, High: [Crystallizers] Mass-fraction conversion is dead code in Batch/Semibatch crystal balances #47, High: [Drying] Energy terms multiply molar drying rates by mass-basis heat quantities #48, High: [Phases] SolidPhase moles use stale mass and miss kg-to-g conversion #50, High: [ThermoModule] Solvent index 0 is treated as missing in concentration conversions #57, Medium: [Distillation] Outlet streams pass mole fractions as mole_conc #58, Medium: [Phases] SolidPhase mass-fraction distribution path is unreachable and dimensionally wrong #64, High: [Phases] VaporPhase constructor maps mole_conc into LiquidPhase mass_conc #73, High: [Containers/Phases] DynamicCollector divides mass-basis enthalpy by mole-basis Cp #84, High: [Phases] Multi-temperature VaporPhase mass-basis latent heat drops component columns #101
• Related documentation-units follow-ups: Audit SimExec documentation and units #129, Document SimExec labor cost assumptions and references #127, Add references for SimExec cost constants #128

Non-goals

• Not adopting Pint (or any units library) inside unit_model/balance RHS code or the solver interface.
• Not changing model equations or numerical results — this is enforcement infrastructure; behavior must be preserved except where a linked bug fix intentionally corrects it.
• Not a mass rename of variables; unit comments remain source annotations per AGENTS.md.