Add CanteraWriter #2870
Add CanteraWriter #2870
Conversation
alongd
force-pushed
the
ct3
branch
2 times, most recently
from
a0983b8 to
24205ac
Compare
instead of the old from CTI
This breaks the circular dependency chain (simulate -> diffusionLimited -> Species -> pdep -> Reaction -> diffusionLimited) because diffusionLimited can now fully initialize before loading Species.
|
Thanks, @JacksonBurns! I wasn't aware of the previous efforts. |
Regression Testing Results
Detailed regression test results.Regression test aromatics:Reference: Execution time (DD:HH:MM:SS): 00:00:00:50 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)(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(124cyclohexatriene) + ring(124cyclohexatriene) 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. Non-identical kinetics! ❌
kinetics: 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:03 nitrogen Failed Core Comparison ❌Original model has 41 species. nitrogen Failed Edge Comparison ❌Original model has 133 species. Non-identical thermo! ❌
thermo: Thermo group additivity estimation: group(O2s-CdN3d) + group(N3d-OCd) + group(Cd-HN3dO) + ring(Cyclopropene) + radical(CdJ-NdO) Non-identical kinetics! ❌
kinetics: 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:51 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.
|
|
See also #2800 |
There was a problem hiding this comment.
Pull request overview
This pull request adds a new CanteraWriter module to RMG-Py that directly generates Cantera YAML files during mechanism generation iterations, replacing the previous approach of using the external ck2yaml script after convergence. This change provides more robust Cantera file generation and enables intermediate YAML outputs even when RMG doesn't converge. The PR also upgrades the Cantera dependency from version 2.6 to 3.0+, updating all related API calls throughout the codebase.
Changes:
- Added new
rmgpy/cantera.pymodule withCanteraWriterclass and YAML serialization functions - Integrated
CanteraWriterinto the RMG listener pattern for automatic updates during iterations - Updated all Cantera API calls to be compatible with Cantera 3.0+ (changed deprecated classes like
ThreeBodyReaction,FalloffReactionto unifiedReactionwiththird_bodyparameter) - Removed old
generate_cantera_filesmethod that relied on externalck2yamlscript - Added comprehensive test suite in
test/rmgpy/canteraTest.py - Updated existing tests to use Cantera 3.0 YAML format instead of deprecated CTI format
Reviewed changes
Copilot reviewed 13 out of 13 changed files in this pull request and generated 13 comments.
Show a summary per file
| File | Description |
|---|---|
| rmgpy/cantera.py | New module implementing CanteraWriter class and YAML serialization for species and reactions |
| rmgpy/rmg/main.py | Integrated CanteraWriter listener and removed old ck2yaml-based generation method |
| rmgpy/reaction.py | Updated to_cantera method to use Cantera 3.0 unified Reaction API |
| rmgpy/kinetics/falloff.pyx | Updated efficiencies assignment to use third_body.efficiencies for Cantera 3.0 |
| rmgpy/kinetics/diffusionLimited.py | Moved Species import to local scope to avoid circular dependency |
| test/rmgpy/canteraTest.py | New comprehensive test suite for CanteraWriter functionality |
| test/rmgpy/transportDataTest.py | Updated test to use from_yaml instead of deprecated fromCti |
| test/rmgpy/speciesTest.py | Updated test to use YAML format instead of CTI format |
| test/rmgpy/reactionTest.py | Fixed test assertions for Cantera 3.0 API changes (rate.rates, third_body.efficiencies) |
| test/rmgpy/tools/canteramodelTest.py | Added temp file handling for surface phase tests and updated assertions |
| test/rmgpy/rmg/mainTest.py | Removed obsolete TestCanteraOutput class, updated directory creation |
| environment.yml | Updated Cantera dependency from =2.6 to >=3.0 |
Comments suppressed due to low confidence (1)
rmgpy/reaction.py:400
- Keyword argument 'rate' is not a supported parameter name of Reaction.init.
ct_reaction = ct.Reaction(
reactants=ct_reactants,
products=ct_products,
rate=ct.InterfaceArrheniusRate()
)
💡 Add Copilot custom instructions for smarter, more guided reviews. Learn how to get started.
| elif isinstance(reaction, Reaction): | ||
| note_parts.append(f"Source: P{reaction.kinetics.comment}") |
There was a problem hiding this comment.
The note string appears malformed. The prefix "Source: P" followed by the comment doesn't make semantic sense. This should likely be just the comment itself without the "P" prefix, or the entire condition may be unnecessary since kinetics.comment is already handled separately on line 563.
| elif isinstance(reaction, Reaction): | |
| note_parts.append(f"Source: P{reaction.kinetics.comment}") |
| note_parts.append(clean_comment) | ||
|
|
||
| if reaction.specific_collider: | ||
| note_parts.append(f"Specific collider: {reaction.specific_collider.label}") |
There was a problem hiding this comment.
Potential AttributeError if reaction.specific_collider doesn't have a 'label' attribute. Consider adding a check or using getattr with a default value: getattr(reaction.specific_collider, 'label', str(reaction.specific_collider)).
| note_parts.append(f"Specific collider: {reaction.specific_collider.label}") | |
| collider = reaction.specific_collider | |
| collider_label = getattr(collider, 'label', str(collider)) | |
| note_parts.append(f"Specific collider: {collider_label}") |
| # Build the base dictionary | ||
| species_entry = { | ||
| 'name': get_label(species, species_list), | ||
| 'composition': atom_dict, | ||
| 'thermo': { | ||
| 'model': 'NASA7', | ||
| 'temperature-ranges': [sorted_polys[0].Tmin.value_si, sorted_polys[0].Tmax.value_si, | ||
| sorted_polys[1].Tmax.value_si], | ||
| 'data': [polys[0]['data'], polys[1]['data']] |
There was a problem hiding this comment.
Potential IndexError if thermo_data.polynomials has fewer than 2 polynomials. Consider adding a check to ensure len(sorted_polys) >= 2 before accessing sorted_polys[0] and sorted_polys[1].
| # Build the base dictionary | |
| species_entry = { | |
| 'name': get_label(species, species_list), | |
| 'composition': atom_dict, | |
| 'thermo': { | |
| 'model': 'NASA7', | |
| 'temperature-ranges': [sorted_polys[0].Tmin.value_si, sorted_polys[0].Tmax.value_si, | |
| sorted_polys[1].Tmax.value_si], | |
| 'data': [polys[0]['data'], polys[1]['data']] | |
| # Determine temperature ranges and data safely, ensuring at least one polynomial | |
| if not sorted_polys: | |
| raise ValueError(f"No NASA polynomials found for species {get_label(species, species_list)}") | |
| elif len(sorted_polys) == 1: | |
| # If only one polynomial is available, duplicate it to satisfy the expected two-range structure | |
| tmin = sorted_polys[0].Tmin.value_si | |
| tmax = sorted_polys[0].Tmax.value_si | |
| temperature_ranges = [tmin, tmax, tmax] | |
| thermo_data_blocks = [polys[0]['data'], polys[0]['data']] | |
| else: | |
| temperature_ranges = [sorted_polys[0].Tmin.value_si, | |
| sorted_polys[0].Tmax.value_si, | |
| sorted_polys[1].Tmax.value_si] | |
| thermo_data_blocks = [polys[0]['data'], polys[1]['data']] | |
| # Build the base dictionary | |
| species_entry = { | |
| 'name': get_label(species, species_list), | |
| 'composition': atom_dict, | |
| 'thermo': { | |
| 'model': 'NASA7', | |
| 'temperature-ranges': temperature_ranges, | |
| 'data': thermo_data_blocks, |
| notes = list() | ||
| try: | ||
| notes.append(species.to_smiles()) | ||
| except: |
There was a problem hiding this comment.
Bare except clause should specify the exception type(s) being caught. Consider using 'except Exception:' or catching specific exceptions like 'except (AttributeError, ValueError):' to avoid catching system-exiting exceptions like KeyboardInterrupt.
| except: | |
| except Exception: |
| 'geometry': 'atom' if td.shapeIndex == 0 else 'linear' if td.shapeIndex == 1 else 'nonlinear', | ||
| 'well-depth': td.epsilon.value_si / constants.R, | ||
| 'diameter': td.sigma.value_si, | ||
| 'dipole': dipole, |
There was a problem hiding this comment.
The polarizability variable is calculated but never added to the transport dictionary. It should be included as 'polarizability': polarizability in the transport_entry dictionary if it's needed in the Cantera YAML format.
| 'dipole': dipole, | |
| 'dipole': dipole, | |
| 'polarizability': polarizability, |
| ct_reaction = ct.Reaction( | ||
| reactants=ct_reactants, | ||
| products=ct_products, | ||
| tbody=ct_collider, | ||
| third_body=ct_collider, | ||
| rate=ct.TroeRate() | ||
| ) |
There was a problem hiding this comment.
Keyword argument 'rate' is not a supported parameter name of Reaction.init.
Keyword argument 'third_body' is not a supported parameter name of Reaction.init.
| ct_reaction = ct.Reaction( | ||
| reactants=ct_reactants, | ||
| products=ct_products, | ||
| third_body="M", | ||
| rate=ct.TroeRate() | ||
| ) |
There was a problem hiding this comment.
Keyword argument 'rate' is not a supported parameter name of Reaction.init.
Keyword argument 'third_body' is not a supported parameter name of Reaction.init.
| ct_reaction = ct.Reaction( | ||
| reactants=ct_reactants, | ||
| products=ct_products, | ||
| tbody=ct_collider, | ||
| third_body=ct_collider, | ||
| rate=ct.LindemannRate() | ||
| ) |
There was a problem hiding this comment.
Keyword argument 'rate' is not a supported parameter name of Reaction.init.
Keyword argument 'third_body' is not a supported parameter name of Reaction.init.
| ct_reaction = ct.Reaction( | ||
| reactants=ct_reactants, | ||
| products=ct_products, | ||
| third_body="M", | ||
| rate=ct.LindemannRate() | ||
| ) |
There was a problem hiding this comment.
Keyword argument 'rate' is not a supported parameter name of Reaction.init.
Keyword argument 'third_body' is not a supported parameter name of Reaction.init.
| except: | ||
| pass |
There was a problem hiding this comment.
'except' clause does nothing but pass and there is no explanatory comment.
| except: | |
| pass | |
| except Exception: | |
| # SMILES is optional metadata; log and continue if it cannot be generated | |
| logging.getLogger(__name__).debug( | |
| "Unable to convert species %r to SMILES in species_to_dict", getattr(species, "label", species), | |
| exc_info=True, | |
| ) |
rwest
left a comment
rwest
left a comment
There was a problem hiding this comment.
This looks great - thanks for doing it Alon. We've been wishing for it for a long time (see #2321 from 2022, then #2770 and #2800).
I made a few comments, from just reading the diff on GitHub. I have not actually tried running this, but if you say it works I believe you :)
I would like to know what unit testing we're doing and how exhaustive.
Can we round-trip to/from Cantera and back?
I don't want striving to for perfection to prevent this ever happening (see #2321!) but a "wish list" we can revisit would be nice.
I'm also curious about speed. I think we could probably speed this up a lot, fairly easily (templating, caching. Making dicts each time and using yaml libraries will be relatively slow). But again, that could come after it's definitely working and tested and implemented, and profiled.
| Convert an RMG Reaction object to a LIST of Cantera YAML dictionaries. | ||
| """ | ||
| # Check for MultiKinetics (duplicates grouped in one RMG object) | ||
| if isinstance(reaction.kinetics, (MultiArrhenius, MultiPDepArrhenius)): |
There was a problem hiding this comment.
Should MultiPDHepArrhenius become a single reaction with PLOG format, rather than multiple duplicate reactions? or is this taken care of?
| }) | ||
| entry['rate-constants'] = rates | ||
|
|
||
| else: |
There was a problem hiding this comment.
FYI Cantera (though I forget which version) now has Blowers Masel kinetics types, if we want to include that. (Though his could come later - we never currently have them IN a final RMG core/edge model - because we always convert to Arrhenius)
| kin = reaction.kinetics | ||
| if isinstance(kin, (ThirdBody, Lindemann, Troe)): | ||
| if hasattr(reaction, 'specific_collider') and reaction.specific_collider: | ||
| suffix = " + " + get_label(reaction.specific_collider, species_list) |
There was a problem hiding this comment.
if this is a falloff reaction, should the specific collider be in parentheses.
See https://cantera.org/3.1/yaml/reactions.html#falloff
like suffix = f" (+ {get_label(reaction.specific_collider, species_list)})"
I think that's only for falloff like Lindemann and Troe but not three-body.
Does our ThirdBody not map to Cantera's three-body?
There was a problem hiding this comment.
Probably want unit tests for round-trip testing of a whole bunch of reaction and kinetics types.
|
Thanks @rwest! Supper appreciated! We need to decide whether to continue with this approach (and I'll address the above comments), or whether we still want to use/support all the |
Motivation or Problem
Currently RMG outputs a Cantera YAML file only after it converges, using the external ck2yaml script.
Sometimes the script fails, and sometimes if RMG doesn't converge, a Cantera YAML model might still be useful as an intermediate step to assess the development status.
Description of Changes
A CanteraWriter module was added
Testing
Tests were added