Fix reactor mode crash paths#107
Conversation
parkyr
left a comment
parkyr
left a comment
There was a problem hiding this comment.
Maintainer review — head 859300d
Reviewed against origin/master @ 1a3ab31. The base has not moved since the branch point (git merge-base origin/master HEAD equals the origin/master tip), so the PR-head tree and the merge-result tree are identical and no separate merge worktree was needed.
All three defects in #70 reproduce on the unchanged base, and all three fixes are real and load-bearing — I checked that the two steady-PFR energy-balance changes are not incidental by applying only the mole_conc one-liner to the base: the PFR test still fails there with CVodeError: 'The right-hand side function failed at the first call.', because CVode swallows the AxisError from .sum(axis=1). That also corrects the triage note on #70, which claimed no further defect was hiding behind the AttributeError; there were two.
That same CVode-swallowing behavior is what the one blocking finding is about.
Blocking
- The coil refusal never reaches a CSTR/Semibatch user. Inline on
PharmaPy/Reactors.py:298.BatchReactor.solve_unithappens to pre-call the RHS atReactors.py:833, so theNotImplementedErrorescapes;CSTR.solve_unitandSemibatchReactor.solve_unitgo straight intoCVode(problem).simulate(...), so the raise is converted to a genericCVodeErrorand the word "coil" never reaches the user. Not a regression (the pre-PRUnboundLocalErrorwas swallowed identically), but #70 acceptance criterion 4 is unmet, and #70 names CSTR as a trigger for exactly this path. - The coil test cannot catch that gap. Inline on
tests/test_reactor_modes.py:98. It asserts at theheat_transfercall boundary, which no real caller crosses directly.
Nonblocking
Four inline: the return_sens=True default leaving CSTR/Semibatch parameter estimation with no working path (Reactors.py:1101); np.sum vs the equivalent-but-shape-strict np.dot plus the stale TODO/commented alternatives (Reactors.py:1566); the evaluate_inputs(0) accessor divergence and per-RHS dict rebuild (Reactors.py:1572); shape-only assertions over a delta_hrxn = 0 fixture (tests/test_reactor_modes.py:78). Plus two here:
Closes #70 is premature while criterion 4 is unmet
Once the solve_unit guard from blocking finding 1 lands, Closes #70 is correct as written. If you would rather land that guard as a follow-up, change the link to Refs #70 so the issue is not auto-closed with a documented acceptance criterion still failing.
a_prime = self.diam / 4 is now on a live code path
Reactors.py:1571 uses self.diam / 4 where the dynamic balance at Reactors.py:1691 uses 4 / self.diam. #70 explicitly carves this out as #33 and not a duplicate, so it is correctly out of scope here — but before this PR solve_steady always crashed before reaching it, and now it executes. Worth a comment pointing at #33 so the next reader does not think the review missed it.
Questions
Should solve_steady silently override the constructor's adiabatic?
solve_steady(self, vol_rxn, adiabatic=False) assigns self.adiabatic = adiabatic unconditionally, so on this head PlugFlowReactor(adiabatic=True).solve_steady(vol) runs non-adiabatic and takes the heat-transfer branch. Separately, solve_steady(vol, adiabatic=True) leaves states_uo == ['mole_conc', 'temp', 'temp'] because nomenclature() already appended 'temp' for the non-isothermal default. Both are pre-existing and currently harmless ('temp' in states_uo is still the only test), and both were unreachable before this PR made solve_steady runnable. Is the constructor argument meant to be the default for the method argument, and is the duplicate append worth cleaning while the function is already being touched?
Tests run
Active environment: pharmapy conda env (no bare python on PATH).
| Command | Result |
|---|---|
pytest tests/test_reactor_modes.py -q -rs |
4 passed |
pytest tests/ -m "not assimulo" -q -rs |
33 passed, 12 deselected, 0 skipped |
pytest tests/ -m assimulo -q -rs |
12 passed, 33 deselected, 0 skipped |
pytest tests/ -m slow -q |
4 passed |
pytest tests/ -m integration -q |
10 passed |
pytest --collect-only -q vs --collect-only -q -m "" |
45 collected both ways; nothing hidden by marker filters |
Additional verification beyond the suite:
- New tests are not vacuously green. In a throwaway worktree at
1a3ab31(imports confirmed to resolve to that tree), all 4 fail with exactly the documentedAttributeError/UnboundLocalError× 3. - Both energy-balance fixes are load-bearing (base +
mole_conconly → stillCVodeError, as above). - Changed-line coverage. Traced
solve_steadyon this head:Reactors.py:1566,1568,1571,1572,1573,1609all execute. The fixture isisothermal: falseandsolve_steadydefaults toadiabatic=False, so the newtemp_htline is genuinely exercised. - Coil propagation.
BatchReactor(isothermal=False, ht_mode='coil').solve_unit(runtime=1)→NotImplementedError(traceback originReactors.py:833→energy_balances→heat_transfer).CSTR(...)andSemibatchReactor(...)with the same settings →assimulo.solvers.sundials.CVodeError: 'The right-hand side function failed at the first call. At time 0.000000.', on both this head and the base. paramest_wrapperfallback. On this head,CSTR(return_sens=True).paramest_wrapper(...)→NotImplementedError;CSTR(return_sens=False).paramest_wrapper(...)→ returnsc_profof shape(5, 3).- Shapes behind
np.sum. In the steady pathdeltah_rxnandratesare both(n_rxns,), sonp.sum(a * b)equalsnp.dot(a, b)exactly.
Warnings in the targeted run are pre-existing DeprecationWarning: invalid escape sequence from unrelated module docstrings and a numpy.matlib PendingDeprecationWarning; none originate in this diff. No generated residue: git status --porcelain is clean, and the __pycache__ / .pytest_cache directories my runs created are gitignored. The temporary base worktree was removed.
GitHub checks (head 859300d)
Latest run set for the reviewed head: Core tests — SUCCESS; Assimulo integration tests — SUCCESS. No skipped or failed checks. Because the base is unchanged since the branch point, these PR-head checks also describe the merge result. mergeStateStatus is CLEAN; master has no branch protection and reviewDecision is empty, so GitHub reports no approval gate.
Merge readiness
The diff is small, correct on the paths it tests, and CI is green — but blocking finding 1 leaves a documented acceptance criterion of the linked issue unmet for two of the three reactors that issue names, and blocking finding 2 means the test suite reports success over that gap. Both are fixed by the same short guard plus one parametrized test.
I would not merge this until the blocking issues above are addressed.
This review is a COMMENT because I authored the PR and GitHub rejects APPROVE/REQUEST_CHANGES from the author; this account cannot provide a formal approval if one is needed. GitHub does not currently report a review requirement on master.
| # Heat transfer area | ||
| if self.ht_mode == 'coil': # Half pipe heat transfer | ||
| pass | ||
| raise NotImplementedError( |
There was a problem hiding this comment.
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.
There was a problem hiding this comment.
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.
There was a problem hiding this comment.
Why do we even allow that ht_mode to be called? Is it implemented elsewhere? Should this become an issue to be implemented?
| check_modeling_objects(self) | ||
|
|
||
| if eval_sens: | ||
| raise NotImplementedError( |
There was a problem hiding this comment.
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")There was a problem hiding this comment.
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.
| # 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 |
There was a problem hiding this comment.
Nonblocking — this new line is exactly the np.dot that sits commented out two lines above it.
I verified the shapes on this head. In the steady path temp is scalar, so getHeatOfRxn collapses its result to 1-D at ThermoModule.py:279-280 (the n_temp == 1 branch), and get_rxn_rates(..., overall_rates=False) returns (n_rxns,). Both operands are (n_rxns,), so np.sum(a * b) and np.dot(a, b) are identical here. The fix is correct.
Prefer np.dot(deltah_rxn, rates) anyway: it is shape-strict, so a later change that makes deltah_rxn 2-D raises instead of silently collapsing the entire array to one scalar. That silent-collapse failure mode is the same class of bug that produced the original .sum(axis=1) crash, and np.sum is the one spelling that cannot detect it.
While here, please delete the now-stale scaffolding directly above:
# 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**3The line is settled and verified; leaving a TODO: Check if this is correct above it tells the next reader the opposite.
There was a problem hiding this comment.
Addressed in 0f28998: steady PFR heat source reduction now uses np.dot(deltah_rxn, rates), and the stale commented alternatives/TODO were removed.
| 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'] |
There was a problem hiding this comment.
Nonblocking — two small things on this line.
Accessor consistency. The sibling dynamic balance reads the same quantity as self.Utility.get_inputs(time)['temp_in'] (Reactors.py:1693). Elsewhere in this file evaluate_inputs(0) is the initial-state idiom only (Reactors.py:1133 and 1329, both for tht_init). Using it inside a balance introduces a second spelling of one concept. Both return temp_in for a static CoolingWater, so this is not a bug today, but the two diverge the moment a DynamicInlet is attached to the utility.
Placement. energy_steady is the CVode RHS, so this rebuilds a dict on every evaluation (93 output points / 121 nonlinear iterations in the new test) to fetch a value that is constant across the whole volume integration. Hoist it into solve_steady and read the attribute here:
# in solve_steady, before Explicit_Problem(...)
self.temp_ht_steady = self.Utility.evaluate_inputs(0)['temp_in']That also gives you a natural place to document why time 0 is the right sample when the independent variable of this solve is volume, not time — which is not obvious from the call site.
There was a problem hiding this comment.
Addressed in 0f28998: solve_steady now samples the steady utility inlet condition once into temp_ht_steady, with a comment explaining the time-zero choice for volume integration; the RHS reads that attribute.
| 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) |
There was a problem hiding this comment.
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.
There was a problem hiding this comment.
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.
|
Addressed the PR #107 review feedback in Changes made:
Tests run locally in the
CI on head
Not changed: the pre-existing |
bernalde
left a comment
bernalde
left a comment
There was a problem hiding this comment.
The changes in the comments also have two major changes to be seen from now on: function docs according to the numpy formatting, and units/magic numbers
| # Heat transfer area | ||
| if self.ht_mode == 'coil': # Half pipe heat transfer | ||
| pass | ||
| raise NotImplementedError( |
There was a problem hiding this comment.
Why do we even allow that ht_mode to be called? Is it implemented elsewhere? Should this become an issue to be implemented?
| # 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 |
There was a problem hiding this comment.
Units formatting and why do we need the 1000 factor?
| 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 |
| return events | ||
|
|
||
| def heat_transfer(self, temp, temp_ht, vol): | ||
| """Return reactor heat transfer duty for supported heat-transfer modes.""" |
There was a problem hiding this comment.
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
2 participants
Summary
Closes #70.
LiquidStream.mole_concinPlugFlowReactor.solve_steadyinstead of the non-existentconcentrattribute.temp_ininput.eval_sens=Truefall-throughs and coil heat-transfer fall-through with explicitNotImplementedErrormessages naming the unsupported modes.Tests
conda run -n pharmapy python -m pytest --collect-only -qconda run -n pharmapy python -m pytest tests/test_reactor_modes.py -qconda run -n pharmapy python -m pytest tests/ -m "not assimulo" -qconda run -n pharmapy python -m pytest tests/ -m assimulo -qgit -c core.whitespace=blank-at-eol,blank-at-eof,space-before-tab,tab-in-indent,cr-at-eol diff --cached --checkBranch Hygiene
masterorigin/masterat1a3ab31origin/mastercontainsupstream/masterat400e1e1