From 859300d15145d2d5044b9b87fd443a90ac1ae646 Mon Sep 17 00:00:00 2001 From: parkyr Date: Thu, 9 Jul 2026 11:12:47 -0400 Subject: [PATCH 1/3] Fix reactor mode crash paths (#70) --- PharmaPy/Reactors.py | 41 +++++++++------- tests/test_reactor_modes.py | 98 +++++++++++++++++++++++++++++++++++++ 2 files changed, 121 insertions(+), 18 deletions(-) create mode 100644 tests/test_reactor_modes.py diff --git a/PharmaPy/Reactors.py b/PharmaPy/Reactors.py index d41bdb5c..7f41dcce 100644 --- a/PharmaPy/Reactors.py +++ b/PharmaPy/Reactors.py @@ -292,9 +292,11 @@ def _eval_state_events(self, time, states, sw): return events def heat_transfer(self, temp, temp_ht, vol): + """Return reactor heat transfer duty for supported heat-transfer modes.""" # Heat transfer area if self.ht_mode == 'coil': # Half pipe heat transfer - pass + raise NotImplementedError( + "heat_transfer with ht_mode='coil' is not supported") else: area_ht = 4 / self.diam * vol + self.area_base # m**2 heat_transf = self.u_ht * area_ht * (temp - temp_ht) @@ -1095,6 +1097,10 @@ def solve_unit(self, runtime=None, time_grid=None, eval_sens=False, check_modeling_objects(self) + if eval_sens: + raise NotImplementedError( + "CSTR sensitivity evaluation is not supported") + self.params_control = params_control self.set_names() @@ -1131,14 +1137,11 @@ def solve_unit(self, runtime=None, time_grid=None, eval_sens=False, # Create problem merged_params = self.Kinetics.concat_params() - if eval_sens: - pass - else: - def fobj(time, states): return self.unit_model( - time, states, merged_params) + def fobj(time, states): return self.unit_model( + time, states, merged_params) - problem = Explicit_Problem(fobj, states_init, - t0=self.elapsed_time) + problem = Explicit_Problem(fobj, states_init, + t0=self.elapsed_time) # Set solver solver = CVode(problem) @@ -1299,6 +1302,10 @@ def solve_unit(self, runtime=None, time_grid=None, eval_sens=False, check_modeling_objects(self) + if eval_sens: + raise NotImplementedError( + "SemibatchReactor sensitivity evaluation is not supported") + self.params_control = params_control self.set_names() @@ -1323,14 +1330,11 @@ def solve_unit(self, runtime=None, time_grid=None, eval_sens=False, states_init = np.append(states_init, tht_init) merged_params = self.Kinetics.concat_params() - if eval_sens: - pass - else: - def fobj(time, states): return self.unit_model( - time, states, merged_params) + def fobj(time, states): return self.unit_model( + time, states, merged_params) - problem = Explicit_Problem(fobj, states_init, - t0=self.elapsed_time) + problem = Explicit_Problem(fobj, states_init, + t0=self.elapsed_time) # Set solver solver = CVode(problem) @@ -1559,13 +1563,14 @@ def energy_steady(self, conc, temp): # source_term = -inner1d(deltah_rxn, rates) * 1000 # W/m**3 # TODO: Check if this is correct # source_term = -np.dot(deltah_rxn, rates) * 1000 # W / m**3 - source_term = -(deltah_rxn * rates).sum(axis=1) * 1000 # W / m**3 + source_term = -np.sum(deltah_rxn * rates) * 1000 # W / m**3 if self.adiabatic: heat_transfer = 0 else: # W/m**3 a_prime = self.diam / 4 # m**2 / m**3 - heat_transfer = self.u_ht * a_prime * (temp - self.Utility.temp) + temp_ht = self.Utility.evaluate_inputs(0)['temp_in'] + heat_transfer = self.u_ht * a_prime * (temp - temp_ht) flow_term = self.Inlet.vol_flow * cp_vol @@ -1601,7 +1606,7 @@ def solve_steady(self, vol_rxn, adiabatic=False): self.isothermal = False self.states_uo.append('temp') - c_inlet = self.Inlet.concentr + c_inlet = self.Inlet.mole_conc self.c_inert = c_inlet[~self.mask_species] c_partic = c_inlet[self.mask_species] diff --git a/tests/test_reactor_modes.py b/tests/test_reactor_modes.py new file mode 100644 index 00000000..eafbb850 --- /dev/null +++ b/tests/test_reactor_modes.py @@ -0,0 +1,98 @@ +# -*- coding: utf-8 -*- + +import copy +import importlib.util +import json + +import numpy as np +import pytest + + +HAS_ASSIMULO = importlib.util.find_spec("assimulo") is not None + +pytestmark = [ + pytest.mark.assimulo, + pytest.mark.integration, + pytest.mark.skipif( + not HAS_ASSIMULO, + reason="assimulo is not installed; solver-backed integration tests skipped", + ), +] + +if HAS_ASSIMULO: + from PharmaPy.Kinetics import RxnKinetics + from PharmaPy.Phases import LiquidPhase + from PharmaPy.Reactors import CSTR, BatchReactor, PlugFlowReactor + from PharmaPy.Reactors import SemibatchReactor + from PharmaPy.Streams import LiquidStream + from PharmaPy.Utilities import CoolingWater + + SENSITIVITY_REACTORS = [CSTR, SemibatchReactor] +else: + SENSITIVITY_REACTORS = [] + + +def _load_pfr_config(data_path): + with open(data_path["integration"] / "pfr_test_constructor_kwargs.json") as f: + config = json.load(f) + + config = copy.deepcopy(config) + tau = config["inlet"].pop("tau") + config["inlet"]["vol_flow"] = config["phase"]["vol"] / tau + + datapath = str(data_path["integration"] / "pfr_test_pure_comp.json") + config["kinetics"]["k_params"] *= 1 / 60 + config["kinetics"]["path"] = datapath + + return config, datapath + + +def _reactor_objects(data_path, reactor): + config, datapath = _load_pfr_config(data_path) + + inlet = LiquidStream(datapath, **config["inlet"]) + phase = LiquidPhase(datapath, **config["phase"]) + kinetics = RxnKinetics(**config["kinetics"]) + utility = CoolingWater(**config["utility"]) + + reactor.Inlet = inlet + reactor.Phases = phase + reactor.Kinetics = kinetics + reactor.Utility = utility + + return reactor + + +def test_pfr_solve_steady_reads_inlet_mole_conc(data_path): + config, _ = _load_pfr_config(data_path) + reactor = _reactor_objects( + data_path, PlugFlowReactor(**config["reactor"]) + ) + + vol_position, states = reactor.solve_steady(reactor.Liquid_1.vol) + + vol_position = np.asarray(vol_position) + + assert vol_position.size > 1 + assert states.shape[0] == vol_position.size + assert reactor.concProfSteady.shape[0] == vol_position.size + + +@pytest.mark.parametrize("reactor_cls", SENSITIVITY_REACTORS) +def test_sensitivity_mode_refuses_unsupported_reactors(data_path, reactor_cls): + if reactor_cls is CSTR: + reactor = reactor_cls() + else: + reactor = reactor_cls(vol_tank=0.002) + + reactor = _reactor_objects(data_path, reactor) + + with pytest.raises(NotImplementedError, match="sensitivity.*not supported"): + reactor.solve_unit(runtime=1, eval_sens=True, verbose=False) + + +def test_coil_ht_mode_refuses_unsupported_heat_transfer(): + reactor = BatchReactor(isothermal=False, ht_mode="coil") + + with pytest.raises(NotImplementedError, match="coil.*not supported"): + reactor.heat_transfer(np.array([300.0]), np.array([290.0]), 0.002) From 0f289984f5a2c6be5ea04445aee627100c2d3504 Mon Sep 17 00:00:00 2001 From: parkyr Date: Thu, 9 Jul 2026 14:47:11 -0400 Subject: [PATCH 2/3] Address reactor mode review feedback --- PharmaPy/Reactors.py | 38 +++++++++++++++++++++++++++---------- tests/test_reactor_modes.py | 24 ++++++++++++++++++++++- 2 files changed, 51 insertions(+), 11 deletions(-) diff --git a/PharmaPy/Reactors.py b/PharmaPy/Reactors.py index 7f41dcce..256eefce 100644 --- a/PharmaPy/Reactors.py +++ b/PharmaPy/Reactors.py @@ -949,7 +949,8 @@ class CSTR(_BaseReactor): whether or not the paramest_wrapper method should return the sensitivity system along with the concentratio profiles. Use False if you want the parameter estimation platform to - estimate the sensitivity system using finite differences + estimate the sensitivity system using finite differences. + Direct sensitivity evaluation is not implemented for CSTR. """ def __init__(self, mask_params=None, @@ -1099,11 +1100,17 @@ def solve_unit(self, runtime=None, time_grid=None, eval_sens=False, if eval_sens: raise NotImplementedError( - "CSTR sensitivity evaluation is not supported") + "CSTR sensitivity evaluation is not supported; construct " + "with return_sens=False to use finite-difference " + "sensitivities") self.params_control = params_control self.set_names() + if self.ht_mode == 'coil' and not self.isothermal: + raise NotImplementedError( + "CSTR heat transfer with ht_mode='coil' is not supported") + self.num_concentr = len(self.Liquid_1.mole_conc) self.args_inputs = (self, self.num_concentr, 0) @@ -1253,7 +1260,8 @@ class SemibatchReactor(CSTR): whether or not the paramest_wrapper method should return the sensitivity system along with the concentratio profiles. Use False if you want the parameter estimation platform to - estimate the sensitivity system using finite differences + estimate the sensitivity system using finite differences. + Direct sensitivity evaluation is not implemented for SemibatchReactor. """ def __init__(self, vol_tank, @@ -1304,11 +1312,18 @@ def solve_unit(self, runtime=None, time_grid=None, eval_sens=False, if eval_sens: raise NotImplementedError( - "SemibatchReactor sensitivity evaluation is not supported") + "SemibatchReactor sensitivity evaluation is not supported; " + "construct with return_sens=False to use finite-difference " + "sensitivities") self.params_control = params_control self.set_names() + if self.ht_mode == 'coil' and not self.isothermal: + raise NotImplementedError( + "SemibatchReactor heat transfer with ht_mode='coil' is " + "not supported") + if runtime is not None: final_time = runtime + self.elapsed_time @@ -1560,17 +1575,15 @@ def energy_steady(self, conc, temp): delta_hrxn=deltah_rxn) # ---------- Balance terms (W) - # source_term = -inner1d(deltah_rxn, rates) * 1000 # W/m**3 - # TODO: Check if this is correct - # source_term = -np.dot(deltah_rxn, rates) * 1000 # W / m**3 - source_term = -np.sum(deltah_rxn * rates) * 1000 # W / m**3 + source_term = -np.dot(deltah_rxn, rates) * 1000 # W / m**3 if self.adiabatic: heat_transfer = 0 else: # W/m**3 + # The steady-PFR area formula is a pre-existing issue tracked in #33. a_prime = self.diam / 4 # m**2 / m**3 - temp_ht = self.Utility.evaluate_inputs(0)['temp_in'] - heat_transfer = self.u_ht * a_prime * (temp - temp_ht) + heat_transfer = self.u_ht * a_prime * ( + temp - self.temp_ht_steady) flow_term = self.Inlet.vol_flow * cp_vol @@ -1618,6 +1631,11 @@ def solve_steady(self, vol_rxn, adiabatic=False): if 'temp' in self.states_uo: states_init = np.append(states_init, self.Inlet.temp) + if 'temp' in self.states_uo and not self.adiabatic: + # The steady solve integrates over volume, not time, so use the + # inlet utility condition at the start of the volume profile. + self.temp_ht_steady = self.Utility.evaluate_inputs(0)['temp_in'] + problem = Explicit_Problem(self.unit_steady, states_init, t0=0) solver = CVode(problem) diff --git a/tests/test_reactor_modes.py b/tests/test_reactor_modes.py index eafbb850..4ddef872 100644 --- a/tests/test_reactor_modes.py +++ b/tests/test_reactor_modes.py @@ -41,7 +41,12 @@ def _load_pfr_config(data_path): config["inlet"]["vol_flow"] = config["phase"]["vol"] / tau datapath = str(data_path["integration"] / "pfr_test_pure_comp.json") - config["kinetics"]["k_params"] *= 1 / 60 + config["kinetics"].update({ + "stoich_matrix": [[-1, -1, 1], [0, -1, 1]], + "k_params": [40 / 60, 10 / 60], + "ea_params": [2e3, 1e3], + "delta_hrxn": [-5e3, -2.5e3], + }) config["kinetics"]["path"] = datapath return config, datapath @@ -76,6 +81,9 @@ def test_pfr_solve_steady_reads_inlet_mole_conc(data_path): assert vol_position.size > 1 assert states.shape[0] == vol_position.size assert reactor.concProfSteady.shape[0] == vol_position.size + assert reactor.Kinetics.num_rxns == 2 + assert reactor.tempProfSteady[-1] > reactor.Inlet.temp + assert reactor.tempProfSteady[-1] < reactor.temp_ht_steady @pytest.mark.parametrize("reactor_cls", SENSITIVITY_REACTORS) @@ -96,3 +104,17 @@ def test_coil_ht_mode_refuses_unsupported_heat_transfer(): with pytest.raises(NotImplementedError, match="coil.*not supported"): reactor.heat_transfer(np.array([300.0]), np.array([290.0]), 0.002) + + +@pytest.mark.parametrize("reactor_cls", SENSITIVITY_REACTORS) +def test_coil_ht_mode_refuses_through_solve_unit(data_path, reactor_cls): + if reactor_cls is CSTR: + reactor = reactor_cls(isothermal=False, ht_mode="coil") + else: + reactor = reactor_cls( + vol_tank=0.002, isothermal=False, ht_mode="coil") + + reactor = _reactor_objects(data_path, reactor) + + with pytest.raises(NotImplementedError, match="coil.*not supported"): + reactor.solve_unit(runtime=1, verbose=False) From 6a0fb16dbfa42cebcc06adf6d42b4ee213c95478 Mon Sep 17 00:00:00 2001 From: parkyr Date: Thu, 6 Aug 2026 16:56:37 -0400 Subject: [PATCH 3/3] Fix steady-PFR heat-transfer area and document changed reactor methods Addresses the second review round on PR #107. Refs #33: PlugFlowReactor.energy_steady computed the specific heat-transfer area as `diam / 4`, the reciprocal of the `4 / diam` the transient balance uses. For a cylindrical tube the wetted area per unit volume is (pi*D*L)/(pi*D**2/4*L) = 4/D, so the old form was wrong by 16/D**2 -- a factor of ~24,800 for the test fixture's 0.0254 m tube. The line was unreachable before this PR fixed the crash in solve_steady, so this is the first release in which it would have produced numbers. Add test_steady_pfr_specific_area_matches_tube_geometry, which recovers a_prime from a single energy_steady call at zero reactant concentration (where the source term vanishes exactly) and pins it to 4/D. Reintroducing `diam / 4` makes it fail with 0.00635 against the expected 157.48, while every other test in the module stays green -- confirming the pre-existing temperature assertions never guarded this. Tighten those temperature assertions accordingly: with the corrected area the tube equilibrates with the utility, so the former `< temp_ht_steady` bound held by ~4e-12 K, inside the solver tolerance. Assert equilibration instead. Annotate units through the steady energy balance in [unit] form and correct cp_vol, which was labelled `W/K` but is `J/m**3/K`. Document the 1000 factors as the L -> m**3 conversions they are. Add Numpy-format docstrings to the methods this PR changed -- heat_transfer, energy_steady, and solve_steady -- covering argument types, units, defaults, returns, and the NotImplementedError raised for ht_mode='coil'. solve_steady now also documents that it overwrites the instance `adiabatic` attribute. Verified in the pharmapy .venv: 48 passed (was 47 collected); lanes "not assimulo" 33, assimulo 15, slow 4, integration 13. Co-Authored-By: Claude Opus 5 --- PharmaPy/Reactors.py | 89 +++++++++++++++++++++++++++++++++---- tests/test_reactor_modes.py | 53 +++++++++++++++++++++- 2 files changed, 132 insertions(+), 10 deletions(-) diff --git a/PharmaPy/Reactors.py b/PharmaPy/Reactors.py index 256eefce..ae7d98da 100644 --- a/PharmaPy/Reactors.py +++ b/PharmaPy/Reactors.py @@ -292,7 +292,29 @@ def _eval_state_events(self, time, states, sw): return events def heat_transfer(self, temp, temp_ht, vol): - """Return reactor heat transfer duty for supported heat-transfer modes.""" + """Return reactor heat transfer duty for supported heat-transfer modes. + + Parameters + ---------- + temp : float or numpy.ndarray + Reactor temperature [K]. + temp_ht : float or numpy.ndarray + Heat-transfer utility temperature [K]. + vol : float + Liquid volume in contact with the heat-transfer surface [m**3]. + + Returns + ------- + heat_transf : float or numpy.ndarray + Heat-transfer duty [W]. Positive when the reactor loses heat to + the utility. + + Raises + ------ + NotImplementedError + If ``ht_mode`` is 'coil', which is documented as an option but + has no implementation. + """ # Heat transfer area if self.ht_mode == 'coil': # Half pipe heat transfer raise NotImplementedError( @@ -1554,14 +1576,29 @@ def material_steady(self, conc, temp): return dconc_dv def energy_steady(self, conc, temp): + """Steady-state energy balance derivative along reactor volume. + + Parameters + ---------- + conc : numpy.ndarray + Molar concentrations of the participating species [mol/L]. + temp : float + Reactor temperature at the current volume coordinate [K]. + + Returns + ------- + dtemp_dv : float + Temperature derivative with respect to reactor volume [K/m**3]. + """ _, cp_j = self.Liquid_1.getCpPure(temp) concentr = np.zeros_like(self.Liquid_1.mole_conc) concentr[self.mask_species] = conc concentr[~self.mask_species] = self.c_inert - # Volumetric heat capacity - cp_vol = np.dot(cp_j, concentr) * 1000 # W/K + # Volumetric heat capacity. cp_j is [J/mol/K] and concentr is [mol/L], + # so the product is [J/L/K]; the 1000 is the L -> m**3 conversion. + cp_vol = np.dot(cp_j, concentr) * 1000 # [J/m**3/K] # Heat of reaction delta_href = self.Kinetics.delta_hrxn @@ -1574,20 +1611,28 @@ def energy_steady(self, conc, temp): rates = self.Kinetics.get_rxn_rates(conc, temp, overall_rates=False, delta_hrxn=deltah_rxn) - # ---------- Balance terms (W) - source_term = -np.dot(deltah_rxn, rates) * 1000 # W / m**3 + # ---------- Balance terms [W/m**3] + # deltah_rxn is [J/mol] and rates is [mol/L/s], so the product is + # [W/L]; the 1000 is the L -> m**3 conversion, matching the transient + # balance in energy_balances. Negative because an exothermic reaction + # carries deltah_rxn < 0 and releases heat. + source_term = -np.dot(deltah_rxn, rates) * 1000 # [W/m**3] if self.adiabatic: - heat_transfer = 0 - else: # W/m**3 - # The steady-PFR area formula is a pre-existing issue tracked in #33. - a_prime = self.diam / 4 # m**2 / m**3 + heat_transfer = 0 # [W/m**3] + else: + # Wetted area per unit reactor volume for a cylindrical tube: + # (pi*D*L) / (pi*D**2/4 * L) = 4/D. Refs #33. + a_prime = 4 / self.diam # [m**2/m**3] + # u_ht [W/m**2/K] * a_prime [m**2/m**3] * dT [K] -> [W/m**3] heat_transfer = self.u_ht * a_prime * ( temp - self.temp_ht_steady) + # vol_flow [m**3/s] * cp_vol [J/m**3/K] -> [W/K] flow_term = self.Inlet.vol_flow * cp_vol # -------- Energy balance + # [W/m**3] / [W/K] -> [K/m**3], integrated over reactor volume dtemp_dv = (source_term - heat_transfer) / flow_term return dtemp_dv @@ -1612,6 +1657,32 @@ def unit_steady(self, time, states, params=None): return deriv def solve_steady(self, vol_rxn, adiabatic=False): + """Integrate the steady-state PFR balances along reactor volume. + + Parameters + ---------- + vol_rxn : float + Reactor volume to integrate over [m**3]. The independent variable + of this solve is volume, not time. + adiabatic : bool (optional, default = False) + Whether to neglect wall heat transfer [-]. When False, the + utility inlet condition is sampled once into ``temp_ht_steady`` + and the energy balance takes the heat-transfer branch. + + Returns + ------- + volPosition : numpy.ndarray + Volume coordinates of the returned profile [m**3]. + states_solver : numpy.ndarray + Solution states at each volume coordinate: participating-species + molar concentrations [mol/L], followed by temperature [K] when + 'temp' is among the unit states. + + Notes + ----- + This method overwrites the instance ``adiabatic`` attribute with the + argument value. + """ self.adiabatic = adiabatic self.set_names() diff --git a/tests/test_reactor_modes.py b/tests/test_reactor_modes.py index 4ddef872..339603de 100644 --- a/tests/test_reactor_modes.py +++ b/tests/test_reactor_modes.py @@ -83,7 +83,58 @@ def test_pfr_solve_steady_reads_inlet_mole_conc(data_path): assert reactor.concProfSteady.shape[0] == vol_position.size assert reactor.Kinetics.num_rxns == 2 assert reactor.tempProfSteady[-1] > reactor.Inlet.temp - assert reactor.tempProfSteady[-1] < reactor.temp_ht_steady + + # With the correct 4/D specific area the 1 inch tube is strongly coupled + # to the utility, so the profile equilibrates to it rather than merely + # staying below it: a bare `< temp_ht_steady` bound holds here by ~4e-12 K, + # which is inside the solver tolerance and would pass for the wrong reason. + assert reactor.tempProfSteady[-1] == pytest.approx( + reactor.temp_ht_steady, abs=1e-6) + assert reactor.tempProfSteady[-1] <= reactor.temp_ht_steady + + +def test_steady_pfr_specific_area_matches_tube_geometry(data_path): + """Pin the steady-PFR specific heat-transfer area to 4/D (Refs #33). + + Probing at zero reactant concentration makes every reaction rate zero, so + the source term vanishes exactly and the steady energy balance reduces to + ``dT/dV = -u_ht * a_prime * (T - T_ht) / (vol_flow * cp_vol)``. That lets + ``a_prime`` [m**2/m**3] be recovered from one call and compared against + the tube geometry 4/D -- the same expression the transient balance uses. + The reciprocal form D/4 fails this by a factor of 16/D**2. + + Zeroing ``delta_hrxn`` would not work here: ``getHeatOfRxn`` applies a + heat-capacity correction between ``tref_hrxn`` and the probe temperature, + so the heat of reaction is nonzero even when the reference value is zero. + """ + config, datapath = _load_pfr_config(data_path) + + reactor = PlugFlowReactor(**config["reactor"]) + reactor.Inlet = LiquidStream(datapath, **config["inlet"]) + reactor.Phases = LiquidPhase(datapath, **config["phase"]) + reactor.Kinetics = RxnKinetics(**config["kinetics"]) + reactor.Utility = CoolingWater(**config["utility"]) + + reactor.solve_steady(reactor.Liquid_1.vol) + + temp_probe = 310.0 # [K], held away from the utility temperature + conc_probe = np.zeros_like(reactor.concProfSteady[0]) # [mol/L], no rates + + dtemp_dv = float(reactor.energy_steady(conc_probe, temp_probe)) # [K/m**3] + + concentr = np.zeros_like(reactor.Liquid_1.mole_conc) # [mol/L] + concentr[reactor.mask_species] = conc_probe + concentr[~reactor.mask_species] = reactor.c_inert + _, cp_j = reactor.Liquid_1.getCpPure(temp_probe) # [J/mol/K] + cp_vol = np.dot(cp_j, concentr) * 1000 # [J/m**3/K] + flow_term = reactor.Inlet.vol_flow * cp_vol # [W/K] + + a_prime = -dtemp_dv * flow_term / ( + reactor.u_ht * (temp_probe - reactor.temp_ht_steady)) # [m**2/m**3] + + assert a_prime == pytest.approx(4 / reactor.diam, rel=1e-8) + # Hand-computed for the fixture's 0.0254 m (1 inch) tube: 4/0.0254. + assert a_prime == pytest.approx(157.4803, rel=1e-4) @pytest.mark.parametrize("reactor_cls", SENSITIVITY_REACTORS)