From 6604443b14e52f6e5ad4888f770f0d3425a3b3d9 Mon Sep 17 00:00:00 2001 From: Gaetanosaure Date: Tue, 1 Sep 2026 15:10:32 +0200 Subject: [PATCH 01/17] modify documentation of "kind" function --- include/cantera/kinetics/ElectronCollisionPlasmaRate.h | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/include/cantera/kinetics/ElectronCollisionPlasmaRate.h b/include/cantera/kinetics/ElectronCollisionPlasmaRate.h index 1cddd48df41..fda8626c664 100644 --- a/include/cantera/kinetics/ElectronCollisionPlasmaRate.h +++ b/include/cantera/kinetics/ElectronCollisionPlasmaRate.h @@ -153,8 +153,10 @@ class ElectronCollisionPlasmaRate : public ReactionRate } //! The kind of the process which will be one of the following: - //! - `"effective"`: A generic effective collision - //! - `"excitation"`: Electronic or vibrational excitation + //! - `"effective"`: Effective momentum-transfer cross section containing + //! elastic and inelastic contributions + //! - `"elastic"`: Elastic momentum-transfer collision + //! - `"excitation"`: Electronic, vibrational or rotational excitation //! - `"ionization"`: Electron-impact ionization //! - `"attachment"`: Electron attachment //! @since New in Cantera 3.2. From 33abcdf8b3751aef715318b3b7f5cabf76e3f6ba Mon Sep 17 00:00:00 2001 From: Gaetanosaure Date: Tue, 1 Sep 2026 15:25:49 +0200 Subject: [PATCH 02/17] update ElectronCollisionPlasmaRate for the new data format --- .../kinetics/ElectronCollisionPlasmaRate.h | 27 ++- src/kinetics/ElectronCollisionPlasmaRate.cpp | 183 ++++++++++++++++-- 2 files changed, 191 insertions(+), 19 deletions(-) diff --git a/include/cantera/kinetics/ElectronCollisionPlasmaRate.h b/include/cantera/kinetics/ElectronCollisionPlasmaRate.h index fda8626c664..31309705682 100644 --- a/include/cantera/kinetics/ElectronCollisionPlasmaRate.h +++ b/include/cantera/kinetics/ElectronCollisionPlasmaRate.h @@ -213,6 +213,19 @@ class ElectronCollisionPlasmaRate : public ReactionRate //! Update the value of #m_crossSectionsInterpolated [m2] void updateInterpolatedCrossSection(span); + //! Name of the electron-collision definition referenced by this rate. + const string& collisionName() const { + return m_collisionName; + } + + //! Return whether tabulated cross-section data have been assigned. + bool hasCrossSectionData() const { + return m_hasCrossSectionData; + } + + //! Assign data from a named entry in the root `electron-collisions` section. + void applyCollisionData(const AnyMap& node); + private: //! The name of the kind of electron collision string m_kind; @@ -224,7 +237,7 @@ class ElectronCollisionPlasmaRate : public ReactionRate string m_product; //! The energy threshold of electron collision - double m_threshold; + double m_threshold = 0.0; //! electron energy levels [eV] vector m_energyLevels; @@ -247,6 +260,18 @@ class ElectronCollisionPlasmaRate : public ReactionRate //! This is used for the calculation of the super-elastic collision reaction //! rate coefficient. Eigen::ArrayXd m_crossSectionsOffset; + + //! Name used to reference the collision data from a reaction. + string m_collisionName; + + //! Whether this rate contains validated tabulated cross-section data. + bool m_hasCrossSectionData = false; + + //! Validate the contents of a named electron-collision definition. + void validateCollisionData(const AnyMap& node) const; + + //! Infer the collision threshold when it is not specified explicitly. + void setDefaultThreshold(); }; } diff --git a/src/kinetics/ElectronCollisionPlasmaRate.cpp b/src/kinetics/ElectronCollisionPlasmaRate.cpp index ce4a1ecbcd6..6312e4f58ed 100644 --- a/src/kinetics/ElectronCollisionPlasmaRate.cpp +++ b/src/kinetics/ElectronCollisionPlasmaRate.cpp @@ -48,32 +48,38 @@ bool ElectronCollisionPlasmaData::update(const ThermoPhase& phase, const Kinetic void ElectronCollisionPlasmaRate::setParameters(const AnyMap& node, const UnitStack& rate_units) { ReactionRate::setParameters(node, rate_units); - if (!node.hasKey("energy-levels") && !node.hasKey("cross-sections")) { - return; + + if (!node.hasKey("collision")) { + throw InputFileError("ElectronCollisionPlasmaRate::setParameters", node, + "Electron-collision reactions require a named 'collision' reference. " + "Tabulated cross-section data must be declared in the root " + "'electron-collisions' section."); } - if (node.hasKey("kind")) { - m_kind = node["kind"].asString(); - } - if (node.hasKey("target")) { - m_target = node["target"].asString(); - } - if (node.hasKey("product")) { - m_product = node["product"].asString(); + m_collisionName = node["collision"].asString(); + + if (m_collisionName.empty()) { + throw InputFileError("ElectronCollisionPlasmaRate::setParameters", node, + "The 'collision' reference cannot be empty."); } - m_energyLevels = node["energy-levels"].asVector(); - m_crossSections = node["cross-sections"].asVector(m_energyLevels.size()); - m_threshold = node.getDouble("threshold", 0.0); + for (const string& key : { + "name", "kind", "target", "product", "threshold", + "energy-levels", "cross-sections" + }) { + if (node.hasKey(key)) { + throw InputFileError("ElectronCollisionPlasmaRate::setParameters", node, + "Electron-collision reaction entries cannot contain '{}'. " + "Collision metadata and tabulated data must be declared in the " + "referenced root 'electron-collisions' entry.", + key); + } + } } void ElectronCollisionPlasmaRate::getParameters(AnyMap& node) const { node["type"] = type(); - node["energy-levels"] = m_energyLevels; - node["cross-sections"] = m_crossSections; - if (!m_kind.empty()) { - node["kind"] = m_kind; - } + node["collision"] = m_collisionName; } void ElectronCollisionPlasmaRate::updateInterpolatedCrossSection( @@ -245,4 +251,145 @@ void ElectronCollisionPlasmaRate::setContext(const Reaction& rxn, const Kinetics } } +void ElectronCollisionPlasmaRate::applyCollisionData(const AnyMap& node) +{ + const string routineName = + "ElectronCollisionPlasmaRate::applyCollisionData"; + + if (!node.hasKey("name")) { + throw InputFileError(routineName, node, + "Electron-collision definitions require a unique 'name'."); + } + + const string name = node["name"].asString(); + + if (name.empty()) { + throw InputFileError(routineName, node, + "Electron-collision definition names cannot be empty."); + } + + if (!m_collisionName.empty() && m_collisionName != name) { + throw InputFileError(routineName, node, + "Reaction references electron collision '{}', but data for '{}' " + "were supplied.", + m_collisionName, name); + } + + m_collisionName = name; + + if (!node.hasKey("kind")) { + throw InputFileError(routineName, node, + "Electron-collision definition '{}' requires 'kind'.", name); + } + + if (!node.hasKey("target")) { + throw InputFileError(routineName, node, + "Electron-collision definition '{}' requires 'target'.", name); + } + + if (!node.hasKey("energy-levels")) { + throw InputFileError(routineName, node, + "Electron-collision definition '{}' requires 'energy-levels'.", + name); + } + + if (!node.hasKey("cross-sections")) { + throw InputFileError(routineName, node, + "Electron-collision definition '{}' requires 'cross-sections'.", + name); + } + + m_kind = node["kind"].asString(); + m_target = node["target"].asString(); + m_product = node.getString("product", ""); + + m_energyLevels = node["energy-levels"].asVector(); + m_crossSections = node["cross-sections"].asVector(m_energyLevels.size()); + + m_threshold = node.getDouble("threshold", 0.0); + setDefaultThreshold(); + validateCollisionData(node); + + // Invalidate all interpolated data after assigning a new table. + m_levelNumber = -3; + m_levelNumberSuperelastic = -2; + m_crossSectionsInterpolated.clear(); + m_crossSectionsOffset.resize(0); + + m_hasCrossSectionData = true; +} + +void ElectronCollisionPlasmaRate::validateCollisionData( + const AnyMap& node) const +{ + const string routineName = "ElectronCollisionPlasmaRate::validateCollisionData"; + + static const set validKinds = { + "effective", "elastic", "excitation", "ionization", "attachment" + }; + + if (!validKinds.count(m_kind)) { + throw InputFileError(routineName, node, + "Unknown electron-collision kind '{}'. Expected one of " + "'effective', 'elastic', 'excitation', 'ionization', or " + "'attachment'.", + m_kind); + } + + if (m_target.empty()) { + throw InputFileError(routineName, node, + "The electron-collision target cannot be empty."); + } + + if (m_energyLevels.size() < 2) { + throw InputFileError(routineName, node, + "Electron-collision data require at least two energy levels."); + } + + if (m_energyLevels.size() != m_crossSections.size()) { + throw InputFileError(routineName, node, + "The 'energy-levels' and 'cross-sections' arrays must have " + "identical lengths."); + } + + for (size_t i = 0; i < m_energyLevels.size(); i++) { + if (!std::isfinite(m_energyLevels[i]) || m_energyLevels[i] < 0.0) { + throw InputFileError(routineName, node, + "Energy levels must be finite and non-negative."); + } + + if (!std::isfinite(m_crossSections[i]) || + m_crossSections[i] < 0.0) { + throw InputFileError(routineName, node, + "Cross sections must be finite and non-negative."); + } + + if (i > 0 && m_energyLevels[i] <= m_energyLevels[i - 1]) { + throw InputFileError(routineName, node, + "Energy levels must be strictly increasing."); + } + } + + if (!std::isfinite(m_threshold) || m_threshold < 0.0) { + throw InputFileError(routineName, node, + "The collision threshold must be finite and non-negative."); + } +} + +void ElectronCollisionPlasmaRate::setDefaultThreshold() +{ + if (m_threshold != 0.0 || + (m_kind != "excitation" && + m_kind != "ionization" && + m_kind != "attachment")) { + return; + } + + for (size_t i = 0; i < m_crossSections.size(); i++) { + if (m_crossSections[i] > 0.0) { + m_threshold = m_energyLevels[i]; + return; + } + } +} } From 1283a3bbc51b04dd76c676b8f42550fddf0491f4 Mon Sep 17 00:00:00 2001 From: Gaetanosaure Date: Tue, 1 Sep 2026 17:36:17 +0200 Subject: [PATCH 03/17] update PlasmaPhase to the new data format --- include/cantera/thermo/PlasmaPhase.h | 14 +++ src/thermo/PlasmaPhase.cpp | 149 +++++++++++++++++++++++---- 2 files changed, 142 insertions(+), 21 deletions(-) diff --git a/include/cantera/thermo/PlasmaPhase.h b/include/cantera/thermo/PlasmaPhase.h index ba8212d8629..4ce464b762a 100644 --- a/include/cantera/thermo/PlasmaPhase.h +++ b/include/cantera/thermo/PlasmaPhase.h @@ -282,6 +282,14 @@ class PlasmaPhase: public IdealGasPhase return m_collisionRates[i]; } + //! Return whether a named electron-collision definition exists. + bool hasElectronCollisionDefinition(const string& name) const { + return m_electronCollisionDefinitions.count(name); + } + + //! Return a named electron-collision definition. + const AnyMap& electronCollisionDefinition(const string& name) const; + //! Update the electron energy distribution. void updateElectronEnergyDistribution(); @@ -1063,6 +1071,12 @@ class PlasmaPhase: public IdealGasPhase //! A function to check that vibrational reservoir species //! are not at risk to hinder phase chemistry. void checkVibrationalReservoirMoleFractions(); + + //! Root-level named electron-collision definitions. + map m_electronCollisionDefinitions; + + //! Names of collision definitions already registered with the EEDF solver. + set m_registeredElectronCollisionNames; }; } diff --git a/src/thermo/PlasmaPhase.cpp b/src/thermo/PlasmaPhase.cpp index 42e337d4fb2..d2a9f91deda 100644 --- a/src/thermo/PlasmaPhase.cpp +++ b/src/thermo/PlasmaPhase.cpp @@ -304,6 +304,19 @@ void PlasmaPhase::setElectronEnergyDistributionParameters(const AnyMap& eedf) } } +const AnyMap& PlasmaPhase::electronCollisionDefinition( + const string& name) const +{ + auto iter = m_electronCollisionDefinitions.find(name); + + if (iter == m_electronCollisionDefinitions.end()) { + throw CanteraError("PlasmaPhase::electronCollisionDefinition", + "Unknown electron-collision definition '{}'.", name); + } + + return iter->second; +} + void PlasmaPhase::setParameters(const AnyMap& phaseNode, const AnyMap& rootNode) { IdealGasPhase::setParameters(phaseNode, rootNode); @@ -318,25 +331,67 @@ void PlasmaPhase::setParameters(const AnyMap& phaseNode, const AnyMap& rootNode) phaseNode["vibrational-reservoir-species-mapping"].asMap(); } + m_electronCollisionDefinitions.clear(); + if (rootNode.hasKey("electron-collisions")) { - for (const auto& item : rootNode["electron-collisions"].asVector()) { - auto rate = make_shared(item); - Composition reactants, products; - reactants[item["target"].asString()] = 1; - reactants[electronSpeciesName()] = 1; - if (item.hasKey("product")) { - products[item["product"].asString()] = 1; + const auto definitions = + rootNode["electron-collisions"].asVector(); + + // Store all definitions before constructing standalone EEDF + // collisions, so reaction references can be resolved independently + // of declaration order. + for (const auto& item : definitions) { + if (!item.hasKey("name")) { + throw InputFileError("PlasmaPhase::setParameters", item, + "Every entry in 'electron-collisions' requires a unique " + "'name'."); + } + + const string name = item["name"].asString(); + + if (name.empty()) { + throw InputFileError("PlasmaPhase::setParameters", item, + "Electron-collision definition names cannot be empty."); + } + + if (m_electronCollisionDefinitions.count(name)) { + throw InputFileError("PlasmaPhase::setParameters", item, + "Duplicate electron-collision definition '{}'.", name); + } + + m_electronCollisionDefinitions.emplace(name, item); + } + + // Every definition contributes to the EEDF, even if no chemical + // reaction references it. + for (const auto& item : definitions) { + auto rate = make_shared(); + rate->applyCollisionData(item); + + Composition reactants; + Composition products; + + reactants[rate->target()] = 1.0; + reactants[electronSpeciesName()] = 1.0; + + if (!rate->product().empty()) { + products[rate->product()] = 1.0; } else { - products[item["target"].asString()] = 1; + products[rate->target()] = 1.0; } - products[electronSpeciesName()] = 1; + + products[electronSpeciesName()] = 1.0; + if (rate->kind() == "ionization") { - products[electronSpeciesName()] += 1; + products[electronSpeciesName()] += 1.0; } else if (rate->kind() == "attachment") { - products[electronSpeciesName()] -= 1; + products.erase(electronSpeciesName()); } - auto R = make_shared(reactants, products, rate); - addCollision(R); + + auto collision = + make_shared(reactants, products, rate); + + addCollision(collision); } } } @@ -580,6 +635,29 @@ void PlasmaPhase::setCollisions() void PlasmaPhase::addCollision(shared_ptr collision) { + // avoid duplications + auto rate = std::dynamic_pointer_cast( + collision->rate()); + + if (!rate) { + throw CanteraError("PlasmaPhase::addCollision", + "Reaction '{}' does not contain an ElectronCollisionPlasmaRate.", + collision->equation()); + } + + if (!rate->hasCrossSectionData()) { + throw CanteraError("PlasmaPhase::addCollision", + "Electron collision '{}' has no tabulated cross-section data.", + rate->collisionName()); + } + + // Root definitions are registered first. A chemical reaction referencing + // the same definition must not add a duplicate EEDF collision. + if (m_registeredElectronCollisionNames.count(rate->collisionName())) { + return; + } + + // Count the number of collisions size_t i = nCollisions(); // setup callback to signal updating the cross-section-related @@ -605,9 +683,16 @@ void PlasmaPhase::addCollision(shared_ptr collision) " collision with equation '{}'", collision->equation()); } + // Check that the target species in the "synthetic" reaction matches the target species + if (target != rate->target()) { + throw CanteraError("PlasmaPhase::addCollision", + "Electron collision '{}' targets '{}', but reaction '{}' uses '{}'.", + rate->collisionName(), rate->target(), + collision->equation(), target); + } + m_collisions.emplace_back(collision); - m_collisionRates.emplace_back( - std::dynamic_pointer_cast(collision->rate())); + m_collisionRates.emplace_back(rate); m_interp_cs_ready.emplace_back(false); // resize parameters @@ -615,8 +700,7 @@ void PlasmaPhase::addCollision(shared_ptr collision) updateInterpolatedCrossSection(i); // Set up data used by Boltzmann solver - auto& rate = *m_collisionRates.back(); - string kind = m_collisionRates.back()->kind(); + string kind = rate->kind(); if ((kind == "effective" || kind == "elastic")) { for (size_t k = 0; k < m_collisions.size() - 1; k++) { @@ -633,9 +717,14 @@ void PlasmaPhase::addCollision(shared_ptr collision) m_kInelastic.push_back(i); } - auto levels = rate.energyLevels(); + auto levels = rate->energyLevels(); m_energyLevels.emplace_back(levels.begin(), levels.end()); - auto sections = rate.crossSections(); + auto sections = rate->crossSections(); + m_crossSections.emplace_back(sections.begin(), sections.end()); + + m_registeredElectronCollisionNames.insert(rate->collisionName()); + m_eedfSolver->setGridCache(); + m_crossSections.emplace_back(sections.begin(), sections.end()); m_eedfSolver->setGridCache(); } @@ -743,10 +832,28 @@ double PlasmaPhase::elasticPowerLoss() updateElasticElectronEnergyLossCoefficients(); // The elastic power loss includes the contributions from inelastic // collisions (inelastic recoil effects). + + vector hasEffectiveCrossSection(nSpecies(), false); + + for (size_t k : m_kElastic) { + if (m_collisionRates[k]->kind() == "effective") { + hasEffectiveCrossSection[m_targetSpeciesIndices[k]] = true; + } + } + double rate = 0.0; for (size_t i = 0; i < nCollisions(); i++) { - rate += concentration(m_targetSpeciesIndices[i]) * - m_elasticElectronEnergyLossCoefficients[i]; + const size_t target = m_targetSpeciesIndices[i]; + const string& kind = m_collisionRates[i]->kind(); + + // The effective cross section already contains all momentum-transfer + // contributions for this target. + if (hasEffectiveCrossSection[target] && kind != "effective") { + continue; + } + + rate += concentration(target) * + m_elasticElectronEnergyLossCoefficients[i]; } const double q_elastic = Avogadro * Avogadro * ElectronCharge * concentration(m_electronSpeciesIndex) * rate; From ed2a91b41315a0dbfff9f82f438d120e882f1760 Mon Sep 17 00:00:00 2001 From: Gaetanosaure Date: Tue, 1 Sep 2026 17:36:45 +0200 Subject: [PATCH 04/17] allow reversible ElectronCollisionPlasma reactions --- src/kinetics/ElectronCollisionPlasmaRate.cpp | 142 ++++++++++++------- 1 file changed, 93 insertions(+), 49 deletions(-) diff --git a/src/kinetics/ElectronCollisionPlasmaRate.cpp b/src/kinetics/ElectronCollisionPlasmaRate.cpp index 6312e4f58ed..35f8215f9cf 100644 --- a/src/kinetics/ElectronCollisionPlasmaRate.cpp +++ b/src/kinetics/ElectronCollisionPlasmaRate.cpp @@ -176,78 +176,122 @@ void ElectronCollisionPlasmaRate::modifyRateConstants( distribution.cwiseProduct(m_crossSectionsOffset)), eps); } -void ElectronCollisionPlasmaRate::setContext(const Reaction& rxn, const Kinetics& kin) +void ElectronCollisionPlasmaRate::setContext( + const Reaction& rxn, const Kinetics& kin) { const ThermoPhase& thermo = kin.thermo(); - // get electron species name - string electronName; - if (thermo.type() == "plasma") { - electronName = dynamic_cast(thermo).electronSpeciesName(); - } else { + + if (thermo.type() != "plasma") { throw CanteraError("ElectronCollisionPlasmaRate::setContext", - "ElectronCollisionPlasmaRate requires plasma phase"); + "ElectronCollisionPlasmaRate requires a plasma phase."); } - // Number of reactants needs to be two + const auto& plasma = dynamic_cast(thermo); + const string electronName = plasma.electronSpeciesName(); + + if (m_collisionName.empty()) { + throw InputFileError("ElectronCollisionPlasmaRate::setContext", rxn.input, + "Electron-collision reaction '{}' does not specify a named " + "'collision' reference.", + rxn.equation()); + } + + if (!plasma.hasElectronCollisionDefinition(m_collisionName)) { + throw InputFileError("ElectronCollisionPlasmaRate::setContext", rxn.input, + "Reaction '{}' references unknown electron collision '{}'.", + rxn.equation(), m_collisionName); + } + + applyCollisionData( + plasma.electronCollisionDefinition(m_collisionName)); + if (rxn.reactants.size() != 2) { throw InputFileError("ElectronCollisionPlasmaRate::setContext", rxn.input, - "ElectronCollisionPlasmaRate requires exactly two reactants"); + "ElectronCollisionPlasmaRate requires exactly two reactants."); } - // Must have only one electron - // @todo add electron-electron collision rate - if (rxn.reactants.at(electronName) != 1) { + auto electronReactant = rxn.reactants.find(electronName); + + if (electronReactant == rxn.reactants.end() || + electronReactant->second != 1.0) { throw InputFileError("ElectronCollisionPlasmaRate::setContext", rxn.input, - "ElectronCollisionPlasmaRate requires one and only one electron"); - } - - // Determine the "kind" of collision if not specified explicitly - if (m_kind.empty()) { - m_kind = "excitation"; // default - if (rxn.reactants == rxn.products) { - m_kind = "effective"; - } else { - for (const auto& [p, stoich] : rxn.products) { - if (p == electronName) { - continue; - } - double q = thermo.charge(thermo.speciesIndex(p, true)); - if (q > 0) { - m_kind = "ionization"; - } else if (q < 0) { - m_kind = "attachment"; - } - } - } + "ElectronCollisionPlasmaRate requires exactly one electron " + "reactant."); } - if (m_threshold == 0.0 && - (m_kind == "excitation" || m_kind == "ionization" || m_kind == "attachment")) - { - for (size_t i = 0; i < m_energyLevels.size(); i++) { - if (m_energyLevels[i] > 0.0) { // Look for first non-zero cross-section - m_threshold = m_energyLevels[i]; - break; - } + string reactionTarget; + + for (const auto& [name, coefficient] : rxn.reactants) { + if (name != electronName && coefficient != 0.0) { + reactionTarget = name; + break; } } + if (reactionTarget != m_target) { + throw InputFileError("ElectronCollisionPlasmaRate::setContext", rxn.input, + "Electron collision '{}' targets '{}', but reaction '{}' uses '{}'.", + m_collisionName, m_target, rxn.equation(), reactionTarget); + } + + double reactantElectrons = electronReactant->second; + double productElectrons = 0.0; + + auto electronProduct = rxn.products.find(electronName); + if (electronProduct != rxn.products.end()) { + productElectrons = electronProduct->second; + } + + string reactionKind; + + if (productElectrons > reactantElectrons) { + reactionKind = "ionization"; + } else if (productElectrons < reactantElectrons) { + reactionKind = "attachment"; + } else if (rxn.reactants == rxn.products) { + reactionKind = "effective"; + } else { + reactionKind = "excitation"; + } + + bool compatibleKind = reactionKind == m_kind; + + // Elastic and effective cross sections both correspond to an unchanged + // electron-target reaction. + if (reactionKind == "effective" && + (m_kind == "effective" || m_kind == "elastic")) { + compatibleKind = true; + } + + // Allow an unresolved inelastic channel whose product is not represented + // as a phase species. The reaction is compositionally unchanged, while the + // collision remains inelastic in the Boltzmann equation. + if (reactionKind == "effective" && m_kind == "excitation") { + compatibleKind = true; + } + + if (!compatibleKind) { + throw InputFileError("ElectronCollisionPlasmaRate::setContext", rxn.input, + "Electron collision '{}' has kind '{}', but reaction '{}' implies " + "kind '{}'.", + m_collisionName, m_kind, rxn.equation(), reactionKind); + } + if (!rxn.reversible) { - return; // end checking of forward reaction + return; } - // For super-elastic collisions if (rxn.products.size() != 2) { throw InputFileError("ElectronCollisionPlasmaRate::setContext", rxn.input, - "ElectronCollisionPlasmaRate requires exactly two products" - " if the reaction is reversible (super-elastic collisions)"); + "ElectronCollisionPlasmaRate requires exactly two products for " + "a reversible excitation reaction."); } - // Must have only one electron - if (rxn.products.at(electronName) != 1) { + if (electronProduct == rxn.products.end() || + electronProduct->second != 1.0) { throw InputFileError("ElectronCollisionPlasmaRate::setContext", rxn.input, - "ElectronCollisionPlasmaRate requires one and only one electron in products" - " if the reaction is reversible (super-elastic collisions)"); + "A reversible electron-collision reaction requires exactly one " + "electron product."); } } From 83267abe3cced1a633ea89e1bceb46737af3da94 Mon Sep 17 00:00:00 2001 From: Gaetanosaure Date: Tue, 1 Sep 2026 17:37:47 +0200 Subject: [PATCH 05/17] fix the elastic/effective incoherence in EEDFTTA --- src/thermo/EEDFTwoTermApproximation.cpp | 100 ++++++++++++++++++++---- 1 file changed, 85 insertions(+), 15 deletions(-) diff --git a/src/thermo/EEDFTwoTermApproximation.cpp b/src/thermo/EEDFTwoTermApproximation.cpp index cbb6afd15b5..72c367c6ee9 100644 --- a/src/thermo/EEDFTwoTermApproximation.cpp +++ b/src/thermo/EEDFTwoTermApproximation.cpp @@ -737,34 +737,104 @@ void EEDFTwoTermApproximation::calculateTotalCrossSection() { m_totalCrossSectionCenter.assign(m_points, 0.0); m_totalCrossSectionEdge.assign(m_points + 1, 0.0); + + vector hasEffectiveCrossSection( + m_phase->nSpecies(), false); + + for (size_t k : m_phase->kElastic()) { + if (m_phase->collisionRate(k)->kind() == "effective") { + hasEffectiveCrossSection[m_phase->targetIndex(k)] = true; + } + } + for (size_t k = 0; k < m_phase->nCollisions(); k++) { - auto x = m_phase->collisionRate(k)->energyLevels(); - auto y = m_phase->collisionRate(k)->crossSections(); + const size_t target = m_phase->targetIndex(k); + const string& kind = m_phase->collisionRate(k)->kind(); + + // An effective cross section already contains the elastic and + // inelastic contributions for its target. Adding the individual + // inelastic channels would count them twice. + if (hasEffectiveCrossSection[target] && kind != "effective") { + continue; + } + + auto levels = m_phase->collisionRate(k)->energyLevels(); + auto sections = m_phase->collisionRate(k)->crossSections(); + const double moleFraction = + m_X_targets[m_klocTargets[k]]; for (size_t i = 0; i < m_points; i++) { - m_totalCrossSectionCenter[i] += m_X_targets[m_klocTargets[k]] * - linearInterp(m_gridCenter[i], x, y); + m_totalCrossSectionCenter[i] += moleFraction * + linearInterp(m_gridCenter[i], levels, sections); } + for (size_t i = 0; i < m_points + 1; i++) { - m_totalCrossSectionEdge[i] += m_X_targets[m_klocTargets[k]] * - linearInterp(m_gridEdge[i], x, y); + m_totalCrossSectionEdge[i] += moleFraction * + linearInterp(m_gridEdge[i], levels, sections); } } } void EEDFTwoTermApproximation::calculateTotalElasticCrossSection() { - m_sigmaElastic.clear(); - m_sigmaElastic.resize(m_points, 0.0); + m_sigmaElastic.assign(m_points, 0.0); + for (size_t k : m_phase->kElastic()) { - auto x = m_phase->collisionRate(k)->energyLevels(); - auto y = m_phase->collisionRate(k)->crossSections(); - // Note: - // moleFraction(m_kTargets[k]) <=> m_X_targets[m_klocTargets[k]] - double mass_ratio = ElectronMass / (m_phase->molecularWeight(m_kTargets[k]) / Avogadro); + auto rate = m_phase->collisionRate(k); + auto levels = rate->energyLevels(); + auto sections = rate->crossSections(); + + const size_t target = m_phase->targetIndex(k); + const double massRatio = + ElectronMass / + (m_phase->molecularWeight(target) / Avogadro); + + const double moleFraction = + m_X_targets[m_klocTargets[k]]; + for (size_t i = 0; i < m_points; i++) { - m_sigmaElastic[i] += 2.0 * mass_ratio * m_X_targets[m_klocTargets[k]] * - linearInterp(m_gridEdge[i], x, y); + const double energy = m_gridEdge[i]; + + double elasticCrossSection = + linearInterp(energy, levels, sections); + + if (rate->kind() == "effective") { + // By definition: + // + // effective = elastic momentum transfer + // + sum(forward inelastic cross sections) + // + // Reconstruct the elastic contribution directly on the EEDF + // grid to avoid interpolation through an intermediate grid. + for (size_t kInelastic : m_phase->kInelastic()) { + if (m_phase->targetIndex(kInelastic) != target) { + continue; + } + + auto inelasticRate = + m_phase->collisionRate(kInelastic); + + elasticCrossSection -= linearInterp( + energy, + inelasticRate->energyLevels(), + inelasticRate->crossSections()); + } + + if (elasticCrossSection < 0.0) { + throw CanteraError( + "EEDFTwoTermApproximation::" + "calculateTotalElasticCrossSection", + "Reconstructed elastic cross section for target '{}' " + "is negative at {} eV: {} m^2. The effective and " + "inelastic cross-section data are inconsistent.", + m_phase->speciesName(target), + energy, elasticCrossSection); + } + } + + m_sigmaElastic[i] += + 2.0 * massRatio * moleFraction * + elasticCrossSection; } } } From fc52f7c6ea383b455d83408b3ff29ed8cd308ffc Mon Sep 17 00:00:00 2001 From: Gaetanosaure Date: Tue, 1 Sep 2026 17:38:08 +0200 Subject: [PATCH 06/17] update format of example data files --- test/data/air-plasma.yaml | 125 ++++++++++++++++++++++++----------- test/data/oxygen-plasma.yaml | 12 +++- 2 files changed, 96 insertions(+), 41 deletions(-) diff --git a/test/data/air-plasma.yaml b/test/data/air-plasma.yaml index 4ff6bb09284..1f746830808 100644 --- a/test/data/air-plasma.yaml +++ b/test/data/air-plasma.yaml @@ -19,6 +19,21 @@ phases: reactions: - equation: O2 + Electron => Electron + Electron + O2+ type: electron-collision-plasma + collision: phelps-O2-ionization + +- equation: N2 + Electron => N2+ + 2 Electron + type: electron-collision-plasma + collision: phelps-N2-ionization + +- equation: O2 + Electron => O2- + type: electron-collision-plasma + collision: phelps-O2-attachment + +electron-collisions: +- name: phelps-O2-ionization + target: O2 + product: O2+ + kind: ionization energy-levels: [12.06, 13.0, 18.0, 28.0, 38.0, 48.0, 58.0, 68.0, 78.0, 88.0, 100.0, 150.0, 200.0, 300.0, 500.0, 700.0, 1000.0, 1500.0, 2000.0, 3000.0, 5000.0, 7000.0, 10000.0] @@ -26,8 +41,11 @@ reactions: 2.5e-20, 2.6e-20, 2.7e-20, 2.7e-20, 2.5e-20, 2.17e-20, 1.66e-20, 1.35e-20, 1.04e-20, 7.6e-21, 6e-21, 4.2e-21, 2.7e-21, 2e-21, 1.4e-21] threshold: 12.06 -- equation: N2 + Electron => N2+ + 2 Electron - type: electron-collision-plasma + +- name: phelps-N2-ionization + target: N2 + product: N2+ + kind: ionization energy-levels: [0.0, 15.6, 16.0, 16.5, 17.0, 17.5, 18.0, 18.5, 19.0, 19.5, 20.0, 21.0, 22.0, 23.0, 25.0, 30.0, 34.0, 45.0, 60.0, 75.0, 100.0, 150.0, 200.0] cross-sections: [0.0, 0.0, 1.95e-22, 4.28e-22, 6.6e-22, 9.11e-22, 1.2e-21, 1.516e-21, @@ -35,8 +53,11 @@ reactions: 9.579e-21, 1.1718e-20, 1.6461e-20, 2.0181e-20, 2.2134e-20, 2.3436e-20, 2.2692e-20, 2.1018e-20] threshold: 15.6 -- equation: O2 + Electron => O2- - type: electron-collision-plasma + +- name: phelps-O2-attachment + target: O2 + product: O2- + kind: attachment energy-levels: [0.0, 0.058, 0.073, 0.083, 0.089, 0.095, 0.103, 0.109, 0.15, 0.17, 0.2, 0.21, 0.23, 0.32, 0.33, 0.35, 0.44, 0.45, 0.47, 0.56, 0.57, 0.59, 0.68, 0.69, 0.71, 0.79, 0.8, 0.82, 0.9, 0.91, 0.93, 1.02, 1.03, 1.05, 1.5, 100.0] @@ -46,25 +67,30 @@ reactions: 0.0, 0.0, 5.50781e-42, 0.0, 0.0, 4.20596e-42, 0.0, 0.0, 0.0] threshold: 0.0 -electron-collisions: -- target: N2 +- name: phelps-N2-effective + target: N2 + kind: effective energy-levels: [0.0, 0.015, 0.03, 0.05, 0.1, 0.15, 0.2, 0.3, 0.4, 0.7, 1.2, 1.5, 1.9, 2.2, 2.8, 3.3, 4.0, 5.0, 7.0, 10.0, 15.0, 20.0, 30.0, 75.0, 150.0] cross-sections: [1.1e-20, 2.55e-20, 3.4e-20, 4.33e-20, 5.95e-20, 7.1e-20, 7.9e-20, 9e-20, 9.7e-20, 1e-19, 1.04e-19, 1.2e-19, 1.96e-19, 2.85e-19, 2.8e-19, 1.72e-19, 1.26e-19, 1.09e-19, 1.01e-19, 1.04e-19, 1.1e-19, 1.02e-19, 9e-20, 6.6e-20, 4.9e-20] - kind: effective -- target: N2 + +- name: phelps-N2-rotational-excitation + target: N2 product: N2(rot) + kind: excitation energy-levels: [0.02, 0.03, 0.4, 0.8, 1.2, 1.6, 1.7, 1.8, 1.9, 2.0, 2.1, 2.2, 2.3, 2.4, 2.5, 2.6, 2.7, 2.8, 2.9, 3.0, 3.1, 3.2, 3.3, 3.6, 5.0] cross-sections: [0.0, 2.5e-22, 2.5e-22, 2.5e-22, 4.7e-22, 8.6e-22, 1.5e-21, 2.35e-21, 1.08e-20, 1.9e-20, 2.03e-20, 2.77e-20, 2.5e-20, 2.19e-20, 2.4e-20, 2.17e-20, 1.62e-20, 1.38e-20, 1.18e-20, 1.03e-20, 8.4e-21, 6.9e-21, 5e-21, 1.7e-21, 0.0] threshold: 0.02 - kind: excitation -- target: N2 + +- name: phelps-N2-v1-excitation + target: N2 product: N2(v1) + kind: excitation energy-levels: [0.29, 0.3, 0.33, 0.4, 0.75, 0.9, 1.0, 1.1, 1.16, 1.2, 1.22, 1.4, 1.5, 1.6, 1.65, 3.6, 4.0, 5.0, 15.0, 18.0, 20.0, 22.0, 23.0, 25.0, 29.0, 32.0, 50.0, 80.0] @@ -72,9 +98,11 @@ electron-collisions: 1.1e-22, 1.25e-22, 1.35e-22, 7e-22, 1e-21, 1.5e-21, 0.0, 0.0, 5.5e-22, 3.5e-22, 3.5e-22, 4e-22, 6.5e-22, 8.5e-22, 8.5e-22, 6e-22, 3e-22, 1.5e-22, 1.2e-22, 0.0] threshold: 0.29 - kind: excitation -- target: N2 + +- name: phelps-N2-v1res-excitation + target: N2 product: N2(v1res) + kind: excitation energy-levels: [0.0, 0.291, 1.6, 1.65, 1.7, 1.8, 1.9, 2.0, 2.1, 2.2, 2.3, 2.4, 2.5, 2.6, 2.7, 2.75, 2.8, 2.9, 3.0, 3.1, 3.2, 3.3, 3.4, 3.5, 3.6, 4.0, 100.0] cross-sections: [0.0, 0.0, 0.0, 2.7e-21, 3.15e-21, 5.4e-21, 1.485e-20, 4.8e-20, @@ -82,27 +110,33 @@ electron-collisions: 1.35e-20, 5.25e-21, 8.7e-21, 1.17e-20, 8.55e-21, 6.6e-21, 6e-21, 5.85e-21, 5.7e-21, 0.0, 0.0] threshold: 0.291 - kind: excitation -- target: N2 + +- name: phelps-N2-v2-excitation + target: N2 product: N2(v2) + kind: excitation energy-levels: [0.0, 0.59, 1.7, 1.8, 1.9, 2.0, 2.1, 2.2, 2.3, 2.4, 2.5, 2.6, 2.7, 2.75, 2.8, 2.9, 3.0, 3.1, 3.2, 3.3, 3.4, 3.5, 3.6, 100.0] cross-sections: [0.0, 0.0, 0.0, 1.5e-22, 6.3e-21, 1.935e-20, 3.3e-20, 1.47e-20, 5.4e-21, 2.115e-20, 3e-20, 5.4e-21, 1.05e-20, 1.725e-20, 1.275e-20, 3.3e-21, 9e-21, 6.45e-21, 3.75e-21, 3.45e-21, 3e-21, 2.13e-21, 0.0, 0.0] threshold: 0.59 - kind: excitation -- target: N2 + +- name: phelps-N2-v3-excitation + target: N2 product: N2(v3) + kind: excitation energy-levels: [0.0, 0.88, 1.9, 2.0, 2.1, 2.2, 2.3, 2.4, 2.5, 2.6, 2.7, 2.75, 2.8, 2.9, 3.0, 3.1, 3.2, 3.3, 3.4, 100.0] cross-sections: [0.0, 0.0, 0.0, 9.6e-21, 2.055e-20, 2.7e-20, 1.695e-20, 7.5e-22, 9.6e-21, 1.47e-20, 4.5e-21, 9.6e-21, 5.4e-21, 8.55e-21, 4.05e-21, 2.82e-21, 2.91e-21, 6.15e-22, 0.0, 0.0] threshold: 0.88 - kind: excitation -- target: N2 + +- name: phelps-N2-C3-excitation + target: N2 product: N2(C3) + kind: excitation energy-levels: [11.03, 11.5, 12.0, 12.5, 13.0, 13.5, 13.8, 14.0, 14.2, 14.5, 15.0, 16.0, 17.0, 18.0, 19.0, 20.0, 22.0, 24.0, 26.0, 28.0, 30.0, 36.0, 40.0, 50.0, 70.0, 100.0, 150.0] @@ -111,52 +145,64 @@ electron-collisions: 1.77e-21, 1.5e-21, 1.28e-21, 1.11e-21, 7.8e-22, 6.3e-22, 3.9e-22, 1.5e-22, 1.5e-23, 0.0] threshold: 11.03 - kind: excitation -- target: N2 + +- name: phelps-N2-B3-excitation + target: N2 product: N2(B3) + kind: excitation energy-levels: [0.0, 7.35, 8.0, 9.0, 10.0, 11.0, 12.0, 13.0, 14.0, 15.0, 16.0, 17.0, 18.0, 20.0, 22.0, 26.0, 30.0, 34.0, 40.0, 50.0, 70.0, 150.0] cross-sections: [0.0, 0.0, 3.62e-22, 9.38e-22, 1.508e-21, 1.863e-21, 2.003e-21, 1.99e-21, 1.816e-21, 1.615e-21, 1.447e-21, 1.307e-21, 1.199e-21, 1.112e-21, 9.51e-22, 8.04e-22, 6.77e-22, 5.63e-22, 4.29e-22, 2.68e-22, 6.7e-23, 0.0] threshold: 7.35 - kind: excitation -- target: N2 + +- name: phelps-N2-a1-excitation + target: N2 product: N2(a1) + kind: excitation energy-levels: [0.0, 8.55, 9.0, 14.0, 15.0, 16.0, 17.0, 18.0, 19.0, 20.0, 24.0, 26.0, 30.0, 40.0, 50.0, 70.0, 100.0, 150.0, 200.0] cross-sections: [0.0, 0.0, 1.27e-22, 1.474e-21, 1.715e-21, 1.916e-21, 2.023e-21, 1.99e-21, 1.923e-21, 1.849e-21, 1.621e-21, 1.528e-21, 1.367e-21, 1.065e-21, 8.51e-22, 6.03e-22, 4.02e-22, 2.68e-22, 2.01e-22] threshold: 8.55 - kind: excitation -- target: N2 + +- name: phelps-N2-w1-excitation + target: N2 product: N2(w1) + kind: excitation energy-levels: [0.0, 8.89, 9.0, 10.0, 11.0, 12.0, 13.0, 14.0, 15.0, 16.0, 17.0, 18.0, 20.0, 22.0, 30.0, 38.0, 50.0, 150.0] cross-sections: [0.0, 0.0, 1.3e-23, 2.61e-22, 4.76e-22, 6.63e-22, 7.84e-22, 7.71e-22, 6.7e-22, 5.43e-22, 4.42e-22, 3.75e-22, 2.88e-22, 2.41e-22, 1.54e-22, 9.4e-23, 4.7e-23, 0.0] threshold: 8.89 - kind: excitation -- target: O2 + +- name: phelps-O2-dissociative-attachment + target: O2 + product: O^-+O + kind: attachment energy-levels: [4.4, 4.9, 5.38, 5.86, 6.1, 6.48, 6.77, 7.05, 7.3, 7.53, 7.77, 8.0, 8.25, 8.73, 9.2, 9.68, 10.15, 11.35, 100.0] cross-sections: [0.0, 0.0, 2.3e-23, 7.2e-23, 1.08e-22, 1.38e-22, 1.52e-22, 1.56e-22, 1.48e-22, 1.31e-22, 1.1e-22, 8.4e-23, 5.4e-23, 2.8e-23, 1.4e-23, 8e-24, 8e-24, 8e-24, 0.0] threshold: 0.0 - product: O^-+O - kind: attachment -- target: O2 + +- name: phelps-O2-effective + target: O2 + kind: effective energy-levels: [0.0, 0.015, 0.03, 0.05, 0.1, 0.15, 0.2, 0.3, 0.4, 0.7, 1.2, 1.5, 1.9, 2.2, 2.8, 3.3, 4.0, 5.0, 7.0, 10.0, 15.0, 20.0, 30.0, 75.0, 150.0, 200.0] cross-sections: [3.5e-21, 8.7e-21, 1.24e-20, 1.6e-20, 2.5e-20, 3.1e-20, 3.6e-20, 4.5e-20, 5.2e-20, 6.1e-20, 7.9e-20, 7.6e-20, 6.9e-20, 6.5e-20, 5.8e-20, 5.5e-20, 5.5e-20, 5.6e-20, 6.6e-20, 8e-20, 8.8e-20, 8.6e-20, 8e-20, 6.8e-20, 6.7e-20, 6e-20] - kind: effective -- target: O2 + +- name: phelps-O2-rotational-excitation + target: O2 product: O2(rot) + kind: excitation energy-levels: [0.07, 0.08, 0.1, 0.2, 0.21, 0.22, 0.32, 0.33, 0.35, 0.44, 0.45, 0.47, 0.56, 0.57, 0.59, 0.68, 0.69, 0.71, 0.79, 0.8, 0.81, 0.9, 0.91, 0.93, 1.02, 1.03, 1.05, 1.13, 1.14, 1.16, 1.22, 1.23, 1.26, 1.34, 1.35, 1.37, 1.44, 1.45, 1.47, 1.54, @@ -166,18 +212,21 @@ electron-collisions: 8.4e-22, 0.0, 0.0, 7.2e-22, 0.0, 0.0, 4.68e-22, 0.0, 0.0, 6e-22, 0.0, 0.0, 3.6e-22, 0.0, 0.0, 2.4e-22, 0.0, 0.0, 1.2e-22, 0.0, 0.0, 4.8e-23, 0.0] threshold: 0.02 - kind: excitation -- target: O2 + +- name: phelps-O2-a1-excitation + target: O2 product: O2(a1) + kind: excitation energy-levels: [0.977, 1.5, 3.5, 5.62, 6.53, 7.89, 13.0, 20.5, 41.0, 100.0] cross-sections: [0.0, 5.8e-23, 4.9e-22, 8.25e-22, 9.08e-22, 8.63e-22, 5.27e-22, 3.24e-22, 1.37e-22, 0.0] threshold: 0.977 - kind: excitation -- target: O2 + +- name: phelps-O2-b1-excitation + target: O2 product: O2(b1) - energy-levels: [1.627, 3.0, 4.0, 7.34, 9.26, 13.0, 17.0, 20.7, 24.0, 35.1, 45.1, 100] + kind: excitation + energy-levels: [1.627, 3.0, 4.0, 7.34, 9.26, 13.0, 17.0, 20.7, 24.0, 35.1, 45.1, 100.0] cross-sections: [0.0, 9.7e-23, 1.49e-22, 1.91e-22, 1.74e-22, 1.3e-22, 1.3e-22, 1.25e-22, 1e-22, 6.3e-23, 5e-24, 0.0] - threshold: 1.627 - kind: excitation + threshold: 1.627 \ No newline at end of file diff --git a/test/data/oxygen-plasma.yaml b/test/data/oxygen-plasma.yaml index 9c172c208bb..f4f8cc5abb5 100644 --- a/test/data/oxygen-plasma.yaml +++ b/test/data/oxygen-plasma.yaml @@ -34,7 +34,7 @@ phases: type: discretized energy-levels: [0.0, 0.1, 1.0, 10.0] distribution: [0.0, 0.2, 0.7, 0.01] - normalize: False + normalize: false species: - name: E @@ -53,7 +53,13 @@ reactions: - equation: O2 + E <=> E + O2 type: electron-collision-plasma - note: This is a electron collision process of plasma + collision: O2-effective + note: This is an electron collision process of plasma. + +electron-collisions: +- name: O2-effective + target: O2 + kind: effective energy-levels: [0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0] cross-sections: [0.0, 5.97e-20, 6.45e-20, 6.74e-20, 6.93e-20, 7.2e-20, - 7.52e-20, 7.86e-20, 8.21e-20, 8.49e-20, 8.8e-20] + 7.52e-20, 7.86e-20, 8.21e-20, 8.49e-20, 8.8e-20] \ No newline at end of file From ca478060ee9abbf7fd1efe629b80986eb9572f4c Mon Sep 17 00:00:00 2001 From: Gaetanosaure Date: Tue, 1 Sep 2026 17:38:31 +0200 Subject: [PATCH 07/17] update kineticsFromYaml test --- test/kinetics/kineticsFromYaml.cpp | 37 +++++++++++++++++++++--------- 1 file changed, 26 insertions(+), 11 deletions(-) diff --git a/test/kinetics/kineticsFromYaml.cpp b/test/kinetics/kineticsFromYaml.cpp index 43f0d9b61e7..202a63940c5 100644 --- a/test/kinetics/kineticsFromYaml.cpp +++ b/test/kinetics/kineticsFromYaml.cpp @@ -643,25 +643,40 @@ TEST(Reaction, TwoTempPlasmaExtendedListFromYaml) TEST(Reaction, ElectronCollisionPlasmaFromYaml) { - auto sol = newSolution("oxygen-plasma.yaml", "discretized-electron-energy-plasma", "none"); + auto sol = newSolution( + "oxygen-plasma.yaml", + "discretized-electron-energy-plasma", + "none" + ); AnyMap rxn = AnyMap::fromYamlString( "{equation: O2 + E => E + O2," " type: electron-collision-plasma," - " energy-levels: [0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0]," - " cross-sections: [0.0, 5.97e-20, 6.45e-20, 6.74e-20, 6.93e-20, 7.2e-20, " - " 7.52e-20, 7.86e-20, 8.21e-20, 8.49e-20, 8.8e-20]}"); - - auto R = newReaction(rxn, *(sol->kinetics())); + " collision: O2-effective}" + ); + auto R = newReaction(rxn, *sol->kinetics()); EXPECT_EQ(R->reactants.at("O2"), 1); EXPECT_EQ(R->reactants.at("E"), 1); EXPECT_EQ(R->products.at("O2"), 1); EXPECT_EQ(R->products.at("E"), 1); - const auto rate = std::dynamic_pointer_cast(R->rate()); - - for (size_t k = 0; k < rate->energyLevels().size(); k++) { - EXPECT_DOUBLE_EQ(rate->energyLevels()[k], rxn["energy-levels"].asVector()[k]); - EXPECT_DOUBLE_EQ(rate->crossSections()[k], rxn["cross-sections"].asVector()[k]); + ASSERT_NE(rate, nullptr); + EXPECT_EQ(rate->collisionName(), "O2-effective"); + EXPECT_EQ(rate->kind(), "effective"); + const vector expectedLevels = { + 0.0, 1.0, 2.0, 3.0, 4.0, 5.0, + 6.0, 7.0, 8.0, 9.0, 10.0 + }; + const vector expectedCrossSections = { + 0.0, 5.97e-20, 6.45e-20, 6.74e-20, 6.93e-20, 7.2e-20, + 7.52e-20, 7.86e-20, 8.21e-20, 8.49e-20, 8.8e-20 + }; + const auto energyLevels = rate->energyLevels(); + const auto crossSections = rate->crossSections(); + ASSERT_EQ(energyLevels.size(), expectedLevels.size()); + ASSERT_EQ(crossSections.size(), expectedCrossSections.size()); + for (size_t k = 0; k < expectedLevels.size(); k++) { + EXPECT_DOUBLE_EQ(energyLevels[k], expectedLevels[k]); + EXPECT_DOUBLE_EQ(crossSections[k], expectedCrossSections[k]); } } From 6e59463d72e8eb74fc8e5b05cba73e073503b4db Mon Sep 17 00:00:00 2001 From: Gaetanosaure Date: Wed, 2 Sep 2026 15:05:14 +0200 Subject: [PATCH 08/17] update round-trip writing with the new format --- include/cantera/thermo/PlasmaPhase.h | 10 +++++++++ src/base/YamlWriter.cpp | 32 ++++++++++++++++++++++++++++ 2 files changed, 42 insertions(+) diff --git a/include/cantera/thermo/PlasmaPhase.h b/include/cantera/thermo/PlasmaPhase.h index 4ce464b762a..5f670bcc0f8 100644 --- a/include/cantera/thermo/PlasmaPhase.h +++ b/include/cantera/thermo/PlasmaPhase.h @@ -155,6 +155,16 @@ class PlasmaPhase: public IdealGasPhase void getParameters(AnyMap& phaseNode) const override; void setParameters(const AnyMap& phaseNode, const AnyMap& rootNode=AnyMap()) override; + + //! Return the named electron-collision definitions used by this phase. + /*! + * The returned map associates each collision name with its root-level + * YAML definition. + */ + const map& electronCollisionDefinitions() const { + return m_electronCollisionDefinitions; + } + //! @} //! @name Electron Species Information //! @{ diff --git a/src/base/YamlWriter.cpp b/src/base/YamlWriter.cpp index ae96c5cd7c5..d8a3e41937e 100644 --- a/src/base/YamlWriter.cpp +++ b/src/base/YamlWriter.cpp @@ -6,6 +6,7 @@ #include "cantera/base/Solution.h" #include "cantera/base/stringUtils.h" #include "cantera/thermo/ThermoPhase.h" +#include "cantera/thermo/PlasmaPhase.h" #include "cantera/thermo/Species.h" #include "cantera/kinetics/Kinetics.h" #include "cantera/kinetics/Reaction.h" @@ -205,6 +206,37 @@ string YamlWriter::toYamlString() const } } + // Build electron-collision definitions used by plasma phases + vector electronCollisionDefs; + std::unordered_map electronCollisionDefIndex; + + for (const auto& phase : m_phases) { + auto plasma = std::dynamic_pointer_cast(phase->thermo()); + if (!plasma) { + continue; + } + + for (const auto& [name, collisionDef] : + plasma->electronCollisionDefinitions()) { + auto iter = electronCollisionDefIndex.find(name); + if (iter == electronCollisionDefIndex.end()) { + electronCollisionDefs.emplace_back(collisionDef); + electronCollisionDefIndex[name] = + electronCollisionDefs.size() - 1; + } else if (electronCollisionDefs[iter->second] != collisionDef) { + throw CanteraError("YamlWriter::toYamlString", + "Multiple electron collisions with name '{}' and different " + "definitions are not supported:\n>>>>>>\n{}\n======\n{}\n<<<<<<\n", + name, collisionDef.toYamlString(), + electronCollisionDefs[iter->second].toYamlString()); + } + } + } + + if (!electronCollisionDefs.empty()) { + output["electron-collisions"] = std::move(electronCollisionDefs); + } + output.setMetadata("precision", AnyValue(m_float_precision)); output.setUnits(m_output_units); return output.toYamlString(); From f938adb9c9cde351efb43490b98478d0a869cfba Mon Sep 17 00:00:00 2001 From: Gaetanosaure Date: Wed, 2 Sep 2026 15:25:40 +0200 Subject: [PATCH 09/17] update lxcat2yaml for compatibility with the old LXCat xml format with respect to the new cross-sections format --- interfaces/cython/cantera/lxcat2yaml.py | 583 ++++++++++++++++++------ 1 file changed, 432 insertions(+), 151 deletions(-) diff --git a/interfaces/cython/cantera/lxcat2yaml.py b/interfaces/cython/cantera/lxcat2yaml.py index eda3b873da0..b810262140a 100644 --- a/interfaces/cython/cantera/lxcat2yaml.py +++ b/interfaces/cython/cantera/lxcat2yaml.py @@ -24,6 +24,7 @@ from __future__ import annotations import argparse +import re import sys import textwrap import xml.etree.ElementTree as etree @@ -47,24 +48,73 @@ BlockMap: type[CommentedMap] = CommentedMap -class Process: - """A class of YAML data for collision of a target species""" - def __init__(self, equation: str, energy_levels: list[float], cross_sections: list[float]) -> None: - self.equation = equation +class ElectronCollision: + """YAML representation of an electron-collision cross section.""" + + def __init__( + self, + name: str, + target: str, + product: str | None, + kind: str, + threshold: float, + energy_levels: list[float], + cross_sections: list[float], + ) -> None: + self.name = name + self.target = target + self.product = product + self.kind = kind + self.threshold = threshold self.energy_levels = energy_levels self.cross_sections = cross_sections + @classmethod + def to_yaml( + cls, + representer: SafeRepresenter, + node: ElectronCollision, + ) -> MappingNode: + out = BlockMap([ + ("name", node.name), + ("target", node.target), + ]) + + if node.product is not None: + out["product"] = node.product + + out["kind"] = node.kind + out["threshold"] = node.threshold + out["energy-levels"] = node.energy_levels + out["cross-sections"] = node.cross_sections + + return representer.represent_dict(out) + + +class Process: + """YAML representation of a reaction referencing collision data.""" + + def __init__(self, equation: str, collision: str) -> None: + self.equation = equation + self.collision = collision + self.duplicate = False + @classmethod def to_yaml(cls, representer: SafeRepresenter, node: Process) -> MappingNode: - out = BlockMap([('equation', node.equation), - ('type', 'electron-collision-plasma'), - ('energy-levels', node.energy_levels), - ('cross-sections', node.cross_sections), - ]) + out = BlockMap([ + ("equation", node.equation), + ("type", "electron-collision-plasma"), + ("collision", node.collision), + ]) + + if node.duplicate: + out["duplicate"] = True + return representer.represent_dict(out) # Define YAML emitter emitter: yaml.YAML = yaml.YAML() +emitter.register_class(ElectronCollision) emitter.register_class(Process) # Return indices of a child name @@ -79,6 +129,93 @@ def Flowlist(*args: Iterable[_VT], **kwargs: _VT) -> list[_VT]: lst.fa.set_flow_style() return cast(list[_VT], lst) +def get_process_kind(process: etree.Element[str]) -> str: + """Return the Cantera collision kind corresponding to an LXCat process.""" + collision_type = process.attrib.get("collisionType", "").lower() + + if collision_type in {"elastic", "effective"}: + return collision_type + + inelastic_type = process.attrib.get("inelasticType", "").lower() + + if "ionization" in inelastic_type: + return "ionization" + + if "attachment" in inelastic_type: + return "attachment" + + # Rotational, vibrational, electronic, and dissociative energy-loss + # processes are all inelastic channels for the EEDF solver. + return "excitation" + + +def make_collision_name( + database: str, + target: str, + kind: str, + product: str | None, + threshold: float, + used_names: set[str], +) -> str: + """Create a deterministic and unique electron-collision name.""" + parts = [database, target, kind] + + if product is not None: + parts.append(product) + + if threshold: + parts.append(f"{threshold:g}-eV") + + normalized = [] + for part in parts: + value = re.sub(r"[^A-Za-z0-9]+", "-", part).strip("-") + if value: + normalized.append(value) + + base_name = "-".join(normalized) or "lxcat-collision" + name = base_name + suffix = 2 + + while name in used_names: + name = f"{base_name}-{suffix}" + suffix += 1 + + used_names.add(name) + return name + + +def reaction_equation(reaction: Process | dict) -> str | None: + """Return the equation stored in a generated or existing reaction.""" + if isinstance(reaction, Process): + return reaction.equation + + if isinstance(reaction, dict): + equation = reaction.get("equation") + if equation is not None: + return str(equation) + + return None + + +def mark_duplicate_reactions(reactions: list) -> None: + """Mark reactions having identical equation strings as duplicates.""" + counts: dict[str, int] = {} + + for reaction in reactions: + equation = reaction_equation(reaction) + if equation is not None: + counts[equation] = counts.get(equation, 0) + 1 + + for reaction in reactions: + equation = reaction_equation(reaction) + if equation is None or counts[equation] <= 1: + continue + + if isinstance(reaction, Process): + reaction.duplicate = True + else: + reaction["duplicate"] = True + class IncorrectXMLNode(LookupError): def __init__(self, message: str = "", node: etree.Element | None = None) -> None: """Error raised when a required node is incorrect in the XML tree. @@ -103,199 +240,343 @@ def __init__(self, message: str = "", node: etree.Element | None = None) -> None def convert( inpfile: str | Path | None = None, database: str | None = None, - mechfile: str | None = None, + mechfile: str | Path | None = None, phase: str | None = None, insert: bool | None = True, outfile: str | Path | None = None, ) -> None: - """Convert an LXCat XML file to a YAML file. - - :param inpfile: - The input LXCat file name. - :param database: - The name of the database. For example, "itikawa". - :param mechfile: - The reaction mechanism file. This option requires using the Cantera library. - :param phase: - The phase name of the mechanism file. This option requires a ``mechfile`` to - also be specified. - :param insert: - The flag of whether to insert the collision reactions or not. - :param outfile: - The output YAML file name. - - All files are assumed to be relative to the current working directory of the Python - process running this script. - """ - if inpfile is not None: - inpfile = Path(inpfile) - lxcat_text = inpfile.read_text().lstrip() - if outfile is None: - outfile = inpfile.with_suffix(".yaml") - else: + """Convert an LXCat XML file to Cantera YAML collision data.""" + if inpfile is None: raise ValueError("'inpfile' must be specified") + inpfile = Path(inpfile) + lxcat_text = inpfile.read_text().lstrip() + + if outfile is None: + outfile = inpfile.with_suffix(".yaml") + + outfile = Path(outfile) + if insert and mechfile is None: raise ValueError("'mech' must be specified if 'insert' is used") + if phase is not None and mechfile is None: + raise ValueError("'mech' must be specified if 'phase' is used.") + gas: OptionalSolutionType = None + source_data = None + if mechfile is not None: if Solution is None: - print("Cantera is not used, so the mechanism file cannot be used.") - sys.exit(1) - elif phase is not None: + raise RuntimeError( + "The Cantera Python module is required when 'mech' is used." + ) + + mechfile = Path(mechfile) + + if phase is not None: gas = Solution(mechfile, phase, transport_model=None) else: gas = Solution(mechfile, transport_model=None) - elif phase is not None: - raise ValueError("'mech' must be specified if 'phase' is used.") + loader = yaml.YAML(typ="rt") + with mechfile.open("r") as mechanism: + source_data = loader.load(mechanism) xml_tree = etree.fromstring(lxcat_text) - # If insert key word is used, create a process list, - # and append all processes together - process_list: list[Process] | None = None - if not insert: - process_list = [] + collision_list: list[ElectronCollision] = [] + process_list: list[Process] = [] + used_names: set[str] = set() + + # Prevent conflicts with collision names already present in a mechanism. + if source_data is not None: + for item in source_data.get("electron-collisions", []): + if isinstance(item, dict) and "name" in item: + used_names.add(str(item["name"])) for database_node in xml_tree: - if database is not None: - if database_node.attrib["id"] != database: + database_id = database_node.attrib.get("id", "lxcat") + + if database is not None and database_id != database: + continue + + groups = get_children(database_node, "groups") + if not groups: + raise IncorrectXMLNode( + "The database requires a 'groups' node.", database_node + ) + + for group in groups[0]: + processes = get_children(group, "processes") + if not processes: continue - # Get groups node - groups_node = get_children(database_node, "groups")[0] + for process in processes[0]: + registerProcess( + process, + collision_list, + process_list, + gas, + database_id, + used_names, + ) - for group in groups_node: - for process in get_children(group, "processes")[0]: - registerProcess(process, process_list, gas) + mark_duplicate_reactions(process_list) if not insert: - # Put process list in collision node - collision_node = {"collisions": process_list} - with Path(outfile).open("w") as output_file: - emitter.dump(collision_node, output_file) - else: - # Get mechanism file unit system - units = None - assert mechfile is not None - with open(mechfile, "r") as mech: - data = yaml.YAML(typ="rt").load(mech) - if "units" in data: - units = data["units"] - assert gas is not None - gas.write_yaml(outfile, units=units) - -def registerProcess(process: etree.Element, - process_list: list[Process] | None, - gas: OptionalSolutionType) -> None: + output = BlockMap([ + ("electron-collisions", collision_list), + ("reactions", process_list), + ]) + + with outfile.open("w") as output_file: + emitter.dump(output, output_file) + + return + + assert gas is not None + assert source_data is not None + assert mechfile is not None + + # Preserve the official converter's behavior of writing a self-contained + # mechanism before adding the newly converted entries. + units = source_data.get("units") + gas.write_yaml(outfile, units=units) + + loader = yaml.YAML(typ="rt") + with outfile.open("r") as mechanism: + output_data = loader.load(mechanism) + + if output_data is None: + output_data = BlockMap() + + if "reactions" not in output_data: + output_data["reactions"] = CommentedSeq() + + if not isinstance(output_data["reactions"], list): + raise ValueError("The top-level 'reactions' entry must be a sequence.") + + output_data["reactions"].extend(process_list) + mark_duplicate_reactions(output_data["reactions"]) + + if "electron-collisions" not in output_data: + output_data["electron-collisions"] = CommentedSeq() + + if not isinstance(output_data["electron-collisions"], list): + raise ValueError( + "The top-level 'electron-collisions' entry must be a sequence." + ) + + output_data["electron-collisions"].extend(collision_list) + + # Keep collision definitions after reactions, following the usual + # mechanism-file organization. + collisions = output_data.pop("electron-collisions") + output_data["electron-collisions"] = collisions + + with outfile.open("w") as output_file: + emitter.dump(output_data, output_file) + +def registerProcess( + process: etree.Element[str], + collision_list: list[ElectronCollision], + process_list: list[Process], + gas: OptionalSolutionType, + database: str, + used_names: set[str], +) -> None: """ - Add a collision process (electron collision reaction) to process_list - and gas object if it exists. - - :param process: - The collision process (electron collision reaction) - :param process_list: - The list of collision processes - :param gas: - The Cantera Solution object + Convert one LXCat process to an electron-collision definition and, + when appropriate, to a chemical reaction referencing that definition. """ - # Get electron specie name electron_name = gas.electron_species_name if gas is not None else "e" - # Parse the threshold + # Read the threshold energy. Other parameters, such as the electron-to- + # target mass ratio, are intentionally ignored here. threshold = 0.0 - parameters_node = get_children(process, "parameters")[0] - if len(get_children(parameters_node, "parameter")) == 1: - parameter = get_children(parameters_node, "parameter")[0] - if parameter.attrib["name"] == 'E': - assert parameter.text is not None - threshold = float(parameter.text) - - # Parse the equation + parameter_nodes = get_children(process, "parameters") + + if parameter_nodes: + for parameter in get_children(parameter_nodes[0], "parameter"): + if parameter.attrib.get("name") == "E": + if parameter.text is None: + raise IncorrectXMLNode( + "The threshold parameter requires a value.", parameter + ) + threshold = float(parameter.text) + break + + # Read the target species. + reactant_nodes = get_children(process, "reactants") + if not reactant_nodes: + raise IncorrectXMLNode( + "The 'process' node requires a 'reactants' node.", process + ) + + target: str | None = None + for reactant_node in reactant_nodes[0]: + if reactant_node.tag.find("molecule") != -1: + if reactant_node.text is None: + raise IncorrectXMLNode( + "A molecular reactant requires a species name.", + reactant_node, + ) + target = reactant_node.text + break + + if target is None: + raise IncorrectXMLNode( + "The electron-collision process requires a target molecule.", + process, + ) + + # A collision whose target is absent from the selected phase cannot be + # used by that phase. + if gas is not None and target not in gas.species_names: + return + + # Read the products and build the chemical equation. product_array: list[str] = [] + molecule_products: list[str] = [] + products_available = True - products: list[etree.Element[str]] = get_children(process, "products") - if products: - for product_node in products[0]: + product_nodes = get_children(process, "products") + if product_nodes: + for product_node in product_nodes[0]: if product_node.tag.find("electron") != -1: product_array.append(electron_name) + continue - if product_node.tag.find("molecule") != -1: - assert product_node.text is not None - product_name = product_node.text - if "state" in product_node.attrib: - state = product_node.attrib["state"].replace(" ","-") - # State is appended in a parenthesis - product_name += f"({state})" - if "charge" in product_node.attrib: - charge = int(product_node.attrib["charge"]) - if charge > 0: - product_name += charge*"+" - else: - product_name += -charge*"-" - - # Filter the collision based on the existed species in the mechanism file - if gas is not None and not product_name in gas.species_names: - return - - product_array.append(product_name) - - for reactant_node in get_children(process, "reactants")[0]: - if reactant_node.tag.find("molecule") != -1: - reactant = reactant_node.text - # Filter the collision based on the existed species in the mechanism file - if gas is not None and not reactant in gas.species_names: - return + if product_node.tag.find("molecule") == -1: + continue + + if product_node.text is None: + raise IncorrectXMLNode( + "A molecular product requires a species name.", + product_node, + ) + + product_name = product_node.text + + if "state" in product_node.attrib: + state = product_node.attrib["state"].replace(" ", "-") + product_name += f"({state})" + + if "charge" in product_node.attrib: + charge = int(product_node.attrib["charge"]) + if charge > 0: + product_name += charge * "+" + elif charge < 0: + product_name += -charge * "-" + + molecule_products.append(product_name) + product_array.append(product_name) + + if gas is not None and product_name not in gas.species_names: + products_available = False - if product_array: # not empty + if product_array: products_string = " + ".join(product_array) else: - # No product is identified. Use the reactant as the product. - products_string = f"{reactant} + {electron_name}" + products_string = f"{target} + {electron_name}" - equation = f"{reactant} + {electron_name} => {products_string}" + equation = f"{target} + {electron_name} => {products_string}" - # Parse the cross-section data - data_x_node = get_children(process, "data_x")[0] - if data_x_node is None: - raise IncorrectXMLNode("The 'process' node requires the 'data_x' node.", process) + # Read the tabulated cross-section data. + data_x_nodes = get_children(process, "data_x") + if not data_x_nodes: + raise IncorrectXMLNode( + "The 'process' node requires the 'data_x' node.", process + ) - data_y_node = get_children(process, "data_y")[0] - if data_y_node is None: - raise IncorrectXMLNode("The 'process' node requires the 'data_y' node.", process) + data_y_nodes = get_children(process, "data_y") + if not data_y_nodes: + raise IncorrectXMLNode( + "The 'process' node requires the 'data_y' node.", process + ) + + data_x_node = data_x_nodes[0] + data_y_node = data_y_nodes[0] + + if data_x_node.text is None or data_y_node.text is None: + raise IncorrectXMLNode( + "Cross-section data nodes cannot be empty.", process + ) - assert data_x_node.text is not None - assert data_y_node.text is not None energy_levels = Flowlist(map(float, data_x_node.text.split())) cross_sections = Flowlist(map(float, data_y_node.text.split())) - # Edit energy levels and cross section if len(energy_levels) != len(cross_sections): - raise IncorrectXMLNode("Energy levels (data_x) and cross section " - "(data_y) must have the same length.", process) - + raise IncorrectXMLNode( + "Energy levels (data_x) and cross section data (data_y) " + "must have the same length.", + process, + ) + + # Ensure that the tabulation starts at the threshold with zero cross + # section, preserving the behavior of the original converter. if energy_levels[0] > threshold: - # Use Flowlist again to ensure correct YAML format energy_levels = Flowlist([threshold, *energy_levels]) cross_sections = Flowlist([0.0, *cross_sections]) else: cross_sections[0] = 0.0 - # If insert mode is on, add the process as a reaction to the gas object. - if gas is not None: - R = ct.Reaction( - equation=equation, - rate=ct.ElectronCollisionPlasmaRate(energy_levels=energy_levels, - cross_sections=cross_sections)) - gas.add_reaction(R) + if len(energy_levels) < 2: + raise IncorrectXMLNode( + "An electron collision requires at least two energy levels.", + process, + ) + + kind = get_process_kind(process) + + # The collision format has a single optional product field. Keep it only + # when the LXCat process has one molecular product. + product = molecule_products[0] if len(molecule_products) == 1 else None + + # When writing a complete mechanism, an unavailable product must not be + # used to construct a synthetic standalone collision reaction. + collision_product = product + if gas is not None and not products_available: + collision_product = None + + collision_name = make_collision_name( + database, + target, + kind, + product, + threshold, + used_names, + ) - # If insert mode is off, process_list is used to store the data. - if process_list is not None: - process_list.append(Process(equation=equation, - energy_levels=energy_levels, - cross_sections=cross_sections)) + collision_list.append( + ElectronCollision( + name=collision_name, + target=target, + product=collision_product, + kind=kind, + threshold=threshold, + energy_levels=energy_levels, + cross_sections=cross_sections, + ) + ) + + # Elastic and effective collisions affect the EEDF but do not represent + # chemical source terms. + if kind in {"elastic", "effective"}: + return + + # Keep the collision data when a product species is unavailable, but do + # not create an invalid chemical reaction. + if not products_available: + return + + process_list.append( + Process( + equation=equation, + collision=collision_name, + ) + ) def create_argparser() -> argparse.ArgumentParser: parser = argparse.ArgumentParser( @@ -342,9 +623,9 @@ def create_argparser() -> argparse.ArgumentParser: help=("This specifies the name of the phase in the mechanism file. Optional.")) parser.add_argument( "--insert", action="store_true", default=False, - help=("Enable inserting the electron-collision reactions into the mechanism file." - "Need to use with the argument --mech to provide the mechanism file" - "Optional.")) + help=("Insert the generated electron-collision definitions and their " + "associated chemical reactions into the specified mechanism. " + "This option requires --mech.")) parser.add_argument( "--output", default=None, help=("Specifies the OUTPUT file name. By default, the output file name is the " From 8fd33d7454f2f3026023ed8c069ef87aed00f513 Mon Sep 17 00:00:00 2001 From: Gaetanosaure Date: Wed, 2 Sep 2026 15:44:55 +0200 Subject: [PATCH 10/17] update lxcat2yaml to be able to convert lxcat 1.1 xml format --- interfaces/cython/cantera/lxcat2yaml.py | 289 +++++++++++++++++------- 1 file changed, 203 insertions(+), 86 deletions(-) diff --git a/interfaces/cython/cantera/lxcat2yaml.py b/interfaces/cython/cantera/lxcat2yaml.py index b810262140a..f792c8c80b6 100644 --- a/interfaces/cython/cantera/lxcat2yaml.py +++ b/interfaces/cython/cantera/lxcat2yaml.py @@ -117,9 +117,26 @@ def to_yaml(cls, representer: SafeRepresenter, node: Process) -> MappingNode: emitter.register_class(ElectronCollision) emitter.register_class(Process) +def normalized_tag(tag: str) -> str: + """Return a normalized XML tag name without namespace or separators.""" + if "}" in tag: + tag = tag.rsplit("}", 1)[1] + + return tag.replace("_", "").lower() + # Return indices of a child name -def get_children(parent: etree.Element[str], child_name: str) -> list[etree.Element[str]]: - return [child for child in parent if child.tag.find(child_name) != -1] +def get_children( + parent: etree.Element[str], + child_name: str, +) -> list[etree.Element[str]]: + """Return direct children matching either legacy or LXCat 1.1 tags.""" + expected = normalized_tag(child_name) + + return [ + child + for child in parent + if normalized_tag(child.tag) == expected + ] _VT = TypeVar("_VT") # Value type. @@ -130,25 +147,33 @@ def Flowlist(*args: Iterable[_VT], **kwargs: _VT) -> list[_VT]: return cast(list[_VT], lst) def get_process_kind(process: etree.Element[str]) -> str: - """Return the Cantera collision kind corresponding to an LXCat process.""" - collision_type = process.attrib.get("collisionType", "").lower() + """Return the Cantera collision kind for legacy and LXCat 1.1 XML.""" + process_type = ( + process.attrib.get("type") + or process.attrib.get("collisionType") + or "" + ).strip().lower() + + inelastic_type = process.attrib.get( + "inelasticType", "" + ).strip().lower() - if collision_type in {"elastic", "effective"}: - return collision_type + if process_type == "effective": + return "effective" - inelastic_type = process.attrib.get("inelasticType", "").lower() + if process_type == "elastic": + return "elastic" - if "ionization" in inelastic_type: + if "ionization" in process_type or "ionization" in inelastic_type: return "ionization" - if "attachment" in inelastic_type: + if "attachment" in process_type or "attachment" in inelastic_type: return "attachment" - # Rotational, vibrational, electronic, and dissociative energy-loss - # processes are all inelastic channels for the EEDF solver. + # Excitation, vibrational, rotational, dissociation, and other + # inelastic energy-loss processes belong to the "excitation" EEDF group. return "excitation" - def make_collision_name( database: str, target: str, @@ -295,7 +320,15 @@ def convert( if isinstance(item, dict) and "name" in item: used_names.add(str(item["name"])) - for database_node in xml_tree: + database_nodes = get_children(xml_tree, "database") + + if not database_nodes: + raise IncorrectXMLNode( + "The LXCat file does not contain a database node.", + xml_tree, + ) + + for database_node in database_nodes: database_id = database_node.attrib.get("id", "lxcat") if database is not None and database_id != database: @@ -307,12 +340,12 @@ def convert( "The database requires a 'groups' node.", database_node ) - for group in groups[0]: + for group in get_children(groups[0], "group"): processes = get_children(group, "processes") if not processes: continue - for process in processes[0]: + for process in get_children(processes[0], "process"): registerProcess( process, collision_list, @@ -378,6 +411,122 @@ def convert( with outfile.open("w") as output_file: emitter.dump(output_data, output_file) +def normalize_lxcat_species( + name: str, + electron_name: str, +) -> tuple[str, bool]: + """Normalize an LXCat species name and identify electrons.""" + name = name.strip() + + if name.lower() in {"e", "electron"}: + return electron_name, True + + # Remove whitespace internal to species labels and convert the LXCat + # charge notation Ar^+ / O^- to the usual Cantera notation Ar+ / O-. + name = re.sub(r"\s+", "", name) + name = name.replace("^+", "+") + name = name.replace("^-", "-") + + return name, False + + +def get_process_species( + process: etree.Element[str], + role: str, + electron_name: str, +) -> list[tuple[str, bool]]: + """ + Read reactants or products from legacy and LXCat 1.1 process nodes. + + The returned boolean indicates whether the species is an electron. + """ + legacy_containers = get_children(process, f"{role}s") + + if legacy_containers: + species_nodes = list(legacy_containers[0]) + else: + species_containers = get_children(process, "species") + if not species_containers: + return [] + + species_nodes = get_children(species_containers[0], role) + + species: list[tuple[str, bool]] = [] + + for node in species_nodes: + tag = normalized_tag(node.tag) + + # Legacy format: + if tag == "electron": + species.append((electron_name, True)) + continue + + if node.text is None: + raise IncorrectXMLNode( + f"An LXCat {role} requires a species name.", + node, + ) + + name, is_electron = normalize_lxcat_species( + node.text, + electron_name, + ) + + # Legacy format represents states and charges using attributes on + # the node. + if tag == "molecule": + if "state" in node.attrib: + state = node.attrib["state"].replace(" ", "-") + name += f"({state})" + + if "charge" in node.attrib: + charge = int(node.attrib["charge"]) + + if charge > 0: + name += charge * "+" + elif charge < 0: + name += -charge * "-" + + species.append((name, is_electron)) + + return species + +def get_process_threshold(process: etree.Element[str]) -> float: + """Read the threshold energy from legacy or LXCat 1.1 XML.""" + parameter_nodes = get_children(process, "parameters") + + if not parameter_nodes: + return 0.0 + + for parameter in parameter_nodes[0]: + tag = normalized_tag(parameter.tag) + + is_legacy_threshold = ( + tag == "parameter" + and parameter.attrib.get("name", "").lower() == "e" + ) + is_lxcat_threshold = tag == "e" + + if not is_legacy_threshold and not is_lxcat_threshold: + continue + + if parameter.text is None: + raise IncorrectXMLNode( + "The threshold-energy node requires a value.", + parameter, + ) + + units = parameter.attrib.get("units", "eV") + if units.lower() != "ev": + raise IncorrectXMLNode( + "LXCat threshold energies must be expressed in eV.", + parameter, + ) + + return float(parameter.text) + + return 0.0 + def registerProcess( process: etree.Element[str], collision_list: list[ElectronCollision], @@ -394,87 +543,54 @@ def registerProcess( # Read the threshold energy. Other parameters, such as the electron-to- # target mass ratio, are intentionally ignored here. - threshold = 0.0 - parameter_nodes = get_children(process, "parameters") + electron_name = gas.electron_species_name if gas is not None else "e" + threshold = get_process_threshold(process) + kind = get_process_kind(process) - if parameter_nodes: - for parameter in get_children(parameter_nodes[0], "parameter"): - if parameter.attrib.get("name") == "E": - if parameter.text is None: - raise IncorrectXMLNode( - "The threshold parameter requires a value.", parameter - ) - threshold = float(parameter.text) - break - - # Read the target species. - reactant_nodes = get_children(process, "reactants") - if not reactant_nodes: - raise IncorrectXMLNode( - "The 'process' node requires a 'reactants' node.", process - ) + reactant_species = get_process_species( + process, + "reactant", + electron_name, + ) - target: str | None = None - for reactant_node in reactant_nodes[0]: - if reactant_node.tag.find("molecule") != -1: - if reactant_node.text is None: - raise IncorrectXMLNode( - "A molecular reactant requires a species name.", - reactant_node, - ) - target = reactant_node.text - break + targets = [ + name + for name, is_electron in reactant_species + if not is_electron + ] - if target is None: + if len(targets) != 1: raise IncorrectXMLNode( - "The electron-collision process requires a target molecule.", + "An electron-collision process requires exactly one " + "non-electron target species.", process, ) - # A collision whose target is absent from the selected phase cannot be - # used by that phase. + target = targets[0] + if gas is not None and target not in gas.species_names: return - # Read the products and build the chemical equation. + product_species = get_process_species( + process, + "product", + electron_name, + ) + product_array: list[str] = [] molecule_products: list[str] = [] products_available = True - product_nodes = get_children(process, "products") - if product_nodes: - for product_node in product_nodes[0]: - if product_node.tag.find("electron") != -1: - product_array.append(electron_name) - continue - - if product_node.tag.find("molecule") == -1: - continue - - if product_node.text is None: - raise IncorrectXMLNode( - "A molecular product requires a species name.", - product_node, - ) - - product_name = product_node.text - - if "state" in product_node.attrib: - state = product_node.attrib["state"].replace(" ", "-") - product_name += f"({state})" + for product_name, is_electron in product_species: + product_array.append(product_name) - if "charge" in product_node.attrib: - charge = int(product_node.attrib["charge"]) - if charge > 0: - product_name += charge * "+" - elif charge < 0: - product_name += -charge * "-" + if is_electron: + continue - molecule_products.append(product_name) - product_array.append(product_name) + molecule_products.append(product_name) - if gas is not None and product_name not in gas.species_names: - products_available = False + if gas is not None and product_name not in gas.species_names: + products_available = False if product_array: products_string = " + ".join(product_array) @@ -514,13 +630,14 @@ def registerProcess( process, ) - # Ensure that the tabulation starts at the threshold with zero cross - # section, preserving the behavior of the original converter. - if energy_levels[0] > threshold: - energy_levels = Flowlist([threshold, *energy_levels]) - cross_sections = Flowlist([0.0, *cross_sections]) - else: - cross_sections[0] = 0.0 + # Effective and elastic cross sections may be non-zero at zero energy. + # Their original tabulated values must therefore be preserved. + if kind not in {"effective", "elastic"}: + if energy_levels[0] > threshold: + energy_levels = Flowlist([threshold, *energy_levels]) + cross_sections = Flowlist([0.0, *cross_sections]) + elif energy_levels[0] == threshold: + cross_sections[0] = 0.0 if len(energy_levels) < 2: raise IncorrectXMLNode( From 2aa650838bb901b4bc0f7441693832d43672c7b9 Mon Sep 17 00:00:00 2001 From: Gaetanosaure Date: Wed, 2 Sep 2026 15:51:55 +0200 Subject: [PATCH 11/17] update tests --- test/data/lxcat-test-convert.yaml | 9 +- test/kinetics/kineticsFromYaml.cpp | 56 ++++---- test/python/test_convert.py | 207 ++++++++++++++++++++--------- test/python/test_reaction.py | 77 ++++++++--- test/python/test_thermo.py | 183 ++++++++++++++++++++----- test/thermo/thermoToYaml.cpp | 135 +++++++++++++++++++ 6 files changed, 527 insertions(+), 140 deletions(-) diff --git a/test/data/lxcat-test-convert.yaml b/test/data/lxcat-test-convert.yaml index 47b9815840c..bbf6f8bd194 100644 --- a/test/data/lxcat-test-convert.yaml +++ b/test/data/lxcat-test-convert.yaml @@ -69,6 +69,13 @@ reactions: - equation: O2 + E => E + O2 type: electron-collision-plasma + collision: O2-effective + +electron-collisions: +- name: O2-effective + target: O2 + kind: effective + threshold: 0.0 energy-levels: [0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0] cross-sections: [0.0, 5.97e-20, 6.45e-20, 6.74e-20, 6.93e-20, 7.2e-20, - 7.52e-20, 7.86e-20, 8.21e-20, 8.49e-20, 8.8e-20] + 7.52e-20, 7.86e-20, 8.21e-20, 8.49e-20, 8.8e-20] \ No newline at end of file diff --git a/test/kinetics/kineticsFromYaml.cpp b/test/kinetics/kineticsFromYaml.cpp index 202a63940c5..f3fe32b4ebe 100644 --- a/test/kinetics/kineticsFromYaml.cpp +++ b/test/kinetics/kineticsFromYaml.cpp @@ -643,40 +643,44 @@ TEST(Reaction, TwoTempPlasmaExtendedListFromYaml) TEST(Reaction, ElectronCollisionPlasmaFromYaml) { - auto sol = newSolution( - "oxygen-plasma.yaml", - "discretized-electron-energy-plasma", - "none" - ); - AnyMap rxn = AnyMap::fromYamlString( - "{equation: O2 + E => E + O2," - " type: electron-collision-plasma," - " collision: O2-effective}" - ); - auto R = newReaction(rxn, *sol->kinetics()); - EXPECT_EQ(R->reactants.at("O2"), 1); - EXPECT_EQ(R->reactants.at("E"), 1); - EXPECT_EQ(R->products.at("O2"), 1); - EXPECT_EQ(R->products.at("E"), 1); - const auto rate = std::dynamic_pointer_cast(R->rate()); + auto sol = newSolution("oxygen-plasma.yaml","discretized-electron-energy-plasma","none"); + + ASSERT_GT(sol->kinetics()->nReactions(), 1); + auto R = sol->kinetics()->reaction(1); + + EXPECT_DOUBLE_EQ(R->reactants.at("O2"), 1.0); + EXPECT_DOUBLE_EQ(R->reactants.at("E"), 1.0); + EXPECT_DOUBLE_EQ(R->products.at("O2"), 1.0); + EXPECT_DOUBLE_EQ(R->products.at("E"), 1.0); + + auto rate = std::dynamic_pointer_cast(R->rate()); ASSERT_NE(rate, nullptr); + EXPECT_EQ(rate->collisionName(), "O2-effective"); EXPECT_EQ(rate->kind(), "effective"); - const vector expectedLevels = { + + const vector expectedEnergyLevels = { 0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0 }; const vector expectedCrossSections = { - 0.0, 5.97e-20, 6.45e-20, 6.74e-20, 6.93e-20, 7.2e-20, - 7.52e-20, 7.86e-20, 8.21e-20, 8.49e-20, 8.8e-20 + 0.0, 5.97e-20, 6.45e-20, 6.74e-20, 6.93e-20, + 7.2e-20, 7.52e-20, 7.86e-20, 8.21e-20, + 8.49e-20, 8.8e-20 }; - const auto energyLevels = rate->energyLevels(); - const auto crossSections = rate->crossSections(); - ASSERT_EQ(energyLevels.size(), expectedLevels.size()); - ASSERT_EQ(crossSections.size(), expectedCrossSections.size()); - for (size_t k = 0; k < expectedLevels.size(); k++) { - EXPECT_DOUBLE_EQ(energyLevels[k], expectedLevels[k]); - EXPECT_DOUBLE_EQ(crossSections[k], expectedCrossSections[k]); + + ASSERT_EQ(rate->energyLevels().size(), expectedEnergyLevels.size()); + ASSERT_EQ(rate->crossSections().size(), expectedCrossSections.size()); + + for (size_t k = 0; k < expectedEnergyLevels.size(); k++) { + EXPECT_DOUBLE_EQ( + rate->energyLevels()[k], + expectedEnergyLevels[k] + ); + EXPECT_DOUBLE_EQ( + rate->crossSections()[k], + expectedCrossSections[k] + ); } } diff --git a/test/python/test_convert.py b/test/python/test_convert.py index ac9d5ceada1..51cc334d021 100644 --- a/test/python/test_convert.py +++ b/test/python/test_convert.py @@ -1686,76 +1686,153 @@ def convert(self, inputFile=None, database=None, mechFile=None, phase=None, mechFile = self.test_data_path / mechFile if output is None: output = Path(inputFile).stem # strip '.xml' - # output to work dir - output = self.test_work_path / output + output = self.test_work_path / output lxcat2yaml.convert(inputFile, database, mechFile, phase, insert, output) return output + @staticmethod + def get_collision(data, target, kind): + """Return the unique collision matching a target and collision kind.""" + matches = [ + collision + for collision in data["electron-collisions"] + if collision["target"] == target + and collision["kind"] == kind + ] + assert len(matches) == 1 + return matches[0] + def test_mechanism_with_lxcat(self): - # get Solution from the mechanism file phase = "isotropic-electron-energy-plasma" - mechFile = "lxcat-test-convert.yaml" - gas1 = ct.Solution(self.test_data_path / mechFile, - phase=phase, transport_model=None) - - # get a stand-alone collisions - standAloneFile = "stand-alone-lxcat.yaml" - self.convert(inputFile='lxcat-test-convert.xml', database="test", - mechFile=mechFile, insert=False, - output=standAloneFile) - - # add collisions to the reaction list - rxn_list = ct.Reaction.list_from_file(self.test_work_path / standAloneFile, - gas1, section="collisions") - for R in rxn_list: - gas1.add_reaction(R) - - # get Solution from the output file - output = "output-lxcat.yaml" - self.convert(inputFile="lxcat-test-convert.xml", database="test", - mechFile=mechFile, insert=True, - output=output) - gas2 = ct.Solution(self.test_work_path / output, - phase=phase, transport_model=None) - - # check number of reactions - assert gas1.n_reactions == gas2.n_reactions == 4 - for i in range(1, gas1.n_reactions): - assert (gas1.reaction(i).rate.energy_levels - == approx(gas2.reaction(i).rate.energy_levels)) - assert (gas1.reaction(i).rate.cross_sections - == approx(gas2.reaction(i).rate.cross_sections)) + mech_file = "lxcat-test-convert.yaml" + + # Verify that the initial mechanism already uses the named-collision + # format and can be loaded by Cantera. + gas1 = ct.Solution(self.test_data_path / mech_file, phase=phase, transport_model=None) + assert gas1.n_reactions == 2 + + # Convert the LXCat data while using the mechanism only as a species + # filter. CO2 collisions must be excluded because CO2 is not present + # in the selected phase. + stand_alone_file = self.convert( + inputFile="lxcat-test-convert.xml", + database="test", + mechFile=mech_file, + phase=phase, + insert=False, + output="stand-alone-lxcat.yaml", + ) + stand_alone = load_yaml(stand_alone_file) + assert len(stand_alone["electron-collisions"]) == 2 + assert len(stand_alone["reactions"]) == 2 + assert {collision["target"] for collision in stand_alone["electron-collisions"]} == {"O2"} + assert {collision["kind"] for collision in stand_alone["electron-collisions"]} == {"ionization", "attachment"} + collision_names = {collision["name"] for collision in stand_alone["electron-collisions"]} + + for reaction in stand_alone["reactions"]: + assert reaction["type"] == "electron-collision-plasma" + assert reaction["collision"] in collision_names + + # Insert the same converted data into the complete mechanism. + output = self.convert( + inputFile="lxcat-test-convert.xml", + database="test", + mechFile=mech_file, + phase=phase, + insert=True, + output="output-lxcat.yaml", + ) + output_data = load_yaml(output) + # One effective collision was already present, and the converter adds + # ionization and attachment. + assert len(output_data["electron-collisions"]) == 3 + names = [collision["name"] for collision in output_data["electron-collisions"]] + assert len(names) == len(set(names)) + + electron_reactions = [reaction for reaction in output_data["reactions"] + if reaction.get("type") == "electron-collision-plasma"] + assert len(electron_reactions) == 3 + for reaction in electron_reactions: + assert "collision" in reaction + assert reaction["collision"] in names + + # Loading the generated mechanism verifies that all references can be + # resolved and that the collision data reach the reaction rates. + gas2 = ct.Solution(output, phase=phase, transport_model=None) + assert gas2.n_reactions == gas1.n_reactions + 2 + assert gas2.n_reactions == 4 + plasma_reactions = [gas2.reaction(i) for i in range(gas2.n_reactions) + if gas2.reaction(i).reaction_type == "electron-collision-plasma"] + assert len(plasma_reactions) == 3 + ionization = next(reaction for reaction in plasma_reactions + if "O2(Total-Ionization)+" in reaction.products) + attachment = next(reaction for reaction in plasma_reactions + if "O2-" in reaction.products) + assert ionization.rate.energy_levels == approx([15.0, 20.0]) + assert ionization.rate.cross_sections == approx([0.0, 5.5e-22]) + assert attachment.rate.energy_levels == approx([0.0, 1.0]) + assert attachment.rate.cross_sections == approx([0.0, 1.0e-22]) def test_stand_alone_lxcat(self): - outfile = "stand-alone-lxcat-without-mech.yaml" - self.convert(inputFile='lxcat-test-convert.xml', - database="test", insert=False, - output=outfile) - - # get Solution from the input file - phase = "isotropic-electron-energy-plasma" - inputFile = "lxcat-test-convert-species.yaml" - gas = ct.Solution(self.test_data_path / inputFile, - phase=phase, transport_model=None) - - # add collisions to the reaction list - rxn_list = ct.Reaction.list_from_file(self.test_work_path / outfile, - gas, section="collisions") - - # verify the data - assert len(rxn_list) == 3 - assert rxn_list[0].equation == "CO2 + e => CO2 + e" - assert rxn_list[0].reaction_type == "electron-collision-plasma" - assert rxn_list[0].rate.energy_levels == approx([0.0, 1.0]) - assert rxn_list[0].rate.cross_sections == approx([0.0, 1.0e-22]) - - assert rxn_list[1].equation == "O2 + e => O2(Total-Ionization)+ + 2 e" - assert rxn_list[1].reaction_type == "electron-collision-plasma" - assert rxn_list[1].rate.energy_levels == approx([15., 20.]) - assert rxn_list[1].rate.cross_sections == approx([0.0, 5.5e-22]) - - assert rxn_list[2].equation == "O2 + e => O2-" - assert rxn_list[2].reaction_type == "electron-collision-plasma" - assert rxn_list[2].rate.energy_levels == approx([0.0, 1.0]) - assert rxn_list[2].rate.cross_sections == approx([0.0, 1.0e-22]) + outfile = self.convert( + inputFile="lxcat-test-convert.xml", + database="test", + insert=False, + output="stand-alone-lxcat-without-mech.yaml", + ) + data = load_yaml(outfile) + assert "electron-collisions" in data + assert "reactions" in data + collisions = data["electron-collisions"] + reactions = data["reactions"] + + # The XML contains four collision data sets: + # - CO2 elastic + # - CO2 ionization + # - O2 ionization + # - O2 attachment + assert len(collisions) == 4 + + # Only inelastic chemical processes generate reactions. + assert len(reactions) == 3 + + names = [collision["name"] for collision in collisions] + assert len(names) == len(set(names)) + + referenced_names = {reaction["collision"] for reaction in reactions} + inelastic_names = {collision["name"] for collision in collisions + if collision["kind"] not in {"elastic", "effective"}} + + # Every generated reaction references exactly one named inelastic + # collision. Elastic/effective collisions remain EEDF-only data. + assert referenced_names == inelastic_names + + for reaction in reactions: + assert reaction["type"] == "electron-collision-plasma" + assert "energy-levels" not in reaction + assert "cross-sections" not in reaction + + co2_elastic = self.get_collision(data, "CO2", "elastic") + assert co2_elastic["threshold"] == approx(0.0) + assert co2_elastic["energy-levels"] == approx([0.0, 1.0]) + assert co2_elastic["cross-sections"] == approx([0.0, 1.0e-22]) + assert co2_elastic["name"] not in referenced_names + + co2_ionization = self.get_collision(data, "CO2", "ionization") + assert co2_ionization["product"] == "CO2(Total-Ionization)+" + assert co2_ionization["threshold"] == approx(15.0) + assert co2_ionization["energy-levels"] == approx([15.0, 20.0]) + assert co2_ionization["cross-sections"] == approx([0.0, 5.5e-22]) + + o2_ionization = self.get_collision(data, "O2", "ionization") + assert o2_ionization["product"] == "O2(Total-Ionization)+" + assert o2_ionization["threshold"] == approx(15.0) + assert o2_ionization["energy-levels"] == approx([15.0, 20.0]) + assert o2_ionization["cross-sections"] == approx([0.0, 5.5e-22]) + + o2_attachment = self.get_collision(data, "O2", "attachment") + assert o2_attachment["product"] == "O2-" + assert o2_attachment["threshold"] == approx(0.0) + assert o2_attachment["energy-levels"] == approx([0.0, 1.0]) + assert o2_attachment["cross-sections"] == approx([0.0, 1.0e-22]) \ No newline at end of file diff --git a/test/python/test_reaction.py b/test/python/test_reaction.py index 9ae3c4c9cdc..f32c2cbc51c 100644 --- a/test/python/test_reaction.py +++ b/test/python/test_reaction.py @@ -2500,27 +2500,22 @@ def electron_reaction_data(request, setup_electron_reaction_tests): @pytest.mark.usefixtures("electron_reaction_data") class TestElectronCollisionPlasmaReaction(ReactionTests): - # This test only test the data input and output but not evaluating the reaction - # rate. The rate evaluation is tested in kineticsFromYaml.cpp because plasma - # reaction rate is much complicated and depends on electron energy distribution - # function. + # Electron-collision rates reference cross-section data stored in the + # root-level electron-collisions section. The tests below verify rate + # construction, serialization, and evaluation in a plasma phase that + # provides the referenced collision definition. _rate_cls = ct.ElectronCollisionPlasmaRate _equation = "O2 + E <=> E + O2" _rate = { - "equation": "O2 + E <=> E + O2", "type": "electron-collision-plasma", - "energy-levels": [0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0], - "cross-sections": [0.0, 5.97e-20, 6.45e-20, 6.74e-20, 6.93e-20, 7.2e-20, - 7.52e-20, 7.86e-20, 8.21e-20, 8.49e-20, 8.8e-20] - } + "collision": "O2-effective", + } _index = 1 _rate_type = "electron-collision-plasma" _yaml = """ equation: O2 + E <=> E + O2 type: electron-collision-plasma - energy-levels: [0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0] - cross-sections: [0.0, 5.97e-20, 6.45e-20, 6.74e-20, 6.93e-20, 7.2e-20, - 7.52e-20, 7.86e-20, 8.21e-20, 8.49e-20, 8.8e-20] + collision: O2-effective """ _phase_def = """ phases: @@ -2535,23 +2530,73 @@ class TestElectronCollisionPlasmaReaction(ReactionTests): energy-levels: [0.0, 0.1, 1.0, 10.0] distribution: [0.0, 0.2, 0.7, 0.01] normalize: False + + electron-collisions: + - name: O2-effective + target: O2 + kind: effective + energy-levels: [0.0, 1.0, 2.0, 3.0, 4.0, 5.0, + 6.0, 7.0, 8.0, 9.0, 10.0] + cross-sections: [0.0, 5.97e-20, 6.45e-20, 6.74e-20, + 6.93e-20, 7.2e-20, 7.52e-20, 7.86e-20, + 8.21e-20, 8.49e-20, 8.8e-20] """ _rc_units = ct.Units("m^3 / kmol / s") - def eval_rate(self, rate): + def _solution_with_reaction(self, reaction): + # A collision reference can only be resolved by a phase containing + # the corresponding root-level electron-collision definition. gas = ct.Solution(yaml=self._phase_def) gas.TDY = self.soln.TDY - gas.add_reaction(ct.Reaction(equation=self._equation, rate=rate)) + gas.add_reaction(reaction) + return gas + + def eval_rate(self, rate): + reaction = ct.Reaction(equation=self._equation, rate=rate) + gas = self._solution_with_reaction(reaction) return gas.forward_rate_constants[0] + def check_rxn(self, rxn): + # Override the generic implementation, which constructs a phase from + # species and reactions only and therefore cannot provide root-level + # electron-collision definitions. + original = self.soln.reaction(self._index) + assert rxn.reactants == original.reactants + assert rxn.products == original.products + + gas = self._solution_with_reaction(rxn) + self.check_solution(gas) + + def test_add_rxn(self): + # The generic test creates a phase without an electron-collisions + # section. Use the complete plasma definition required to resolve + # the collision reference. + rxn = self.from_yaml() + gas = self._solution_with_reaction(rxn) + self.check_solution(gas) + @pytest.mark.skip("No rate is not supported") def test_no_rate(self): pass + def test_replace_rate(self): + # An electron-collision rate cannot be constructed without a named + # collision reference. Test replacement using two valid rate objects. + rxn = self.from_yaml() + rxn.rate = ct.ReactionRate.from_dict(self._rate) + self.check_rxn(rxn) + def test_roundtrip(self): - # check round-trip instantiation via input_data + # Check that serialization retains the reference and does not restore + # the deprecated inline cross-section representation. rxn = self.from_rate(self._rate_obj) rate_input_data = dict(rxn.rate.input_data) + + assert rate_input_data["type"] == "electron-collision-plasma" + assert rate_input_data["collision"] == "O2-effective" + assert "energy-levels" not in rate_input_data + assert "cross-sections" not in rate_input_data + rate_obj = rxn.rate.__class__(input_data=rate_input_data) rxn2 = self.from_rate(rate_obj) - self.check_rxn(rxn2) + self.check_rxn(rxn2) \ No newline at end of file diff --git a/test/python/test_thermo.py b/test/python/test_thermo.py index 1ae6493ec96..39c00960421 100644 --- a/test/python/test_thermo.py +++ b/test/python/test_thermo.py @@ -1378,17 +1378,6 @@ def phase(self): phase.isotropic_shape_factor = 1.0 return phase - @property - def collision_data(self): - return { - "equation": "O2 + E => E + O2", - "type": "electron-collision-plasma", - "energy-levels": [0.0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0], - "cross-sections": [0.0, 3.83e-20, 4.47e-20, 4.79e-20, 5.07e-20, 5.31e-20, - 5.49e-20, 5.64e-20, 5.77e-20, 5.87e-20, 5.97e-20], - "duplicate": True, - } - def test_converting_electron_energy_to_temperature(self, phase): phase.mean_electron_energy = 1.0 Te = 2.0 / 3.0 * ct.electron_charge / ct.boltzmann @@ -1456,28 +1445,35 @@ def test_elastic_power_loss_high_T(self, phase): corr_factor = (phase.mean_temperature / phase.T )**2 assert phase.elastic_power_loss * corr_factor == approx(2865540) - def test_elastic_power_loss_replace_rate(self, phase): - phase.TPX = 1000, ct.one_atm, "O2:1, E:1e-5" - rate = ct.ReactionRate.from_dict(self.collision_data) - phase.reaction(1).rate = rate - # Since `elasticPowerLoss` involves two multiplications with the concentration - # of species, and since now it is computed as C=P/(R * ), we need to apply a - # correction factor to get back the value that was obtained with C=P/(R * T), - # i.e. C_old = C_new * (/T), and q_ela_old = q_ela_new * (C_old/C_new)^2. - corr_factor = (phase.mean_temperature / phase.T )**2 - assert phase.elastic_power_loss * corr_factor == approx(11765800095) + @staticmethod + def electron_collision_reaction(phase): + for index in range(phase.n_reactions): + reaction = phase.reaction(index) + if reaction.rate.type == "electron-collision-plasma": + return index, reaction + raise AssertionError("No electron-collision-plasma reaction found") - def test_elastic_power_loss_add_reaction(self, phase): - phase2 = ct.Solution(thermo="plasma", kinetics="bulk", - species=phase.species(), reactions=[]) - phase.TPX = 1000, ct.one_atm, "O2:1, E:1e-5" - phase.add_reaction(ct.Reaction.from_dict(self.collision_data, phase)) - # Since `elasticPowerLoss` involves two multiplications with the concentration - # of species, and since now it is computed as C=P/(R * ), we need to apply a - # correction factor to get back the value that was obtained with C=P/(R * T), - # i.e. C_old = C_new * (/T), and q_ela_old = q_ela_new * (C_old/C_new)^2. - corr_factor = (phase.mean_temperature / phase.T )**2 - assert phase.elastic_power_loss * corr_factor == approx(18612132428) + def test_electron_collision_reference(self, phase): + _, reaction = self.electron_collision_reaction(phase) + rate_data = reaction.rate.input_data + + assert rate_data["type"] == "electron-collision-plasma" + assert isinstance(rate_data["collision"], str) + assert rate_data["collision"] + assert "energy-levels" not in rate_data + assert "cross-sections" not in rate_data + + def test_adding_collision_reference_does_not_duplicate_collision(self, phase): + phase.TPX = 1000.0, ct.one_atm, "O2:1, E:1e-5" + expected = phase.elastic_power_loss + + _, reaction = self.electron_collision_reaction(phase) + reaction_data = reaction.input_data + reaction_data["duplicate"] = True + + phase.add_reaction(ct.Reaction.from_dict(reaction_data, phase)) + + assert phase.elastic_power_loss == approx(expected) def test_elastic_power_loss_change_levels(self, phase): phase.electron_energy_levels = np.linspace(0,10,101) @@ -1699,6 +1695,129 @@ def test_vibrational_reservoir_invalid_mapping(self, phase_name, message): with pytest.raises(ct.CanteraError, match=message): ct.Solution("vibrational-relaxation.yaml", phase_name) + @staticmethod + def collision_test_yaml(momentum_kind): + levels = ", ".join(str(value) for value in range(21)) + elastic = ", ".join(["2.0e-20"] * 21) + inelastic = ", ".join(["0.0", "0.0"] + ["1.0e-20"] * 19) + + if momentum_kind == "effective": + momentum = ", ".join(["2.0e-20", "2.0e-20"] + ["3.0e-20"] * 19) + elif momentum_kind == "elastic": + momentum = elastic + else: + raise ValueError(f"Unexpected collision kind: {momentum_kind}") + + return f""" + phases: + - name: plasma + thermo: plasma + kinetics: bulk + elements: [O, E] + species: [E, O2, O2_excited] + reactions: all + electron-energy-distribution: + type: Boltzmann-two-term + energy-levels: [{levels}] + + species: + - name: E + composition: {{E: 1}} + thermo: + model: constant-cp + h0: 0.0 J/kmol + s0: 0.0 J/kmol/K + cp0: 20786.0 J/kmol/K + + - name: O2 + composition: {{O: 2}} + thermo: + model: constant-cp + h0: 0.0 J/kmol + s0: 0.0 J/kmol/K + cp0: 30000.0 J/kmol/K + + - name: O2_excited + composition: {{O: 2}} + thermo: + model: constant-cp + h0: 0.0 J/kmol + s0: 0.0 J/kmol/K + cp0: 30000.0 J/kmol/K + + reactions: + - equation: O2 + E => O2_excited + E + type: electron-collision-plasma + collision: O2-excitation + + electron-collisions: + - name: O2-momentum-transfer + target: O2 + kind: {momentum_kind} + energy-levels: [{levels}] + cross-sections: [{momentum}] + + - name: O2-excitation + target: O2 + product: O2_excited + kind: excitation + threshold: 2.0 + energy-levels: [{levels}] + cross-sections: [{inelastic}] + """ + + @staticmethod + def solve_collision_test(yaml): + phase = ct.Solution(yaml=yaml, transport_model=None) + phase.TPX = ( + 300.0, + ct.one_atm, + "O2:1, O2_excited:0, E:1e-10", + ) + phase.reduced_electric_field = 200.0e-21 + phase.update_electron_energy_distribution() + + return (np.array(phase.electron_energy_levels), + np.array(phase.electron_energy_distribution), + phase.elastic_power_loss) + + def test_effective_cross_section_equivalent_to_decomposition(self): + effective = self.solve_collision_test(self.collision_test_yaml("effective")) + decomposed = self.solve_collision_test(self.collision_test_yaml("elastic")) + + effective_grid, effective_eedf, effective_loss = effective + elastic_grid, elastic_eedf, elastic_loss = decomposed + assert effective_grid == approx(elastic_grid) + assert effective_eedf == approx(elastic_eedf, rel=1e-8, abs=1e-12) + assert effective_loss == approx(elastic_loss, rel=1e-8) + + def test_unknown_electron_collision_reference(self): + yaml = self.collision_test_yaml("effective").replace( + "collision: O2-excitation", + "collision: undefined-collision") + with pytest.raises(ct.CanteraError, match="undefined-collision"): + ct.Solution(yaml=yaml, transport_model=None) + + def test_duplicate_electron_collision_name(self): + yaml = self.collision_test_yaml("effective") + """ + - name: O2-momentum-transfer + target: O2 + kind: elastic + energy-levels: [0.0, 1.0] + cross-sections: [1.0e-20, 1.0e-20] + """ + with pytest.raises(ct.CanteraError, match="Duplicate electron-collision definition"): + ct.Solution(yaml=yaml, transport_model=None) + + def test_missing_electron_collision_name(self): + yaml = self.collision_test_yaml("effective").replace( + "- name: O2-momentum-transfer\n" + " target: O2", + "- target: O2", + ) + with pytest.raises(ct.CanteraError, match="name"): + ct.Solution(yaml=yaml, transport_model=None) + class TestImport: """ Tests the various ways of creating a Solution object diff --git a/test/thermo/thermoToYaml.cpp b/test/thermo/thermoToYaml.cpp index 087fb1feb6f..9e8c2660e33 100644 --- a/test/thermo/thermoToYaml.cpp +++ b/test/thermo/thermoToYaml.cpp @@ -5,6 +5,9 @@ #include "cantera/base/YamlWriter.h" #include "cantera/thermo/Species.h" #include "cantera/thermo/PlasmaPhase.h" +#include "cantera/kinetics/Kinetics.h" +#include "cantera/kinetics/Reaction.h" +#include "cantera/kinetics/ElectronCollisionPlasmaRate.h" using namespace Cantera; typedef vector strvec; @@ -325,6 +328,138 @@ TEST_F(ThermoToYaml, DiscretizedElectronEnergyPlasma) EXPECT_DOUBLE_EQ(dist[3], 0.01); } +// TEST_F(ThermoToYaml, ElectronCollisionReferences) +// { +// const string phaseName = "discretized-electron-energy-plasma"; + +// auto original = newSolution("oxygen-plasma.yaml", phaseName, "none"); +// auto second = newSolution( +// "oxygen-plasma.yaml", +// "isotropic-electron-energy-plasma", +// "none" +// ); + +// YamlWriter writer; +// writer.addPhase(original); +// writer.addPhase(second); +// writer.skipUserDefined(); + +// AnyMap root = AnyMap::fromYamlString(writer.toYamlString()); + +// ASSERT_TRUE(root.hasKey("electron-collisions")); +// const auto& collisionDefs = +// root["electron-collisions"].asVector(); + +// // Both plasma phases use the same definition, which should only be +// // written once. +// ASSERT_EQ(collisionDefs.size(), 1u); + +// const AnyMap& collision = root["electron-collisions"].getMapWhere( +// "name", "test_O2_effective_O2_0" +// ); + +// EXPECT_EQ(collision["target"], "O2"); +// EXPECT_EQ(collision["product"], "O2"); +// EXPECT_EQ(collision["kind"], "effective"); + +// const auto& levels = +// collision["energy-levels"].asVector(); +// const auto& crossSections = +// collision["cross-sections"].asVector(); + +// ASSERT_EQ(levels.size(), 11u); +// ASSERT_EQ(crossSections.size(), levels.size()); +// EXPECT_DOUBLE_EQ(levels[1], 1.0); +// EXPECT_DOUBLE_EQ(crossSections[1], 5.97e-20); + +// ASSERT_TRUE(root.hasKey("reactions")); +// const AnyMap& reaction = root["reactions"].getMapWhere( +// "type", "electron-collision-plasma" +// ); + +// EXPECT_EQ(reaction["collision"], "test_O2_effective_O2_0"); +// EXPECT_FALSE(reaction.hasKey("energy-levels")); +// EXPECT_FALSE(reaction.hasKey("cross-sections")); + +// // Verify that the serialized document can be used to reconstruct a +// // complete Solution and that the collision reference is resolved. +// auto duplicate = newSolution( +// root["phases"].getMapWhere("name", phaseName), +// root, +// "none" +// ); + +// ASSERT_EQ(duplicate->kinetics()->nReactions(), 2u); + +// auto rate = std::dynamic_pointer_cast( +// duplicate->kinetics()->reaction(1)->rate() +// ); +// ASSERT_TRUE(rate); + +// EXPECT_EQ(rate->collisionName(), "test_O2_effective_O2_0"); +// EXPECT_EQ(rate->kind(), "effective"); +// ASSERT_EQ(rate->energyLevels().size(), 11u); +// ASSERT_EQ(rate->crossSections().size(), 11u); +// EXPECT_DOUBLE_EQ(rate->energyLevels()[1], 1.0); +// EXPECT_DOUBLE_EQ(rate->crossSections()[1], 5.97e-20); +// } + +TEST_F(ThermoToYaml, ElectronCollisionReferences) +{ + const string phaseName = "isotropic-electron-energy-plasma"; + auto original = newSolution("oxygen-plasma.yaml", phaseName, "none"); + + YamlWriter writer; + writer.addPhase(original); + writer.skipUserDefined(); + + AnyMap root = AnyMap::fromYamlString(writer.toYamlString()); + + ASSERT_TRUE(root.hasKey("electron-collisions")); + const auto collisions = + root["electron-collisions"].asVector(); + + ASSERT_FALSE(collisions.empty()); + + for (const auto& collision : collisions) { + EXPECT_TRUE(collision.hasKey("name")); + EXPECT_FALSE(collision["name"].asString().empty()); + EXPECT_TRUE(collision.hasKey("target")); + EXPECT_TRUE(collision.hasKey("kind")); + EXPECT_TRUE(collision.hasKey("energy-levels")); + EXPECT_TRUE(collision.hasKey("cross-sections")); + } + + ASSERT_TRUE(root.hasKey("reactions")); + const auto reactions = root["reactions"].asVector(); + + bool foundElectronCollision = false; + for (const auto& reaction : reactions) { + if (reaction.getString("type", "") != + "electron-collision-plasma") { + continue; + } + + foundElectronCollision = true; + EXPECT_TRUE(reaction.hasKey("collision")); + EXPECT_FALSE(reaction.hasKey("energy-levels")); + EXPECT_FALSE(reaction.hasKey("cross-sections")); + } + + EXPECT_TRUE(foundElectronCollision); + + auto duplicate = newSolution( + root["phases"].getMapWhere("name", phaseName), + root, + "none" + ); + + ASSERT_NE(duplicate, nullptr); + ASSERT_NE(duplicate->kinetics(), nullptr); + EXPECT_EQ(original->kinetics()->nReactions(), + duplicate->kinetics()->nReactions()); +} + class ThermoYamlRoundTrip : public testing::Test { From 92b45ca521a4c9337c2690dc4d6aa8fc6a47a01e Mon Sep 17 00:00:00 2001 From: Gaetanosaure Date: Wed, 2 Sep 2026 15:52:50 +0200 Subject: [PATCH 12/17] remove unnecessary lxcat-test-convert-species.yaml --- test/data/lxcat-test-convert-species.yaml | 23 ----------------------- 1 file changed, 23 deletions(-) delete mode 100644 test/data/lxcat-test-convert-species.yaml diff --git a/test/data/lxcat-test-convert-species.yaml b/test/data/lxcat-test-convert-species.yaml deleted file mode 100644 index 48e129f8bde..00000000000 --- a/test/data/lxcat-test-convert-species.yaml +++ /dev/null @@ -1,23 +0,0 @@ -description: |- - This file is for testing purpose. The species data is needed to construct the - reactions in the generated list of the plasma collision reactions from the lxcat file. - Reactions which contains species not in the species list will be ignored. - -units: {length: cm, quantity: molec, activation-energy: K} - -phases: -- name: isotropic-electron-energy-plasma - thermo: plasma - elements: [O, C, E] - species: - - nasa_gas.yaml/species: [O2, O2-, CO2] - - lxcat-test-convert.yaml/species: [E, O2(Total-Ionization)+] - - kinetics: gas - reactions: none - transport: Ion - electron-energy-distribution: - type: isotropic - shape-factor: 2.0 - mean-electron-energy: 1.0 eV - energy-levels: [0.0, 0.1, 1.0, 10.0] From 746802efe776481b9c1b04a877ea02b1dad0545a Mon Sep 17 00:00:00 2001 From: Gaetanosaure Date: Thu, 3 Sep 2026 11:35:40 +0200 Subject: [PATCH 13/17] Make the warning readable when inconsistent cross-section dataset are loaded --- .../cantera/thermo/EEDFTwoTermApproximation.h | 8 ++ src/thermo/EEDFTwoTermApproximation.cpp | 115 ++++++++++++++++-- test/python/test_thermo.py | 59 ++++++++- 3 files changed, 171 insertions(+), 11 deletions(-) diff --git a/include/cantera/thermo/EEDFTwoTermApproximation.h b/include/cantera/thermo/EEDFTwoTermApproximation.h index ad89b4098b6..a8f7d465c9e 100644 --- a/include/cantera/thermo/EEDFTwoTermApproximation.h +++ b/include/cantera/thermo/EEDFTwoTermApproximation.h @@ -369,6 +369,14 @@ class EEDFTwoTermApproximation //! First call to calculateDistributionFunction bool m_first_call; + //! Flags indicating whether a negative reconstructed elastic cross-section warning + //! has been emitted for each species during the first EEDF calculation. + vector m_negativeElasticCrossSectionWarningsIssued; + + //! Enable negative reconstructed elastic cross-section warnings during the first + //! EEDF calculation only. + bool m_negativeElasticCrossSectionWarningsEnabled = true; + //! Energy grid spacing type. Can be linear, quadratic or geometric. string m_gridType = "linear"; diff --git a/src/thermo/EEDFTwoTermApproximation.cpp b/src/thermo/EEDFTwoTermApproximation.cpp index 72c367c6ee9..99a6f1b6c3a 100644 --- a/src/thermo/EEDFTwoTermApproximation.cpp +++ b/src/thermo/EEDFTwoTermApproximation.cpp @@ -13,7 +13,9 @@ #include "cantera/numerics/funcs.h" #include "cantera/thermo/PlasmaPhase.h" #include "cantera/kinetics/ElectronCollisionPlasmaRate.h" +#include #include +#include namespace Cantera { @@ -177,6 +179,7 @@ int EEDFTwoTermApproximation::calculateDistributionFunction() { if (m_first_call) { initSpeciesIndexCrossSections(); + m_negativeElasticCrossSectionWarningsIssued.assign(m_phase->nSpecies(), false); m_first_call = false; } @@ -212,6 +215,11 @@ int EEDFTwoTermApproximation::calculateDistributionFunction() // update electron mobility m_electronMobility = electronMobility(m_f0); + + // Warnings are diagnostic information for the initial calculation only. Keep them + // disabled for every subsequent user-requested EEDF update. + m_negativeElasticCrossSectionWarningsEnabled = false; + return 0; } @@ -779,6 +787,16 @@ void EEDFTwoTermApproximation::calculateTotalElasticCrossSection() { m_sigmaElastic.assign(m_points, 0.0); + struct NegativeCrossSectionSample { + size_t gridIndex; + double value; + }; + + // Negative reconstructed values are collected first so that diagnostics can be + // emitted once per target instead to avoid the spam once per energy-grid point. + vector> negativeSamples( + m_phase->nSpecies()); + for (size_t k : m_phase->kElastic()) { auto rate = m_phase->collisionRate(k); auto levels = rate->energyLevels(); @@ -821,14 +839,7 @@ void EEDFTwoTermApproximation::calculateTotalElasticCrossSection() } if (elasticCrossSection < 0.0) { - throw CanteraError( - "EEDFTwoTermApproximation::" - "calculateTotalElasticCrossSection", - "Reconstructed elastic cross section for target '{}' " - "is negative at {} eV: {} m^2. The effective and " - "inelastic cross-section data are inconsistent.", - m_phase->speciesName(target), - energy, elasticCrossSection); + negativeSamples[target].push_back({i, elasticCrossSection}); } } @@ -837,6 +848,94 @@ void EEDFTwoTermApproximation::calculateTotalElasticCrossSection() elasticCrossSection; } } + // the function already performed its main job, the following is to emit + // readable diagnostics to the user instead of spamming the log + // with one warning per negqtive grid point. + for (size_t target = 0; target < negativeSamples.size(); target++) { + auto& samples = negativeSamples[target]; + if (samples.empty() || !m_negativeElasticCrossSectionWarningsEnabled + || m_negativeElasticCrossSectionWarningsIssued[target]) { + continue; + } + + // More than one effective data set for the same target is normally rejected + // by PlasmaPhase. Sort and merge duplicate grid indices defensively, retaining + // the most negative reconstructed value at each point. + std::sort(samples.begin(), samples.end(), + [](const auto& left, const auto& right) { + return left.gridIndex < right.gridIndex; + }); + + vector uniqueSamples; + uniqueSamples.reserve(samples.size()); + + for (const auto& sample : samples) { + if (!uniqueSamples.empty() + && uniqueSamples.back().gridIndex == sample.gridIndex) { + uniqueSamples.back().value = std::min( + uniqueSamples.back().value, sample.value); + } else { + uniqueSamples.push_back(sample); + } + } + + std::ostringstream details; + details << std::setprecision(8); + + size_t first = 0; + size_t regionCount = 0; + + while (first < uniqueSamples.size()) { + size_t last = first; + size_t minimum = first; + + // Consecutive negative grid points form one negative interval. + while (last + 1 < uniqueSamples.size() + && uniqueSamples[last + 1].gridIndex + == uniqueSamples[last].gridIndex + 1) { + last++; + + if (uniqueSamples[last].value < uniqueSamples[minimum].value) { + minimum = last; + } + } + + if (regionCount > 0) { + details << "; "; + } + + const size_t firstIndex = uniqueSamples[first].gridIndex; + const size_t lastIndex = uniqueSamples[last].gridIndex; + const size_t minimumIndex = uniqueSamples[minimum].gridIndex; + + if (first == last) { + // A single isolated negative grid point. + details << "point " << m_gridEdge[firstIndex] << " eV" + << " (sigma_el = " << uniqueSamples[first].value + << " m^2)"; + } else { + // Two or more consecutive negative points. + details << "interval [" << m_gridEdge[firstIndex] << ", " + << m_gridEdge[lastIndex] << "] eV" + << " (" << last - first + 1 + << " grid points; minimum sigma_el = " + << uniqueSamples[minimum].value + << " m^2 at " << m_gridEdge[minimumIndex] << " eV)"; + } + + regionCount++; + first = last + 1; + } + + warn_user("EEDFTwoTermApproximation::" + "calculateTotalElasticCrossSection", + "Reconstructed elastic cross section for target '{}' is negative at {}. " + "The effective cross section is smaller than the sum of the inelastic " + "cross sections over these energies; the input data are inconsistent " + "with the LXCat definition of an effective cross section.", + m_phase->speciesName(target), details.str()); + m_negativeElasticCrossSectionWarningsIssued[target] = true; + } } void EEDFTwoTermApproximation::setGridCache() diff --git a/test/python/test_thermo.py b/test/python/test_thermo.py index 39c00960421..1fde3a2f479 100644 --- a/test/python/test_thermo.py +++ b/test/python/test_thermo.py @@ -1696,18 +1696,35 @@ def test_vibrational_reservoir_invalid_mapping(self, phase_name, message): ct.Solution("vibrational-relaxation.yaml", phase_name) @staticmethod - def collision_test_yaml(momentum_kind): + def collision_test_yaml(momentum_kind, *, negative_elastic=False, adapt_grid=False): levels = ", ".join(str(value) for value in range(21)) elastic = ", ".join(["2.0e-20"] * 21) inelastic = ", ".join(["0.0", "0.0"] + ["1.0e-20"] * 19) if momentum_kind == "effective": - momentum = ", ".join(["2.0e-20", "2.0e-20"] + ["3.0e-20"] * 19) + high_energy_value = "9.9e-21" if negative_elastic else "3.0e-20" + momentum = ", ".join( + ["2.0e-20", "2.0e-20"] + [high_energy_value] * 19) elif momentum_kind == "elastic": momentum = elastic else: raise ValueError(f"Unexpected collision kind: {momentum_kind}") + if adapt_grid: + grid_definition = ( + " initial-max-energy-level: 20.0\n" + " grid-cell-count: 20\n" + " energy-levels-spacing: linear\n" + " energy-grid-adaptation:\n" + " enabled: true\n" + " min-decay-decades: 1000.0\n" + " max-decay-decades: 1001.0\n" + " update-factor: 0.25\n" + " max-iterations: 3" + ) + else: + grid_definition = f" energy-levels: [{levels}]" + return f""" phases: - name: plasma @@ -1718,7 +1735,7 @@ def collision_test_yaml(momentum_kind): reactions: all electron-energy-distribution: type: Boltzmann-two-term - energy-levels: [{levels}] +{grid_definition} species: - name: E @@ -1791,6 +1808,40 @@ def test_effective_cross_section_equivalent_to_decomposition(self): assert effective_eedf == approx(elastic_eedf, rel=1e-8, abs=1e-12) assert effective_loss == approx(elastic_loss, rel=1e-8) + def test_negative_elastic_warning_only_on_first_solve(self, recwarn): + yaml = self.collision_test_yaml( + "effective", negative_elastic=True, adapt_grid=True) + phase = ct.Solution(yaml=yaml, transport_model=None) + phase.TPX = 300.0, ct.one_atm, "O2:1, O2_excited:0, E:1e-10" + phase.reduced_electric_field = 200.0e-21 + + warning_text = "Reconstructed elastic cross section for target 'O2'" + with pytest.warns(UserWarning, match=warning_text) as caught: + phase.update_electron_energy_distribution() + + matching = [ + item for item in caught + if warning_text in str(item.message) + ] + assert len(matching) == 1 + assert "interval [" in str(matching[0].message) + + # The deliberately unreachable decay target forces all three internal + # grid-adaptation iterations during the first user-requested solve. + assert phase.electron_energy_levels[-1] == approx(39.0625) + first_max_energy = phase.electron_energy_levels[-1] + + # Change the field to ensure that this is a real second EEDF solve. The + # negative reconstructed cross section must no longer issue a warning. + recwarn.clear() + phase.reduced_electric_field = 300.0e-21 + phase.update_electron_energy_distribution() + + assert phase.electron_energy_levels[-1] > first_max_energy + assert not any( + warning_text in str(item.message) for item in recwarn + ) + def test_unknown_electron_collision_reference(self): yaml = self.collision_test_yaml("effective").replace( "collision: O2-excitation", @@ -1818,6 +1869,8 @@ def test_missing_electron_collision_name(self): with pytest.raises(ct.CanteraError, match="name"): ct.Solution(yaml=yaml, transport_model=None) + + class TestImport: """ Tests the various ways of creating a Solution object From 0fd57c88b94126bf4ea45411c88b3b1ed69a819a Mon Sep 17 00:00:00 2001 From: Gaetanosaure Date: Thu, 3 Sep 2026 12:55:18 +0200 Subject: [PATCH 14/17] force cross section extrapolation to zero outside of their definition --- .../cantera/thermo/EEDFTwoTermApproximation.h | 6 ++++++ src/thermo/EEDFTwoTermApproximation.cpp | 19 +++++++++++++------ test/python/test_thermo.py | 10 +++++----- 3 files changed, 24 insertions(+), 11 deletions(-) diff --git a/include/cantera/thermo/EEDFTwoTermApproximation.h b/include/cantera/thermo/EEDFTwoTermApproximation.h index a8f7d465c9e..6edee711ee3 100644 --- a/include/cantera/thermo/EEDFTwoTermApproximation.h +++ b/include/cantera/thermo/EEDFTwoTermApproximation.h @@ -426,6 +426,12 @@ class EEDFTwoTermApproximation span fpts, double below_value, double above_value); + //! Linearly interpolate a cross section, returning zero outside its + //! tabulated energy range. + double interpolateCrossSection(double energy, + span energyLevels, + span crossSections); + //! The threshold in reduced electric field [townsend, Td] below which no EEDF will //! be computed, but a Maxwellian at the gas temperature will be imposed instead. double m_thresholdToMaxwellian = 1; diff --git a/src/thermo/EEDFTwoTermApproximation.cpp b/src/thermo/EEDFTwoTermApproximation.cpp index 99a6f1b6c3a..bccbff6ff21 100644 --- a/src/thermo/EEDFTwoTermApproximation.cpp +++ b/src/thermo/EEDFTwoTermApproximation.cpp @@ -243,6 +243,14 @@ double EEDFTwoTermApproximation::linearInterpBounded( return linearInterp(x, xpts, fpts); } +double EEDFTwoTermApproximation::interpolateCrossSection( + double energy, span energyLevels, + span crossSections) +{ + return linearInterpBounded( + energy, energyLevels, crossSections, 0.0, 0.0); +} + void EEDFTwoTermApproximation::projectPreviousEEDFOnCurrentGrid( const Eigen::VectorXd& oldGridCenter, const Eigen::VectorXd& oldF0) { @@ -773,12 +781,12 @@ void EEDFTwoTermApproximation::calculateTotalCrossSection() for (size_t i = 0; i < m_points; i++) { m_totalCrossSectionCenter[i] += moleFraction * - linearInterp(m_gridCenter[i], levels, sections); + interpolateCrossSection(m_gridCenter[i], levels, sections); } for (size_t i = 0; i < m_points + 1; i++) { m_totalCrossSectionEdge[i] += moleFraction * - linearInterp(m_gridEdge[i], levels, sections); + interpolateCrossSection(m_gridEdge[i], levels, sections); } } } @@ -813,8 +821,7 @@ void EEDFTwoTermApproximation::calculateTotalElasticCrossSection() for (size_t i = 0; i < m_points; i++) { const double energy = m_gridEdge[i]; - double elasticCrossSection = - linearInterp(energy, levels, sections); + double elasticCrossSection = interpolateCrossSection(energy, levels, sections); if (rate->kind() == "effective") { // By definition: @@ -832,7 +839,7 @@ void EEDFTwoTermApproximation::calculateTotalElasticCrossSection() auto inelasticRate = m_phase->collisionRate(kInelastic); - elasticCrossSection -= linearInterp( + elasticCrossSection -= interpolateCrossSection( energy, inelasticRate->energyLevels(), inelasticRate->crossSections()); @@ -976,7 +983,7 @@ void EEDFTwoTermApproximation::setGridCache() nodes.resize(std::distance(nodes.begin(), last)); vector sigma0(nodes.size()); for (size_t i = 0; i < nodes.size(); i++) { - sigma0[i] = linearInterp(nodes[i], x, y); + sigma0[i] = interpolateCrossSection(nodes[i], x, y); } // search position of cell j diff --git a/test/python/test_thermo.py b/test/python/test_thermo.py index 1fde3a2f479..0b7a1b05e22 100644 --- a/test/python/test_thermo.py +++ b/test/python/test_thermo.py @@ -1717,8 +1717,8 @@ def collision_test_yaml(momentum_kind, *, negative_elastic=False, adapt_grid=Fal " energy-levels-spacing: linear\n" " energy-grid-adaptation:\n" " enabled: true\n" - " min-decay-decades: 1000.0\n" - " max-decay-decades: 1001.0\n" + " min-decay-decades: 1e-12\n" + " max-decay-decades: 2e-12\n" " update-factor: 0.25\n" " max-iterations: 3" ) @@ -1826,9 +1826,9 @@ def test_negative_elastic_warning_only_on_first_solve(self, recwarn): assert len(matching) == 1 assert "interval [" in str(matching[0].message) - # The deliberately unreachable decay target forces all three internal + # The deliberately tiny decay window forces all three internal # grid-adaptation iterations during the first user-requested solve. - assert phase.electron_energy_levels[-1] == approx(39.0625) + assert phase.electron_energy_levels[-1] == approx(10.24) first_max_energy = phase.electron_energy_levels[-1] # Change the field to ensure that this is a real second EEDF solve. The @@ -1837,7 +1837,7 @@ def test_negative_elastic_warning_only_on_first_solve(self, recwarn): phase.reduced_electric_field = 300.0e-21 phase.update_electron_energy_distribution() - assert phase.electron_energy_levels[-1] > first_max_energy + assert phase.electron_energy_levels[-1] < first_max_energy assert not any( warning_text in str(item.message) for item in recwarn ) From dd891a7c7a34d3c05638b9c724f6b4b0da74e4e9 Mon Sep 17 00:00:00 2001 From: Gaetanosaure Date: Thu, 3 Sep 2026 16:43:06 +0200 Subject: [PATCH 15/17] correct elecron temperature synchronisation with EEDF updates --- src/thermo/PlasmaPhase.cpp | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/src/thermo/PlasmaPhase.cpp b/src/thermo/PlasmaPhase.cpp index d2a9f91deda..5743c04362c 100644 --- a/src/thermo/PlasmaPhase.cpp +++ b/src/thermo/PlasmaPhase.cpp @@ -414,6 +414,10 @@ void PlasmaPhase::updateElectronEnergyDistribution() m_electronEnergyDist = asVectorXd(m_eedfSolver->getEEDFEdge()); m_nPoints = m_electronEnergyLevels.size(); electronEnergyLevelChanged(); + + // Keep the electron temperature consistent with the EEDF returned by + // the Boltzmann solver. + updateElectronTemperatureFromEnergyDist(); } else { throw CanteraError("PlasmaPhase::updateElectronEnergyDistribution", "Call to calculateDistributionFunction failed."); From 7647a48ae03adca84c4f4150fa8a26311349026f Mon Sep 17 00:00:00 2001 From: Gaetanosaure Date: Fri, 4 Sep 2026 09:28:17 +0200 Subject: [PATCH 16/17] update documentation --- doc/doxygen/cantera.bib | 12 ++ .../reference/kinetics/rate-constants.md | 56 +++++++- doc/sphinx/yaml/lxcat2yaml.md | 106 ++++++++++++++- doc/sphinx/yaml/reactions.md | 125 +++++++++++++----- 4 files changed, 261 insertions(+), 38 deletions(-) diff --git a/doc/doxygen/cantera.bib b/doc/doxygen/cantera.bib index 8f536d151d8..99950856065 100644 --- a/doc/doxygen/cantera.bib +++ b/doc/doxygen/cantera.bib @@ -381,6 +381,18 @@ @article{nishioka1996 url = {https://doi.org/10.1016/0010-2180(95)00132-8}, doi = {10.1016/0010-2180(95)00132-8}, year = {1996}} +@article{pancheshnyi2012, + author = {S.~Pancheshnyi and S.~Biagi and M.~C.~Bordage and + G.~J.~M.~Hagelaar and W.~L.~Morgan and A.~V.~Phelps and + L.~C.~Pitchford}, + journal = {Chemical Physics}, + pages = {148--153}, + title = {The {LXCat} project: Electron scattering cross sections and swarm + parameters for low temperature plasma modeling}, + doi = {10.1016/j.chemphys.2011.04.020}, + url = {https://doi.org/10.1016/j.chemphys.2011.04.020}, + volume = {398}, + year = {2012}} @article{pedersen1993, author = {T.~Pedersen and R.~.C.~Brown}, title = {Simulation of electric field effects in premixed methane flames}, diff --git a/doc/sphinx/reference/kinetics/rate-constants.md b/doc/sphinx/reference/kinetics/rate-constants.md index e4ac45a5023..e49fa870468 100644 --- a/doc/sphinx/reference/kinetics/rate-constants.md +++ b/doc/sphinx/reference/kinetics/rate-constants.md @@ -687,12 +687,12 @@ different representations are the responsibility of the mechanism author. ## Electron Collision Plasma Reactions The electron collision plasma reaction rate uses the electron collision data and the -electron energy distribution to calculate the reaction rate. Hagelaar and Pitchford +electron energy distribution to calculate the reaction rate. {cite:t}`hagelaar2005` define the reaction rate coefficient (Eqn. 63) as, $$ k = \gamma \int_0^{\infty} \epsilon \sigma F_0 d\epsilon $$ -where $\gamma = \sqrt{2/m_e}$ (Eqn.4 in Hagelaar {cite:t}`hagelaar2015`), $m_e$ [kg] is +where $\gamma = \sqrt{2/m_e}$ (Eqn.4 in {cite:t}`hagelaar2015`), $m_e$ [kg] is the electron mass, $\epsilon$ [J] is the electron energy, $\sigma(\epsilon)$ [m²] is the reaction collision cross section, $F_0(\epsilon)$ [$\t{J^{-3/2}}$] is the normalized electron energy distribution function, and $k$ has units of [m³/s]. @@ -707,13 +707,59 @@ $$ where $e$ is the elementary charge [C] and $N_A$ is the Avogadro constant [$\t{kmol^{-1}}$]. +### Effective and elastic cross sections + +For a given target species, the momentum-transfer cross section may be provided +either as an `elastic` cross section or as an `effective` cross section. As detailed in +{cite:t}`pancheshnyi2012`, an effective cross section already +includes the contribution of the inelastic collision processes for the same target: + +$$ +\sigma_\mathrm{eff}(\epsilon) += +\sigma_\mathrm{el}(\epsilon) ++ +\sum_j \sigma_{\mathrm{inel},j}(\epsilon). +$$ + +When an effective cross section is provided, Cantera uses it directly in the total +momentum-transfer cross section and reconstructs the elastic contribution as + +$$ +\sigma_\mathrm{el}(\epsilon) += +\sigma_\mathrm{eff}(\epsilon) +- +\sum_j \sigma_{\mathrm{inel},j}(\epsilon). +$$ + +When an effective cross section is available, the individual inelastic cross sections +are therefore not added again to the total momentum-transfer cross section. This avoids +counting their contribution twice in the electron energy distribution calculation. + +Cross sections are linearly interpolated over their tabulated energy range and are +taken to be zero outside this range. If the reconstructed elastic cross section is +negative, its value is retained and a warning identifies the affected target species +and energy intervals during the first calculation of the electron energy distribution. +Negative reconstructed values may indicate that the effective and inelastic +cross-section data are not mutually consistent. + ```{versionadded} 3.1 ``` +:::{versionchanged} 4.0 +Cross-section data are now stored in the top-level `electron-collisions` section. +An `electron-collision-plasma` reaction references one of these definitions using +the `collision` field. Cross-section data can no longer be specified directly in +the reaction entry. +::: + :::{admonition} YAML Usage :class: tip -Electron collision reactions can be defined in the YAML format by specifying -[`electron-collision-plasma`](sec-yaml-electron-collision-plasma) as the reaction `type` -and providing lists with the `cross-sections` and corresponding `energy-levels`. +Electron collision reactions are defined by specifying +[`electron-collision-plasma`](sec-yaml-electron-collision-plasma) as the reaction +`type`. The reaction references a named cross-section dataset from the top-level +[`electron-collisions`](sec-yaml-electron-collisions) section using the `collision` +field. ::: diff --git a/doc/sphinx/yaml/lxcat2yaml.md b/doc/sphinx/yaml/lxcat2yaml.md index 3f941948e2a..6df275af9c2 100644 --- a/doc/sphinx/yaml/lxcat2yaml.md +++ b/doc/sphinx/yaml/lxcat2yaml.md @@ -1,9 +1,113 @@ (sec-lxcat2yaml)= + # LXCat to YAML conversion +The `lxcat2yaml` utility converts electron-collision cross-section data +from an LXCat XML file to Cantera's YAML format. Both the legacy LXCat XML +format and LXCat XML version 1.1 are supported. + +The generated YAML separates collision data from chemical reactions: + +- Cross-section data and collision metadata are stored in the root-level + `electron-collisions` section. +- Electron-collision reactions refer to these data using the `collision` + field. +- Elastic and effective collisions are included in `electron-collisions` + for calculation of the electron energy distribution function (EEDF), but + do not generate chemical reactions. + +## Basic conversion + +To convert an LXCat XML file without inserting the results into an existing +mechanism, run: + +```bash +lxcat2yaml --input=my-cross-sections.xml --output=my-cross-sections.yaml +``` + +If the `lxcat2yaml` executable is not available on the system path, the +converter can instead be invoked as a Python module: + +```bash +python -m cantera.lxcat2yaml \ + --input=my-cross-sections.xml \ + --output=my-cross-sections.yaml +``` + +If `--output` is omitted, the output filename is generated by replacing the +input filename extension with `.yaml`. + +The resulting file contains collision definitions and any associated +chemical reactions. For example: + +```yaml +electron-collisions: +- name: Phelps-Ar-ionization-Ar-15-8-eV + target: Ar + product: Ar+ + kind: ionization + threshold: 15.8 + energy-levels: [15.8, 16.0, 17.0] + cross-sections: [0.0, 2.02e-22, 1.34e-21] + +reactions: +- equation: Ar + e => e + e + Ar+ + type: electron-collision-plasma + collision: Phelps-Ar-ionization-Ar-15-8-eV +``` + +Collision names are generated deterministically from the database identifier, +target species, collision kind, product, and threshold. A numerical suffix is +added when needed to keep names unique. + +## Inserting data into a mechanism + +The converted data can be inserted into an existing Cantera mechanism using +the `--insert` option: + +```bash +lxcat2yaml \ + --input=my-cross-sections.xml \ + --database=Phelps \ + --mech=plasma-mechanism.yaml \ + --phase=plasma \ + --insert \ + --output=plasma-mechanism-with-cross-sections.yaml +``` + +The `--mech` option is required when `--insert` is used. If the mechanism +contains multiple phases, `--phase` selects the phase whose species are used +when filtering the LXCat processes. + +The value supplied to `--database` is matched against the `id` attribute of +the LXCat `Database` element. If this option is omitted, processes from all +databases in the XML file are considered. + +When a mechanism is supplied: + +- Processes whose target species is absent from the selected phase are + omitted. +- A collision definition is retained when its target exists but one or more + product species are unavailable. +- A chemical reaction is generated only when all its product species exist + in the selected phase. +- Reactions with identical equations are marked as duplicates. +- The generated mechanism is loaded by Cantera after conversion to verify + that it is valid. + +:::{note} +Electron energies and cross sections in the LXCat input are expected to be +expressed in eV and m{sup}`2`, respectively. Species names in the LXCat file +must correspond to species in the selected Cantera phase. Common electron +names and LXCat charge notation, such as `e`, `Electron`, `Ar^+`, and `O^-`, +are normalized automatically. +::: + +## Command-line options + ```{eval-rst} .. argparse:: :module: cantera.lxcat2yaml :func: create_argparser :prog: lxcat2yaml -``` +``` \ No newline at end of file diff --git a/doc/sphinx/yaml/reactions.md b/doc/sphinx/yaml/reactions.md index 669df73db2d..fe7f5bc20e9 100644 --- a/doc/sphinx/yaml/reactions.md +++ b/doc/sphinx/yaml/reactions.md @@ -454,74 +454,135 @@ Example: ``` (sec-yaml-electron-collision-plasma)= + ### `electron-collision-plasma` -Electron collision plasma reactions involve an electron as one of the reactants, and are -parameterized by the collision cross section as a function of the electron energy. The -rate calculation is [described here](sec-electron-collision-plasma-rate). The rate -parameters are specified using the following additional fields in the reaction entry: +Electron collision plasma reactions involve an electron and a target species as +reactants. Their rate coefficient is calculated from the electron energy distribution +and a named collision cross-section dataset, as +[described here](sec-electron-collision-plasma-rate). -`energy-levels` -: A list of electron energy levels [eV] +The reaction entry uses the following additional field: -`cross-sections` -: A list of collision cross sections [m²] for the reaction at the specified energy - levels. +`collision` +: The name of an entry in the top-level + [`electron-collisions`](sec-yaml-electron-collisions) section. The referenced entry + provides the collision kind, target species, energy levels, and cross sections. + +The target specified by the collision definition must match the non-electron reactant +of the reaction. The collision `kind` must also be consistent with the reaction +stoichiometry. Example: ```yaml -- equation: O2 + e => e + e + O2+ +- equation: O2 + E => E + E + O2+ type: electron-collision-plasma - energy-levels: [13.0, 15.5, 18, 23] - cross-sections: [1.17e-22, 7.3e-22, 1.64e-21, 3.66e-21] + collision: O2-ionization ``` :::{versionadded} 3.1 ::: +:::{versionchanged} 4.0 +Cross-section data can no longer be specified directly in an +`electron-collision-plasma` reaction entry. The reaction must instead reference a +named entry in the top-level `electron-collisions` section using the `collision` +field. +::: + (sec-yaml-electron-collisions)= + ### `electron-collisions` -The `electron-collisions` field defines a list of cross-section datasets for -electron-impact processes that are used in plasma-phase simulations. These entries -are not formal reactions (they are not added to `Kinetics` objects), but serve -as data inputs for computing the electron energy distribution function. +The top-level `electron-collisions` section contains named electron collision +cross-section datasets used by plasma phases. Every entry is included when calculating +the electron energy distribution, whether or not it is referenced by a chemical +reaction. + +Entries in this section are not added to the phase's `Kinetics` object. To make a +collision process contribute to chemical source terms, define an +[`electron-collision-plasma`](sec-yaml-electron-collision-plasma) reaction whose +`collision` field references the corresponding dataset. + +Cross-section data in this format can be generated from XML files downloaded from the +[LXCat website](https://nl.lxcat.net/home/news.php) using the +[`lxcat2yaml`](sec-lxcat2yaml) conversion tool. + +Each entry uses the following fields: -Each entry includes: +`name` +: A unique, non-empty name identifying the collision dataset. `target` -: The name of the species that is the target of the collision +: The name of the species targeted by the electron collision. The target species must + be present in the plasma phase. + +`kind` +: The type of electron collision process. Supported values are: + + - `effective`: An effective momentum-transfer cross section containing the elastic + and inelastic contributions for the target. + - `elastic`: An elastic momentum-transfer cross section. + - `excitation`: An electronic, vibrational, or rotational excitation cross section. + - `ionization`: An electron-impact ionization cross section. + - `attachment`: An electron attachment cross section. + + At most one `effective` or `elastic` cross-section dataset may be defined for each + target species. `energy-levels` -: A list of electron energy values [eV] at which the cross-section is provided +: A list of at least two electron energy values [eV]. Values must be finite, + non-negative, and strictly increasing. `cross-sections` -: Corresponding cross-section values [m²] for each energy level +: A list of collision cross sections [m²] corresponding to `energy-levels`. The two + lists must have the same length. Cross-section values must be finite and + non-negative. -`kind` -: A string indicating the process type. Options include: - - `"effective"` – lumped or total effect of several channels - - `"excitation"` – electronic excitation - - `"ionization"` – electron-impact ionization - - `"attachment"` – electron attachment processes +`product` +: An optional description of the product or excited state produced by the collision. + This value does not need to correspond to a distinct species in the phase and does + not define the products of a chemical reaction. + +`threshold` +: An optional non-negative collision threshold [eV]. If omitted or set to zero for an + `excitation`, `ionization`, or `attachment` process, the threshold is inferred from + the first energy level whose cross section is greater than zero. Example: ```yaml electron-collisions: -- target: N2 - energy-levels: [0.0, 0.015, 0.03, 0.05, 0.1, 0.15, 0.2, 0.3, 0.4, 0.7, 1.2, 1.5, 1.9, - 2.2, 2.8, 3.3, 4.0, 5.0, 7.0, 10.0, 15.0, 20.0, 30.0, 75.0, 150.0] - cross-sections: [1.1e-20, 2.55e-20, 3.4e-20, 4.33e-20, 5.95e-20, 7.1e-20, 7.9e-20, - 9e-20, 9.7e-20, 1e-19, 1.04e-19, 1.2e-19, 1.96e-19, 2.85e-19, 2.8e-19, 1.72e-19, - 1.26e-19, 1.09e-19, 1.01e-19, 1.04e-19, 1.1e-19, 1.02e-19, 9e-20, 6.6e-20, 4.9e-20] +- name: O2-effective + target: O2 kind: effective + energy-levels: [0.0, 1.0, 2.0, 3.0] + cross-sections: [3.5e-21, 7.9e-20, 6.5e-20, 5.5e-20] + +- name: O2-ionization + target: O2 + product: O2+ + kind: ionization + threshold: 12.06 + energy-levels: [12.06, 13.0, 18.0, 28.0] + cross-sections: [0.0, 2.3e-22, 2.0e-21, 7.4e-21] ``` +The `O2-effective` entry is used by the electron energy distribution solver without +requiring a corresponding chemical reaction. The `O2-ionization` entry is additionally +connected to the chemical mechanism by the `electron-collision-plasma` reaction shown +above. + :::{versionadded} 3.2 ::: +:::{versionchanged} 4.0 +Each electron collision definition now requires a unique `name`. Electron collision +reactions reference these definitions using the `collision` field, and tabulated +cross-section data are stored exclusively in this section. +::: + (sec-yaml-falloff)= ### `falloff` From 7eea1eb3652ae320e59e4b1aa46ddf5e6891a800 Mon Sep 17 00:00:00 2001 From: Gaetanosaure Date: Fri, 4 Sep 2026 10:30:21 +0200 Subject: [PATCH 17/17] solve formatting issues --- include/cantera/thermo/PlasmaPhase.h | 2 +- src/kinetics/ElectronCollisionPlasmaRate.cpp | 2 +- test/data/soot-therm.dat | 2 +- test/python/test_thermo.py | 1 - 4 files changed, 3 insertions(+), 4 deletions(-) diff --git a/include/cantera/thermo/PlasmaPhase.h b/include/cantera/thermo/PlasmaPhase.h index 5f670bcc0f8..c7be030e1b0 100644 --- a/include/cantera/thermo/PlasmaPhase.h +++ b/include/cantera/thermo/PlasmaPhase.h @@ -155,7 +155,7 @@ class PlasmaPhase: public IdealGasPhase void getParameters(AnyMap& phaseNode) const override; void setParameters(const AnyMap& phaseNode, const AnyMap& rootNode=AnyMap()) override; - + //! Return the named electron-collision definitions used by this phase. /*! * The returned map associates each collision name with its root-level diff --git a/src/kinetics/ElectronCollisionPlasmaRate.cpp b/src/kinetics/ElectronCollisionPlasmaRate.cpp index 35f8215f9cf..fa75a61a977 100644 --- a/src/kinetics/ElectronCollisionPlasmaRate.cpp +++ b/src/kinetics/ElectronCollisionPlasmaRate.cpp @@ -48,7 +48,7 @@ bool ElectronCollisionPlasmaData::update(const ThermoPhase& phase, const Kinetic void ElectronCollisionPlasmaRate::setParameters(const AnyMap& node, const UnitStack& rate_units) { ReactionRate::setParameters(node, rate_units); - + if (!node.hasKey("collision")) { throw InputFileError("ElectronCollisionPlasmaRate::setParameters", node, "Electron-collision reactions require a named 'collision' reference. " diff --git a/test/data/soot-therm.dat b/test/data/soot-therm.dat index e68ec73f66d..7b5b0dcd3b4 100644 --- a/test/data/soot-therm.dat +++ b/test/data/soot-therm.dat @@ -47,5 +47,5 @@ H 264 3.65839677E+01 3.36764102E-02-1.16783938E-05 1.83077466E-09-1.06963777E-13 2 9.29809483E+03-1.81272070E+02-1.29758980E+01 1.63790064E-01-1.43851166E-04 3 6.31057915E-08-1.09568047E-11 2.48866399E+04 7.94950474E+01 4 - + END diff --git a/test/python/test_thermo.py b/test/python/test_thermo.py index 0b7a1b05e22..5ab7cfdff48 100644 --- a/test/python/test_thermo.py +++ b/test/python/test_thermo.py @@ -1869,7 +1869,6 @@ def test_missing_electron_collision_name(self): with pytest.raises(ct.CanteraError, match="name"): ct.Solution(yaml=yaml, transport_model=None) - class TestImport: """