Fix reversible reaction Jacobian#114
Conversation
![chatgpt-codex-connector[bot]](https://avatars.githubusercontent.com/in/1144995?s=60&v=4){kind=small}
There was a problem hiding this comment.
💡 Codex Review
Here are some automated review suggestions for this pull request.
Reviewed commit: 22c62c9250
ℹ️ About Codex in GitHub
Your team has set up Codex to review pull requests in this repo. Reviews are triggered when you
- Open a pull request for review
- Mark a draft as ready
- Comment "@codex review".
If Codex has suggestions, it will comment; otherwise it will react with 👍.
Codex can also answer questions or update the PR. Try commenting "@codex address that feedback".
| conc = np.asarray(conc) | ||
| conc_correc = np.maximum(eps, conc) | ||
|
|
||
| keq_temp = self.equil_term(temp, self.delta_hrxn) |
There was a problem hiding this comment.
For reversible calls where the rate path supplies a temperature-corrected delta_hrxn (e.g. Reactors.material_balances computes deltah_rxn with Liquid_1.getHeatOfRxn(...) and passes it to get_rxn_rates), this new reverse derivative always recomputes Keq from self.delta_hrxn instead. When heat-capacity corrections make deltah_rxn differ from the reference value, the Jacobian is no longer the derivative of the residual used by the solver, so non-isothermal reversible reactors get the wrong reverse-term concentration columns; derivatives/get_rxn_rates(..., jac=True) needs a way to use the same per-call ΔH as the rate evaluation.
Useful? React with 👍 / 👎.
There was a problem hiding this comment.
Addressed in 6421028: runtime delta_hrxn is now accepted by get_rxn_rates(..., jac=True) / derivatives, _reverse_df_dstates uses it for Keq, and _BaseReactor.get_jacobians passes the reactor's corrected heat of reaction into the Jacobian path. Added tests for the runtime-DeltaH mismatch and reactor handoff.
parkyr
left a comment
parkyr
left a comment
There was a problem hiding this comment.
Maintainer review — COMMENT (posted by the PR author)
The reverse-term math is correct and the fix is a clear improvement: the added _reverse_df_dstates computes ∂/∂C_k [∏ C_prod^β / Keq] = β_k · r_term / C_k, derivatives now returns df − dr (verified against a finite-difference Jacobian for the isothermal A↔B case — matches within FD tolerance). The forward-only irreversible path is preserved, batched 2-D concentration/time Jacobians are supported via einsum, and the elem_df_dstates rewrite incidentally fixes a latent 2-D broadcasting bug (the old / (conc + eps) would misalign for conc.ndim == 2). CI on the head is green (Core tests + Assimulo integration).
One correctness gap remains, described below. It does not regress master (the base Jacobian omitted the reverse term entirely), so I am not treating it as blocking — but it is the natural follow-up and the PR body's "match a finite-difference Jacobian of the reversible net species rates" claim only holds for constant ΔH.
Nonblocking: reverse-term Keq uses reference ΔH, not the runtime ΔH the solver integrates
_reverse_df_dstates recomputes Keq from self.equil_term(temp, self.delta_hrxn) — the reference ΔH stored at construction. But the reactor residual passes a temperature-corrected ΔH into the rate path: Reactors.material_balances / energy_balances compute deltah_rxn = self.Liquid_1.getHeatOfRxn(..., self.Kinetics.delta_hrxn, self.Kinetics.tref_hrxn) and call get_rxn_rates(..., delta_hrxn=deltah_rxn), while get_jacobians calls self.Kinetics.derivatives(conc, temp) with no ΔH. When Δcp ≠ 0 (non-isothermal reversible reactor) the Jacobian is no longer the derivative of the residual CVode integrates. This matters most for the sensitivity system (dsens_dt = jac_states.dot(sens) + jac_params in Reactors.get_jacobians), which uses jac_states directly, and secondarily degrades Newton convergence.
Demonstration (same fix, keq_params=[2], delta_hrxn=[-5e4], tref_hrxn=298.15, temp=350; residual fed a runtime delta_hrxn=[-4e4] as getHeatOfRxn would):
Jacobian (uses self.delta_hrxn = -5e4): B-column = [ 9.92, -9.92]
FD of residual (uses runtime = -4e4): B-column = [ 5.46, -5.46]
max abs diff over (A,B) block: 4.47
Concrete fix: thread the per-call ΔH through get_rxn_rates(..., jac=True, delta_hrxn=...) → derivatives(conc, temp, delta_hrxn=...) → _reverse_df_dstates, and pass the reactor's already-computed deltah_rxn from get_jacobians (mirroring material_balances). If that plumbing is out of scope for this PR, document the limitation in the docstring and open a follow-up, since derivatives never accepted runtime ΔH — this PR is the first to make the Jacobian ΔH-dependent.
Nonblocking: test coverage does not exercise the runtime-ΔH / Keq(T) path
All three new tests use delta_hrxn=[0.0], so Keq is constant and the reference-vs-runtime ΔH question above is invisible to them. Add a case with nonzero ΔH at temp ≠ tref_hrxn that compares derivatives against a finite difference of the same residual the reactor integrates (i.e. get_rxn_rates(..., delta_hrxn=deltah_runtime)), so the gap above is either guarded or explicitly documented as deferred.
Review summary
- Blocking: none.
- Nonblocking: (1) reverse-term Keq uses reference ΔH instead of the runtime ΔH the residual uses — wrong Jacobian/sensitivities for non-isothermal reversible reactors with Δcp ≠ 0 (strict improvement over base, not a regression); (2) tests only cover constant ΔH.
- Questions: none.
- Tests run (repo
.venv, Python 3.11.15, at head22c62c9):pytest tests/test_reversible_jacobian_reverse_term.py tests/test_cryst_kinetics.py→ 8 passed;pytest tests/ -m "not assimulo"→ 36 passed, 8 deselected;--collect-only→ 44 collected. Independently reproduced the ΔH mismatch shown above. - Merge-readiness: No blocking issues; head CI green (Core tests, Assimulo integration). This review is a
COMMENTbecause this account authored the PR and GitHub rejects authorAPPROVE; an eligible maintainer's approval is required only if branch protection reports a review gate.
|
Addressed the runtime heat-of-reaction review feedback in Commits pushed:
Changes made:
Tests run:
Intentionally not addressed: none. Remaining notes:
|
|
|
||
| df_dconc = f_term[..., np.newaxis] * self.params_f / (conc + eps) | ||
| if conc.ndim == 1: | ||
| conc_term = conc + eps |
There was a problem hiding this comment.
Is this eps a constant across the whole code? What’s its value? Is it only added to avoid division by zero is a concentration is small? I’m not very satisfied with this from a numerical point of view
| @@ -532,12 +532,55 @@ def equilibrium_model(self, conc, temp, deltah_rxn): | |||
|
|
|||
| def elem_df_dstates(self, conc): | |||
|
|
|||
There was a problem hiding this comment.
Let’s use the chance to document this and all new and exisiting functions in this file according to the numpy format
|
|
||
| return jac_states | ||
| else: # --------------- wrt parameters | ||
| dk_dphi = self.dk_dkparams(temp).T |
There was a problem hiding this comment.
We need units for all these, espcially since we are computing derivatives
| @@ -0,0 +1,268 @@ | |||
| # -*- coding: utf-8 -*- | |||
| """Regression test for issue #23. | |||
There was a problem hiding this comment.
I really like this format for the tests. Commenting each test function sounds unnecessary but this is a nice summary. We should keep applying it
2 participants
Closes #23.
Summary
RxnKineticsstate Jacobians.Triage / Acceptance Criteria
Uses the maintainer triage artifact at #23 (comment).
Covered:
derivativesandget_rxn_rates(..., jac=True)match a finite-difference Jacobian of the reversible net species rates.Related Issues
Tests
.venv/bin/python -m pytest tests/test_reversible_jacobian_reverse_term.py -qfailed on unchangedorigin/masterwith the missing reverse-termC_Bcolumn..venv/bin/python -m pytest tests/test_reversible_jacobian_reverse_term.py -q.venv/bin/python -m pytest tests/test_reversible_jacobian_reverse_term.py tests/test_cryst_kinetics.py -q.venv/bin/python -m pytest tests/ -m "not assimulo" -q.venv/bin/python -m pytest --collect-only -qgit -c core.whitespace=blank-at-eol,blank-at-eof,space-before-tab,tab-in-indent,cr-at-eol diff --cached --checkLocal environment: repository
.venvwith Python 3.11.15. The changed kinetics code is pure NumPy and does not require the assimulo backend.Branch Hygiene
masterfix/issue-23-reversible-jacobianorigin/masterat1a3ab3155216f07e9a18c60687a7b6b1c5959202upstream/masterat400e1e13ad4ecd3ba202d5fae8c70ec73bdf9c3eis already an ancestor oforigin/master.PharmaPy/Kinetics.py, but it addresses issue High: [Kinetics] Custom rate-law setup and vectorized equilibrium reverse term crash #71. No open PR touchestests/test_reversible_jacobian_reverse_term.py.