Implement sensitivity for SurfaceReactor#2935
Implement sensitivity for SurfaceReactor#2935sevyharris wants to merge 1 commit intoReactionMechanismGenerator:mainfrom
Conversation
There was a problem hiding this comment.
Pull request overview
Note
Copilot was unable to run its full agentic suite in this review.
Implements analytical Jacobian support for SurfaceReactor so solver sensitivities can be computed, and adjusts base solver logic to avoid accessing P for constant-volume systems (where SurfaceReactor tracks P_initial only).
Changes:
- Enables sensitivity residual evaluation for
SurfaceReactor. - Adds an analytical
SurfaceReactor.jacobian()implementation (ported fromSimpleReactor). - Avoids computing
R*T/Pin the base reactor path whenconstant_volumeis enabled.
Reviewed changes
Copilot reviewed 2 out of 2 changed files in this pull request and generated 10 comments.
| File | Description |
|---|---|
| rmgpy/solver/surface.pyx | Enables sensitivities and adds an analytical Jacobian implementation for SurfaceReactor. |
| rmgpy/solver/base.pyx | Guards self.P access when running sensitivity logic under constant-volume operation. |
| V = self.V # constant volume | ||
|
|
||
| Ctot = self.P.value_si / (constants.R * self.T.value_si) |
|
|
||
|
|
||
| elif ir[j, 2] == -1: # only two reactants | ||
| corr = - k * C[ir[j, 0]] * C[ir[j, 1]] / Ctot |
| for i in range(num_core_species): | ||
| pd[ir[j, 0], i] -= corr | ||
| pd[ir[j, 1], i] -= corr | ||
|
|
||
| pd[ip[j, 0], ir[j, 1]] += deriv | ||
| for i in range(num_core_species): | ||
| pd[ip[j, 0], i] += corr | ||
| if ip[j, 1] != -1: | ||
| pd[ip[j, 1], ir[j, 1]] += deriv | ||
| for i in range(num_core_species): | ||
| pd[ip[j, 1], i] += corr | ||
| if ip[j, 2] != -1: | ||
| pd[ip[j, 2], ir[j, 1]] += deriv | ||
| for i in range(num_core_species): | ||
| pd[ip[j, 2], i] += corr | ||
|
|
||
|
|
||
| else: # three reactants | ||
| corr = - 2 * k * C[ir[j, 0]] * C[ir[j, 1]] * C[ir[j, 2]] / Ctot |
|
|
||
|
|
||
| elif ip[j, 2] == -1: # only two reactants | ||
| corr = -k * C[ip[j, 0]] * C[ip[j, 1]] / Ctot |
| for i in range(num_core_species): | ||
| pd[ip[j, 0], i] -= corr | ||
| pd[ip[j, 1], i] -= corr | ||
|
|
||
| pd[ir[j, 0], ip[j, 1]] += deriv | ||
| for i in range(num_core_species): | ||
| pd[ir[j, 0], i] += corr | ||
| if ir[j, 1] != -1: | ||
| pd[ir[j, 1], ip[j, 1]] += deriv | ||
| for i in range(num_core_species): | ||
| pd[ir[j, 1], i] += corr | ||
| if ir[j, 2] != -1: | ||
| pd[ir[j, 2], ip[j, 1]] += deriv | ||
| for i in range(num_core_species): | ||
| pd[ir[j, 2], i] += corr | ||
|
|
||
|
|
||
| else: # three reactants | ||
| corr = - 2 * k * C[ip[j, 0]] * C[ip[j, 1]] * C[ip[j, 2]] / Ctot | ||
| if (ip[j, 0] == ip[j, 1] & ip[j, 0] == ip[j, 2]): | ||
| deriv = 3 * k * C[ip[j, 0]] * C[ip[j, 0]] | ||
| pd[ip[j, 0], ip[j, 0]] -= 3 * deriv | ||
| for i in range(num_core_species): | ||
| pd[ip[j, 0], i] -= 3 * corr | ||
|
|
||
| pd[ir[j, 0], ip[j, 0]] += deriv | ||
| for i in range(num_core_species): | ||
| pd[ir[j, 0], i] += corr | ||
| if ir[j, 1] != -1: | ||
| pd[ir[j, 1], ip[j, 0]] += deriv | ||
| for i in range(num_core_species): | ||
| pd[ir[j, 1], i] += corr | ||
| if ir[j, 2] != -1: | ||
| pd[ir[j, 2], ip[j, 0]] += deriv | ||
| for i in range(num_core_species): | ||
| pd[ir[j, 2], i] += corr |
|
|
||
| else: # three reactants | ||
| corr = - 2 * k * C[ir[j, 0]] * C[ir[j, 1]] * C[ir[j, 2]] / Ctot | ||
| if (ir[j, 0] == ir[j, 1] & ir[j, 0] == ir[j, 2]): |
|
|
||
| else: # three reactants | ||
| corr = - 2 * k * C[ip[j, 0]] * C[ip[j, 1]] * C[ip[j, 2]] / Ctot | ||
| if (ip[j, 0] == ip[j, 1] & ip[j, 0] == ip[j, 2]): |
| def jacobian(self, double t, np.ndarray[np.float64_t, ndim=1] y, np.ndarray[np.float64_t, ndim=1] dydt, | ||
| double cj, np.ndarray[np.float64_t, ndim=1] senpar = np.zeros(1, float)): |
| @cython.boundscheck(False) | ||
| def jacobian(self, double t, np.ndarray[np.float64_t, ndim=1] y, np.ndarray[np.float64_t, ndim=1] dydt, | ||
| double cj, np.ndarray[np.float64_t, ndim=1] senpar = np.zeros(1, float)): | ||
| """ | ||
| Return the analytical Jacobian for the reaction system. | ||
| """ |
| @cython.boundscheck(False) | ||
| def jacobian(self, double t, np.ndarray[np.float64_t, ndim=1] y, np.ndarray[np.float64_t, ndim=1] dydt, | ||
| double cj, np.ndarray[np.float64_t, ndim=1] senpar = np.zeros(1, float)): | ||
| """ | ||
| Return the analytical Jacobian for the reaction system. | ||
| """ |
cd789c1 to
9fb6750
Compare
Regression Testing Results
Detailed regression test results.Regression test aromatics:Reference: Execution time (DD:HH:MM:SS): 00:00:00:55 aromatics Passed Core Comparison ✅Original model has 15 species. aromatics Failed Edge Comparison ❌Original model has 106 species. Non-identical thermo! ❌
thermo: Thermo group additivity estimation: group(Cs-(Cds-Cds)CsCsH) + group(Cs-(Cds-Cds)CsCsH) + group(Cs-(Cds-Cds)(Cds-Cds)CsH) + group(Cs-(Cds-Cds)CsHH) + group(Cds-CdsCsCs) + group(Cds-CdsCsH) + group(Cds-CdsCsH) + group(Cds-CdsCsH) + polycyclic(s3_4_5_ene_3) + polycyclic(s2_4_5_diene_1_5) + polycyclic(s2_5_5_diene_1_5) - ring(Cyclobutene) - ring(Cyclopentene) - ring(Cyclopentene) + radical(cyclopentene-allyl) Non-identical thermo! ❌
thermo: Thermo group additivity estimation: group(Cs-(Cds-Cds)CsCsH) + group(Cs-(Cds-Cds)(Cds-Cds)CsH) + group(Cs-(Cds-Cds)(Cds-Cds)CsH) + group(Cs-CsCsHH) + group(Cds-CdsCsCs) + group(Cds-CdsCsH) + group(Cds-CdsCsH) + group(Cds-CdsCsH) + Estimated bicyclic component: polycyclic(s2_3_5_ane) - ring(Cyclopropane) - ring(Cyclopentane) + ring(Cyclopentene) + ring(Cyclopropane) + polycyclic(s2_3_6_diene_0_3) + Estimated bicyclic component: polycyclic(s3_5_6_ane) - ring(Cyclopentane) - ring(Cyclohexane) + ring(Cyclopentene) + ring(1,4-Cyclohexadiene) - ring(Cyclopropane) - ring(Cyclopentene) - ring(1,4-Cyclohexadiene) + radical(cyclopentene-4) Non-identical thermo! ❌
thermo: Thermo group additivity estimation: group(Cs-(Cds-Cds)(Cds-Cds)CsCs) + group(Cs-(Cds-Cds)CsCsH) + group(Cs-(Cds-Cds)CsCsH) + group(Cs-CsCsHH) + group(Cds-CdsCsH) + group(Cds-CdsCsH) + group(Cds-CdsCsH) + group(Cds-CdsCsH) + polycyclic(s2_4_4_ene_1) + polycyclic(s1_4_5_diene_1_6) + polycyclic(s3_4_5_ene_1) - ring(Cyclobutene) - ring(Cyclobutane) - ring(Cyclopentene) + radical(bicyclo[2.1.1]hex-2-ene-C5) Non-identical thermo! ❌
thermo: Thermo group additivity estimation: group(Cs-(Cds-Cds)(Cds-Cds)(Cds-Cds)H) + group(Cds-Cds(Cds-Cds)(Cds-Cds)) + group(Cds-CdsCsH) + group(Cds-CdsCsH) + group(Cds-Cds(Cds-Cds)H) + group(Cds-Cds(Cds-Cds)H) + group(Cds-CdsCsH) + group(Cdd-CdsCds) + Estimated bicyclic component: polycyclic(s4_6_6_ane) - ring(Cyclohexane) - ring(Cyclohexane) + ring(1,4-Cyclohexadiene) + ring(124cyclohexatriene) Non-identical thermo! ❌
thermo: Thermo group additivity estimation: group(Cs-(Cds-Cds)(Cds-Cds)CsH) + group(Cs-(Cds-Cds)CsHH) + group(Cds-Cds(Cds-Cds)(Cds-Cds)) + group(Cds-CdsCsH) + group(Cds-CdsCsH) + group(Cds-CdsCsH) + group(Cds-Cds(Cds-Cds)H) + group(Cds-Cds(Cds-Cds)H) + Estimated bicyclic component: polycyclic(s4_6_6_ane) - ring(Cyclohexane) - ring(Cyclohexane) + ring(1,3-Cyclohexadiene) + ring(1,3-Cyclohexadiene) Non-identical thermo! ❌
thermo: Thermo group additivity estimation: group(Cs-(Cds-Cds)(Cds-Cds)CsH) + group(Cs-(Cds-Cds)(Cds-Cds)CsH) + group(Cds-Cds(Cds-Cds)(Cds-Cds)) + group(Cds- CdsCsH) + group(Cds-CdsCsH) + group(Cds-CdsCsH) + group(Cds-CdsCsH) + group(Cds-Cds(Cds-Cds)H) + group(Cds-Cds(Cds-Cds)H) + group(Cds-CdsHH) + Estimated bicyclic component: polycyclic(s4_6_6_ane) - ring(Cyclohexane) - ring(Cyclohexane) + ring(1,3-Cyclohexadiene) + ring(1,3-Cyclohexadiene) + radical(Cds_P) Non-identical thermo! ❌
thermo: Thermo group additivity estimation: group(Cs-(Cds-Cds)(Cds-Cds)CsH) + group(Cs-(Cds-Cds)(Cds-Cds)CsH) + group(Cds-Cds(Cds-Cds)(Cds-Cds)) + group(Cds- CdsCsH) + group(Cds-CdsCsH) + group(Cds-CdsCsH) + group(Cds-CdsCsH) + group(Cds-Cds(Cds-Cds)H) + group(Cds-Cds(Cds-Cds)H) + group(Cds-CdsHH) + Estimated bicyclic component: polycyclic(s4_6_6_ane) - ring(Cyclohexane) - ring(Cyclohexane) + ring(1,3-Cyclohexadiene) + ring(1,3-Cyclohexadiene) Non-identical kinetics! ❌
kinetics: Non-identical kinetics! ❌
kinetics: Non-identical kinetics! ❌
kinetics: Non-identical kinetics! ❌
kinetics: Non-identical kinetics! ❌
kinetics: Non-identical kinetics! ❌
kinetics: Non-identical kinetics! ❌
kinetics: Non-identical kinetics! ❌
kinetics: Non-identical kinetics! ❌
kinetics: Non-identical kinetics! ❌
kinetics: Non-identical kinetics! ❌
kinetics: Non-identical kinetics! ❌
kinetics: Non-identical kinetics! ❌
kinetics: Non-identical kinetics! ❌
kinetics: Non-identical kinetics! ❌
kinetics: Non-identical kinetics! ❌
kinetics: Non-identical kinetics! ❌
kinetics: Non-identical kinetics! ❌
kinetics: Non-identical kinetics! ❌
kinetics: Non-identical kinetics! ❌
kinetics: Non-identical kinetics! ❌
kinetics: DetailsObservables Test Case: Aromatics Comparison✅ All Observables varied by less than 0.500 on average between old model and new model in all conditions! aromatics Passed Observable Testing ✅Regression test liquid_oxidation:Reference: Execution time (DD:HH:MM:SS): 00:00:02:01 liquid_oxidation Passed Core Comparison ✅Original model has 37 species. liquid_oxidation Failed Edge Comparison ❌Original model has 214 species. DetailsObservables Test Case: liquid_oxidation Comparison✅ All Observables varied by less than 0.100 on average between old model and new model in all conditions! liquid_oxidation Passed Observable Testing ✅Regression test nitrogen:Reference: Execution time (DD:HH:MM:SS): 00:00:01:04 nitrogen Failed Core Comparison ❌Original model has 41 species. nitrogen Failed Edge Comparison ❌Original model has 133 species. DetailsObservables Test Case: NC Comparison✅ All Observables varied by less than 0.200 on average between old model and new model in all conditions! nitrogen Passed Observable Testing ✅Regression test oxidation:Reference: Execution time (DD:HH:MM:SS): 00:00:01:47 oxidation Passed Core Comparison ✅Original model has 59 species. oxidation Passed Edge Comparison ✅Original model has 230 species. DetailsObservables Test Case: Oxidation Comparison✅ All Observables varied by less than 0.500 on average between old model and new model in all conditions! oxidation Passed Observable Testing ✅Errors occurred during observable testing
WARNING:root:Initial mole fractions do not sum to one; normalizing.
|
This adds a Jacobian function (the same one copied over from SimpleReactor) so that we can run get sensitivities from a SurfaceReactor
9fb6750 to
fa79adf
Compare
Motivation or Problem
The sensitivity implementation was never completed for the SurfaceReactor, so this adds a Jacobian function, (the same one copied over from SimpleReactor) so that we can run get sensitivities from a SurfaceReactor
Description of Changes
Testing
Demo notebook showing that surface sensitivities line up across RMG and Cantera
TODO: annotated version of the Jacobian function to convince myself it's fine to use SimpleReactor's version of the function
TODO: Test over a much wider range of concentrations, temperatures, and surface volume ratios
Reviewer Tips
Reminder to rebuild RMG to get this branch working.
Be careful if you're comparing surface mechanisms across RMG and Cantera. There are still a few edge cases where the kinetics don't match perfectly.
Here's an example mechanism where the kinetics have been closely vetted to match across RMG and Cantera:
chem_annotated.yaml.txt
chem_annotated-gas.inp.txt
chem_annotated-surface.inp.txt
species_dictionary.txt