Vibrational-relaxation from plasma treatment enhancements by Gaetanosaure · Pull Request #2145 · Cantera/cantera

Gaetanosaure · 2026-06-24T11:58:35Z

Changes proposed in this pull request

This PR proposes to implement the point 6 from the former pull request plasma treatment enhancements. I write it here again for reference:

Added a reaction class for vibrational relaxation reactions, including detailed relaxation and relaxation modelled using the vibrational energy equation. Four possible models are currently available : one for detailed vibration named “multi-state-resolved” (taking reactions having a base rate constant in the form of the reactions given in the Zhong mechanism and a scaling using the harmonic approximation of the SSH theory), and three options implementing the relaxation rates of the mean vibrational energy for mean vibrational energy models: “constant”, “castela” and “starikovskiy”. More details about these models can be found in the documentation of the header of the new class VibrationalRelaxationRate.h. The main point is that all of these models implement the following equation for the reaction rates:

However, depending on the chosen model, all parameters are not required and some default values are applied.

This PR also incorporates some modifications to TwoTempPlasmaRate that were done in the previous PR without having a particular number.

Closes #2144

AI Statement (required)

Limited use of generative AI.
Standard or boilerplate code snippets were generated with AI and manually reviewed;
all design, logic, and implementation decisions were made by the contributor.
Examples: IDE code-completions or brief LLM queries for common patterns.

codecov · 2026-06-25T22:24:30Z

Codecov Report

❌ Patch coverage is 5.85366% with 386 lines in your changes missing coverage. Please review.
✅ Project coverage is 77.62%. Comparing base (ddb114a) to head (b50eb2f).
⚠️ Report is 59 commits behind head on main.

Files with missing lines	Patch %	Lines
src/kinetics/VibrationalRelaxationRate.cpp	0.00%	341 Missing ⚠️
src/thermo/PlasmaPhase.cpp	39.53%	23 Missing and 3 partials ⚠️
...clude/cantera/kinetics/VibrationalRelaxationRate.h	0.00%	14 Missing ⚠️
src/kinetics/TwoTempPlasmaRate.cpp	50.00%	4 Missing ⚠️
src/kinetics/ReactionRateFactory.cpp	0.00%	1 Missing ⚠️

Additional details and impacted files

@@            Coverage Diff             @@
##             main    #2145      +/-   ##
==========================================
- Coverage   77.80%   77.62%   -0.19%     
==========================================
  Files         452      455       +3     
  Lines       54332    55838    +1506     
  Branches     9111     9174      +63     
==========================================
+ Hits        42275    43343    +1068     
- Misses       9016     9443     +427     
- Partials     3041     3052      +11

☔ View full report in Codecov by Harness.
📢 Have feedback on the report? Share it here.

🚀 New features to boost your workflow:

❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.

speth · 2026-07-09T13:38:46Z

@Gaetanosaure - before you go any further on making changes in this branch, let's get this rebased onto the current main branch and resolve the merge conflicts. If the conflicts are too tricky, let me know and I can do it for you.

Gaetanosaure · 2026-07-09T14:08:28Z

@speth ok let me try and do the rebase now, if I struggle too much I'll let you know!

Add the VibrationalRelaxationRate class for vibrational relaxation reactions, registered under the `vibrational-relaxation` rate type. It supports four models: `multi-state-resolved` (detailed V-V/V-T relaxation with an SSH harmonic-oscillator scaling of a base rate constant) and the mean-vibrational-energy V-T relaxation models `constant`, `castela`, and `starikovskiy`. See the class documentation for the governing rate expression and references. Also generalize TwoTempPlasmaRate with an optional gas-temperature exponent `b-gas` (default 0, which recovers the Kossyi form); a non-zero value applies the generalized k_f = A T^{b_g} T_e^b exp(-Ea_g/RT) exp(Ea_e(Te-T)/(R T Te)) expression.

…ional reservoir species Detect fictive vibrational-reservoir species (named "X(v)") and warn when the corresponding base species "X" is missing, or when the reservoir mole fraction grows into a non-negligible fraction of the base species -- a sign that the chosen reservoir energy is too low and is perturbing the chemistry.

…ate.h

Gaetanosaure · 2026-07-09T15:18:28Z

@speth the rebase should be clean now. Let me know if you notice anything strange but I think i managed to do it properly.

….cpp

…mmit

…pression completely. The minus sign will be born by YAML coefficients directly.

…ous commits

…iny declaration error

…tch / case

…dable

…type

Gaetanosaure · 2026-07-09T22:58:35Z

Hi @speth , after going through your comments on the code of the previous PR I tried to imagine the comments you would have done on the code of the present PR, and thus tried to enhance it as much as possible. Hopefully that will reduce your workload while processing this part of the integration (although it is very likely far from perfect). I believe I have gone as far as possible without your input, please tell me if you have any questions and let me know if I should do any further code modifications.

I would just like to emphasise that I have thought of this class as a container for all possible vibration models that could be later included in addition to the four already present. This is the reason for the parameter vibration-model: avoiding declaring a new reaction class every time a new vibration relaxation implementation is required. I thus tried to make it as modular as possible in the way it is organised.

The 'constant' option is redundant with a simple Arrhenius, but I find it interesting for the user to still be able to classify it as vibrational relaxation. If you think it should not be included we can get rid of it.

BangShiuh · 2026-07-11T11:04:50Z

@Gaetanosaure This looks clean and easier to review! Can you clarify if number 5. Implemented a new format for the YAML files ... is still in this PR or will be in this PR?

Gaetanosaure · 2026-07-13T06:44:03Z

Hi @BangShiuh ! Thanks for going through this bit of code as well.

Yes it will be but for more clarity and to make the review process easier we broke down the first big pull request (#2121) into three smaller ones.

The first chunk was chunk integrated last week: #2144

The second chunk is the present one on vibrational treatment, and then the last points (the main being taking into account the effective / elastic difference and changing the YAML format) will be done on a third one not yet opened because it depends on the prior implementation of the two other PRs.

…spellings and added getParameters function

Gaetanosaure · 2026-07-13T15:43:40Z

Hi @speth and @BangShiuh ,

I added a second extension to the TwoTemperaturePlasmaRate class that allows to take into account more plasma reactions via the parameter T-inv. This can serve for example in the case of R40 (Capitelli book Table 10.7):
N4+ + N2 => N2+ + N2 + N2
which is necessary in the simulation of the Rusterholtz test case (https://iopscience.iop.org/article/10.1088/0022-3727/46/46/464010).

BangShiuh · 2026-07-14T22:53:14Z

@Gaetanosaure Extending the TwoTemperaturePlasmaRate class is fine to me. I prefer a separate PR for it because this is a distinguish change. The better we document the change the less problems in the future. Also I see you use images showing the code / format in the PR description. This is not idea because it can hinder finding / maintaining the information. (also AI read text better.) I don't see code change for 5. so I assume you will do something about it?

Gaetanosaure · 2026-07-15T07:08:52Z

Hi @BangShiuh , yes indeed as I said in my answer to your first message point 5 will be integrated in a later PR separate from the present one. If you are interested in this point specifically I can notify you once it is open (that I, once the present one will be resolved).

I think there was a slight mix up between the commit splitting @speth did and the way I understood them, as a result the description of this PR is misleading in the sense that we are actually integrating the changed to the kinetics classes here: vibrational relaxation and TwoTempPlasmaRate. Handling of effective/elastic cross sections and accompaniying changes will be handled in the next PR building on this one.

Gaetanosaure · 2026-07-15T07:15:52Z

I have updated the PR description to make it match the changes actually in the commits

speth

Thanks for splitting this off, from the other changes, @Gaetanosaure. I have a number of suggestions for improving this implementation. Besides the specific in-line comments, to broader notes:

The main description of the new rates (formulas, concept, references) should be added in doc/sphinx/reference/kinetics/rate-constants.md. I'd suggest that much of what's in the C++ class documentation would be more accessible as part of this scientific reference documentation than buried in the Doxygen documentation, with the focus for Doxygen more on any implementation details. The YAML documentation in doc/sphinx/yaml/reactions.md should link to the section in the "reference" docs and focus on just describing the input fields, though repeating the equation is okay.
The major missing piece in this PR currently is any tests whatsoever. I'm not quite sure how this interacts with the other upcoming changes to the input format for electron collisions, but rather than deferring this, I think it would be worth adding some bare-bones inputs that can be used now in test/data. To parallel the existing test structure for reaction rates, there should be new tests of both pure-C++ construction (kineticsFromScratch.cpp), YAML-based construction (kineticsFromYaml.cpp) and the Python-based pathways (test_reaction.py). If you have any issues with the tests (the Python suite can be a bit tricky to debug) just give me a shout.

speth · 2026-07-14T16:59:22Z

+(sec-vibrational-relaxation-constant-rate)=
+
+### Vibrational relaxation: constant
+
+The `constant` vibrational relaxation rate expression is


I think it would be better to describe all of these variants further down, with the vibrational-relaxation reaction type. The Arrhenius and Blowers-Masel description were separated out because those parameterizations can be used in multiple different reaction types, whereas this is only used for the one reaction type. Probably having the two-temperature plasma rate expression explained up here is also a bit of organizational mistake.

speth · 2026-07-14T17:00:46Z

+No other rate parameters are allowed for this model. In particular, `b`, `n`, `B`,
+`C`, `D`, `m`, `E`, `z`, and `scaling` are not permitted.


I think the second statement here is unnecessary. Besides the redundancy, none of these other parameters have been introduced, so there's really nothing that would lead a user to imagining that they could be used. Likewise the corresponding statements for the other model options.

speth · 2026-07-14T17:05:14Z

+`scaling`
+: Dimensionless multiplicative scaling factor applied to the rate coefficient.
+Defaults to 1.0.


I don't understand the purpose of this parameter, if we already have A.

speth · 2026-07-14T17:06:53Z

+ C T^{-1/3}
+ D T^{-2/3}
+\right)
+$$


This is a bit excessive for splitting an equation onto multiple lines.

speth · 2026-07-14T17:09:57Z

+
+### Vibrational relaxation: castela
+
+The `castela` vibrational relaxation rate expression is based on the relaxation time


We capitalize proper nouns in YAML names:

Suggested change

The `castela` vibrational relaxation rate expression is based on the relaxation time

The `Castela` vibrational relaxation rate expression is based on the relaxation time

(so, likewise for Starikovskiy)

speth · 2026-07-17T19:44:35Z

+const string ModelConstant = "constant";
+const string ModelMultiState = "multi-state-resolved";
+const string ModelStarikovskiy = "starikovskiy";
+const string ModelCastela = "castela";


This is less clear than just writing these names in the handful of places where they are used.

speth · 2026-07-17T19:47:34Z

+enum class VibModel {
+    Constant,
+    MultiStateResolved,
+    Starikovskiy,
+    Castela
+};


The pattern where this enum is used is needlessly complex. Since it's not in the hot path of the reaction rate evaluation, you would be better off just having the sequential

if (model == "constant") { ... } else if (model == "multi-state-resolved") { ... } ...

block in each case. The only one that needs to handle the exceptional case of an undefined model is setParameters.

speth · 2026-07-17T19:50:25Z

+const AnyMap& getRateConstantMap(const AnyMap& node)
+{
+    if (!node.hasKey("rate-constant")) {
+        throw InputFileError(WhereSetParameters, node,
+            "A vibrational-relaxation reaction requires a 'rate-constant' "
+            "mapping.");
+    }
+
+    const auto& rate = node["rate-constant"];
+
+    if (!rate.is<AnyMap>()) {
+        throw InputFileError(WhereSetParameters, node,
+            "The 'rate-constant' field must be a mapping.");
+    }
+
+    return rate.as<AnyMap>();
+}
+


What's the exception message given if the rate-constant node is missing normally? I'm not sure it's so bad that all this extra checking is worthwhile.

speth · 2026-07-17T22:26:45Z

+    switch (parseVibrationModel(m_vibration_model, node, WhereGetParameters)) {
+    case VibModel::Constant:
+        getConstantParameters(node, rateNode);
+        break;
+    case VibModel::MultiStateResolved:
+        getMultiStateParameters(rateNode);
+        break;
+    case VibModel::Starikovskiy:
+        getStarikovskiyParameters(rateNode);
+        break;
+    case VibModel::Castela:
+        getCastelaParameters(node, rateNode);
+        break;
+    }


These hidden helper functions for different variations on the reaction type are telling me that this should probably be implemented as a family of four different derived classes, using just a single type string to distinguish each one, instead of having both a type and vibration-model as a "sub-type". If this was only true in one place, I'd say we could leave it, but this pattern repeats several times -- here for getParameters, then again for setParameters and setContext (validation). You can leave all the common infrastructure implemented in the intermediate VibrationalRelaxationRate class, at which point the specific classes should just have implementations of these few specific methods.

speth · 2026-07-17T22:29:15Z

+    if (std::abs(m_b - 1.0) > VibTolerance
+        || std::abs(m_D) > VibTolerance
+        || std::abs(m_E) > VibTolerance
+        || std::abs(m_scaling - 1.0) > VibTolerance)
+    {
+        throw InputFileError(WhereGetParameters, node,
+            "Cannot serialize this rate as 'castela': the internal "
+            "parameters are not consistent with the Castela form. "
+            "Expected b = 1, D = 0, E = 0, and scaling = 1.");
+    }


It really shouldn't be necessary to do validation at the output stage -- there shouldn't be a code path by which these parameters take on invalid values. If these parameters had incorrect values before this, then there must have been rate calculations that were being done incorrectly.

Gaetanosaure changed the title ~~Vibrational relaxation~~ Vibrational-relaxation from plasma treatment enhancements Jun 24, 2026

ischoegl added Input Input parsing and conversion (for example, ck2yaml) Kinetics labels Jun 26, 2026

Gaetanosaure added 10 commits July 9, 2026 16:15

Updated TwoTempPlasmaRate class comments and formatting

abd1c09

accept only b-gas as entry for yaml

927ab86

update two-temperature-plasma yaml documentation

0ba7ea0

simplify twoTempPlasmaRate anymap constructor code

085147d

updated cantera.bib

6260a5a

update documentation content and format in VibrationalRelaxationRate.h

f31351d

line length formatting + further formatting of VibrationalRelaxationR…

ad4f41d

…ate.h

modify basic formatting in vibrationalRelaxationRate.cpp

54d996b

Gaetanosaure force-pushed the vibrational-relaxation branch from 9d525fc to 54d996b Compare July 9, 2026 15:14

Gaetanosaure added 12 commits July 9, 2026 18:00

Code formatting and commenting in Plasma phase .h and .cpp

1165179

Add comments and simplify some code bits in vibrationalRelaxationRate…

dc3580d

….cpp

correct mistake made on vibration-model assignment during previous co…

0437f84

…mmit

drop starikovskiy model in-house sign convention to generalise the ex…

c639bb2

…pression completely. The minus sign will be born by YAML coefficients directly.

replace all YAML parameters' underscores with hyphens

53a1605

centralise string handling for often used strings to avoid typos on them

07807c9

make required and forbidden key tests more compact

624c0fc

make data input in setParameters more readable and compact

32cb13b

make sameComposition more compact and readable

133d934

correct typos and finish corrections that were forgotten during previ…

06a051b

…ous commits

Factor vibrational relaxation rate-constant map parsing and correct T…

72d98c6

…iny declaration error

Split vibrational relaxation parameter treatment by model and use swi…

5b2e735

…tch / case

Gaetanosaure added 3 commits July 10, 2026 00:12

split getParameter by model as well

e773816

validateDetailedRelaxationReaction made a little more compact and rea…

162e01e

…dable

update reactions documentation for the vibration-relaxation reaction …

cc10c0c

…type

BangShiuh reviewed Jul 11, 2026

View reviewed changes

Comment thread src/thermo/PlasmaPhase.cpp Outdated

Gaetanosaure added 2 commits July 13, 2026 08:47

correct mispelling in code comments

b50eb2f

added the T-inv optional parameter in TwoTempPlasmaRate, corrected mi…

5338645

…spellings and added getParameters function

Gaetanosaure mentioned this pull request Jul 16, 2026

Plasma treatment enhancements #2121

Closed

speth requested changes Jul 17, 2026

View reviewed changes

		No other rate parameters are allowed for this model. In particular, `b`, `n`, `B`,
		`C`, `D`, `m`, `E`, `z`, and `scaling` are not permitted.


		### Vibrational relaxation: castela

		The `castela` vibrational relaxation rate expression is based on the relaxation time

Uh oh!

Uh oh!

Conversation

Gaetanosaure commented Jun 24, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

codecov Bot commented Jun 25, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Codecov Report

Uh oh!

speth commented Jul 9, 2026

Uh oh!

Gaetanosaure commented Jul 9, 2026

Uh oh!

Gaetanosaure commented Jul 9, 2026

Uh oh!

Gaetanosaure commented Jul 9, 2026

Uh oh!

Uh oh!

BangShiuh commented Jul 11, 2026

Uh oh!

Gaetanosaure commented Jul 13, 2026

Uh oh!

Gaetanosaure commented Jul 13, 2026

Uh oh!

BangShiuh commented Jul 14, 2026

Uh oh!

Gaetanosaure commented Jul 15, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

Gaetanosaure commented Jul 15, 2026

Uh oh!

speth left a comment

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

4 participants

Gaetanosaure commented Jun 24, 2026 •

edited

Loading

codecov Bot commented Jun 25, 2026 •

edited

Loading

Gaetanosaure commented Jul 15, 2026 •

edited

Loading