Skip to content

Fix reactor mode crash paths #107

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
parkyr wants to merge 2 commits into
base: master
Choose a base branch
from
Open

Fix reactor mode crash paths #107

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
859300d
Fix reactor mode crash paths (#70)
parkyr Jul 9, 2026
0f28998
Address reactor mode review feedback
parkyr Jul 9, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
69 changes: 46 additions & 23 deletions PharmaPy/Reactors.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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."""

@bernalde bernalde Jul 10, 2026

Copy link
Copy Markdown

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

For this, and any upcoming PR, we need to have a docs template that also discusses the arguments, their type, units (in [unit] format), and possible default value. I propose to use the Numpy format. This would also help us with the future type hinting

# Heat transfer area
if self.ht_mode == 'coil': # Half pipe heat transfer
pass
raise NotImplementedError(

@parkyr parkyr Jul 9, 2026

Copy link
Copy Markdown
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Blocking — this refusal never reaches a CSTR/Semibatch user.

BatchReactor.solve_unit pre-calls the RHS at Reactors.py:833 (self.derivatives = call_fn(...)) before constructing the solver, so this NotImplementedError escapes with its message intact. CSTR.solve_unit and SemibatchReactor.solve_unit go straight from Explicit_Problem to CVode(problem).simulate(...), so the raise happens inside the SUNDIALS RHS callback and is converted to:

assimulo.solvers.sundials.CVodeError: 'The right-hand side function failed at the first call. At time 0.000000.'

I reproduced this on this head for both CSTR(isothermal=False, ht_mode='coil').solve_unit(runtime=1) and the SemibatchReactor equivalent. The word "coil" never appears. This is not a regression — the pre-PR UnboundLocalError was swallowed the same way — but it leaves acceptance criterion 4 of #70 unmet ("the refusal is raised at call time and names the unsupported mode, so paramest_wrapper and the non-isothermal energy balances fail loudly rather than at an unrelated line"), and #70 names CSTR (1060/1066) as a trigger for exactly this path.

Fix: mirror the eval_sens guard this PR already adds. At the top of CSTR.solve_unit and SemibatchReactor.solve_unit, after self.set_names() and before CVode is constructed:

if self.ht_mode == 'coil' and not self.isothermal:
    raise NotImplementedError(
        "CSTR heat transfer with ht_mode='coil' is not supported")

Scoping on not self.isothermal is safe: an isothermal CSTR never calls heat_transfer. unit_model skips energy_balances entirely, and the heat_prof=True branch at Reactors.py:1059-1061 takes ht_term = -(source_term + flow_term). Keep this raise as defense in depth.

@parkyr parkyr Jul 9, 2026

Copy link
Copy Markdown
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Addressed in 0f28998: CSTR.solve_unit and SemibatchReactor.solve_unit now raise call-time NotImplementedError for non-isothermal ht_mode='coil', before CVode can swallow it.

@bernalde bernalde Jul 10, 2026

Copy link
Copy Markdown

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Why do we even allow that ht_mode to be called? Is it implemented elsewhere? Should this become an issue to be implemented?

"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)
Expand Down Expand Up @@ -947,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,
Expand Down Expand Up @@ -1095,9 +1098,19 @@ def solve_unit(self, runtime=None, time_grid=None, eval_sens=False,

check_modeling_objects(self)

if eval_sens:
raise NotImplementedError(

@parkyr parkyr Jul 9, 2026

Copy link
Copy Markdown
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Nonblocking — this improves the error but still leaves no working path.

return_sens defaults to True for CSTR (Reactors.py:958) and SemibatchReactor (Reactors.py:1260), and paramest_wrapper branches on it (Reactors.py:476-479). So parameter estimation on a default-constructed CSTR now raises NotImplementedError instead of UnboundLocalError — clearer, but equally unusable.

Meanwhile paramest_wrapper's else branch (Reactors.py:490-496) already implements the finite-difference fallback the class docstring advertises: "Use False if you want the parameter estimation platform to estimate the sensitivity system using finite differences." I checked both on this head:

CSTR(return_sens=True).paramest_wrapper(...)   -> NotImplementedError
CSTR(return_sens=False).paramest_wrapper(...)  -> OK, c_prof shape (5, 3)

Flipping the CSTR/SemibatchReactor default to return_sens=False would make parameter estimation actually work rather than trade one exception for another, and it is a one-word change in each constructor.

If you would rather not change a public default in this PR, please at least name the escape hatch in the message and fix the docstring at Reactors.py:948-952, which currently implies the default is usable:

raise NotImplementedError(
    "CSTR sensitivity evaluation is not supported; construct with "
    "return_sens=False to use finite-difference sensitivities")

@parkyr parkyr Jul 9, 2026

Copy link
Copy Markdown
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Addressed the lower-risk option in 0f28998: the CSTR/Semibatch sensitivity errors and constructor docs now name return_sens=False as the finite-difference escape hatch. I left the public default unchanged in this PR.

"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)

Expand Down Expand Up @@ -1131,14 +1144,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)
Expand Down Expand Up @@ -1250,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,
Expand Down Expand Up @@ -1299,9 +1310,20 @@ 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; "
"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

Expand All @@ -1323,14 +1345,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)
Expand Down Expand Up @@ -1556,16 +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 = -(deltah_rxn * rates).sum(axis=1) * 1000 # W / m**3
source_term = -np.dot(deltah_rxn, rates) * 1000 # W / m**3

@bernalde bernalde Jul 10, 2026

Copy link
Copy Markdown

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Units formatting and why do we need the 1000 factor?


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

@bernalde bernalde Jul 10, 2026

Copy link
Copy Markdown

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

What is 4 here?

heat_transfer = self.u_ht * a_prime * (temp - self.Utility.temp)
heat_transfer = self.u_ht * a_prime * (
temp - self.temp_ht_steady)

flow_term = self.Inlet.vol_flow * cp_vol

Expand Down Expand Up @@ -1601,7 +1619,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]
Expand All @@ -1613,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)

Expand Down
120 changes: 120 additions & 0 deletions tests/test_reactor_modes.py
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,120 @@
# -*- 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"].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


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
Comment thread
bernalde marked this conversation as resolved.
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)
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)

@parkyr parkyr Jul 9, 2026

Copy link
Copy Markdown
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Blocking (pairs with my comment on PharmaPy/Reactors.py:298) — this asserts the refusal at the heat_transfer call boundary, which no real caller crosses directly.

Every user reaches heat_transfer through solve_unit, and for CSTR/SemibatchReactor the NotImplementedError is swallowed by CVode into a generic CVodeError. So this test passes today while the user-visible behavior still violates #70's acceptance criterion 4 — the test gives false confidence over exactly the gap it is meant to close.

Please add a solve_unit-level regression for the reactors whose solve_unit does not pre-call the RHS:

@pytest.mark.parametrize("reactor_cls", SENSITIVITY_REACTORS)
def test_coil_ht_mode_refuses_through_solve_unit(data_path, reactor_cls):
    kwargs = {} if reactor_cls is CSTR else {"vol_tank": 0.002}
    reactor = _reactor_objects(
        data_path, reactor_cls(isothermal=False, ht_mode="coil", **kwargs))

    with pytest.raises(NotImplementedError, match="coil.*not supported"):
        reactor.solve_unit(runtime=1, verbose=False)

That test is red on this head (it gets CVodeError) and green once the solve_unit guard is added. Keeping the current direct-call test alongside it is fine.

@parkyr parkyr Jul 9, 2026

Copy link
Copy Markdown
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Addressed in 0f28998: added a solve_unit-level parametrized regression for CSTR and Semibatch coil mode; it now passes with the new call-time guards.



@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)