From c00109b6b429068cb2d799499c26c9d2a0e3da9e Mon Sep 17 00:00:00 2001 From: Gaetanosaure Date: Wed, 6 May 2026 11:02:58 +0200 Subject: [PATCH 01/29] check_that_commits_work --- src/thermo/PlasmaPhase.cpp | 1 + 1 file changed, 1 insertion(+) diff --git a/src/thermo/PlasmaPhase.cpp b/src/thermo/PlasmaPhase.cpp index 26f83aeb469..992dbbdf915 100644 --- a/src/thermo/PlasmaPhase.cpp +++ b/src/thermo/PlasmaPhase.cpp @@ -777,6 +777,7 @@ double PlasmaPhase::jouleHeatingPower() const return sigma * E * E; // W/m^3 } +// This is a cool function double PlasmaPhase::intrinsicHeating() { // Joule heating: sigma * E^2 [W/m^3] From a1f940162957c1e4c4f3c6476e7649166f69ccc2 Mon Sep 17 00:00:00 2001 From: Gaetanosaure Date: Wed, 6 May 2026 12:30:45 +0200 Subject: [PATCH 02/29] [YAML] Grid specification: Linear, Quadratic and Geometric now available, in addition to the already present custom grid --- .../cantera/thermo/EEDFTwoTermApproximation.h | 6 ++ src/thermo/EEDFTwoTermApproximation.cpp | 93 +++++++++++++++++++ src/thermo/PlasmaPhase.cpp | 76 +++++++++++++++ 3 files changed, 175 insertions(+) diff --git a/include/cantera/thermo/EEDFTwoTermApproximation.h b/include/cantera/thermo/EEDFTwoTermApproximation.h index fddc47d5448..7acdb9d74f7 100644 --- a/include/cantera/thermo/EEDFTwoTermApproximation.h +++ b/include/cantera/thermo/EEDFTwoTermApproximation.h @@ -69,6 +69,12 @@ class EEDFTwoTermApproximation int calculateDistributionFunction(); void setLinearGrid(double& kTe_max, size_t& ncell); + + void setQuadraticGrid(double& kTe_max, size_t& ncell); + + void setGeometricGrid(double& kTe_max, size_t& ncell); + + void setCustomGrid(span levels); void setGridCache(); diff --git a/src/thermo/EEDFTwoTermApproximation.cpp b/src/thermo/EEDFTwoTermApproximation.cpp index 7fa38c916e5..a5a654f996e 100644 --- a/src/thermo/EEDFTwoTermApproximation.cpp +++ b/src/thermo/EEDFTwoTermApproximation.cpp @@ -43,6 +43,99 @@ void EEDFTwoTermApproximation::setLinearGrid(double& kTe_max, size_t& ncell) setGridCache(); } + +void EEDFTwoTermApproximation::setQuadraticGrid(double& kTe_max, size_t& ncell) +{ + m_points = ncell; + + m_gridCenter.resize(m_points); + m_gridEdge.resize(m_points + 1); + m_f0.resize(m_points); + m_f0_edge.resize(m_points + 1); + + double n = static_cast(m_points); + + for (size_t j = 0; j <= m_points; j++) { + double x = static_cast(j); + m_gridEdge[j] = kTe_max * x * (x + 1.0) / (n * (n + 1.0)); + } + + for (size_t j = 0; j < m_points; j++) { + m_gridCenter[j] = 0.5 * (m_gridEdge[j] + m_gridEdge[j + 1]); + } + + setGridCache(); +} + +void EEDFTwoTermApproximation::setGeometricGrid(double& kTe_max, size_t& ncell) +{ + m_points = ncell; + + m_gridCenter.resize(m_points); + m_gridEdge.resize(m_points + 1); + m_f0.resize(m_points); + m_f0_edge.resize(m_points + 1); + + double ratio = 1.05; + + double denominator = std::pow(ratio, m_points) - 1.0; + + if (std::abs(denominator) < 1e-14) { + throw CanteraError("EEDFTwoTermApproximation::setGeometricGrid", + "Invalid geometric-grid ratio."); + } + + for (size_t j = 0; j < m_points; j++) { + m_gridEdge[j] = kTe_max * (std::pow(ratio, j) - 1.0) / denominator; + + m_gridCenter[j] = kTe_max + * (std::pow(ratio, j + 0.5) - 1.0) + / denominator; + } + + m_gridEdge[m_points] = kTe_max; + + setGridCache(); +} + +void EEDFTwoTermApproximation::setCustomGrid(span levels) +{ + if (levels.size() < 2) { + throw CanteraError("EEDFTwoTermApproximation::setCustomGrid", + "Energy grid must contain at least two edge points."); + } + + m_points = levels.size() - 1; + + m_gridCenter.resize(m_points); + m_gridEdge.resize(m_points + 1); + m_f0.resize(m_points); + m_f0_edge.resize(m_points + 1); + + for (size_t j = 0; j < m_points + 1; j++) { + if (!std::isfinite(levels[j])) { + throw CanteraError("EEDFTwoTermApproximation::setCustomGrid", + "Energy grid contains a non-finite value."); + } + if (levels[j] < 0.0) { + throw CanteraError("EEDFTwoTermApproximation::setCustomGrid", + "Energy grid values must be non-negative."); + } + if (j > 0 && levels[j] <= levels[j - 1]) { + throw CanteraError("EEDFTwoTermApproximation::setCustomGrid", + "Energy grid values must be strictly increasing."); + } + + m_gridEdge[j] = levels[j]; + } + + for (size_t j = 0; j < m_points; j++) { + m_gridCenter[j] = 0.5 * (m_gridEdge[j] + m_gridEdge[j + 1]); + } + + setGridCache(); +} + int EEDFTwoTermApproximation::calculateDistributionFunction() { if (m_first_call) { diff --git a/src/thermo/PlasmaPhase.cpp b/src/thermo/PlasmaPhase.cpp index 992dbbdf915..193266cd9bc 100644 --- a/src/thermo/PlasmaPhase.cpp +++ b/src/thermo/PlasmaPhase.cpp @@ -318,7 +318,83 @@ void PlasmaPhase::setParameters(const AnyMap& phaseNode, const AnyMap& rootNode) auto levels = eedf["energy-levels"].asVector(); auto distribution = eedf["distribution"].asVector(levels.size()); setDiscretizedElectronEnergyDist(levels, distribution); + } else if (m_distributionType == "Boltzmann-two-term") { + m_eedfSolver = make_unique(this); + + if (eedf.hasKey("energy-levels")) { + // Mode A: user-provided energy grid edges + auto levels = eedf["energy-levels"].asVector(); + + m_eedfSolver->setCustomGrid(levels); + + m_nPoints = levels.size(); + + } else { + // Mode B: generated grid from readable YAML parameters + if (!eedf.hasKey("initial_max_energy_level")) { + throw CanteraError("PlasmaPhase::setParameters", + "Boltzmann-two-term requires either 'energy-levels' or " + "'initial_max_energy_level'."); + } + + if (!eedf.hasKey("initial_number_of_energy_grid_cells")) { + throw CanteraError("PlasmaPhase::setParameters", + "Boltzmann-two-term requires either 'energy-levels' or " + "'initial_number_of_energy_grid_cells'."); + } + + double kTe_max = eedf["initial_max_energy_level"].asDouble(); + size_t nGridCells = static_cast( + eedf["initial_number_of_energy_grid_cells"].asInt()); + + if (kTe_max <= 0.0) { + throw CanteraError("PlasmaPhase::setParameters", + "initial_max_energy_level must be greater than zero."); + } + + if (nGridCells == 0) { + throw CanteraError("PlasmaPhase::setParameters", + "initial_number_of_energy_grid_cells must be greater than zero."); + } + + string energyLevelsDistribution = "Linear"; + if (eedf.hasKey("energy_levels_distribution")) { + energyLevelsDistribution = + eedf["energy_levels_distribution"].asString(); + } else { + writelog("No energy_levels_distribution key found in the input file. " + "Defaulting to linear grid.\n"); + } + + if (energyLevelsDistribution == "Linear") { + m_eedfSolver->setLinearGrid(kTe_max, nGridCells); + } else if (energyLevelsDistribution == "Quadratic") { + m_eedfSolver->setQuadraticGrid(kTe_max, nGridCells); + } else if (energyLevelsDistribution == "Geometric") { + m_eedfSolver->setGeometricGrid(kTe_max, nGridCells); + } else { + throw CanteraError("PlasmaPhase::setParameters", + "energy_levels_distribution should be Linear, Quadratic " + "or Geometric."); + } + + // In PlasmaPhase, m_nPoints is the number of grid edges + m_nPoints = nGridCells + 1; + } + + m_electronEnergyLevels = Eigen::Map( + m_eedfSolver->getGridEdge().data(), m_nPoints); + + m_electronEnergyDist.resize(m_nPoints); + m_electronEnergyDist.setZero(); + + checkElectronEnergyLevels(); + electronEnergyLevelChanged(); + } else { + throw CanteraError("PlasmaPhase::setParameters", + "Unknown type for electron energy distribution."); } + } if (rootNode.hasKey("electron-collisions")) { From a8deddf6dd5378986eb725e856fd72fd81614fe2 Mon Sep 17 00:00:00 2001 From: Gaetanosaure Date: Wed, 6 May 2026 14:54:55 +0200 Subject: [PATCH 03/29] [thermo] Correct the instrinsic heating expression in PlasmaPhase: the energy deposited by the plasma corresponds to Joules heating, elastic power losses are a fraction of it and were thus counted twice in the previous expression --- include/cantera/thermo/PlasmaPhase.h | 2 ++ src/thermo/PlasmaPhase.cpp | 14 +++++++++++--- 2 files changed, 13 insertions(+), 3 deletions(-) diff --git a/include/cantera/thermo/PlasmaPhase.h b/include/cantera/thermo/PlasmaPhase.h index 61bb85c9338..53b3062a417 100644 --- a/include/cantera/thermo/PlasmaPhase.h +++ b/include/cantera/thermo/PlasmaPhase.h @@ -430,6 +430,8 @@ class PlasmaPhase: public IdealGasPhase double intrinsicHeating() override; + double inelasticPower(); + protected: void updateThermo() const override; diff --git a/src/thermo/PlasmaPhase.cpp b/src/thermo/PlasmaPhase.cpp index 193266cd9bc..9c833cd65e2 100644 --- a/src/thermo/PlasmaPhase.cpp +++ b/src/thermo/PlasmaPhase.cpp @@ -853,8 +853,7 @@ double PlasmaPhase::jouleHeatingPower() const return sigma * E * E; // W/m^3 } -// This is a cool function -double PlasmaPhase::intrinsicHeating() +double PlasmaPhase::inelasticPower() { // Joule heating: sigma * E^2 [W/m^3] const double qJ = jouleHeatingPower(); @@ -864,7 +863,16 @@ double PlasmaPhase::intrinsicHeating() double qElastic = elasticPowerLoss(); checkFinite(qElastic); - return qJ + qElastic; + return qJ - qElastic; +} + +double PlasmaPhase::intrinsicHeating() +{ + // Joule heating: sigma * E^2 [W/m^3] + const double qJ = jouleHeatingPower(); + checkFinite(qJ); + + return qJ; } From b96f091fe422875cc176801944e54226bd7a3e9c Mon Sep 17 00:00:00 2001 From: Gaetanosaure Date: Tue, 12 May 2026 14:34:09 +0200 Subject: [PATCH 04/29] [thermo] include the option to have the eedf grid adapting it's maximum maximum value automatically to match the problem's requirements --- .../cantera/thermo/EEDFTwoTermApproximation.h | 28 ++++ src/thermo/EEDFTwoTermApproximation.cpp | 150 ++++++++++++++++++ src/thermo/PlasmaPhase.cpp | 99 ++++++++---- 3 files changed, 249 insertions(+), 28 deletions(-) diff --git a/include/cantera/thermo/EEDFTwoTermApproximation.h b/include/cantera/thermo/EEDFTwoTermApproximation.h index 7acdb9d74f7..16309ee6be3 100644 --- a/include/cantera/thermo/EEDFTwoTermApproximation.h +++ b/include/cantera/thermo/EEDFTwoTermApproximation.h @@ -78,6 +78,18 @@ class EEDFTwoTermApproximation void setGridCache(); + void setGridType(const string& gridType); + + void setInitialGridParameters(double initialMaxEnergy, size_t nGridCells); + + void enableGridAdaptation(bool enabled); + + void setGridAdaptationParameters(bool enabled, + double minDecayDecades, + double maxDecayDecades, + double updateFactor, + size_t maxIterations); + vector getGridEdge() const { return m_gridEdge; } @@ -287,6 +299,22 @@ class EEDFTwoTermApproximation //! First call to calculateDistributionFunction bool m_first_call; + + + string m_gridType = "Linear"; + + double m_kTeMax = 60.0; + size_t m_initialGridCells = 301; + + bool m_adaptGrid = false; + double m_minEedfDecay = 8.0; + double m_maxEedfDecay = 14.0; + double m_gridUpdateFactor = 0.25; + size_t m_maxGridAdaptIterations = 5; + + void updateGrid(double maxEnergy); + + }; // end of class EEDFTwoTermApproximation } // end of namespace Cantera diff --git a/src/thermo/EEDFTwoTermApproximation.cpp b/src/thermo/EEDFTwoTermApproximation.cpp index a5a654f996e..a9f1896d04a 100644 --- a/src/thermo/EEDFTwoTermApproximation.cpp +++ b/src/thermo/EEDFTwoTermApproximation.cpp @@ -160,6 +160,68 @@ int EEDFTwoTermApproximation::calculateDistributionFunction() converge(m_f0); + if (m_adaptGrid) { + const double fFloor = 1e-300; + + for (size_t n = 0; n < m_maxGridAdaptIterations; n++) { + double fLeft = std::max(std::abs(m_f0(0)), fFloor); + double fRight = std::max(std::abs(m_f0(m_points - 1)), fFloor); + + double decades = std::log10(fLeft) - std::log10(fRight); + + if (!std::isfinite(decades)) { + throw CanteraError("EEDFTwoTermApproximation::calculateDistributionFunction", + "Non-finite EEDF decay detected during grid adaptation."); + } + + if (decades < m_minEedfDecay) { + // The right boundary is too low: the tail has not decayed enough. + double newMaxEnergy = m_kTeMax * (1.0 + m_gridUpdateFactor); + + updateGrid(newMaxEnergy); + + if (m_firstguess == "maxwell") { + for (size_t j = 0; j < m_points; j++) { + m_f0(j) = 2.0 * std::sqrt(1.0 / Pi) * + std::pow(m_init_kTe, -1.5) * + std::exp(-m_gridCenter[j] / m_init_kTe); + } + m_f0 /= norm(m_f0, m_gridCenter); + } else { + throw CanteraError("EEDFTwoTermApproximation::calculateDistributionFunction", + "Unknown EEDF first guess."); + } + + updateCrossSections(); + converge(m_f0); + + } else if (decades > m_maxEedfDecay) { + // The right boundary is unnecessarily high. + double newMaxEnergy = m_kTeMax / (1.0 + m_gridUpdateFactor); + + updateGrid(newMaxEnergy); + + if (m_firstguess == "maxwell") { + for (size_t j = 0; j < m_points; j++) { + m_f0(j) = 2.0 * std::sqrt(1.0 / Pi) * + std::pow(m_init_kTe, -1.5) * + std::exp(-m_gridCenter[j] / m_init_kTe); + } + m_f0 /= norm(m_f0, m_gridCenter); + } else { + throw CanteraError("EEDFTwoTermApproximation::calculateDistributionFunction", + "Unknown EEDF first guess."); + } + + updateCrossSections(); + converge(m_f0); + + } else { + break; + } + } + } + // write the EEDF at grid edges vector f(m_f0.data(), m_f0.data() + m_f0.rows() * m_f0.cols()); vector x(m_gridCenter.data(), m_gridCenter.data() + m_gridCenter.rows() * m_gridCenter.cols()); @@ -689,4 +751,92 @@ double EEDFTwoTermApproximation::norm(const Eigen::VectorXd& f, const Eigen::Vec return numericalQuadrature(m_quadratureMethod, p, grid); } + +void EEDFTwoTermApproximation::setGridType(const string& gridType) +{ + if (gridType != "Linear" && + gridType != "Quadratic" && + gridType != "Geometric") { + throw CanteraError("EEDFTwoTermApproximation::setGridType", + "Unknown energy grid type '{}'. Expected Linear, Quadratic or Geometric.", + gridType); + } + + m_gridType = gridType; +} + +void EEDFTwoTermApproximation::setInitialGridParameters(double initialMaxEnergy, + size_t nGridCells) +{ + if (!std::isfinite(initialMaxEnergy) || initialMaxEnergy <= 0.0) { + throw CanteraError("EEDFTwoTermApproximation::setInitialGridParameters", + "initialMaxEnergy must be finite and greater than zero."); + } + + if (nGridCells == 0) { + throw CanteraError("EEDFTwoTermApproximation::setInitialGridParameters", + "nGridCells must be greater than zero."); + } + + m_kTeMax = initialMaxEnergy; + m_initialGridCells = nGridCells; +} + +void EEDFTwoTermApproximation::enableGridAdaptation(bool enabled) +{ + m_adaptGrid = enabled; +} + +void EEDFTwoTermApproximation::setGridAdaptationParameters(bool enabled, + double minDecayDecades, + double maxDecayDecades, + double updateFactor, + size_t maxIterations) +{ + if (!std::isfinite(minDecayDecades) || !std::isfinite(maxDecayDecades) || + minDecayDecades <= 0.0 || maxDecayDecades <= minDecayDecades) { + throw CanteraError("EEDFTwoTermApproximation::setGridAdaptationParameters", + "Require 0 < min_decay_decades < max_decay_decades."); + } + + if (!std::isfinite(updateFactor) || updateFactor <= 0.0) { + throw CanteraError("EEDFTwoTermApproximation::setGridAdaptationParameters", + "update_factor must be finite and greater than zero."); + } + + if (maxIterations == 0) { + throw CanteraError("EEDFTwoTermApproximation::setGridAdaptationParameters", + "max_iterations must be greater than zero."); + } + + m_adaptGrid = enabled; + m_minEedfDecay = minDecayDecades; + m_maxEedfDecay = maxDecayDecades; + m_gridUpdateFactor = updateFactor; + m_maxGridAdaptIterations = maxIterations; +} + +void EEDFTwoTermApproximation::updateGrid(double maxEnergy) +{ + if (!std::isfinite(maxEnergy) || maxEnergy <= 0.0) { + throw CanteraError("EEDFTwoTermApproximation::updateGrid", + "Maximum grid energy must be finite and greater than zero."); + } + + m_kTeMax = maxEnergy; + + if (m_gridType == "Linear") { + setLinearGrid(m_kTeMax, m_initialGridCells); + } else if (m_gridType == "Quadratic") { + setQuadraticGrid(m_kTeMax, m_initialGridCells); + } else if (m_gridType == "Geometric") { + setGeometricGrid(m_kTeMax, m_initialGridCells); + } else { + throw CanteraError("EEDFTwoTermApproximation::updateGrid", + "Unknown energy grid type '{}'.", m_gridType); + } + + m_has_EEDF = false; +} + } diff --git a/src/thermo/PlasmaPhase.cpp b/src/thermo/PlasmaPhase.cpp index 9c833cd65e2..c83e30a2e6c 100644 --- a/src/thermo/PlasmaPhase.cpp +++ b/src/thermo/PlasmaPhase.cpp @@ -62,24 +62,27 @@ void PlasmaPhase::updateElectronEnergyDistribution() } else if (m_distributionType == "Boltzmann-two-term") { auto ierr = m_eedfSolver->calculateDistributionFunction(); if (ierr == 0) { + auto levels = m_eedfSolver->getGridEdge(); auto y = m_eedfSolver->getEEDFEdge(); - m_electronEnergyDist = Eigen::Map(y.data(), m_nPoints); + + if (levels.size() != y.size()) { + throw CanteraError("PlasmaPhase::updateElectronEnergyDistribution", + "Inconsistent EEDF solver output: grid edge size and EEDF edge size differ."); + } + + m_nPoints = levels.size(); + + m_electronEnergyLevels = Eigen::Map( + levels.data(), m_nPoints); + + m_electronEnergyDist = Eigen::Map( + y.data(), m_nPoints); + + electronEnergyLevelChanged(); } else { throw CanteraError("PlasmaPhase::updateElectronEnergyDistribution", "Call to calculateDistributionFunction failed."); } - bool validEEDF = ( - static_cast(m_electronEnergyDist.size()) == m_nPoints && - m_electronEnergyDist.allFinite() && - m_electronEnergyDist.maxCoeff() > 0.0 && - m_electronEnergyDist.sum() > 0.0 - ); - - if (validEEDF) { - updateElectronTemperatureFromEnergyDist(); - } else { - writelog("Skipping Te update: EEDF is empty, non-finite, or unnormalized.\n"); - } } else { throw CanteraError("PlasmaPhase::updateElectronEnergyDistribution", "Unknown method '{}' for determining EEDF", m_distributionType); @@ -322,15 +325,17 @@ void PlasmaPhase::setParameters(const AnyMap& phaseNode, const AnyMap& rootNode) m_eedfSolver = make_unique(this); if (eedf.hasKey("energy-levels")) { - // Mode A: user-provided energy grid edges + // Mode A: user-provided grid edges. + // In this mode, the grid is considered explicit and fixed by default. auto levels = eedf["energy-levels"].asVector(); m_eedfSolver->setCustomGrid(levels); + m_eedfSolver->enableGridAdaptation(false); m_nPoints = levels.size(); } else { - // Mode B: generated grid from readable YAML parameters + // Mode B: generated initial grid. if (!eedf.hasKey("initial_max_energy_level")) { throw CanteraError("PlasmaPhase::setParameters", "Boltzmann-two-term requires either 'energy-levels' or " @@ -343,13 +348,13 @@ void PlasmaPhase::setParameters(const AnyMap& phaseNode, const AnyMap& rootNode) "'initial_number_of_energy_grid_cells'."); } - double kTe_max = eedf["initial_max_energy_level"].asDouble(); + double initialMaxEnergy = eedf["initial_max_energy_level"].asDouble(); size_t nGridCells = static_cast( eedf["initial_number_of_energy_grid_cells"].asInt()); - if (kTe_max <= 0.0) { + if (!std::isfinite(initialMaxEnergy) || initialMaxEnergy <= 0.0) { throw CanteraError("PlasmaPhase::setParameters", - "initial_max_energy_level must be greater than zero."); + "initial_max_energy_level must be finite and greater than zero."); } if (nGridCells == 0) { @@ -366,35 +371,73 @@ void PlasmaPhase::setParameters(const AnyMap& phaseNode, const AnyMap& rootNode) "Defaulting to linear grid.\n"); } + m_eedfSolver->setGridType(energyLevelsDistribution); + m_eedfSolver->setInitialGridParameters(initialMaxEnergy, nGridCells); + if (energyLevelsDistribution == "Linear") { - m_eedfSolver->setLinearGrid(kTe_max, nGridCells); + m_eedfSolver->setLinearGrid(initialMaxEnergy, nGridCells); } else if (energyLevelsDistribution == "Quadratic") { - m_eedfSolver->setQuadraticGrid(kTe_max, nGridCells); + m_eedfSolver->setQuadraticGrid(initialMaxEnergy, nGridCells); } else if (energyLevelsDistribution == "Geometric") { - m_eedfSolver->setGeometricGrid(kTe_max, nGridCells); + m_eedfSolver->setGeometricGrid(initialMaxEnergy, nGridCells); } else { throw CanteraError("PlasmaPhase::setParameters", - "energy_levels_distribution should be Linear, Quadratic " - "or Geometric."); + "energy_levels_distribution should be Linear, Quadratic or Geometric."); + } + + if (eedf.hasKey("energy_grid_adaptation")) { + const AnyMap adapt = eedf["energy_grid_adaptation"].as(); + + bool enabled = true; + if (adapt.hasKey("enabled")) { + enabled = adapt["enabled"].asBool(); + } + + double minDecayDecades = 8.0; + if (adapt.hasKey("min_decay_decades")) { + minDecayDecades = adapt["min_decay_decades"].asDouble(); + } + + double maxDecayDecades = 14.0; + if (adapt.hasKey("max_decay_decades")) { + maxDecayDecades = adapt["max_decay_decades"].asDouble(); + } + + double updateFactor = 0.25; + if (adapt.hasKey("update_factor")) { + updateFactor = adapt["update_factor"].asDouble(); + } + + size_t maxIterations = 5; + if (adapt.hasKey("max_iterations")) { + maxIterations = static_cast( + adapt["max_iterations"].asInt()); + } + + m_eedfSolver->setGridAdaptationParameters( + enabled, minDecayDecades, maxDecayDecades, + updateFactor, maxIterations); + } else { + m_eedfSolver->enableGridAdaptation(false); } - // In PlasmaPhase, m_nPoints is the number of grid edges + // In PlasmaPhase, m_nPoints is the number of grid edges. m_nPoints = nGridCells + 1; } m_electronEnergyLevels = Eigen::Map( - m_eedfSolver->getGridEdge().data(), m_nPoints); + m_eedfSolver->getGridEdge().data(), m_eedfSolver->getGridEdge().size()); + + m_nPoints = m_eedfSolver->getGridEdge().size(); m_electronEnergyDist.resize(m_nPoints); m_electronEnergyDist.setZero(); checkElectronEnergyLevels(); electronEnergyLevelChanged(); - } else { - throw CanteraError("PlasmaPhase::setParameters", - "Unknown type for electron energy distribution."); } + } if (rootNode.hasKey("electron-collisions")) { From 5d4e0665096737fbd0e81b626070c1aada5745c2 Mon Sep 17 00:00:00 2001 From: Gaetanosaure Date: Tue, 12 May 2026 15:13:31 +0200 Subject: [PATCH 05/29] [thermo] Verify plasma EEDF input data are acceptable and correct --- include/cantera/thermo/PlasmaPhase.h | 29 +++++- src/thermo/PlasmaPhase.cpp | 132 +++++++++++++++++++++------ 2 files changed, 130 insertions(+), 31 deletions(-) diff --git a/include/cantera/thermo/PlasmaPhase.h b/include/cantera/thermo/PlasmaPhase.h index 53b3062a417..35934128624 100644 --- a/include/cantera/thermo/PlasmaPhase.h +++ b/include/cantera/thermo/PlasmaPhase.h @@ -113,7 +113,10 @@ class PlasmaPhase: public IdealGasPhase //! Get electron energy levels. //! @param levels The vector of electron energy levels (eV). Length: #m_nPoints void getElectronEnergyLevels(span levels) const { - Eigen::Map(levels.data(), levels.size()) = m_electronEnergyLevels; + checkArraySize("PlasmaPhase::getElectronEnergyLevels", + levels.size(), m_nPoints); + Eigen::Map(levels.data(), levels.size()) = + m_electronEnergyLevels; } //! Set discretized electron energy distribution. @@ -128,7 +131,10 @@ class PlasmaPhase: public IdealGasPhase //! @param distrb The vector of electron energy distribution. //! Length: #m_nPoints. void getElectronEnergyDistribution(span distrb) const { - Eigen::Map(distrb.data(), distrb.size()) = m_electronEnergyDist; + checkArraySize("PlasmaPhase::getElectronEnergyDistribution", + distrb.size(), m_nPoints); + Eigen::Map(distrb.data(), distrb.size()) = + m_electronEnergyDist; } //! Set the shape factor of isotropic electron energy distribution. @@ -368,6 +374,10 @@ class PlasmaPhase: public IdealGasPhase //! Set the absolute electric field strength [V/m] void setElectricField(double E) { + if (!std::isfinite(E) || E < 0.0) { + throw CanteraError("PlasmaPhase::setElectricField", + "Electric field must be finite and non-negative."); + } m_electricField = E; } @@ -379,8 +389,19 @@ class PlasmaPhase: public IdealGasPhase //} //! Get the reduced electric field strength [V·m²] - double reducedElectricField() const { - return m_electricField / (molarDensity() * Avogadro); + void setReducedElectricField(double EN) { + if (!std::isfinite(EN) || EN < 0.0) { + throw CanteraError("PlasmaPhase::setReducedElectricField", + "Reduced electric field must be finite and non-negative."); + } + + const double nDensity = molarDensity() * Avogadro; + if (!std::isfinite(nDensity) || nDensity <= 0.0) { + throw CanteraError("PlasmaPhase::setReducedElectricField", + "Cannot set reduced electric field with non-positive number density."); + } + + m_electricField = EN * nDensity; } //! Set reduced electric field given in [V·m²] diff --git a/src/thermo/PlasmaPhase.cpp b/src/thermo/PlasmaPhase.cpp index c83e30a2e6c..403b4a88428 100644 --- a/src/thermo/PlasmaPhase.cpp +++ b/src/thermo/PlasmaPhase.cpp @@ -91,16 +91,22 @@ void PlasmaPhase::updateElectronEnergyDistribution() electronEnergyDistributionChanged(); } -void PlasmaPhase::normalizeElectronEnergyDistribution() { - Eigen::ArrayXd eps32 = m_electronEnergyLevels.pow(3./2.); - double norm = 2./3. * numericalQuadrature(m_quadratureMethod, - m_electronEnergyDist, eps32); - if (norm < 0.0) { +void PlasmaPhase::normalizeElectronEnergyDistribution() +{ + checkElectronEnergyLevels(); + checkElectronEnergyDistribution(); + + Eigen::ArrayXd eps32 = m_electronEnergyLevels.pow(3.0 / 2.0); + double norm = 2.0 / 3.0 * numericalQuadrature( + m_quadratureMethod, m_electronEnergyDist, eps32); + + if (!std::isfinite(norm) || norm <= 0.0) { throw CanteraError("PlasmaPhase::normalizeElectronEnergyDistribution", - "The norm is negative. This might be caused by bad " - "electron energy distribution"); + "Electron energy distribution has invalid norm: {}.", norm); } + m_electronEnergyDist /= norm; + checkElectronEnergyDistribution(); } void PlasmaPhase::setElectronEnergyDistributionType(const string& type) @@ -130,7 +136,13 @@ void PlasmaPhase::setIsotropicElectronEnergyDistribution() checkElectronEnergyDistribution(); } -void PlasmaPhase::setElectronTemperature(const double Te) { +void PlasmaPhase::setElectronTemperature(const double Te) +{ + if (!std::isfinite(Te) || Te <= 0.0) { + throw CanteraError("PlasmaPhase::setElectronTemperature", + "Electron temperature must be finite and positive."); + } + m_electronTemp = Te; updateElectronEnergyDistribution(); } @@ -158,7 +170,13 @@ void PlasmaPhase::endEquilibrate() ThermoPhase::endEquilibrate(); } -void PlasmaPhase::setMeanElectronEnergy(double energy) { +void PlasmaPhase::setMeanElectronEnergy(double energy) +{ + if (!std::isfinite(energy) || energy <= 0.0) { + throw CanteraError("PlasmaPhase::setMeanElectronEnergy", + "Mean electron energy must be finite and positive."); + } + m_electronTemp = 2.0 / 3.0 * energy * ElectronCharge / Boltzmann; updateElectronEnergyDistribution(); } @@ -192,39 +210,92 @@ void PlasmaPhase::electronEnergyLevelChanged() void PlasmaPhase::checkElectronEnergyLevels() const { - Eigen::ArrayXd h = m_electronEnergyLevels.tail(m_nPoints - 1) - - m_electronEnergyLevels.head(m_nPoints - 1); - if (m_electronEnergyLevels[0] < 0.0 || (h <= 0.0).any()) { + if (m_nPoints < 2 || + static_cast(m_electronEnergyLevels.size()) != m_nPoints) { throw CanteraError("PlasmaPhase::checkElectronEnergyLevels", - "Values of electron energy levels need to be positive and " - "monotonically increasing."); + "Electron energy levels must contain at least two points and match m_nPoints."); + } + + for (size_t i = 0; i < m_nPoints; i++) { + const double eps = m_electronEnergyLevels[i]; + if (!std::isfinite(eps)) { + throw CanteraError("PlasmaPhase::checkElectronEnergyLevels", + "Electron energy level {} is non-finite.", i); + } + if (eps < 0.0) { + throw CanteraError("PlasmaPhase::checkElectronEnergyLevels", + "Electron energy levels must be non-negative."); + } + if (i > 0 && eps <= m_electronEnergyLevels[i - 1]) { + throw CanteraError("PlasmaPhase::checkElectronEnergyLevels", + "Electron energy levels must be strictly increasing."); + } } } void PlasmaPhase::checkElectronEnergyDistribution() const { - Eigen::ArrayXd h = m_electronEnergyLevels.tail(m_nPoints - 1) - - m_electronEnergyLevels.head(m_nPoints - 1); - if ((m_electronEnergyDist < 0.0).any()) { + if (m_nPoints < 2 || + static_cast(m_electronEnergyDist.size()) != m_nPoints || + static_cast(m_electronEnergyLevels.size()) != m_nPoints) { throw CanteraError("PlasmaPhase::checkElectronEnergyDistribution", - "Values of electron energy distribution cannot be negative."); + "Electron energy distribution and energy levels must have matching " + "sizes of at least two points."); } + + double sum = 0.0; + double maxVal = -BigNumber; + + for (size_t i = 0; i < m_nPoints; i++) { + const double f = m_electronEnergyDist[i]; + if (!std::isfinite(f)) { + throw CanteraError("PlasmaPhase::checkElectronEnergyDistribution", + "Electron energy distribution contains a non-finite value at index {}.", + i); + } + if (f < 0.0) { + throw CanteraError("PlasmaPhase::checkElectronEnergyDistribution", + "Electron energy distribution cannot contain negative values."); + } + sum += f; + maxVal = std::max(maxVal, f); + } + + if (!(sum > 0.0) || !(maxVal > 0.0)) { + throw CanteraError("PlasmaPhase::checkElectronEnergyDistribution", + "Electron energy distribution must contain at least one positive value."); + } + if (m_electronEnergyDist[m_nPoints - 1] > 0.01) { warn_user("PlasmaPhase::checkElectronEnergyDistribution", - "The value of the last element of electron energy distribution exceed 0.01. " - "This indicates that the value of electron energy level is not high enough " - "to contain the isotropic distribution at mean electron energy of " - "{} eV", meanElectronEnergy()); + "The last value of the electron energy distribution exceeds 0.01. " + "The electron energy grid may be too short for a mean electron energy " + "of {} eV.", meanElectronEnergy()); } } + void PlasmaPhase::setDiscretizedElectronEnergyDist(span levels, span dist) { + if (levels.size() != dist.size()) { + throw CanteraError("PlasmaPhase::setDiscretizedElectronEnergyDist", + "Energy levels and electron energy distribution must have the same size. " + "Got {} levels and {} distribution values.", levels.size(), dist.size()); + } + + if (levels.size() < 2) { + throw CanteraError("PlasmaPhase::setDiscretizedElectronEnergyDist", + "A discretized electron energy distribution requires at least two points."); + } + m_distributionType = "discretized"; m_nPoints = levels.size(); - m_electronEnergyLevels = asVectorXd(levels); - m_electronEnergyDist = asVectorXd(dist); + m_electronEnergyLevels = + Eigen::Map(levels.data(), m_nPoints); + m_electronEnergyDist = + Eigen::Map(dist.data(), m_nPoints); + checkElectronEnergyLevels(); if (m_do_normalizeElectronEnergyDist) { normalizeElectronEnergyDistribution(); @@ -248,15 +319,22 @@ void PlasmaPhase::updateElectronTemperatureFromEnergyDist() "trapezoidal", m_electronEnergyDist, eps52); } - if (epsilon_m < 0.0) { + if (!std::isfinite(epsilon_m) || epsilon_m <= 0.0) { throw CanteraError("PlasmaPhase::updateElectronTemperatureFromEnergyDist", - "The electron energy distribution produces negative electron temperature."); + "The electron energy distribution produces an invalid mean electron energy: {}.", + epsilon_m); } m_electronTemp = 2.0 / 3.0 * epsilon_m * ElectronCharge / Boltzmann; } -void PlasmaPhase::setIsotropicShapeFactor(double x) { +void PlasmaPhase::setIsotropicShapeFactor(double x) +{ + if (!std::isfinite(x) || x <= 0.0) { + throw CanteraError("PlasmaPhase::setIsotropicShapeFactor", + "The isotropic shape factor must be finite and positive."); + } + m_isotropicShapeFactor = x; updateElectronEnergyDistribution(); } From f9ddf93a8910aba063f381c3eed6b9aba3b48c0e Mon Sep 17 00:00:00 2001 From: Gaetanosaure Date: Tue, 12 May 2026 15:19:19 +0200 Subject: [PATCH 06/29] [thermo] Initialize plasma EEDF solver before loading input (PlasmaPhase constructor slightly modified) --- src/thermo/PlasmaPhase.cpp | 35 +++++++++++++++++++++++++---------- 1 file changed, 25 insertions(+), 10 deletions(-) diff --git a/src/thermo/PlasmaPhase.cpp b/src/thermo/PlasmaPhase.cpp index 403b4a88428..9e1d1de9010 100644 --- a/src/thermo/PlasmaPhase.cpp +++ b/src/thermo/PlasmaPhase.cpp @@ -23,21 +23,36 @@ namespace { PlasmaPhase::PlasmaPhase(const string& inputFile, const string& id_) { - initThermoFile(inputFile, id_); - - // initial electron temperature + // Initialize electron temperature before setParameters() can trigger EEDF updates. m_electronTemp = temperature(); - // Initialize the Boltzmann Solver + // The EEDF solver must exist before initThermoFile(), because setParameters() + // may add electron collisions and addCollision() updates the solver cache. m_eedfSolver = make_unique(this); - // Set Energy Grid (Hardcoded Defaults for Now) - double kTe_max = 60; - size_t nGridCells = 301; - m_nPoints = nGridCells + 1; + // Safe default grid used before / unless the input file overrides it. + double kTe_max = 100.0; + size_t nGridCells = 1001; + m_eedfSolver->setInitialGridParameters(kTe_max, nGridCells); m_eedfSolver->setLinearGrid(kTe_max, nGridCells); - m_electronEnergyLevels = asVectorXd(m_eedfSolver->getGridEdge()); - m_electronEnergyDist.setZero(m_nPoints); + + auto levels = m_eedfSolver->getGridEdge(); + m_nPoints = levels.size(); + m_electronEnergyLevels = Eigen::Map( + levels.data(), m_nPoints); + m_electronEnergyDist.resize(m_nPoints); + m_electronEnergyDist.setZero(); + + if (!inputFile.empty()) { + initThermoFile(inputFile, id_); + } + + // If no EEDF was supplied by input, initialize a valid default isotropic EEDF. + if (m_distributionType == "isotropic" && + static_cast(m_electronEnergyDist.size()) == m_nPoints && + m_electronEnergyDist.sum() == 0.0) { + updateElectronEnergyDistribution(); + } } PlasmaPhase::~PlasmaPhase() From d6fd205aea35de19ce4f4d89b7169c6f668af9c6 Mon Sep 17 00:00:00 2001 From: Gaetanosaure Date: Tue, 12 May 2026 15:26:47 +0200 Subject: [PATCH 07/29] [thermo] Update electron temperature after solving Boltzmann EEDF if the EEDF is correct --- src/thermo/PlasmaPhase.cpp | 13 +++++++++++++ 1 file changed, 13 insertions(+) diff --git a/src/thermo/PlasmaPhase.cpp b/src/thermo/PlasmaPhase.cpp index 9e1d1de9010..2203efcfbd1 100644 --- a/src/thermo/PlasmaPhase.cpp +++ b/src/thermo/PlasmaPhase.cpp @@ -94,6 +94,19 @@ void PlasmaPhase::updateElectronEnergyDistribution() y.data(), m_nPoints); electronEnergyLevelChanged(); + + bool validEEDF = ( + static_cast(m_electronEnergyDist.size()) == m_nPoints && + m_electronEnergyDist.allFinite() && + m_electronEnergyDist.maxCoeff() > 0.0 && + m_electronEnergyDist.sum() > 0.0 + ); + + if (validEEDF) { + updateElectronTemperatureFromEnergyDist(); + } else { + writelog("Skipping Te update: EEDF is empty, non-finite, or unnormalized.\n"); + } } else { throw CanteraError("PlasmaPhase::updateElectronEnergyDistribution", "Call to calculateDistributionFunction failed."); From 49e19bc68e879242d5076ff05d5ab7d2f9cba6f5 Mon Sep 17 00:00:00 2001 From: Gaetanosaure Date: Tue, 12 May 2026 15:37:39 +0200 Subject: [PATCH 08/29] [thermo] Impose Maxwellian distribution with T_e = T_gas at low E/N (for now hard-coded at 1 Td) --- .../cantera/thermo/EEDFTwoTermApproximation.h | 2 + include/cantera/thermo/PlasmaPhase.h | 10 +- src/thermo/EEDFTwoTermApproximation.cpp | 138 +++++++++++------- 3 files changed, 91 insertions(+), 59 deletions(-) diff --git a/include/cantera/thermo/EEDFTwoTermApproximation.h b/include/cantera/thermo/EEDFTwoTermApproximation.h index 16309ee6be3..0cc90fa190d 100644 --- a/include/cantera/thermo/EEDFTwoTermApproximation.h +++ b/include/cantera/thermo/EEDFTwoTermApproximation.h @@ -314,6 +314,8 @@ class EEDFTwoTermApproximation void updateGrid(double maxEnergy); + const double EN_min = 1e-21; // Reduced electric field below which the EEDF is not computed. + // Instead, a Maxwellian at the gas temperature is imposed. }; // end of class EEDFTwoTermApproximation diff --git a/include/cantera/thermo/PlasmaPhase.h b/include/cantera/thermo/PlasmaPhase.h index 35934128624..1d704c2ad4b 100644 --- a/include/cantera/thermo/PlasmaPhase.h +++ b/include/cantera/thermo/PlasmaPhase.h @@ -388,6 +388,11 @@ class PlasmaPhase: public IdealGasPhase // return ne / n_total; //} + //! Get the reduced electric field strength [V·m²] + double reducedElectricField() const { + return m_electricField / (molarDensity() * Avogadro); + } + //! Get the reduced electric field strength [V·m²] void setReducedElectricField(double EN) { if (!std::isfinite(EN) || EN < 0.0) { @@ -404,11 +409,6 @@ class PlasmaPhase: public IdealGasPhase m_electricField = EN * nDensity; } - //! Set reduced electric field given in [V·m²] - void setReducedElectricField(double EN) { - m_electricField = EN * molarDensity() * Avogadro; // [V/m] - } - virtual void setSolution(std::weak_ptr soln) override; /** diff --git a/src/thermo/EEDFTwoTermApproximation.cpp b/src/thermo/EEDFTwoTermApproximation.cpp index a9f1896d04a..815419d4915 100644 --- a/src/thermo/EEDFTwoTermApproximation.cpp +++ b/src/thermo/EEDFTwoTermApproximation.cpp @@ -146,93 +146,123 @@ int EEDFTwoTermApproximation::calculateDistributionFunction() updateMoleFractions(); updateCrossSections(); - if (!m_has_EEDF) { - if (m_firstguess == "maxwell") { - for (size_t j = 0; j < m_points; j++) { - m_f0(j) = 2.0 * pow(1.0 / Pi, 0.5) * pow(m_init_kTe, -3. / 2.) * - exp(-m_gridCenter[j] / m_init_kTe); - } - } else { + const double EN = m_phase->reducedElectricField(); + + auto setMaxwellian = [&](double kTe_eV) { + if (!std::isfinite(kTe_eV) || kTe_eV <= 0.0) { throw CanteraError("EEDFTwoTermApproximation::calculateDistributionFunction", - " unknown EEDF first guess"); + "Invalid kTe value for Maxwellian first guess."); } - } - converge(m_f0); - - if (m_adaptGrid) { const double fFloor = 1e-300; - for (size_t n = 0; n < m_maxGridAdaptIterations; n++) { - double fLeft = std::max(std::abs(m_f0(0)), fFloor); - double fRight = std::max(std::abs(m_f0(m_points - 1)), fFloor); + for (size_t j = 0; j < m_points; j++) { + double arg = -m_gridCenter[j] / kTe_eV; + + if (arg < -700.0) { + m_f0(j) = fFloor; + } else { + m_f0(j) = std::max(fFloor, + 2.0 * std::sqrt(1.0 / Pi) * std::pow(kTe_eV, -1.5) + * std::exp(arg)); + } + } + + double fnorm = norm(m_f0, m_gridCenter); + + if (!std::isfinite(fnorm) || fnorm <= 0.0) { + throw CanteraError("EEDFTwoTermApproximation::calculateDistributionFunction", + "Invalid norm for Maxwellian initialization."); + } - double decades = std::log10(fLeft) - std::log10(fRight); + m_f0 /= fnorm; + }; - if (!std::isfinite(decades)) { + // At very low reduced electric field, force a Maxwellian at the gas + // temperature and skip the Boltzmann convergence. + if (EN <= EN_min) { + setMaxwellian(Boltzmann * m_phase->temperature() / ElectronCharge); + } else { + // At non-zero reduced electric field, use a hot Maxwellian first guess + // only when no previously converged EEDF is available. + if (!m_has_EEDF) { + if (m_firstguess == "maxwell") { + setMaxwellian(m_init_kTe); + } else { throw CanteraError("EEDFTwoTermApproximation::calculateDistributionFunction", - "Non-finite EEDF decay detected during grid adaptation."); + "Unknown EEDF first guess '{}'.", m_firstguess); } + } - if (decades < m_minEedfDecay) { - // The right boundary is too low: the tail has not decayed enough. - double newMaxEnergy = m_kTeMax * (1.0 + m_gridUpdateFactor); + converge(m_f0); - updateGrid(newMaxEnergy); + if (m_adaptGrid) { + const double fFloor = 1e-300; - if (m_firstguess == "maxwell") { - for (size_t j = 0; j < m_points; j++) { - m_f0(j) = 2.0 * std::sqrt(1.0 / Pi) * - std::pow(m_init_kTe, -1.5) * - std::exp(-m_gridCenter[j] / m_init_kTe); - } - m_f0 /= norm(m_f0, m_gridCenter); - } else { + for (size_t n = 0; n < m_maxGridAdaptIterations; n++) { + double fLeft = std::max(std::abs(m_f0(0)), fFloor); + double fRight = std::max(std::abs(m_f0(m_points - 1)), fFloor); + + double decades = std::log10(fLeft) - std::log10(fRight); + + if (!std::isfinite(decades)) { throw CanteraError("EEDFTwoTermApproximation::calculateDistributionFunction", - "Unknown EEDF first guess."); + "Non-finite EEDF decay detected during grid adaptation."); } - updateCrossSections(); - converge(m_f0); + if (decades < m_minEedfDecay) { + // The right boundary is too low: the tail has not decayed enough. + double newMaxEnergy = m_kTeMax * (1.0 + m_gridUpdateFactor); + + updateGrid(newMaxEnergy); + + if (m_firstguess == "maxwell") { + setMaxwellian(m_init_kTe); + } else { + throw CanteraError("EEDFTwoTermApproximation::calculateDistributionFunction", + "Unknown EEDF first guess '{}'.", m_firstguess); + } + + updateCrossSections(); + converge(m_f0); - } else if (decades > m_maxEedfDecay) { - // The right boundary is unnecessarily high. - double newMaxEnergy = m_kTeMax / (1.0 + m_gridUpdateFactor); + } else if (decades > m_maxEedfDecay) { + // The right boundary is unnecessarily high. + double newMaxEnergy = m_kTeMax / (1.0 + m_gridUpdateFactor); - updateGrid(newMaxEnergy); + updateGrid(newMaxEnergy); - if (m_firstguess == "maxwell") { - for (size_t j = 0; j < m_points; j++) { - m_f0(j) = 2.0 * std::sqrt(1.0 / Pi) * - std::pow(m_init_kTe, -1.5) * - std::exp(-m_gridCenter[j] / m_init_kTe); + if (m_firstguess == "maxwell") { + setMaxwellian(m_init_kTe); + } else { + throw CanteraError("EEDFTwoTermApproximation::calculateDistributionFunction", + "Unknown EEDF first guess '{}'.", m_firstguess); } - m_f0 /= norm(m_f0, m_gridCenter); - } else { - throw CanteraError("EEDFTwoTermApproximation::calculateDistributionFunction", - "Unknown EEDF first guess."); - } - updateCrossSections(); - converge(m_f0); + updateCrossSections(); + converge(m_f0); - } else { - break; + } else { + break; + } } } } - // write the EEDF at grid edges + // Write the EEDF at grid edges. vector f(m_f0.data(), m_f0.data() + m_f0.rows() * m_f0.cols()); - vector x(m_gridCenter.data(), m_gridCenter.data() + m_gridCenter.rows() * m_gridCenter.cols()); + vector x(m_gridCenter.data(), + m_gridCenter.data() + m_gridCenter.rows() * m_gridCenter.cols()); + for (size_t i = 0; i < m_points + 1; i++) { m_f0_edge[i] = linearInterp(m_gridEdge[i], x, f); } m_has_EEDF = true; - // update electron mobility + // Update electron mobility. m_electronMobility = electronMobility(m_f0); + return 0; } From fc569813fb53a0a602598c0aaffbdfa4cb5483b2 Mon Sep 17 00:00:00 2001 From: Gaetanosaure Date: Tue, 12 May 2026 15:52:46 +0200 Subject: [PATCH 09/29] [thermo] check electron mobility value before returning it --- include/cantera/thermo/EEDFTwoTermApproximation.h | 4 +--- src/thermo/EEDFTwoTermApproximation.cpp | 9 +++++++++ 2 files changed, 10 insertions(+), 3 deletions(-) diff --git a/include/cantera/thermo/EEDFTwoTermApproximation.h b/include/cantera/thermo/EEDFTwoTermApproximation.h index 0cc90fa190d..13644713cb1 100644 --- a/include/cantera/thermo/EEDFTwoTermApproximation.h +++ b/include/cantera/thermo/EEDFTwoTermApproximation.h @@ -98,9 +98,7 @@ class EEDFTwoTermApproximation return m_f0_edge; } - double getElectronMobility() const { - return m_electronMobility; - } + double getElectronMobility() const; protected: diff --git a/src/thermo/EEDFTwoTermApproximation.cpp b/src/thermo/EEDFTwoTermApproximation.cpp index 815419d4915..1c451ee5bbc 100644 --- a/src/thermo/EEDFTwoTermApproximation.cpp +++ b/src/thermo/EEDFTwoTermApproximation.cpp @@ -869,4 +869,13 @@ void EEDFTwoTermApproximation::updateGrid(double maxEnergy) m_has_EEDF = false; } +double EEDFTwoTermApproximation::getElectronMobility() const +{ + if (!m_has_EEDF || !std::isfinite(m_electronMobility)) { + throw CanteraError("EEDFTwoTermApproximation::getElectronMobility", + "Electron mobility is not available before a valid EEDF has been computed."); + } + return m_electronMobility; +} + } From 63904875b328b002a0ac19e3ca8943f96ff4634a Mon Sep 17 00:00:00 2001 From: Gaetanosaure Date: Tue, 12 May 2026 16:00:14 +0200 Subject: [PATCH 10/29] [thermo] Exclude attachment collisions from elastic power loss --- src/thermo/PlasmaPhase.cpp | 18 +++++++++++++++++- 1 file changed, 17 insertions(+), 1 deletion(-) diff --git a/src/thermo/PlasmaPhase.cpp b/src/thermo/PlasmaPhase.cpp index 2203efcfbd1..13885824e80 100644 --- a/src/thermo/PlasmaPhase.cpp +++ b/src/thermo/PlasmaPhase.cpp @@ -778,7 +778,17 @@ void PlasmaPhase::updateElasticElectronEnergyLossCoefficients() void PlasmaPhase::updateElasticElectronEnergyLossCoefficient(size_t i) { - // @todo exclude attachment collisions + if (m_elasticElectronEnergyLossCoefficients.size() != nCollisions()) { + m_elasticElectronEnergyLossCoefficients.resize(nCollisions(), 0.0); + } + + const string kind = m_collisionRates[i]->kind(); + + if (kind == "attachment") { + m_elasticElectronEnergyLossCoefficients[i] = 0.0; + return; + } + size_t k = m_targetSpeciesIndices[i]; // Map cross sections to Eigen::ArrayXd @@ -812,6 +822,12 @@ double PlasmaPhase::elasticPowerLoss() // collisions (inelastic recoil effects). double rate = 0.0; for (size_t i = 0; i < nCollisions(); i++) { + const string kind = m_collisionRates[i]->kind(); + + if (kind == "attachment") { + continue; + } + rate += concentration(m_targetSpeciesIndices[i]) * m_elasticElectronEnergyLossCoefficients[i]; } From 66f82a35c74729d3ee5a699c8965dc068b426077 Mon Sep 17 00:00:00 2001 From: Gaetanosaure Date: Tue, 12 May 2026 16:03:23 +0200 Subject: [PATCH 11/29] [thermo] change vector handling in electronEnergyLevelChanged --- src/thermo/PlasmaPhase.cpp | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/src/thermo/PlasmaPhase.cpp b/src/thermo/PlasmaPhase.cpp index 13885824e80..5d6545cfb4b 100644 --- a/src/thermo/PlasmaPhase.cpp +++ b/src/thermo/PlasmaPhase.cpp @@ -228,10 +228,12 @@ void PlasmaPhase::electronEnergyLevelChanged() m_levelNum++; // Cross sections are interpolated on the energy levels if (m_collisions.size() > 0) { + vector energyLevels(m_nPoints); + MappedVector(energyLevels.data(), m_nPoints) = m_electronEnergyLevels; for (shared_ptr collision : m_collisions) { const auto& rate = boost::polymorphic_pointer_downcast (collision->rate()); - rate->updateInterpolatedCrossSection(asSpan(m_electronEnergyLevels)); + rate->updateInterpolatedCrossSection(energyLevels); } } } From 23df3c388e7de688d3709bce8acbf52d802050f3 Mon Sep 17 00:00:00 2001 From: Gaetanosaure Date: Wed, 13 May 2026 14:21:30 +0200 Subject: [PATCH 12/29] [YAML] Change the YAML input file format. Compatibility is retained with the older YAML file format but it throws a deprecation warning. --- .../kinetics/ElectronCollisionPlasmaRate.h | 23 +- .../cantera/thermo/EEDFTwoTermApproximation.h | 8 +- include/cantera/thermo/PlasmaPhase.h | 14 +- src/kinetics/ElectronCollisionPlasmaRate.cpp | 253 +++++++++++++-- src/thermo/EEDFTwoTermApproximation.cpp | 6 +- src/thermo/PlasmaPhase.cpp | 297 ++++++++++++++++-- 6 files changed, 533 insertions(+), 68 deletions(-) diff --git a/include/cantera/kinetics/ElectronCollisionPlasmaRate.h b/include/cantera/kinetics/ElectronCollisionPlasmaRate.h index 1cddd48df41..2a822adc70d 100644 --- a/include/cantera/kinetics/ElectronCollisionPlasmaRate.h +++ b/include/cantera/kinetics/ElectronCollisionPlasmaRate.h @@ -211,6 +211,18 @@ class ElectronCollisionPlasmaRate : public ReactionRate //! Update the value of #m_crossSectionsInterpolated [m2] void updateInterpolatedCrossSection(span); + //! Returns the name of the collision linking the reaction to the data stored in electron-collision + const string& collisionName() const; + + //! Returns the information whether the considered reaction node has cross-sections data + bool hasCrossSectionData() const; + + //! Enters the collision data found in the YAML when it is given in the old format. + void applyCollisionData(const AnyMap& node); + + //! Checks the validity of the YAML entry. + void validateCollisionData(const AnyMap& node) const; + private: //! The name of the kind of electron collision string m_kind; @@ -222,7 +234,7 @@ class ElectronCollisionPlasmaRate : public ReactionRate string m_product; //! The energy threshold of electron collision - double m_threshold; + double m_threshold = 0; //! electron energy levels [eV] vector m_energyLevels; @@ -245,6 +257,15 @@ class ElectronCollisionPlasmaRate : public ReactionRate //! This is used for the calculation of the super-elastic collision reaction //! rate coefficient. Eigen::ArrayXd m_crossSectionsOffset; + + //! The name of the collision field in the new YAML PlasmaPhase implementation + string m_collisionName; + + //! Check whether a yaml node entry offers some cross section data + bool m_hasCrossSectionData = false; + + void setDefaultThreshold(); + }; } diff --git a/include/cantera/thermo/EEDFTwoTermApproximation.h b/include/cantera/thermo/EEDFTwoTermApproximation.h index 13644713cb1..eef2128c5ca 100644 --- a/include/cantera/thermo/EEDFTwoTermApproximation.h +++ b/include/cantera/thermo/EEDFTwoTermApproximation.h @@ -69,7 +69,7 @@ class EEDFTwoTermApproximation int calculateDistributionFunction(); void setLinearGrid(double& kTe_max, size_t& ncell); - + void setQuadraticGrid(double& kTe_max, size_t& ncell); void setGeometricGrid(double& kTe_max, size_t& ncell); @@ -79,7 +79,6 @@ class EEDFTwoTermApproximation void setGridCache(); void setGridType(const string& gridType); - void setInitialGridParameters(double initialMaxEnergy, size_t nGridCells); void enableGridAdaptation(bool enabled); @@ -298,7 +297,6 @@ class EEDFTwoTermApproximation //! First call to calculateDistributionFunction bool m_first_call; - string m_gridType = "Linear"; double m_kTeMax = 60.0; @@ -312,8 +310,8 @@ class EEDFTwoTermApproximation void updateGrid(double maxEnergy); - const double EN_min = 1e-21; // Reduced electric field below which the EEDF is not computed. - // Instead, a Maxwellian at the gas temperature is imposed. + const double EN_min = 1e-21; // the reduced electric field below which the EEDF id not computed. Instead, a maxwellian at the + // gas temperature is imposed. }; // end of class EEDFTwoTermApproximation diff --git a/include/cantera/thermo/PlasmaPhase.h b/include/cantera/thermo/PlasmaPhase.h index 1d704c2ad4b..169e39870e3 100644 --- a/include/cantera/thermo/PlasmaPhase.h +++ b/include/cantera/thermo/PlasmaPhase.h @@ -393,7 +393,7 @@ class PlasmaPhase: public IdealGasPhase return m_electricField / (molarDensity() * Avogadro); } - //! Get the reduced electric field strength [V·m²] + //! Set reduced electric field given in [V·m²] void setReducedElectricField(double EN) { if (!std::isfinite(EN) || EN < 0.0) { throw CanteraError("PlasmaPhase::setReducedElectricField", @@ -583,6 +583,11 @@ class PlasmaPhase: public IdealGasPhase */ void updateElasticElectronEnergyLossCoefficients(); + //! Add an electron collision with respect to the old formats of electron collision. Function made for compatibility. + void addStandaloneElectronCollision(const AnyMap& item); + + string inferElectronCollisionKind(const shared_ptr& collision) const; + private: //! Solver used to calculate the EEDF based on electron collision rates @@ -623,6 +628,13 @@ class PlasmaPhase: public IdealGasPhase //! Work array mutable std::vector m_work; + + //! The array containing the names of the electron collisions + map m_electronCollisionDefinitions; + + //! The array containing the references of each electron collision. + set m_referencedElectronCollisions; + }; } diff --git a/src/kinetics/ElectronCollisionPlasmaRate.cpp b/src/kinetics/ElectronCollisionPlasmaRate.cpp index ce4a1ecbcd6..fa4502f8ce0 100644 --- a/src/kinetics/ElectronCollisionPlasmaRate.cpp +++ b/src/kinetics/ElectronCollisionPlasmaRate.cpp @@ -45,35 +45,90 @@ bool ElectronCollisionPlasmaData::update(const ThermoPhase& phase, const Kinetic return true; } -void ElectronCollisionPlasmaRate::setParameters(const AnyMap& node, const UnitStack& rate_units) +void ElectronCollisionPlasmaRate::setParameters(const AnyMap& node, + const UnitStack& rate_units) { + + // Diagnostics de débuggage: + + // writelog("[plasma-collision-debug] ElectronCollisionPlasmaRate::setParameters: start\n"); + + // if (node.hasKey("collision")) { + // writelog("[plasma-collision-debug] found collision reference: '{}'\n", + // node["collision"].asString()); + // } + + // bool hasEnergyLevels = node.hasKey("energy-levels"); + // bool hasCrossSections = node.hasKey("cross-sections"); + + // writelog("[plasma-collision-debug] has energy-levels: {}, has cross-sections: {}\n", + // hasEnergyLevels, hasCrossSections); + + // if (node.hasKey("kind")) { + // writelog("[plasma-collision-debug] YAML kind: '{}'\n", + // node["kind"].asString()); + // } + + // if (node.hasKey("threshold")) { + // writelog("[plasma-collision-debug] YAML threshold: {}\n", + // node["threshold"].asDouble()); + // } + + // fin des diagnostics + + ReactionRate::setParameters(node, rate_units); - if (!node.hasKey("energy-levels") && !node.hasKey("cross-sections")) { - return; - } - if (node.hasKey("kind")) { - m_kind = node["kind"].asString(); + if (node.hasKey("collision")) { + m_collisionName = node["collision"].asString(); } - if (node.hasKey("target")) { - m_target = node["target"].asString(); - } - if (node.hasKey("product")) { - m_product = node["product"].asString(); + + bool hasInlineCrossSectionData = + node.hasKey("energy-levels") && node.hasKey("cross-sections"); + + bool isNamedElectronCollision = node.hasKey("name"); + bool isCollisionReference = node.hasKey("collision"); + + if (hasInlineCrossSectionData) { + if (!isNamedElectronCollision && !isCollisionReference) { + writelog("CAREFUL! Inline electron-collision cross-section data without a " + "'name' entry are deprecated. Please move these data to a named " + "'electron-collisions' entry and reference it using 'collision'.\n"); + } + + applyCollisionData(node); } - m_energyLevels = node["energy-levels"].asVector(); - m_crossSections = node["cross-sections"].asVector(m_energyLevels.size()); - m_threshold = node.getDouble("threshold", 0.0); + + // un peu plus de débuggage. + // writelog("[plasma-collision-debug] after setParameters: collisionName='{}', kind='{}', " + // "target='{}', product='{}', threshold={}, hasCrossSectionData={}, nLevels={}, nSections={}\n", + // m_collisionName, m_kind, m_target, m_product, m_threshold, + // m_hasCrossSectionData, m_energyLevels.size(), m_crossSections.size()); } -void ElectronCollisionPlasmaRate::getParameters(AnyMap& node) const { +void ElectronCollisionPlasmaRate::getParameters(AnyMap& node) const +{ node["type"] = type(); + node["energy-levels"] = m_energyLevels; node["cross-sections"] = m_crossSections; + + if (m_threshold != 0.0) { + node["threshold"] = m_threshold; + } + if (!m_kind.empty()) { node["kind"] = m_kind; } + + if (!m_target.empty()) { + node["target"] = m_target; + } + + if (!m_product.empty()) { + node["product"] = m_product; + } } void ElectronCollisionPlasmaRate::updateInterpolatedCrossSection( @@ -172,6 +227,8 @@ void ElectronCollisionPlasmaRate::modifyRateConstants( void ElectronCollisionPlasmaRate::setContext(const Reaction& rxn, const Kinetics& kin) { + // writelog("[plasma-collision-debug] ElectronCollisionPlasmaRate::setContext: equation='{}'\n", rxn.equation()); + const ThermoPhase& thermo = kin.thermo(); // get electron species name string electronName; @@ -215,16 +272,12 @@ void ElectronCollisionPlasmaRate::setContext(const Reaction& rxn, const Kinetics } } - 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; - } - } - } + // writelog("[plasma-collision-debug] setContext result: collisionName='{}', kind='{}', " + // "threshold={}, hasCrossSectionData={}, nLevels={}\n", + // m_collisionName, m_kind, m_threshold, + // m_hasCrossSectionData, m_energyLevels.size()); + + setDefaultThreshold(); if (!rxn.reversible) { return; // end checking of forward reaction @@ -245,4 +298,154 @@ void ElectronCollisionPlasmaRate::setContext(const Reaction& rxn, const Kinetics } } +void ElectronCollisionPlasmaRate::applyCollisionData(const AnyMap& node) +{ + + // logging de debug + // writelog("[plasma-collision-debug] ElectronCollisionPlasmaRate::applyCollisionData: start\n"); + + // if (node.hasKey("name")) { + // writelog("[plasma-collision-debug] collision entry name: '{}'\n", + // node["name"].asString()); + // } + + // if (node.hasKey("target")) { + // writelog("[plasma-collision-debug] collision target: '{}'\n", + // node["target"].asString()); + // } + + // if (node.hasKey("product")) { + // writelog("[plasma-collision-debug] collision product: '{}'\n", + // node["product"].asString()); + // } + + // if (node.hasKey("kind")) { + // writelog("[plasma-collision-debug] collision kind: '{}'\n", + // node["kind"].asString()); + // } + + // fin du logging de debug + + if (node.hasKey("kind")) { + string collisionKind = node["kind"].asString(); + + if (!m_kind.empty() && m_kind != collisionKind) { + string collisionName = node.hasKey("name") ? node["name"].asString() : m_collisionName; + + throw InputFileError("applyCollisionData", node, + "Electron collision '{}' has kind '{}', but the reaction was inferred as '{}'.", + collisionName, collisionKind, m_kind); + } + + m_kind = collisionKind; + } + + + if (node.hasKey("target")) { + m_target = node["target"].asString(); + } + + if (node.hasKey("product")) { + m_product = node["product"].asString(); + } + + if (!node.hasKey("energy-levels")) { + throw InputFileError("applyCollisionData", node, "Missing 'energy-levels'"); + } + + if (!node.hasKey("cross-sections")) { + throw InputFileError("applyCollisionData", node, "Missing 'cross-sections'"); + } + + m_energyLevels = node["energy-levels"].asVector(); + m_crossSections = node["cross-sections"].asVector(m_energyLevels.size()); + m_threshold = node.getDouble("threshold", 0.0); + + // un peu plus de logging de debug + + // writelog("[plasma-collision-debug] loaded collision data: kind='{}', target='{}', " + // "product='{}', threshold={}, nLevels={}, nSections={}\n", + // m_kind, m_target, m_product, m_threshold, + // m_energyLevels.size(), m_crossSections.size()); + + // if (!m_energyLevels.empty()) { + // writelog("[plasma-collision-debug] energy range: [{}, {}]\n", + // m_energyLevels.front(), m_energyLevels.back()); + // } + + // if (!m_crossSections.empty()) { + // writelog("[plasma-collision-debug] cross-section range values: first={}, last={}\n", + // m_crossSections.front(), m_crossSections.back()); + // } + + // le debug se stop ici + + setDefaultThreshold(); + + // encore du debug + // writelog("[plasma-collision-debug] final threshold after default logic: {}\n", + // m_threshold); + + validateCollisionData(node); + m_hasCrossSectionData = true; + + +} + +void ElectronCollisionPlasmaRate::validateCollisionData(const AnyMap& node) const +{ + if (m_energyLevels.size() < 2) { + throw InputFileError("validateCollisionData" , node, "Need at least two energy levels."); + } + + if (m_energyLevels.size() != m_crossSections.size()) { + throw InputFileError("validateCollisionData" , node, "energy-levels and cross-sections size mismatch."); + } + + for (size_t i = 0; i < m_energyLevels.size(); i++) { + if (!std::isfinite(m_energyLevels[i]) || m_energyLevels[i] < 0.0) { + throw InputFileError("validateCollisionData" , node, "Inifnite or negative energy level value"); + } + if (!std::isfinite(m_crossSections[i]) || m_crossSections[i] < 0.0) { + throw InputFileError("validateCollisionData" , node, "Inifnite or negative cross-section value"); + } + if (i > 0 && m_energyLevels[i] <= m_energyLevels[i - 1]) { + throw InputFileError("validateCollisionData" , node, "energy-levels must be strictly increasing."); + } + } + + if (!std::isfinite(m_threshold) || m_threshold < 0.0) { + throw InputFileError("validateCollisionData" , node, "Inifnite or negative threshold value"); + } +} + +void ElectronCollisionPlasmaRate::setDefaultThreshold() +{ + if (m_threshold != 0.0 || m_energyLevels.empty()) { + return; + } + + if (m_kind != "excitation" && m_kind != "ionization" && m_kind != "attachment") { + return; + } + + for (double level : m_energyLevels) { + if (level > 0.0) { + m_threshold = level; + break; + } + } +} + +const string& ElectronCollisionPlasmaRate::collisionName() const +{ + return m_collisionName; +} + +bool ElectronCollisionPlasmaRate::hasCrossSectionData() const +{ + return m_hasCrossSectionData; +} + + } diff --git a/src/thermo/EEDFTwoTermApproximation.cpp b/src/thermo/EEDFTwoTermApproximation.cpp index 1c451ee5bbc..6c8abb91eae 100644 --- a/src/thermo/EEDFTwoTermApproximation.cpp +++ b/src/thermo/EEDFTwoTermApproximation.cpp @@ -43,7 +43,6 @@ void EEDFTwoTermApproximation::setLinearGrid(double& kTe_max, size_t& ncell) setGridCache(); } - void EEDFTwoTermApproximation::setQuadraticGrid(double& kTe_max, size_t& ncell) { m_points = ncell; @@ -592,7 +591,9 @@ double EEDFTwoTermApproximation::electronMobility(const Eigen::VectorXd& f0) } } double nDensity = m_phase->molarDensity() * Avogadro; - return -1./3. * m_gamma * simpson(asVectorXd(y), asVectorXd(m_gridEdge)) / nDensity; + auto f = ConstMappedVector(y.data(), y.size()); + auto x = ConstMappedVector(m_gridEdge.data(), m_gridEdge.size()); + return -1./3. * m_gamma * simpson(f, x) / nDensity; } void EEDFTwoTermApproximation::initSpeciesIndexCrossSections() @@ -781,7 +782,6 @@ double EEDFTwoTermApproximation::norm(const Eigen::VectorXd& f, const Eigen::Vec return numericalQuadrature(m_quadratureMethod, p, grid); } - void EEDFTwoTermApproximation::setGridType(const string& gridType) { if (gridType != "Linear" && diff --git a/src/thermo/PlasmaPhase.cpp b/src/thermo/PlasmaPhase.cpp index 5d6545cfb4b..a58f114929b 100644 --- a/src/thermo/PlasmaPhase.cpp +++ b/src/thermo/PlasmaPhase.cpp @@ -21,6 +21,27 @@ namespace { const double gamma = sqrt(2 * ElectronCharge / ElectronMass); } +// Initial constructor for PlasmaPhase. Might have generated some crashes so was changed by the implementation below: +// PlasmaPhase::PlasmaPhase(const string& inputFile, const string& id_) +// { +// initThermoFile(inputFile, id_); + +// // initial electron temperature +// m_electronTemp = temperature(); + +// // Initialize the Boltzmann Solver +// m_eedfSolver = make_unique(this); + +// // Set Energy Grid (Hardcoded Defaults for Now) +// double kTe_max = 60; +// size_t nGridCells = 301; +// m_nPoints = nGridCells + 1; +// m_eedfSolver->setLinearGrid(kTe_max, nGridCells); +// m_electronEnergyLevels = MappedVector(m_eedfSolver->getGridEdge().data(), m_nPoints); +// m_electronEnergyDist.setZero(m_nPoints); +// } + +// New PlasmaPhase constructor PlasmaPhase::PlasmaPhase(const string& inputFile, const string& id_) { // Initialize electron temperature before setParameters() can trigger EEDF updates. @@ -94,23 +115,22 @@ void PlasmaPhase::updateElectronEnergyDistribution() y.data(), m_nPoints); electronEnergyLevelChanged(); - - bool validEEDF = ( - static_cast(m_electronEnergyDist.size()) == m_nPoints && - m_electronEnergyDist.allFinite() && - m_electronEnergyDist.maxCoeff() > 0.0 && - m_electronEnergyDist.sum() > 0.0 - ); - - if (validEEDF) { - updateElectronTemperatureFromEnergyDist(); - } else { - writelog("Skipping Te update: EEDF is empty, non-finite, or unnormalized.\n"); - } } else { throw CanteraError("PlasmaPhase::updateElectronEnergyDistribution", "Call to calculateDistributionFunction failed."); } + bool validEEDF = ( + static_cast(m_electronEnergyDist.size()) == m_nPoints && + m_electronEnergyDist.allFinite() && + m_electronEnergyDist.maxCoeff() > 0.0 && + m_electronEnergyDist.sum() > 0.0 + ); + + if (validEEDF) { + updateElectronTemperatureFromEnergyDist(); + } else { + writelog("Skipping Te update: EEDF is empty, non-finite, or unnormalized.\n"); + } } else { throw CanteraError("PlasmaPhase::updateElectronEnergyDistribution", "Unknown method '{}' for determining EEDF", m_distributionType); @@ -263,6 +283,7 @@ void PlasmaPhase::checkElectronEnergyLevels() const } } + void PlasmaPhase::checkElectronEnergyDistribution() const { if (m_nPoints < 2 || @@ -304,7 +325,6 @@ void PlasmaPhase::checkElectronEnergyDistribution() const } } - void PlasmaPhase::setDiscretizedElectronEnergyDist(span levels, span dist) { @@ -318,14 +338,12 @@ void PlasmaPhase::setDiscretizedElectronEnergyDist(span levels, throw CanteraError("PlasmaPhase::setDiscretizedElectronEnergyDist", "A discretized electron energy distribution requires at least two points."); } - m_distributionType = "discretized"; m_nPoints = levels.size(); m_electronEnergyLevels = Eigen::Map(levels.data(), m_nPoints); m_electronEnergyDist = Eigen::Map(dist.data(), m_nPoints); - checkElectronEnergyLevels(); if (m_do_normalizeElectronEnergyDist) { normalizeElectronEnergyDistribution(); @@ -548,29 +566,136 @@ void PlasmaPhase::setParameters(const AnyMap& phaseNode, const AnyMap& rootNode) } + // writelog("[plasma-collision-debug] PlasmaPhase::setParameters: scanning electron-collisions\n"); 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; + if (item.hasKey("name")) { + string name = item["name"].asString(); + if (name.empty()) { + throw InputFileError("setParameters", rootNode, "Empty electron-collision name."); + } + if (m_electronCollisionDefinitions.count(name)) { + throw InputFileError("setParameters", rootNode, "Duplicate electron-collision name '{}'.", name); + } + m_electronCollisionDefinitions[name] = item; + } + // if (item.hasKey("name")) { + // writelog("[plasma-collision-debug] found named electron-collision: '{}'\n", + // item["name"].asString()); + // } else { + // writelog("[plasma-collision-debug] found anonymous electron-collision entry\n"); + // } + + // if (item.hasKey("target")) { + // writelog("[plasma-collision-debug] target: '{}'\n", + // item["target"].asString()); + // } + + // if (item.hasKey("kind")) { + // writelog("[plasma-collision-debug] kind: '{}'\n", + // item["kind"].asString()); + // } + } + // writelog("[plasma-collision-debug] PlasmaPhase::setParameters: registered {} named electron-collision definitions\n", + // m_electronCollisionDefinitions.size()); + } + + // in the case of a wrong combination of the two entry formats, ensure that each collision is only loaded once + // writelog("[plasma-collision-debug] PlasmaPhase::setParameters: scanning reactions for collision references\n"); + if (rootNode.hasKey("reactions")) { + for (const auto& rxnNode : rootNode["reactions"].asVector()) { + if (rxnNode.hasKey("type") + && rxnNode["type"].asString() == "electron-collision-plasma" + && rxnNode.hasKey("collision")) { + string name = rxnNode["collision"].asString(); + // writelog("[plasma-collision-debug] reaction references collision '{}'\n", name); + + if (!m_electronCollisionDefinitions.count(name)) { + // writelog("[plasma-collision-debug] reference '{}' found in electron-collisions\n", + // name); + throw InputFileError("setParameters", rootNode, + "Reaction references unknown electron collision '{}'.", name); + } + m_referencedElectronCollisions.insert(name); + } + } + } + + // writelog("[plasma-collision-debug] PlasmaPhase::setParameters: {} electron-collisions are referenced by reactions\n", + // m_referencedElectronCollisions.size()); + + if (rootNode.hasKey("electron-collisions")) { + for (const auto& item : rootNode["electron-collisions"].asVector()) { + if (item.hasKey("name")) { + string name = item["name"].asString(); + + if (m_referencedElectronCollisions.count(name)) { + // writelog("[plasma-collision-debug] skipping electron-collision '{}' because it is referenced by a reaction\n", + // name); + continue; + } + + // writelog("[plasma-collision-debug] adding named standalone electron-collision '{}'\n", + // name); } else { - products[item["target"].asString()] = 1; + // writelog("[plasma-collision-debug] adding anonymous standalone electron-collision\n"); } - products[electronSpeciesName()] = 1; - if (rate->kind() == "ionization") { - products[electronSpeciesName()] += 1; - } else if (rate->kind() == "attachment") { - products[electronSpeciesName()] -= 1; + bool hasName = item.hasKey("name"); + + if (hasName) { + string name = item["name"].asString(); + + // Nouveau format : si la collision est référencée par une réaction, + // ne pas créer de réaction synthétique ici. + if (m_referencedElectronCollisions.count(name)) { + continue; + } } - auto R = make_shared(reactants, products, rate); - addCollision(R); + + // Ancien format anonyme, ou nouvelle collision nommée non référencée : + // elle reste une collision EEDF autonome. + addStandaloneElectronCollision(item); + } } } + +void PlasmaPhase::addStandaloneElectronCollision(const AnyMap& item) +{ + auto rate = make_shared(item); + + if (!rate->hasCrossSectionData()) { + throw InputFileError("addStandaloneElectronCollision", item, + "Cross-sections are required in the reaction if you are using the deprecated input format!"); + } + + Composition reactants; + Composition products; + + string target = item["target"].asString(); + + reactants[target] = 1; + reactants[electronSpeciesName()] = 1; + + if (item.hasKey("product")) { + products[item["product"].asString()] = 1; + } else { + products[target] = 1; + } + + products[electronSpeciesName()] = 1; + + if (rate->kind() == "ionization") { + products[electronSpeciesName()] += 1; + } else if (rate->kind() == "attachment") { + products[electronSpeciesName()] -= 1; + } + + auto R = make_shared(reactants, products, rate); + addCollision(R); +} + bool PlasmaPhase::addSpecies(shared_ptr spec) { bool added = IdealGasPhase::addSpecies(spec); @@ -649,6 +774,7 @@ void PlasmaPhase::setCollisions() void PlasmaPhase::addCollision(shared_ptr collision) { + // writelog("[plasma-collision-debug] Adding electron collision '{}'\n", collision->equation()); size_t i = nCollisions(); // setup callback to signal updating the cross-section-related @@ -673,7 +799,72 @@ void PlasmaPhase::addCollision(shared_ptr collision) throw CanteraError("PlasmaPhase::addCollision", "Error identifying target for" " collision with equation '{}'", collision->equation()); } + // writelog("[plasma-collision-debug] identified target species: '{}'\n", target); + + // management of the new data file format + + auto ratePtr = std::dynamic_pointer_cast(collision->rate()); + // writelog("[plasma-collision-debug] rate before resolution: collisionName='{}', kind='{}', " + // "target='{}', product='{}', threshold={}, hasCrossSectionData={}, nLevels={}, nSections={}\n", + // ratePtr->collisionName(), ratePtr->kind(), ratePtr->target(), ratePtr->product(), + // ratePtr->threshold(), ratePtr->hasCrossSectionData(), + // ratePtr->energyLevels().size(), ratePtr->crossSections().size()); + + if (!ratePtr) { + throw CanteraError("addCollision", "The rate is not initialised"); + } + + if (!ratePtr->collisionName().empty() && !ratePtr->hasCrossSectionData()) { + auto it = m_electronCollisionDefinitions.find(ratePtr->collisionName()); + if (it == m_electronCollisionDefinitions.end()) { + throw CanteraError("addCollision", "Unknown electron collision '{}'.", ratePtr->collisionName()); + } + + // writelog("[plasma-collision-debug] resolving collision reference '{}'\n", + // ratePtr->collisionName()); + ratePtr->applyCollisionData(it->second); + // writelog("[plasma-collision-debug] rate after resolution: collisionName='{}', kind='{}', " + // "target='{}', product='{}', threshold={}, hasCrossSectionData={}, nLevels={}, nSections={}\n", + // ratePtr->collisionName(), ratePtr->kind(), ratePtr->target(), ratePtr->product(), + // ratePtr->threshold(), ratePtr->hasCrossSectionData(), + // ratePtr->energyLevels().size(), ratePtr->crossSections().size()); + } + + if (!ratePtr->hasCrossSectionData()) { + throw CanteraError("addCollision", + "ElectronCollisionPlasmaRate requires either inline cross-section data " + "or a valid 'collision' reference."); + } + + if (!ratePtr->target().empty() && ratePtr->target() != target) { + throw CanteraError("PlasmaPhase::addCollision", + "Electron collision '{}' targets '{}', but reaction '{}' uses target '{}'.", + ratePtr->collisionName(), ratePtr->target(), + collision->equation(), target); + } + // writelog("[plasma-collision-debug] target validation OK: reaction target='{}', collision target='{}'\n", + // target, ratePtr->target()); + + string kindFromReaction = inferElectronCollisionKind(collision); + string kindFromCollision = ratePtr->kind(); + + bool compatibleKind = kindFromReaction == kindFromCollision; + + if ((kindFromReaction == "elastic" || kindFromReaction == "effective") && + (kindFromCollision == "elastic" || kindFromCollision == "effective")) { + compatibleKind = true; + } + if (!compatibleKind) { + throw CanteraError("PlasmaPhase::addCollision", + "Electron collision '{}' has kind '{}', but reaction '{}' is inferred as '{}'.", + ratePtr->collisionName(), kindFromCollision, + collision->equation(), kindFromReaction); + } + // writelog("[plasma-collision-debug] kind validation OK: reaction inferred kind='{}', collision kind='{}'\n", + // kindFromReaction, ratePtr->kind()); + + // writelog("[plasma-collision-debug] adding electron collision at index {}\n", i); m_collisions.emplace_back(collision); m_collisionRates.emplace_back( std::dynamic_pointer_cast(collision->rate())); @@ -684,7 +875,6 @@ 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(); if ((kind == "effective" || kind == "elastic")) { @@ -701,12 +891,53 @@ void PlasmaPhase::addCollision(shared_ptr collision) } else { m_kInelastic.push_back(i); } +// if ((kind == "effective" || kind == "elastic")) { +// writelog("[plasma-collision-debug] classified as elastic/effective collision, index={}\n", i); +// } else { +// writelog("[plasma-collision-debug] classified as inelastic collision, index={}\n", i); +// } - auto levels = rate.energyLevels(); + auto levels = ratePtr->energyLevels(); m_energyLevels.emplace_back(levels.begin(), levels.end()); - auto sections = rate.crossSections(); + auto sections = ratePtr->crossSections(); m_crossSections.emplace_back(sections.begin(), sections.end()); m_eedfSolver->setGridCache(); + + // writelog("[plasma-collision-debug] addCollision done: nCollisions={}, nElastic={}, nInelastic={}\n", + // nCollisions(), m_kElastic.size(), m_kInelastic.size()); +} + +string PlasmaPhase::inferElectronCollisionKind( + const shared_ptr& collision) const +{ + const string eName = electronSpeciesName(); + + double nReactantElectrons = 0.0; + double nProductElectrons = 0.0; + + auto reactantElectron = collision->reactants.find(eName); + if (reactantElectron != collision->reactants.end()) { + nReactantElectrons = reactantElectron->second; + } + + auto productElectron = collision->products.find(eName); + if (productElectron != collision->products.end()) { + nProductElectrons = productElectron->second; + } + + if (nProductElectrons > nReactantElectrons) { + return "ionization"; + } + + if (nProductElectrons < nReactantElectrons) { + return "attachment"; + } + + if (collision->reactants == collision->products) { + return "elastic"; + } + + return "excitation"; } bool PlasmaPhase::updateInterpolatedCrossSection(size_t i) From e4b3596dfdc0232a86044e30c815832874159c40 Mon Sep 17 00:00:00 2001 From: Gaetanosaure Date: Wed, 13 May 2026 14:44:13 +0200 Subject: [PATCH 13/29] [YAML] added parameter "reduced_field_threshold_before_maxwellian_Td" for the user to be able to choose the reduced electric field below which they want to impose a Maxwellian --- include/cantera/thermo/EEDFTwoTermApproximation.h | 4 +++- src/thermo/EEDFTwoTermApproximation.cpp | 10 ++++++++++ src/thermo/PlasmaPhase.cpp | 9 +++++++++ 3 files changed, 22 insertions(+), 1 deletion(-) diff --git a/include/cantera/thermo/EEDFTwoTermApproximation.h b/include/cantera/thermo/EEDFTwoTermApproximation.h index eef2128c5ca..9091ad77082 100644 --- a/include/cantera/thermo/EEDFTwoTermApproximation.h +++ b/include/cantera/thermo/EEDFTwoTermApproximation.h @@ -99,6 +99,8 @@ class EEDFTwoTermApproximation double getElectronMobility() const; + void setReducedFieldThresholdBeforeMaxwellianTd(double threshold); + protected: //! Formerly options for the EEDF solver @@ -310,7 +312,7 @@ class EEDFTwoTermApproximation void updateGrid(double maxEnergy); - const double EN_min = 1e-21; // the reduced electric field below which the EEDF id not computed. Instead, a maxwellian at the + double EN_min = 1e-21; // the reduced electric field below which the EEDF id not computed. Instead, a maxwellian at the // gas temperature is imposed. }; // end of class EEDFTwoTermApproximation diff --git a/src/thermo/EEDFTwoTermApproximation.cpp b/src/thermo/EEDFTwoTermApproximation.cpp index 6c8abb91eae..6f659d67520 100644 --- a/src/thermo/EEDFTwoTermApproximation.cpp +++ b/src/thermo/EEDFTwoTermApproximation.cpp @@ -878,4 +878,14 @@ double EEDFTwoTermApproximation::getElectronMobility() const return m_electronMobility; } +// Set the reduced electric field threshold below which the EEDF is forced to be Maxwellian at the gas temperature. +// The input to this function is expected to be in Townsend. +void EEDFTwoTermApproximation::setReducedFieldThresholdBeforeMaxwellianTd(double threshold){ + if (!std::isfinite(threshold) || threshold < 0.0) { + throw CanteraError("EEDFTwoTermApproximation::setReducedFieldThresholdBeforeMaxwellianTd", + "Reduced field threshold must be finite and non-negative."); + } + EN_min = threshold*1e-21; +} + } diff --git a/src/thermo/PlasmaPhase.cpp b/src/thermo/PlasmaPhase.cpp index a58f114929b..2ddf1430c89 100644 --- a/src/thermo/PlasmaPhase.cpp +++ b/src/thermo/PlasmaPhase.cpp @@ -511,6 +511,15 @@ void PlasmaPhase::setParameters(const AnyMap& phaseNode, const AnyMap& rootNode) "energy_levels_distribution should be Linear, Quadratic or Geometric."); } + if (eedf.hasKey("reduced_field_threshold_before_maxwellian_Td")){ + double maxwellian_threshold = eedf["reduced_field_threshold_before_maxwellian_Td"].asDouble(); + if (!std::isfinite(maxwellian_threshold) || maxwellian_threshold < 0.0) { + throw CanteraError("PlasmaPhase::setParameters", + "reduced_field_threshold_before_maxwellian_Td must be finite and non-negative."); + } + m_eedfSolver->setReducedFieldThresholdBeforeMaxwellianTd(maxwellian_threshold); // The input to this function is expected to be in Townsend. + } + if (eedf.hasKey("energy_grid_adaptation")) { const AnyMap adapt = eedf["energy_grid_adaptation"].as(); From 93d515be140d73c56f6b5423502b1c047f97b22d Mon Sep 17 00:00:00 2001 From: Gaetanosaure Date: Wed, 13 May 2026 14:54:02 +0200 Subject: [PATCH 14/29] [thermo] changed the name of electronMobility to computeElectronMobility for more clarity --- include/cantera/thermo/EEDFTwoTermApproximation.h | 2 +- src/thermo/EEDFTwoTermApproximation.cpp | 8 ++++---- 2 files changed, 5 insertions(+), 5 deletions(-) diff --git a/include/cantera/thermo/EEDFTwoTermApproximation.h b/include/cantera/thermo/EEDFTwoTermApproximation.h index 9091ad77082..e1c257fad24 100644 --- a/include/cantera/thermo/EEDFTwoTermApproximation.h +++ b/include/cantera/thermo/EEDFTwoTermApproximation.h @@ -209,7 +209,7 @@ class EEDFTwoTermApproximation double electronDiffusivity(const Eigen::VectorXd& f0); //! Mobility - double electronMobility(const Eigen::VectorXd& f0); + double computeElectronMobility(const Eigen::VectorXd& f0); //! Initialize species indices associated with cross-section data void initSpeciesIndexCrossSections(); diff --git a/src/thermo/EEDFTwoTermApproximation.cpp b/src/thermo/EEDFTwoTermApproximation.cpp index 6f659d67520..a502d5ae430 100644 --- a/src/thermo/EEDFTwoTermApproximation.cpp +++ b/src/thermo/EEDFTwoTermApproximation.cpp @@ -260,7 +260,7 @@ int EEDFTwoTermApproximation::calculateDistributionFunction() m_has_EEDF = true; // Update electron mobility. - m_electronMobility = electronMobility(m_f0); + m_electronMobility = computeElectronMobility(m_f0); return 0; } @@ -462,7 +462,7 @@ SparseMat EEDFTwoTermApproximation::matrix_A(const Eigen::VectorXd& f0) double alpha; double E = m_phase->electricField(); if (m_growth == "spatial") { - double mu = electronMobility(f0); + double mu = computeElectronMobility(f0); double D = electronDiffusivity(f0); alpha = (mu * E - sqrt(pow(mu * E, 2) - 4 * D * nu * nDensity)) / 2.0 / D / nDensity; } else { @@ -578,7 +578,7 @@ double EEDFTwoTermApproximation::electronDiffusivity(const Eigen::VectorXd& f0) return 1./3. * m_gamma * simpson(f, x) / nDensity; } -double EEDFTwoTermApproximation::electronMobility(const Eigen::VectorXd& f0) +double EEDFTwoTermApproximation::computeElectronMobility(const Eigen::VectorXd& f0) { double nu = netProductionFrequency(f0); vector y(m_points + 1, 0.0); @@ -592,7 +592,7 @@ double EEDFTwoTermApproximation::electronMobility(const Eigen::VectorXd& f0) } double nDensity = m_phase->molarDensity() * Avogadro; auto f = ConstMappedVector(y.data(), y.size()); - auto x = ConstMappedVector(m_gridEdge.data(), m_gridEdge.size()); + auto x = ConstMappedVector(m_gridEdge.data(), m_gridEdge.size()); return -1./3. * m_gamma * simpson(f, x) / nDensity; } From a5665bc02fa62253238a3e0f7ca326e5f0b5e0bd Mon Sep 17 00:00:00 2001 From: Gaetanosaure Date: Wed, 13 May 2026 17:04:21 +0200 Subject: [PATCH 15/29] [general] ElectronCollisionPlasmaRate, EEDFTwoTermApproximation and PlasmaPhase code cleaning and commenting --- .../kinetics/ElectronCollisionPlasmaRate.h | 13 +- .../cantera/thermo/EEDFTwoTermApproximation.h | 115 +++++++++++++++--- include/cantera/thermo/PlasmaPhase.h | 15 ++- src/kinetics/ElectronCollisionPlasmaRate.cpp | 98 +-------------- src/thermo/EEDFTwoTermApproximation.cpp | 21 +++- src/thermo/PlasmaPhase.cpp | 101 +++------------ 6 files changed, 150 insertions(+), 213 deletions(-) diff --git a/include/cantera/kinetics/ElectronCollisionPlasmaRate.h b/include/cantera/kinetics/ElectronCollisionPlasmaRate.h index 2a822adc70d..7b8508f6c00 100644 --- a/include/cantera/kinetics/ElectronCollisionPlasmaRate.h +++ b/include/cantera/kinetics/ElectronCollisionPlasmaRate.h @@ -214,10 +214,10 @@ class ElectronCollisionPlasmaRate : public ReactionRate //! Returns the name of the collision linking the reaction to the data stored in electron-collision const string& collisionName() const; - //! Returns the information whether the considered reaction node has cross-sections data + //! Return whether this rate object has tabulated cross-section data. bool hasCrossSectionData() const; - //! Enters the collision data found in the YAML when it is given in the old format. + //! Enters the collision data found in the YAML when it is given in the old (still accepted but deprecated) format. void applyCollisionData(const AnyMap& node); //! Checks the validity of the YAML entry. @@ -252,18 +252,17 @@ class ElectronCollisionPlasmaRate : public ReactionRate //! collision cross sections [m2] after interpolation vector m_crossSectionsInterpolated; - //! collision cross section [m2] interpolated on #m_energyLevels offset by the - //! threshold energy (the first energy level). - //! This is used for the calculation of the super-elastic collision reaction - //! rate coefficient. + //! Super-elastic cross sections [m²] interpolated on the shared EEDF grid + //! after shifting the tabulated energy levels by #m_threshold. Eigen::ArrayXd m_crossSectionsOffset; //! The name of the collision field in the new YAML PlasmaPhase implementation string m_collisionName; - //! Check whether a yaml node entry offers some cross section data + //! Check whether a YAML node entry offers some cross section data bool m_hasCrossSectionData = false; + //! Get the energy threshold for the reaction from the energy levels and cross sections data if threshold is not given in the reaction void setDefaultThreshold(); }; diff --git a/include/cantera/thermo/EEDFTwoTermApproximation.h b/include/cantera/thermo/EEDFTwoTermApproximation.h index e1c257fad24..666bbfca4d6 100644 --- a/include/cantera/thermo/EEDFTwoTermApproximation.h +++ b/include/cantera/thermo/EEDFTwoTermApproximation.h @@ -61,44 +61,101 @@ class EEDFTwoTermApproximation virtual ~EEDFTwoTermApproximation() = default; - // compute the EEDF given an electric field - // CQM The solver will take the species to consider and the set of cross-sections - // from the PlasmaPhase object. - // It will write the EEDF and its grid into the PlasmaPhase object. - // Successful returns are indicated by a return value of 0. + //! Compute the electron energy distribution function for the current plasma state. + /*! + * The solver uses the electric field, species mole fractions, and electron + * collision cross sections provided by the associated PlasmaPhase object. + * On success, the internal EEDF and energy grid are updated. + * + * @return 0 if the EEDF calculation succeeds. + */ int calculateDistributionFunction(); + //! Sets a linear energy grid for the EEDF solver, defined by the maximum energy and the number of grid cells. void setLinearGrid(double& kTe_max, size_t& ncell); + //! Sets a quadratic energy grid for the EEDF solver, defined by the maximum energy and the number of grid cells. void setQuadraticGrid(double& kTe_max, size_t& ncell); + //! Sets a geometric energy grid for the EEDF solver, defined by the maximum energy and the number of grid cells. + //! @todo the geometric ratio is currently fixed to 1.05 but this parameter could be chosen by the user in the yaml file. void setGeometricGrid(double& kTe_max, size_t& ncell); + //! Sets a custom energy grid for the EEDF solver, defined by the user-provided vector of energy levels. void setCustomGrid(span levels); + //! Build or rebuild the grid-dependent cache used for scattering matrices. void setGridCache(); + //! Set the type of generated electron energy grid. + /*! + * Supported values are "Linear", "Quadratic", and "Geometric". + * + * @param gridType Type of grid spacing to use when generating or adapting + * the electron energy grid. + */ void setGridType(const string& gridType); + + //! Set the initial grid parameters used by generated EEDF grids. + /*! + * These values are used when creating Linear, Quadratic, or Geometric grids, + * and are also reused during grid adaptation. + * + * @param initialMaxEnergy Maximum electron energy of the initial grid [eV]. + * @param nGridCells Number of grid cells. The number of grid edges is + * nGridCells + 1. + */ void setInitialGridParameters(double initialMaxEnergy, size_t nGridCells); + //! This funcion enables or disables the grid adaptation based on the EEDF decay at its tail. void enableGridAdaptation(bool enabled); + //! Set parameters controlling automatic adaptation of the EEDF energy grid. + /*! + * Grid adaptation adjusts the maximum grid energy based on the number of + * decades by which the EEDF decays between the low- and high-energy ends of + * the grid. The number of grid cells is kept fixed. + * + * @param enabled True to enable grid adaptation. + * @param minDecayDecades Minimum acceptable EEDF tail decay in decades. If + * the decay is smaller, the maximum grid energy is + * increased. + * @param maxDecayDecades Maximum acceptable EEDF tail decay in decades. If + * the decay is larger, the maximum grid energy is + * decreased. + * @param updateFactor Relative factor used to increase or decrease the + * maximum grid energy during adaptation. + * @param maxIterations Maximum number of grid adaptation iterations per + * EEDF solve. + */ void setGridAdaptationParameters(bool enabled, double minDecayDecades, double maxDecayDecades, double updateFactor, size_t maxIterations); + //! Return the electron energy grid edges [eV]. + /*! + * The returned vector contains m_points + 1 values corresponding to cell + * boundaries. + */ vector getGridEdge() const { return m_gridEdge; } + //! Return the EEDF values interpolated at the electron energy grid edges. + /*! + * These values are copied back to PlasmaPhase after a successful + * Boltzmann-two-term EEDF solve. + */ vector getEEDFEdge() const { return m_f0_edge; } + //! Return the latest value of the computed electron mobility computed from the EEDF double getElectronMobility() const; + //! Controls the threshold below which the EEDF solver does not solve for the EEDF but imposes a Maxwellian distribution. void setReducedFieldThresholdBeforeMaxwellianTd(double threshold); protected: @@ -232,7 +289,7 @@ class EEDFTwoTermApproximation //! @param grid Vector representing the energy grid corresponding to f double norm(const Eigen::VectorXd& f, const Eigen::VectorXd& grid); - //! Electron mobility [m²/V·s] + //! Electron mobility computed from the most recent valid EEDF [m²/V/s]. double m_electronMobility; //! Grid of electron energy (cell center) [eV] @@ -253,7 +310,7 @@ class EEDFTwoTermApproximation //! The energy boundaries of the overlap of cell i and j vector>> m_eps; - //! normalized electron energy distribution function + //! Normalized electron energy distribution function Eigen::VectorXd m_f0; //! EEDF at grid edges (cell boundaries) @@ -266,13 +323,13 @@ class EEDFTwoTermApproximation //! energy grid vector m_totalCrossSectionEdge; - //! vector of total elastic cross section weighted with mass ratio + //! Vector of total elastic cross section weighted with mass ratio vector m_sigmaElastic; - //! list of target species indices in global Cantera numbering (1 index per cs) + //! List of target species indices in global Cantera numbering (1 index per cs) vector m_kTargets; - //! list of target species indices in local X EEDF numbering (1 index per cs) + //! List of target species indices in local X EEDF numbering (1 index per cs) vector m_klocTargets; //! Indices of species which has no cross-section data @@ -287,10 +344,14 @@ class EEDFTwoTermApproximation //! Previous mole fraction of targets used to compute eedf vector m_X_targets_prev; - //! in factor. This is used for calculating the Q matrix of - //! scattering-in processes. + //! Multiplicity factor for scattering-in terms in matrix_Q(). + /*! + * Values are 2 for ionization, 0 for attachment, and 1 for other collision + * types. + */ vector m_inFactor; + //! Conversion factor sqrt(2 e / m_e) used in electron energy-space transport terms. double m_gamma; //! flag of having an EEDF @@ -299,21 +360,41 @@ class EEDFTwoTermApproximation //! First call to calculateDistributionFunction bool m_first_call; + //! The grid type used by the EEDF solver. Supported values are "Linear", "Quadratic", and "Geometric". string m_gridType = "Linear"; + //! Initial maximum energy of the grid [eV] used when generating EEDF grids. Can be changed by user input of by grid adaptation. double m_kTeMax = 60.0; + + //! Initial / default number of grid cells. Can be changed by user input. size_t m_initialGridCells = 301; + //! Flag indicating whether grid adaptation based on EEDF tail decay is enabled. bool m_adaptGrid = false; - double m_minEedfDecay = 8.0; - double m_maxEedfDecay = 14.0; + + // Parameters controlling grid adaptation based on EEDF tail decay. They will be used if no parameters are provided by the user in the YAML file. + + //! Minimum acceptable EEDF tail decay in decades; lower values trigger grid expansion. + double m_minEedfDecay = 10.0; + + //! Maximum acceptable EEDF tail decay in decades; higher values trigger grid contraction. + double m_maxEedfDecay = 12.0; + + //! Relative factor used to increase or decrease the maximum grid energy. double m_gridUpdateFactor = 0.25; - size_t m_maxGridAdaptIterations = 5; + //! Maximum number of grid adaptation iterations per EEDF solve. + size_t m_maxGridAdaptIterations = 100; + + //! Sets a new grid for the EEDF solver once the maximum energy is updated. This function is only called in the case of grid adaptation. void updateGrid(double maxEnergy); - double EN_min = 1e-21; // the reduced electric field below which the EEDF id not computed. Instead, a maxwellian at the - // gas temperature is imposed. + //! Reduced electric field threshold below which the EEDF is forced to a Maxwellian. + /*! + * Stored internally in SI units [V m²]. User-facing input is provided in + * Townsend by setReducedFieldThresholdBeforeMaxwellianTd(). + */ + double EN_min = 1e-21; }; // end of class EEDFTwoTermApproximation diff --git a/include/cantera/thermo/PlasmaPhase.h b/include/cantera/thermo/PlasmaPhase.h index 169e39870e3..1a6f5a085c8 100644 --- a/include/cantera/thermo/PlasmaPhase.h +++ b/include/cantera/thermo/PlasmaPhase.h @@ -22,7 +22,7 @@ class ElectronCollisionPlasmaRate; //! Base class for handling plasma properties, specifically focusing on the //! electron energy distribution. /*! - * This class provides functionality to manage the the electron energy distribution + * This class provides functionality to manage the electron energy distribution * using two primary methods for defining the electron distribution and electron * temperature. * @@ -192,6 +192,7 @@ class PlasmaPhase: public IdealGasPhase return m_do_normalizeElectronEnergyDist; } + //! adds a species to the phase. Override from IdealGasPhase to take into account electrons. bool addSpecies(shared_ptr spec) override; //! Electron Temperature (K) @@ -343,7 +344,7 @@ class PlasmaPhase: public IdealGasPhase return m_levelNum; } - //! Get the indicies for inelastic electron collisions + //! Get the indices for inelastic electron collisions //! @since New in %Cantera 3.2. const vector& kInelastic() const { return m_kInelastic; @@ -381,7 +382,8 @@ class PlasmaPhase: public IdealGasPhase m_electricField = E; } - //! Calculate the degree of ionization + // @todo Add the method to compute the degree of ionization of the plasma. + //!Calculate the degree of ionization //double ionDegree() const { // double ne = concentration(m_electronSpeciesIndex); // [kmol/m³] // double n_total = molarDensity(); // [kmol/m³] @@ -449,8 +451,10 @@ class PlasmaPhase: public IdealGasPhase void endEquilibrate() override; + //! Calculates the intrinsic heating power (W/m³) of the plasma. double intrinsicHeating() override; + //! Calculates the power losses (W/m³) of the plasma electrons through inelastic collisions. double inelasticPower(); protected: @@ -461,7 +465,7 @@ class PlasmaPhase: public IdealGasPhase void electronEnergyDistributionChanged(); //! When electron energy level changed, plasma properties such as - //! electron-collision reaction rates need to be re-evaluate. + //! electron-collision reaction rates need to be re-evaluated. //! In addition, the cross-sections need to be interpolated at //! the new level. void electronEnergyLevelChanged(); @@ -475,7 +479,7 @@ class PlasmaPhase: public IdealGasPhase //! Check the electron energy distribution /*! - * This method check the electron energy distribution for the criteria + * This method checks the electron energy distribution for the criteria * below. * * 1. The values of electron energy distribution cannot be negative. @@ -550,6 +554,7 @@ class PlasmaPhase: public IdealGasPhase //! where i is the specific process, j is the index of vector. Unit: [eV] vector> m_energyLevels; + // @todo Add a variable to track the ionization degree of the plasma. //! ionization degree for the electron-electron collisions (tmp is the previous one) //double m_ionDegree = 0.0; diff --git a/src/kinetics/ElectronCollisionPlasmaRate.cpp b/src/kinetics/ElectronCollisionPlasmaRate.cpp index fa4502f8ce0..6223ba64aac 100644 --- a/src/kinetics/ElectronCollisionPlasmaRate.cpp +++ b/src/kinetics/ElectronCollisionPlasmaRate.cpp @@ -31,12 +31,12 @@ bool ElectronCollisionPlasmaData::update(const ThermoPhase& phase, const Kinetic return false; } - // Update distribution + // Update the cached eedf from the plasma phase. m_dist_number = pp.distributionNumber(); distribution.resize(pp.nElectronEnergyLevels()); pp.getElectronEnergyDistribution(distribution); - // Update energy levels + // Update the cached energy grid only when the phase grid revision changes. if (pp.levelNumber() != levelNumber || energyLevels.empty()) { levelNumber = pp.levelNumber(); energyLevels.resize(pp.nElectronEnergyLevels()); @@ -49,34 +49,6 @@ void ElectronCollisionPlasmaRate::setParameters(const AnyMap& node, const UnitStack& rate_units) { - // Diagnostics de débuggage: - - // writelog("[plasma-collision-debug] ElectronCollisionPlasmaRate::setParameters: start\n"); - - // if (node.hasKey("collision")) { - // writelog("[plasma-collision-debug] found collision reference: '{}'\n", - // node["collision"].asString()); - // } - - // bool hasEnergyLevels = node.hasKey("energy-levels"); - // bool hasCrossSections = node.hasKey("cross-sections"); - - // writelog("[plasma-collision-debug] has energy-levels: {}, has cross-sections: {}\n", - // hasEnergyLevels, hasCrossSections); - - // if (node.hasKey("kind")) { - // writelog("[plasma-collision-debug] YAML kind: '{}'\n", - // node["kind"].asString()); - // } - - // if (node.hasKey("threshold")) { - // writelog("[plasma-collision-debug] YAML threshold: {}\n", - // node["threshold"].asDouble()); - // } - - // fin des diagnostics - - ReactionRate::setParameters(node, rate_units); if (node.hasKey("collision")) { @@ -99,12 +71,6 @@ void ElectronCollisionPlasmaRate::setParameters(const AnyMap& node, applyCollisionData(node); } - - // un peu plus de débuggage. - // writelog("[plasma-collision-debug] after setParameters: collisionName='{}', kind='{}', " - // "target='{}', product='{}', threshold={}, hasCrossSectionData={}, nLevels={}, nSections={}\n", - // m_collisionName, m_kind, m_target, m_product, m_threshold, - // m_hasCrossSectionData, m_energyLevels.size(), m_crossSections.size()); } void ElectronCollisionPlasmaRate::getParameters(AnyMap& node) const @@ -197,7 +163,7 @@ void ElectronCollisionPlasmaRate::modifyRateConstants( m_crossSectionsOffset.resize(shared_data.energyLevels.size()); for (size_t i = 1; i < m_energyLevels.size(); i++) { // The energy levels are offset by the first energy level (threshold) - superElasticEnergyLevels.push_back(m_energyLevels[i] - m_energyLevels[0]); + superElasticEnergyLevels.push_back(m_energyLevels[i] - m_threshold); } for (size_t i = 0; i < shared_data.energyLevels.size(); i++) { // The interpolated super-elastic cross section is evaluated @@ -214,20 +180,19 @@ void ElectronCollisionPlasmaRate::modifyRateConstants( shared_data.energyLevels.data(), shared_data.energyLevels.size() ); - // Map energyLevels in Eigen::ArrayXd + // Map the electron energy distribution to Eigen::ArrayXd. auto distribution = Eigen::Map( shared_data.distribution.data(), shared_data.distribution.size() ); // unit in kmol/m3/s kr = pow(2.0 * ElectronCharge / ElectronMass, 0.5) * Avogadro * - simpson((eps + m_energyLevels[0]).cwiseProduct( + simpson((eps + m_threshold).cwiseProduct( distribution.cwiseProduct(m_crossSectionsOffset)), eps); } void ElectronCollisionPlasmaRate::setContext(const Reaction& rxn, const Kinetics& kin) { - // writelog("[plasma-collision-debug] ElectronCollisionPlasmaRate::setContext: equation='{}'\n", rxn.equation()); const ThermoPhase& thermo = kin.thermo(); // get electron species name @@ -272,11 +237,6 @@ void ElectronCollisionPlasmaRate::setContext(const Reaction& rxn, const Kinetics } } - // writelog("[plasma-collision-debug] setContext result: collisionName='{}', kind='{}', " - // "threshold={}, hasCrossSectionData={}, nLevels={}\n", - // m_collisionName, m_kind, m_threshold, - // m_hasCrossSectionData, m_energyLevels.size()); - setDefaultThreshold(); if (!rxn.reversible) { @@ -301,31 +261,6 @@ void ElectronCollisionPlasmaRate::setContext(const Reaction& rxn, const Kinetics void ElectronCollisionPlasmaRate::applyCollisionData(const AnyMap& node) { - // logging de debug - // writelog("[plasma-collision-debug] ElectronCollisionPlasmaRate::applyCollisionData: start\n"); - - // if (node.hasKey("name")) { - // writelog("[plasma-collision-debug] collision entry name: '{}'\n", - // node["name"].asString()); - // } - - // if (node.hasKey("target")) { - // writelog("[plasma-collision-debug] collision target: '{}'\n", - // node["target"].asString()); - // } - - // if (node.hasKey("product")) { - // writelog("[plasma-collision-debug] collision product: '{}'\n", - // node["product"].asString()); - // } - - // if (node.hasKey("kind")) { - // writelog("[plasma-collision-debug] collision kind: '{}'\n", - // node["kind"].asString()); - // } - - // fin du logging de debug - if (node.hasKey("kind")) { string collisionKind = node["kind"].asString(); @@ -361,31 +296,8 @@ void ElectronCollisionPlasmaRate::applyCollisionData(const AnyMap& node) m_crossSections = node["cross-sections"].asVector(m_energyLevels.size()); m_threshold = node.getDouble("threshold", 0.0); - // un peu plus de logging de debug - - // writelog("[plasma-collision-debug] loaded collision data: kind='{}', target='{}', " - // "product='{}', threshold={}, nLevels={}, nSections={}\n", - // m_kind, m_target, m_product, m_threshold, - // m_energyLevels.size(), m_crossSections.size()); - - // if (!m_energyLevels.empty()) { - // writelog("[plasma-collision-debug] energy range: [{}, {}]\n", - // m_energyLevels.front(), m_energyLevels.back()); - // } - - // if (!m_crossSections.empty()) { - // writelog("[plasma-collision-debug] cross-section range values: first={}, last={}\n", - // m_crossSections.front(), m_crossSections.back()); - // } - - // le debug se stop ici - setDefaultThreshold(); - // encore du debug - // writelog("[plasma-collision-debug] final threshold after default logic: {}\n", - // m_threshold); - validateCollisionData(node); m_hasCrossSectionData = true; diff --git a/src/thermo/EEDFTwoTermApproximation.cpp b/src/thermo/EEDFTwoTermApproximation.cpp index a502d5ae430..615c9f1e93c 100644 --- a/src/thermo/EEDFTwoTermApproximation.cpp +++ b/src/thermo/EEDFTwoTermApproximation.cpp @@ -21,7 +21,7 @@ typedef Eigen::SparseMatrix SparseMat; EEDFTwoTermApproximation::EEDFTwoTermApproximation(PlasmaPhase* s) { - // store a pointer to s. + // Store the PlasmaPhase context used by the solver (pointer to s). m_phase = s; m_first_call = true; m_has_EEDF = false; @@ -75,8 +75,8 @@ void EEDFTwoTermApproximation::setGeometricGrid(double& kTe_max, size_t& ncell) m_f0.resize(m_points); m_f0_edge.resize(m_points + 1); + // @todo Make the geometric-grid ratio configurable from input. double ratio = 1.05; - double denominator = std::pow(ratio, m_points) - 1.0; if (std::abs(denominator) < 1e-14) { @@ -147,6 +147,7 @@ int EEDFTwoTermApproximation::calculateDistributionFunction() const double EN = m_phase->reducedElectricField(); + // Helper used to impose a Maxwellian EEDF. auto setMaxwellian = [&](double kTe_eV) { if (!std::isfinite(kTe_eV) || kTe_eV <= 0.0) { throw CanteraError("EEDFTwoTermApproximation::calculateDistributionFunction", @@ -178,7 +179,8 @@ int EEDFTwoTermApproximation::calculateDistributionFunction() }; // At very low reduced electric field, force a Maxwellian at the gas - // temperature and skip the Boltzmann convergence. + // temperature and skip the Boltzmann convergence to avoid wrong EEDF convergence and numerical instabilities + // when integrating over time. if (EN <= EN_min) { setMaxwellian(Boltzmann * m_phase->temperature() / ElectronCharge); } else { @@ -195,6 +197,11 @@ int EEDFTwoTermApproximation::calculateDistributionFunction() converge(m_f0); + // Grid adaptation based on EEDF tail decay. If enabled, this will iteratively adjust the grid + // until the EEDF tail decays within the specified bounds given in the YAML file. + // @todo Implement a more robust version which also varies the number of grid points if necessary. + // The current version only adjusts the maximum energy of the grid. + if (m_adaptGrid) { const double fFloor = 1e-300; @@ -299,7 +306,7 @@ void EEDFTwoTermApproximation::converge(Eigen::VectorXd& f0) if (err1 < m_rtol) { break; } else if (n == m_maxn - 1) { - throw CanteraError("WeaklyIonizedGas::converge", "Convergence failed"); + throw CanteraError("EEDFTwoTermApproximation::converge", "Convergence failed"); } } } @@ -316,7 +323,7 @@ Eigen::VectorXd EEDFTwoTermApproximation::iterate(const Eigen::VectorXd& f0, dou for (size_t k : m_phase->kInelastic()) { SparseMat Q_k = matrix_Q(g, k); SparseMat P_k = matrix_P(g, k); - PQ += (matrix_Q(g, k) - matrix_P(g, k)) * m_X_targets[m_klocTargets[k]]; + PQ += (Q_k - P_k) * m_X_targets[m_klocTargets[k]]; } SparseMat A = matrix_A(f0); @@ -526,7 +533,7 @@ SparseMat EEDFTwoTermApproximation::matrix_A(const Eigen::VectorXd& f0) SparseMat A(m_points, m_points); A.setFromTriplets(tripletList.begin(), tripletList.end()); - //plus G + // plus G SparseMat G(m_points, m_points); if (m_growth == "temporal") { for (size_t i = 0; i < m_points; i++) { @@ -650,6 +657,7 @@ void EEDFTwoTermApproximation::updateCrossSections() } // Update the species mole fractions used for EEDF computation +// Renormalize over species with electron-collision cross-section data. void EEDFTwoTermApproximation::updateMoleFractions() { double tmp_sum = 0.0; @@ -866,6 +874,7 @@ void EEDFTwoTermApproximation::updateGrid(double maxEnergy) "Unknown energy grid type '{}'.", m_gridType); } + // put back the flag at false because since the grid has changed the EEDF must be computed again. m_has_EEDF = false; } diff --git a/src/thermo/PlasmaPhase.cpp b/src/thermo/PlasmaPhase.cpp index 2ddf1430c89..58143707a15 100644 --- a/src/thermo/PlasmaPhase.cpp +++ b/src/thermo/PlasmaPhase.cpp @@ -21,27 +21,6 @@ namespace { const double gamma = sqrt(2 * ElectronCharge / ElectronMass); } -// Initial constructor for PlasmaPhase. Might have generated some crashes so was changed by the implementation below: -// PlasmaPhase::PlasmaPhase(const string& inputFile, const string& id_) -// { -// initThermoFile(inputFile, id_); - -// // initial electron temperature -// m_electronTemp = temperature(); - -// // Initialize the Boltzmann Solver -// m_eedfSolver = make_unique(this); - -// // Set Energy Grid (Hardcoded Defaults for Now) -// double kTe_max = 60; -// size_t nGridCells = 301; -// m_nPoints = nGridCells + 1; -// m_eedfSolver->setLinearGrid(kTe_max, nGridCells); -// m_electronEnergyLevels = MappedVector(m_eedfSolver->getGridEdge().data(), m_nPoints); -// m_electronEnergyDist.setZero(m_nPoints); -// } - -// New PlasmaPhase constructor PlasmaPhase::PlasmaPhase(const string& inputFile, const string& id_) { // Initialize electron temperature before setParameters() can trigger EEDF updates. @@ -96,6 +75,9 @@ void PlasmaPhase::updateElectronEnergyDistribution() } else if (m_distributionType == "isotropic") { setIsotropicElectronEnergyDistribution(); } else if (m_distributionType == "Boltzmann-two-term") { + // The two-term Boltzmann solver may adapt the electron energy grid, so both + // the grid and the distribution are copied back before invalidating dependent + // cross-section data. auto ierr = m_eedfSolver->calculateDistributionFunction(); if (ierr == 0) { auto levels = m_eedfSolver->getGridEdge(); @@ -129,6 +111,9 @@ void PlasmaPhase::updateElectronEnergyDistribution() if (validEEDF) { updateElectronTemperatureFromEnergyDist(); } else { + // Keep the previous electron temperature if the solver did not return a usable + // distribution. This avoids replacing Te with a value derived from invalid data. + // The user will be warned of this behaviour. writelog("Skipping Te update: EEDF is empty, non-finite, or unnormalized.\n"); } } else { @@ -575,7 +560,6 @@ void PlasmaPhase::setParameters(const AnyMap& phaseNode, const AnyMap& rootNode) } - // writelog("[plasma-collision-debug] PlasmaPhase::setParameters: scanning electron-collisions\n"); if (rootNode.hasKey("electron-collisions")) { for (const auto& item : rootNode["electron-collisions"].asVector()) { if (item.hasKey("name")) { @@ -588,40 +572,18 @@ void PlasmaPhase::setParameters(const AnyMap& phaseNode, const AnyMap& rootNode) } m_electronCollisionDefinitions[name] = item; } - // if (item.hasKey("name")) { - // writelog("[plasma-collision-debug] found named electron-collision: '{}'\n", - // item["name"].asString()); - // } else { - // writelog("[plasma-collision-debug] found anonymous electron-collision entry\n"); - // } - - // if (item.hasKey("target")) { - // writelog("[plasma-collision-debug] target: '{}'\n", - // item["target"].asString()); - // } - - // if (item.hasKey("kind")) { - // writelog("[plasma-collision-debug] kind: '{}'\n", - // item["kind"].asString()); - // } } - // writelog("[plasma-collision-debug] PlasmaPhase::setParameters: registered {} named electron-collision definitions\n", - // m_electronCollisionDefinitions.size()); } - // in the case of a wrong combination of the two entry formats, ensure that each collision is only loaded once - // writelog("[plasma-collision-debug] PlasmaPhase::setParameters: scanning reactions for collision references\n"); + // in the case of a wrong combination of the two entry formats (should the user mix new and old formats somehow), ensure that each collision is only loaded once if (rootNode.hasKey("reactions")) { for (const auto& rxnNode : rootNode["reactions"].asVector()) { if (rxnNode.hasKey("type") && rxnNode["type"].asString() == "electron-collision-plasma" && rxnNode.hasKey("collision")) { string name = rxnNode["collision"].asString(); - // writelog("[plasma-collision-debug] reaction references collision '{}'\n", name); if (!m_electronCollisionDefinitions.count(name)) { - // writelog("[plasma-collision-debug] reference '{}' found in electron-collisions\n", - // name); throw InputFileError("setParameters", rootNode, "Reaction references unknown electron collision '{}'.", name); } @@ -630,39 +592,28 @@ void PlasmaPhase::setParameters(const AnyMap& phaseNode, const AnyMap& rootNode) } } - // writelog("[plasma-collision-debug] PlasmaPhase::setParameters: {} electron-collisions are referenced by reactions\n", - // m_referencedElectronCollisions.size()); - if (rootNode.hasKey("electron-collisions")) { for (const auto& item : rootNode["electron-collisions"].asVector()) { if (item.hasKey("name")) { string name = item["name"].asString(); if (m_referencedElectronCollisions.count(name)) { - // writelog("[plasma-collision-debug] skipping electron-collision '{}' because it is referenced by a reaction\n", - // name); continue; } - // writelog("[plasma-collision-debug] adding named standalone electron-collision '{}'\n", - // name); - } else { - // writelog("[plasma-collision-debug] adding anonymous standalone electron-collision\n"); } bool hasName = item.hasKey("name"); if (hasName) { string name = item["name"].asString(); - // Nouveau format : si la collision est référencée par une réaction, - // ne pas créer de réaction synthétique ici. + // New YAML format: if the collision is referenced by a reaction, do not create a synthetic reaction here. if (m_referencedElectronCollisions.count(name)) { continue; } } - // Ancien format anonyme, ou nouvelle collision nommée non référencée : - // elle reste une collision EEDF autonome. + // Old YAML format for anonymous collisions or unreferenced named collisions: they remain as standalone EEDF collisions. addStandaloneElectronCollision(item); } @@ -783,7 +734,7 @@ void PlasmaPhase::setCollisions() void PlasmaPhase::addCollision(shared_ptr collision) { - // writelog("[plasma-collision-debug] Adding electron collision '{}'\n", collision->equation()); + size_t i = nCollisions(); // setup callback to signal updating the cross-section-related @@ -808,16 +759,10 @@ void PlasmaPhase::addCollision(shared_ptr collision) throw CanteraError("PlasmaPhase::addCollision", "Error identifying target for" " collision with equation '{}'", collision->equation()); } - // writelog("[plasma-collision-debug] identified target species: '{}'\n", target); // management of the new data file format auto ratePtr = std::dynamic_pointer_cast(collision->rate()); - // writelog("[plasma-collision-debug] rate before resolution: collisionName='{}', kind='{}', " - // "target='{}', product='{}', threshold={}, hasCrossSectionData={}, nLevels={}, nSections={}\n", - // ratePtr->collisionName(), ratePtr->kind(), ratePtr->target(), ratePtr->product(), - // ratePtr->threshold(), ratePtr->hasCrossSectionData(), - // ratePtr->energyLevels().size(), ratePtr->crossSections().size()); if (!ratePtr) { throw CanteraError("addCollision", "The rate is not initialised"); @@ -829,14 +774,7 @@ void PlasmaPhase::addCollision(shared_ptr collision) throw CanteraError("addCollision", "Unknown electron collision '{}'.", ratePtr->collisionName()); } - // writelog("[plasma-collision-debug] resolving collision reference '{}'\n", - // ratePtr->collisionName()); ratePtr->applyCollisionData(it->second); - // writelog("[plasma-collision-debug] rate after resolution: collisionName='{}', kind='{}', " - // "target='{}', product='{}', threshold={}, hasCrossSectionData={}, nLevels={}, nSections={}\n", - // ratePtr->collisionName(), ratePtr->kind(), ratePtr->target(), ratePtr->product(), - // ratePtr->threshold(), ratePtr->hasCrossSectionData(), - // ratePtr->energyLevels().size(), ratePtr->crossSections().size()); } if (!ratePtr->hasCrossSectionData()) { @@ -851,8 +789,6 @@ void PlasmaPhase::addCollision(shared_ptr collision) ratePtr->collisionName(), ratePtr->target(), collision->equation(), target); } - // writelog("[plasma-collision-debug] target validation OK: reaction target='{}', collision target='{}'\n", - // target, ratePtr->target()); string kindFromReaction = inferElectronCollisionKind(collision); string kindFromCollision = ratePtr->kind(); @@ -870,10 +806,7 @@ void PlasmaPhase::addCollision(shared_ptr collision) ratePtr->collisionName(), kindFromCollision, collision->equation(), kindFromReaction); } - // writelog("[plasma-collision-debug] kind validation OK: reaction inferred kind='{}', collision kind='{}'\n", - // kindFromReaction, ratePtr->kind()); - // writelog("[plasma-collision-debug] adding electron collision at index {}\n", i); m_collisions.emplace_back(collision); m_collisionRates.emplace_back( std::dynamic_pointer_cast(collision->rate())); @@ -900,11 +833,7 @@ void PlasmaPhase::addCollision(shared_ptr collision) } else { m_kInelastic.push_back(i); } -// if ((kind == "effective" || kind == "elastic")) { -// writelog("[plasma-collision-debug] classified as elastic/effective collision, index={}\n", i); -// } else { -// writelog("[plasma-collision-debug] classified as inelastic collision, index={}\n", i); -// } + auto levels = ratePtr->energyLevels(); m_energyLevels.emplace_back(levels.begin(), levels.end()); @@ -912,10 +841,12 @@ void PlasmaPhase::addCollision(shared_ptr collision) m_crossSections.emplace_back(sections.begin(), sections.end()); m_eedfSolver->setGridCache(); - // writelog("[plasma-collision-debug] addCollision done: nCollisions={}, nElastic={}, nInelastic={}\n", - // nCollisions(), m_kElastic.size(), m_kInelastic.size()); } +// Infer the collision kind from the electron balance: producing electrons +// indicates ionization, consuming electrons indicates attachment, unchanged +// reactants/products indicate elastic scattering, and the remaining case is +// treated as excitation. string PlasmaPhase::inferElectronCollisionKind( const shared_ptr& collision) const { @@ -993,7 +924,7 @@ void PlasmaPhase::updateElasticElectronEnergyLossCoefficients() static const int cacheId = m_cache.getId(); CachedScalar last_stateNum = m_cache.getScalar(cacheId); - // combine the distribution and energy level number + // combine the distribution and energy level number to get a single cache variable. int stateNum = m_distNum + m_levelNum; vector interpChanged(m_collisions.size()); From d042c3789a67c01e86e189dd1bf2314e165b3796 Mon Sep 17 00:00:00 2001 From: Gaetanosaure Date: Thu, 14 May 2026 13:19:07 +0200 Subject: [PATCH 16/29] [Kinetics] Added the class DetailedVVVTRate to be able to compute the reaction rates from detailed vibrational-translational relaxation and vibrational-vibrational relaxation when vibrational species are involved in the phase --- include/cantera/kinetics/DetailedVVVTRate.h | 208 ++++++++++++++++++++ src/kinetics/DetailedVVVTRate.cpp | 194 ++++++++++++++++++ src/kinetics/ReactionRateFactory.cpp | 6 + test/kinetics/kineticsFromScratch.cpp | 1 + test/kinetics/kineticsFromYaml.cpp | 1 + 5 files changed, 410 insertions(+) create mode 100644 include/cantera/kinetics/DetailedVVVTRate.h create mode 100644 src/kinetics/DetailedVVVTRate.cpp diff --git a/include/cantera/kinetics/DetailedVVVTRate.h b/include/cantera/kinetics/DetailedVVVTRate.h new file mode 100644 index 00000000000..43cba7d2e20 --- /dev/null +++ b/include/cantera/kinetics/DetailedVVVTRate.h @@ -0,0 +1,208 @@ +//! @file DetailedVVVTRate.h +//! Header for plasma reaction rates corresponding to detailed vibration handling. +//! +//! This file is part of Cantera. See License.txt in the top-level directory or +//! at https://cantera.org/license.txt for license and copyright information. + +#ifndef CT_DETAILEDVVVTRATE_H +#define CT_DETAILEDVVVTRATE_H + +#include "Arrhenius.h" + +#include + +namespace Cantera +{ + +//! Data container holding shared data specific to DetailedVVVTRate. +struct DetailedVibData : public ReactionData +{ + //! Update cached temperature-dependent data. + /*! + * @param phase Thermodynamic phase used to retrieve the gas temperature. + * @param kin Kinetics object. Not used here, but required by the + * `ReactionData` interface. + * + * @return `true` if the temperature has changed and rates need to be + * recomputed; `false` otherwise. + */ + bool update(const ThermoPhase& phase, const Kinetics& kin) override; + + using ReactionData::update; +}; + + +//! Detailed VV-VT reaction rate. +/*! + * This reaction rate implements the following temperature-dependent + * expression: + * + * @f[ + * k_f = + * scaling \, A \, + * \exp \left( + * b \ln T + * + B + * + C T^{-1/3} + * + D T^{-2/3} + * \right) + * @f] + * + * where `T` is the gas temperature in K. + * + * Unit conventions: + * + * - `A` uses the standard Cantera rate coefficient units. Its units depend + * on the reaction order and dimensionality, and are converted by + * `ArrheniusBase`. + * - `b` is dimensionless. + * - `B` is dimensionless. + * - `C` is interpreted as being expressed in K^(1/3), assuming `T` is in K. + * - `D` is interpreted as being expressed in K^(2/3), assuming `T` is in K. + * - `scaling` is dimensionless. + * + * Important: `B`, `C`, `D`, and `scaling` are read as raw floating-point + * numbers from YAML. They are not converted by Cantera's unit system. + * + * @ingroup arrheniusGroup + */ +class DetailedVVVTRate : public ArrheniusBase +{ +public: + //! Default constructor. + DetailedVVVTRate(); + + //! Constructor. + /*! + * @param A Pre-exponential factor. The unit system is (kmol, m, s); + * actual units depend on the reaction order and the + * dimensionality, surface or bulk. + * @param B Dimensionless constant in the exponential. + * @param C Constant multiplying T^(-1/3), interpreted in K^(1/3). + * @param D Constant multiplying T^(-2/3), interpreted in K^(2/3). + * @param b Dimensionless temperature exponent. + * @param scaling Dimensionless scaling factor from the harmonic oscillator + * theory. + */ + DetailedVVVTRate(double A, double B, double C, double D, + double b, double scaling = 1.0); + + //! Constructor based on AnyMap content. + /*! + * Expected YAML form: + * + * @code{.yaml} + * type: detailed-vv-vt + * rate-constant: + * A: ... + * b: ... + * B: ... + * C: ... + * D: ... + * scaling: ... + * @endcode + * + * The `scaling` entry is optional and defaults to 1.0. + */ + explicit DetailedVVVTRate(const AnyMap& node, + const UnitStack& rate_units = {}); + + //! Set rate parameters from an AnyMap. + /*! + * `A` and `b` are handled by `ArrheniusBase`. + * `B`, `C`, `D`, and `scaling` are handled by this class. + */ + void setParameters(const AnyMap& node, + const UnitStack& rate_units) override; + + //! Get rate parameters for YAML serialization. + /*! + * This method is needed so that Cantera can reconstruct or export the + * full custom rate, including `B`, `C`, `D`, and `scaling`. + */ + void getParameters(AnyMap& node) const override; + + //! Create a rate evaluator for this reaction rate type. + unique_ptr newMultiRate() const override { + return make_unique>(); + } + + //! String identifying this reaction rate specialization. + const string type() const override { + return "detailed-vv-vt"; + } + + //! Set context of reaction rate evaluation. + /*! + * This rate is intended for irreversible non-equilibrium plasma reactions. + * Reversible reactions are rejected because the reverse rate cannot be + * obtained from conventional thermochemistry for this model. + */ + void setContext(const Reaction& rxn, const Kinetics& kin) override; + + //! Evaluate reaction rate. + /*! + * @param shared_data Data shared by all reactions of this type. + * + * @return Forward rate coefficient. + */ + double evalFromStruct(const DetailedVibData& shared_data) const { + const double recipT13 = std::cbrt(shared_data.recipT); + + return m_scaling * m_A * std::exp( + m_b * shared_data.logT + + m_B + + m_C * recipT13 + + m_D * recipT13 * recipT13 + ); + } + + //! Evaluate derivative of reaction rate with respect to temperature, + //! divided by the reaction rate. + /*! + * This returns: + * + * @f[ + * \frac{1}{k_f} \frac{d k_f}{dT} + * = + * \frac{b}{T} + * - \frac{C}{3} T^{-4/3} + * - \frac{2D}{3} T^{-5/3} + * @f] + * + * This derivative is consistent with `evalFromStruct` and does not modify + * the reaction rate expression itself. + * + * @param shared_data Data shared by all reactions of this type. + */ + double ddTScaledFromStruct(const DetailedVibData& shared_data) const; + +private: + //! Dimensionless constant in the exponential. + double m_B = 0.0; + + //! Constant multiplying T^(-1/3), interpreted in K^(1/3). + double m_C = 0.0; + + //! Constant multiplying T^(-2/3), interpreted in K^(2/3). + double m_D = 0.0; + + //! Dimensionless scaling factor. + double m_scaling = 1.0; + + //! YAML key for `B`. + string m_B_str = "B"; + + //! YAML key for `C`. + string m_C_str = "C"; + + //! YAML key for `D`. + string m_D_str = "D"; + + //! YAML key for `scaling`. + string m_scaling_str = "scaling"; +}; + +} + +#endif \ No newline at end of file diff --git a/src/kinetics/DetailedVVVTRate.cpp b/src/kinetics/DetailedVVVTRate.cpp new file mode 100644 index 00000000000..90ea8aef863 --- /dev/null +++ b/src/kinetics/DetailedVVVTRate.cpp @@ -0,0 +1,194 @@ +//! @file DetailedVVVTRate.cpp + +// This file is part of Cantera. See License.txt in the top-level directory or +// at https://cantera.org/license.txt for license and copyright information. + +#include "cantera/kinetics/DetailedVVVTRate.h" +#include "cantera/kinetics/Reaction.h" +#include "cantera/thermo/ThermoPhase.h" + +#include +#include + +namespace Cantera +{ + +bool DetailedVibData::update(const ThermoPhase& phase, const Kinetics& kin) +{ + double T = phase.temperature(); + + if (T == temperature) { + return false; + } + + ReactionData::update(T); + return true; +} + + +DetailedVVVTRate::DetailedVVVTRate() +{ +} + + +DetailedVVVTRate::DetailedVVVTRate(double A, double B, double C, double D, + double b, double scaling) + : ArrheniusBase(A, b, 0.0) + , m_B(B) + , m_C(C) + , m_D(D) + , m_scaling(scaling) +{ +} + + +DetailedVVVTRate::DetailedVVVTRate(const AnyMap& node, + const UnitStack& rate_units) + : DetailedVVVTRate() +{ + setParameters(node, rate_units); +} + + +void DetailedVVVTRate::setParameters(const AnyMap& node, + const UnitStack& rate_units) +{ + // Let ArrheniusBase handle: + // + // - storage of the input AnyMap + // - negative-A handling + // - conversion of A using the normal unit system + // - reading of b + // + // The Ea parameter, if present, is parsed by ArrheniusBase but is not used + // by DetailedVVVTRate. The rate expression implemented here uses B, C, D + // instead of an Arrhenius activation energy. + ArrheniusBase::setParameters(node, rate_units); + + if (!node.hasKey("rate-constant")) { + return; + } + + const auto& rate = node["rate-constant"]; + + if (!rate.is()) { + throw InputFileError("DetailedVVVTRate::setParameters", node, + "The 'rate-constant' field for a 'detailed-vv-vt' reaction " + "must be a mapping containing A, b, B, C, D, and optionally " + "scaling."); + } + + const auto& rate_map = rate.as(); + + // B is dimensionless. + if (rate_map.hasKey(m_B_str)) { + m_B = rate_map[m_B_str].asDouble(); + } + + // C is interpreted as K^(1/3), assuming T is in K. + // It is intentionally read as a raw number and is not converted by the + // Cantera unit system. + if (rate_map.hasKey(m_C_str)) { + m_C = rate_map[m_C_str].asDouble(); + } + + // D is interpreted as K^(2/3), assuming T is in K. + // It is intentionally read as a raw number and is not converted by the + // Cantera unit system. + if (rate_map.hasKey(m_D_str)) { + m_D = rate_map[m_D_str].asDouble(); + } + + // scaling is dimensionless. + if (rate_map.hasKey(m_scaling_str)) { + m_scaling = rate_map[m_scaling_str].asDouble(); + } +} + + +void DetailedVVVTRate::getParameters(AnyMap& node) const +{ + if (!valid()) { + return; + } + + if (allowNegativePreExponentialFactor()) { + node["negative-A"] = true; + } + + AnyMap rateNode; + + // Store A using the same convention as ArrheniusBase::getRateParameters. + // When the reaction has been associated with a Kinetics object, Cantera + // knows the conversion units for the leading rate coefficient. + if (conversionUnits().factor() != 0.0) { + rateNode[m_A_str].setQuantity(m_A, conversionUnits()); + } else { + // This case can occur when the rate was created outside the context of + // a Kinetics object and therefore the reaction-order-dependent units + // are not known. + rateNode[m_A_str] = m_A; + rateNode["__unconvertible__"] = true; + } + + // b is dimensionless. + rateNode[m_b_str] = m_b; + + // Custom DetailedVVVTRate parameters. + // + // B is dimensionless. + // C is interpreted as K^(1/3), assuming T is in K. + // D is interpreted as K^(2/3), assuming T is in K. + // scaling is dimensionless. + // + // These values are deliberately serialized as raw floating-point numbers. + rateNode[m_B_str] = m_B; + rateNode[m_C_str] = m_C; + rateNode[m_D_str] = m_D; + rateNode[m_scaling_str] = m_scaling; + + rateNode.setFlowStyle(); + + node["rate-constant"] = std::move(rateNode); +} + + +double DetailedVVVTRate::ddTScaledFromStruct( + const DetailedVibData& shared_data) const +{ + const double invT = shared_data.recipT; + const double invT13 = std::cbrt(invT); + + // For: + // + // k = scaling * A * exp(b*log(T) + B + C*T^(-1/3) + D*T^(-2/3)) + // + // the logarithmic derivative is: + // + // (1/k) * dk/dT + // = b/T + // - (C/3) * T^(-4/3) + // - (2D/3) * T^(-5/3) + // + // Using invT = 1/T and invT13 = T^(-1/3): + // + // T^(-4/3) = invT13 * invT + // T^(-5/3) = invT13 * invT13 * invT + return m_b * invT + - (m_C / 3.0) * invT13 * invT + - (2.0 * m_D / 3.0) * invT13 * invT13 * invT; +} + + +void DetailedVVVTRate::setContext(const Reaction& rxn, const Kinetics& kin) +{ + // DetailedVVVTRate is intended for non-equilibrium plasma kinetics. + // The reverse rate cannot be calculated from conventional thermochemistry + // without an additional non-equilibrium model. + if (rxn.reversible) { + throw InputFileError("DetailedVVVTRate::setContext", rxn.input, + "DetailedVVVTRate does not support reversible reactions."); + } +} + +} \ No newline at end of file diff --git a/src/kinetics/ReactionRateFactory.cpp b/src/kinetics/ReactionRateFactory.cpp index 1019fe4f4ed..b8eb268bf70 100644 --- a/src/kinetics/ReactionRateFactory.cpp +++ b/src/kinetics/ReactionRateFactory.cpp @@ -18,6 +18,7 @@ #include "cantera/kinetics/InterfaceRate.h" #include "cantera/kinetics/PlogRate.h" #include "cantera/kinetics/TwoTempPlasmaRate.h" +#include "cantera/kinetics/DetailedVVVTRate.h" namespace Cantera { @@ -40,6 +41,11 @@ ReactionRateFactory::ReactionRateFactory() return new TwoTempPlasmaRate(node, rate_units); }); + // DetailedVVVTRate evaluator + reg("detailed-vv-vt", [](const AnyMap& node, const UnitStack& rate_units) { + return new DetailedVVVTRate(node, rate_units); + }); + // ElectronCollisionPlasmaRate evaluator reg("electron-collision-plasma", [](const AnyMap& node, const UnitStack& rate_units) { return new ElectronCollisionPlasmaRate(node, rate_units); diff --git a/test/kinetics/kineticsFromScratch.cpp b/test/kinetics/kineticsFromScratch.cpp index b4d120099db..bcd4e1d92f9 100644 --- a/test/kinetics/kineticsFromScratch.cpp +++ b/test/kinetics/kineticsFromScratch.cpp @@ -14,6 +14,7 @@ #include "cantera/kinetics/InterfaceRate.h" #include "cantera/kinetics/PlogRate.h" #include "cantera/kinetics/TwoTempPlasmaRate.h" +#include "cantera/kinetics/DetailedVVVTRate.h" #include "cantera/base/Array.h" #include "cantera/base/stringUtils.h" diff --git a/test/kinetics/kineticsFromYaml.cpp b/test/kinetics/kineticsFromYaml.cpp index fbe699b3889..c34a8633f8a 100644 --- a/test/kinetics/kineticsFromYaml.cpp +++ b/test/kinetics/kineticsFromYaml.cpp @@ -13,6 +13,7 @@ #include "cantera/kinetics/InterfaceRate.h" #include "cantera/kinetics/PlogRate.h" #include "cantera/kinetics/TwoTempPlasmaRate.h" +#include "cantera/kinetics/DetailedVVVTRate.h" #include "cantera/thermo/SurfPhase.h" #include "cantera/thermo/ThermoFactory.h" #include "cantera/base/Array.h" From 1527744b63f5503afefcb941d4f21c64ba07f1fc Mon Sep 17 00:00:00 2001 From: Gaetanosaure Date: Thu, 14 May 2026 15:13:15 +0200 Subject: [PATCH 17/29] [python] Adding the access to reading the electronMobility of the gas --- interfaces/cython/cantera/thermo.pxd | 1 + interfaces/cython/cantera/thermo.pyx | 11 +++++++++++ 2 files changed, 12 insertions(+) diff --git a/interfaces/cython/cantera/thermo.pxd b/interfaces/cython/cantera/thermo.pxd index c80bed26316..53eb23113e2 100644 --- a/interfaces/cython/cantera/thermo.pxd +++ b/interfaces/cython/cantera/thermo.pxd @@ -213,6 +213,7 @@ cdef extern from "cantera/thermo/PlasmaPhase.h": double electricField() void updateElectronEnergyDistribution() except +translate_exception double elasticPowerLoss() except +translate_exception + double electronMobility() cdef extern from "cantera/cython/thermo_utils.h": diff --git a/interfaces/cython/cantera/thermo.pyx b/interfaces/cython/cantera/thermo.pyx index 021f80ca359..b4529abf97f 100644 --- a/interfaces/cython/cantera/thermo.pyx +++ b/interfaces/cython/cantera/thermo.pyx @@ -1991,6 +1991,17 @@ cdef class ThermoPhase(_SolutionBase): raise ThermoModelMethodError(self.thermo_model) return self.plasma.elasticPowerLoss() + property electron_mobility: + """ + Electron mobility (m^2/(V.s)) + + .. versionadded:: 4.0 + """ + def __get__(self): + if not self._enable_plasma: + raise ThermoModelMethodError(self.thermo_model) + return self.plasma.electronMobility() + cdef class InterfacePhase(ThermoPhase): """ A class representing a surface, edge phase """ def __cinit__(self, *args, **kwargs): From b27674b5a2688657be2865e23d0e644b701f3173 Mon Sep 17 00:00:00 2001 From: Gaetanosaure Date: Fri, 15 May 2026 11:53:44 +0200 Subject: [PATCH 18/29] [thermo] Allowing collapsed inelastic channels entries in the YAML --- src/kinetics/ElectronCollisionPlasmaRate.cpp | 24 +++++++++++++++----- src/thermo/PlasmaPhase.cpp | 12 ++++++++++ src/zeroD/ConstPressureReactor.cpp | 14 +++++++++++- 3 files changed, 43 insertions(+), 7 deletions(-) diff --git a/src/kinetics/ElectronCollisionPlasmaRate.cpp b/src/kinetics/ElectronCollisionPlasmaRate.cpp index 6223ba64aac..a21f4fd5e0f 100644 --- a/src/kinetics/ElectronCollisionPlasmaRate.cpp +++ b/src/kinetics/ElectronCollisionPlasmaRate.cpp @@ -260,22 +260,36 @@ void ElectronCollisionPlasmaRate::setContext(const Reaction& rxn, const Kinetics void ElectronCollisionPlasmaRate::applyCollisionData(const AnyMap& node) { - if (node.hasKey("kind")) { string collisionKind = node["kind"].asString(); if (!m_kind.empty() && m_kind != collisionKind) { string collisionName = node.hasKey("name") ? node["name"].asString() : m_collisionName; - throw InputFileError("applyCollisionData", node, - "Electron collision '{}' has kind '{}', but the reaction was inferred as '{}'.", + bool allowCollapsedInelastic = + (m_kind == "effective" || m_kind == "elastic") && + collisionKind == "excitation"; + + if (!allowCollapsedInelastic) { + throw InputFileError("applyCollisionData", node, + "Electron collision '{}' has kind '{}', but the reaction was inferred as '{}'.", + collisionName, collisionKind, m_kind); + } + + warn_user("ElectronCollisionPlasmaRate::applyCollisionData", + "Electron collision '{}' has kind '{}', but the reaction was inferred " + "as '{}'. Treating this as an intentionally collapsed inelastic " + "electron-collision channel. No species source term will be generated " + "for the unresolved product.", collisionName, collisionKind, m_kind); } + // Important: keep the collision classified using the data-file kind. + // For collapsed channels such as N2(rot), this ensures that the + // Boltzmann solver treats the channel as inelastic. m_kind = collisionKind; } - if (node.hasKey("target")) { m_target = node["target"].asString(); } @@ -300,8 +314,6 @@ void ElectronCollisionPlasmaRate::applyCollisionData(const AnyMap& node) validateCollisionData(node); m_hasCrossSectionData = true; - - } void ElectronCollisionPlasmaRate::validateCollisionData(const AnyMap& node) const diff --git a/src/thermo/PlasmaPhase.cpp b/src/thermo/PlasmaPhase.cpp index 58143707a15..9398f03efcc 100644 --- a/src/thermo/PlasmaPhase.cpp +++ b/src/thermo/PlasmaPhase.cpp @@ -800,6 +800,16 @@ void PlasmaPhase::addCollision(shared_ptr collision) compatibleKind = true; } + // Allow collapsed inelastic channels. + // Example: + // reaction: Electron + N2 => Electron + N2 + // collision: kind: excitation, product: N2(rot) + if (!compatibleKind && + (kindFromReaction == "effective" || kindFromReaction == "elastic") && + kindFromCollision == "excitation") { + compatibleKind = true; + } + if (!compatibleKind) { throw CanteraError("PlasmaPhase::addCollision", "Electron collision '{}' has kind '{}', but reaction '{}' is inferred as '{}'.", @@ -1210,6 +1220,8 @@ double PlasmaPhase::intrinsicHeating() const double qJ = jouleHeatingPower(); checkFinite(qJ); + writelog("Calling PlasmaPhase::intrinsicHeating(): qJ = {}\n", qJ); + return qJ; } diff --git a/src/zeroD/ConstPressureReactor.cpp b/src/zeroD/ConstPressureReactor.cpp index 7972806be55..965dba12f58 100644 --- a/src/zeroD/ConstPressureReactor.cpp +++ b/src/zeroD/ConstPressureReactor.cpp @@ -87,8 +87,20 @@ void ConstPressureReactor::eval(double time, span LHS, span RHS) // external heat transfer double dHdt = m_Qdot; + // if (m_energy) { + // dHdt += m_thermo->intrinsicHeating() * m_vol; + // } + if (m_energy) { - dHdt += m_thermo->intrinsicHeating() * m_vol; + double qint = m_thermo->intrinsicHeating(); + + static int count = 0; + if (count++ < 20) { + writelog("intrinsicHeating = {} W/m3, V = {}, contribution = {} W\n", + qint, m_vol, qint * m_vol); + } + + dHdt += qint * m_vol; } // add terms for outlets From 55eba45ef1d16132f205e6fbc257ef60f6b37eab Mon Sep 17 00:00:00 2001 From: Gaetanosaure Date: Fri, 15 May 2026 12:39:51 +0200 Subject: [PATCH 19/29] cleaning the code from debug messages + adding an optional parameter to two-temperature-plasma reactions to extend the formula when needed --- include/cantera/kinetics/TwoTempPlasmaRate.h | 18 +++++++++++++----- src/kinetics/TwoTempPlasmaRate.cpp | 19 ++++++++++++++++++- src/thermo/PlasmaPhase.cpp | 2 -- src/zeroD/ConstPressureReactor.cpp | 14 +------------- 4 files changed, 32 insertions(+), 21 deletions(-) diff --git a/include/cantera/kinetics/TwoTempPlasmaRate.h b/include/cantera/kinetics/TwoTempPlasmaRate.h index 7f17e9de723..20c8ec779da 100644 --- a/include/cantera/kinetics/TwoTempPlasmaRate.h +++ b/include/cantera/kinetics/TwoTempPlasmaRate.h @@ -48,12 +48,14 @@ struct TwoTempPlasmaData : public ReactionData * activation energy for electron is included. * * @f[ - * k_f = A T_e^b \exp (-E_{a,g}/RT) \exp (E_{a,e} (T_e - T)/(R T T_e)) + * k_f = A T^{b_g} T_e^b * exp(-E_{a,g}/RT) * exp(E_{a,e}(T_e - T)/(R T T_e)) * @f] * * where @f$ T_e @f$ is the electron temperature, @f$ E_{a,g} @f$ is the activation * energy for gas, and @f$ E_{a,e} @f$ is the activation energy for electron, see * Kossyi, et al. @cite kossyi1992. + * The optional gas temperature exponent b_g defaults to zero, which stricly corresponds to @cite kossyi1992. + * If b_g is non-zero, a generalisation is used. * * @ingroup arrheniusGroup */ @@ -69,7 +71,10 @@ class TwoTempPlasmaRate : public ArrheniusBase * @param b Temperature exponent (non-dimensional) * @param Ea Activation energy in energy units [J/kmol] * @param EE Activation electron energy in energy units [J/kmol] + * @param bg Gas temperature exponent (non-dimensional). If not specified, defaults to 0. */ + TwoTempPlasmaRate(double A, double b, double Ea=0.0, double EE=0.0, double bg=0.0); + TwoTempPlasmaRate(double A, double b, double Ea=0.0, double EE=0.0); TwoTempPlasmaRate(const AnyMap& node, const UnitStack& rate_units={}); @@ -90,12 +95,12 @@ class TwoTempPlasmaRate : public ArrheniusBase */ double evalFromStruct(const TwoTempPlasmaData& shared_data) const { // m_E4_R is the electron activation (in temperature units) - return m_A * std::exp(m_b * shared_data.logTe - - m_Ea_R * shared_data.recipT + - m_E4_R * (shared_data.electronTemp - shared_data.temperature) - * shared_data.recipTe * shared_data.recipT); + return m_A * std::exp(m_bg * shared_data.logT + m_b * shared_data.logTe + - m_Ea_R * shared_data.recipT + m_E4_R * (shared_data.electronTemp - shared_data.temperature) + * shared_data.recipTe * shared_data.recipT); } + //! Evaluate derivative of reaction rate with respect to temperature //! divided by reaction rate /*! @@ -109,6 +114,9 @@ class TwoTempPlasmaRate : public ArrheniusBase double activationElectronEnergy() const { return m_E4_R * GasConstant; } + +protected: + double m_bg = 0.0; //!< Gas temperature exponent }; } diff --git a/src/kinetics/TwoTempPlasmaRate.cpp b/src/kinetics/TwoTempPlasmaRate.cpp index ce3f4694969..bd7bbe73dbd 100644 --- a/src/kinetics/TwoTempPlasmaRate.cpp +++ b/src/kinetics/TwoTempPlasmaRate.cpp @@ -57,19 +57,36 @@ TwoTempPlasmaRate::TwoTempPlasmaRate(double A, double b, double Ea, double EE) m_Ea_str = "Ea-gas"; m_E4_str = "Ea-electron"; m_E4_R = EE / GasConstant; + m_bg = 0.0; +} + +TwoTempPlasmaRate::TwoTempPlasmaRate(double A, double b, double Ea, double EE, double bg) + : ArrheniusBase(A, b, Ea) +{ + m_Ea_str = "Ea-gas"; + m_E4_str = "Ea-electron"; + m_E4_R = EE / GasConstant; + m_bg = bg; } TwoTempPlasmaRate::TwoTempPlasmaRate(const AnyMap& node, const UnitStack& rate_units) : TwoTempPlasmaRate() { setParameters(node, rate_units); + + if (node.hasKey("b-gas")) { + m_bg = node["b-gas"].asDouble(); + } else if (node.hasKey("b_gas")) { + m_bg = node["b_gas"].asDouble(); + } } double TwoTempPlasmaRate::ddTScaledFromStruct(const TwoTempPlasmaData& shared_data) const { warn_user("TwoTempPlasmaRate::ddTScaledFromStruct", "Temperature derivative does not consider changes of electron temperature."); - return (m_Ea_R - m_E4_R) * shared_data.recipT * shared_data.recipT; + return m_bg * shared_data.recipT + + (m_Ea_R - m_E4_R) * shared_data.recipT * shared_data.recipT; } void TwoTempPlasmaRate::setContext(const Reaction& rxn, const Kinetics& kin) diff --git a/src/thermo/PlasmaPhase.cpp b/src/thermo/PlasmaPhase.cpp index 9398f03efcc..5aa549015c6 100644 --- a/src/thermo/PlasmaPhase.cpp +++ b/src/thermo/PlasmaPhase.cpp @@ -1220,8 +1220,6 @@ double PlasmaPhase::intrinsicHeating() const double qJ = jouleHeatingPower(); checkFinite(qJ); - writelog("Calling PlasmaPhase::intrinsicHeating(): qJ = {}\n", qJ); - return qJ; } diff --git a/src/zeroD/ConstPressureReactor.cpp b/src/zeroD/ConstPressureReactor.cpp index 965dba12f58..7972806be55 100644 --- a/src/zeroD/ConstPressureReactor.cpp +++ b/src/zeroD/ConstPressureReactor.cpp @@ -87,20 +87,8 @@ void ConstPressureReactor::eval(double time, span LHS, span RHS) // external heat transfer double dHdt = m_Qdot; - // if (m_energy) { - // dHdt += m_thermo->intrinsicHeating() * m_vol; - // } - if (m_energy) { - double qint = m_thermo->intrinsicHeating(); - - static int count = 0; - if (count++ < 20) { - writelog("intrinsicHeating = {} W/m3, V = {}, contribution = {} W\n", - qint, m_vol, qint * m_vol); - } - - dHdt += qint * m_vol; + dHdt += m_thermo->intrinsicHeating() * m_vol; } // add terms for outlets From 49cbc703ab1db2b3d6e4c6a0a6053ebb7a3a2628 Mon Sep 17 00:00:00 2001 From: Gaetanosaure Date: Tue, 19 May 2026 00:43:14 +0200 Subject: [PATCH 20/29] [Kinetics] Enlarge DetailedVVVTRate into VibrationalRelaxationRate to be able to handle the vibrational energy equation model with V-T relaxation models constant, Castela and Starikovskiy. This works by declaring fictitious vibrational species in the YAML that act as energy reservoir, hence curbing the need for a new PlasmaReactor class. --- include/cantera/kinetics/DetailedVVVTRate.h | 208 ----- include/cantera/kinetics/TwoTempPlasmaRate.h | 2 +- .../kinetics/VibrationalRelaxationRate.h | 253 ++++++ src/kinetics/DetailedVVVTRate.cpp | 194 ---- src/kinetics/ReactionRateFactory.cpp | 8 +- src/kinetics/TwoTempPlasmaRate.cpp | 4 +- src/kinetics/VibrationalRelaxationRate.cpp | 835 ++++++++++++++++++ test/kinetics/kineticsFromScratch.cpp | 2 +- test/kinetics/kineticsFromYaml.cpp | 2 +- 9 files changed, 1096 insertions(+), 412 deletions(-) delete mode 100644 include/cantera/kinetics/DetailedVVVTRate.h create mode 100644 include/cantera/kinetics/VibrationalRelaxationRate.h delete mode 100644 src/kinetics/DetailedVVVTRate.cpp create mode 100644 src/kinetics/VibrationalRelaxationRate.cpp diff --git a/include/cantera/kinetics/DetailedVVVTRate.h b/include/cantera/kinetics/DetailedVVVTRate.h deleted file mode 100644 index 43cba7d2e20..00000000000 --- a/include/cantera/kinetics/DetailedVVVTRate.h +++ /dev/null @@ -1,208 +0,0 @@ -//! @file DetailedVVVTRate.h -//! Header for plasma reaction rates corresponding to detailed vibration handling. -//! -//! This file is part of Cantera. See License.txt in the top-level directory or -//! at https://cantera.org/license.txt for license and copyright information. - -#ifndef CT_DETAILEDVVVTRATE_H -#define CT_DETAILEDVVVTRATE_H - -#include "Arrhenius.h" - -#include - -namespace Cantera -{ - -//! Data container holding shared data specific to DetailedVVVTRate. -struct DetailedVibData : public ReactionData -{ - //! Update cached temperature-dependent data. - /*! - * @param phase Thermodynamic phase used to retrieve the gas temperature. - * @param kin Kinetics object. Not used here, but required by the - * `ReactionData` interface. - * - * @return `true` if the temperature has changed and rates need to be - * recomputed; `false` otherwise. - */ - bool update(const ThermoPhase& phase, const Kinetics& kin) override; - - using ReactionData::update; -}; - - -//! Detailed VV-VT reaction rate. -/*! - * This reaction rate implements the following temperature-dependent - * expression: - * - * @f[ - * k_f = - * scaling \, A \, - * \exp \left( - * b \ln T - * + B - * + C T^{-1/3} - * + D T^{-2/3} - * \right) - * @f] - * - * where `T` is the gas temperature in K. - * - * Unit conventions: - * - * - `A` uses the standard Cantera rate coefficient units. Its units depend - * on the reaction order and dimensionality, and are converted by - * `ArrheniusBase`. - * - `b` is dimensionless. - * - `B` is dimensionless. - * - `C` is interpreted as being expressed in K^(1/3), assuming `T` is in K. - * - `D` is interpreted as being expressed in K^(2/3), assuming `T` is in K. - * - `scaling` is dimensionless. - * - * Important: `B`, `C`, `D`, and `scaling` are read as raw floating-point - * numbers from YAML. They are not converted by Cantera's unit system. - * - * @ingroup arrheniusGroup - */ -class DetailedVVVTRate : public ArrheniusBase -{ -public: - //! Default constructor. - DetailedVVVTRate(); - - //! Constructor. - /*! - * @param A Pre-exponential factor. The unit system is (kmol, m, s); - * actual units depend on the reaction order and the - * dimensionality, surface or bulk. - * @param B Dimensionless constant in the exponential. - * @param C Constant multiplying T^(-1/3), interpreted in K^(1/3). - * @param D Constant multiplying T^(-2/3), interpreted in K^(2/3). - * @param b Dimensionless temperature exponent. - * @param scaling Dimensionless scaling factor from the harmonic oscillator - * theory. - */ - DetailedVVVTRate(double A, double B, double C, double D, - double b, double scaling = 1.0); - - //! Constructor based on AnyMap content. - /*! - * Expected YAML form: - * - * @code{.yaml} - * type: detailed-vv-vt - * rate-constant: - * A: ... - * b: ... - * B: ... - * C: ... - * D: ... - * scaling: ... - * @endcode - * - * The `scaling` entry is optional and defaults to 1.0. - */ - explicit DetailedVVVTRate(const AnyMap& node, - const UnitStack& rate_units = {}); - - //! Set rate parameters from an AnyMap. - /*! - * `A` and `b` are handled by `ArrheniusBase`. - * `B`, `C`, `D`, and `scaling` are handled by this class. - */ - void setParameters(const AnyMap& node, - const UnitStack& rate_units) override; - - //! Get rate parameters for YAML serialization. - /*! - * This method is needed so that Cantera can reconstruct or export the - * full custom rate, including `B`, `C`, `D`, and `scaling`. - */ - void getParameters(AnyMap& node) const override; - - //! Create a rate evaluator for this reaction rate type. - unique_ptr newMultiRate() const override { - return make_unique>(); - } - - //! String identifying this reaction rate specialization. - const string type() const override { - return "detailed-vv-vt"; - } - - //! Set context of reaction rate evaluation. - /*! - * This rate is intended for irreversible non-equilibrium plasma reactions. - * Reversible reactions are rejected because the reverse rate cannot be - * obtained from conventional thermochemistry for this model. - */ - void setContext(const Reaction& rxn, const Kinetics& kin) override; - - //! Evaluate reaction rate. - /*! - * @param shared_data Data shared by all reactions of this type. - * - * @return Forward rate coefficient. - */ - double evalFromStruct(const DetailedVibData& shared_data) const { - const double recipT13 = std::cbrt(shared_data.recipT); - - return m_scaling * m_A * std::exp( - m_b * shared_data.logT - + m_B - + m_C * recipT13 - + m_D * recipT13 * recipT13 - ); - } - - //! Evaluate derivative of reaction rate with respect to temperature, - //! divided by the reaction rate. - /*! - * This returns: - * - * @f[ - * \frac{1}{k_f} \frac{d k_f}{dT} - * = - * \frac{b}{T} - * - \frac{C}{3} T^{-4/3} - * - \frac{2D}{3} T^{-5/3} - * @f] - * - * This derivative is consistent with `evalFromStruct` and does not modify - * the reaction rate expression itself. - * - * @param shared_data Data shared by all reactions of this type. - */ - double ddTScaledFromStruct(const DetailedVibData& shared_data) const; - -private: - //! Dimensionless constant in the exponential. - double m_B = 0.0; - - //! Constant multiplying T^(-1/3), interpreted in K^(1/3). - double m_C = 0.0; - - //! Constant multiplying T^(-2/3), interpreted in K^(2/3). - double m_D = 0.0; - - //! Dimensionless scaling factor. - double m_scaling = 1.0; - - //! YAML key for `B`. - string m_B_str = "B"; - - //! YAML key for `C`. - string m_C_str = "C"; - - //! YAML key for `D`. - string m_D_str = "D"; - - //! YAML key for `scaling`. - string m_scaling_str = "scaling"; -}; - -} - -#endif \ No newline at end of file diff --git a/include/cantera/kinetics/TwoTempPlasmaRate.h b/include/cantera/kinetics/TwoTempPlasmaRate.h index 20c8ec779da..b91531c01bc 100644 --- a/include/cantera/kinetics/TwoTempPlasmaRate.h +++ b/include/cantera/kinetics/TwoTempPlasmaRate.h @@ -73,7 +73,7 @@ class TwoTempPlasmaRate : public ArrheniusBase * @param EE Activation electron energy in energy units [J/kmol] * @param bg Gas temperature exponent (non-dimensional). If not specified, defaults to 0. */ - TwoTempPlasmaRate(double A, double b, double Ea=0.0, double EE=0.0, double bg=0.0); + TwoTempPlasmaRate(double A, double b, double Ea, double EE, double bg); TwoTempPlasmaRate(double A, double b, double Ea=0.0, double EE=0.0); diff --git a/include/cantera/kinetics/VibrationalRelaxationRate.h b/include/cantera/kinetics/VibrationalRelaxationRate.h new file mode 100644 index 00000000000..59592f6ec27 --- /dev/null +++ b/include/cantera/kinetics/VibrationalRelaxationRate.h @@ -0,0 +1,253 @@ +//! @file VibrationalRelaxationRate.h +//! Header for vibrational relaxation reaction rates in plasma kinetics. +//! +//! This file is part of Cantera. See License.txt in the top-level directory or +//! at https://cantera.org/license.txt for license and copyright information. + +#ifndef CT_VIBRATIONALRELAXATIONRATE_H +#define CT_VIBRATIONALRELAXATIONRATE_H + +#include "Arrhenius.h" + +#include + +namespace Cantera +{ + +//! Shared temperature data for vibrational relaxation rates. +struct DetailedVibData : public ReactionData +{ + //! Update cached temperature-dependent data. + /*! + * @param phase Thermodynamic phase used to retrieve the gas temperature. + * @param kin Kinetics object. Not used here, but required by the + * ReactionData interface. + * + * @return `true` if the temperature has changed and rates need to be + * recomputed; `false` otherwise. + */ + bool update(const ThermoPhase& phase, const Kinetics& kin) override; + + using ReactionData::update; +}; + + +//! Vibrational relaxation reaction rate. +/*! + * This class provides a common implementation for several vibrational + * relaxation models: + * + * - `constant` + * - `detailed-vv-vt` + * - `starikovskiy` + * - `castela` + * + * Internally, all models are mapped to the following generic expression: + * + * @f[ + * k_f = + * scaling \, A \, + * \exp \left( + * b \ln T + * + B + * + C T^{-1/3} + * + D T^{-m} + * + E T^{-z} + * \right) + * @f] + * + * where `T` is the gas temperature in K. + * + * The YAML reaction type is: + * + * @code{.yaml} + * type: vibrational-relaxation + * @endcode + * + * The selected physical model is specified separately using: + * + * @code{.yaml} + * vibration_model: constant + * vibration_model: detailed-vv-vt + * vibration_model: starikovskiy + * vibration_model: castela + * @endcode + * + * Unit conventions: + * + * - `A` uses standard Cantera rate coefficient units. Its units depend on the + * reaction order and are converted by `ArrheniusBase`. + * - `b`, `B`, `m`, `z`, and `scaling` are dimensionless. + * - `C` is interpreted as K^(1/3), assuming `T` is in K. + * - `D` is interpreted as K^m, assuming `T` is in K. + * - `E` is interpreted as K^z, assuming `T` is in K. + * + * The coefficients `B`, `C`, `D`, `E`, `m`, `z`, and `scaling` are read as + * raw floating-point values. They are not converted by Cantera's unit system. + * + * @ingroup arrheniusGroup + */ +class VibrationalRelaxationRate : public ArrheniusBase +{ +public: + //! Default constructor. + VibrationalRelaxationRate(); + + //! Constructor using the internal generic representation. + /*! + * @param A Pre-exponential factor. + * @param B Dimensionless constant in the exponential. + * @param C Coefficient multiplying T^(-1/3). + * @param D Coefficient multiplying T^(-m). + * @param b Dimensionless temperature exponent. + * @param scaling Dimensionless scaling factor. + * @param m Temperature exponent used by the D term. + * @param E Coefficient multiplying T^(-z). + * @param z Temperature exponent used by the E term. + */ + VibrationalRelaxationRate(double A, double B, double C, double D, + double b, double scaling = 1.0, + double m = 2.0 / 3.0, + double E = 0.0, double z = 1.0); + + //! Constructor based on AnyMap content. + explicit VibrationalRelaxationRate(const AnyMap& node, + const UnitStack& rate_units = {}); + + //! Set rate parameters from an AnyMap. + void setParameters(const AnyMap& node, + const UnitStack& rate_units) override; + + //! Get rate parameters for YAML serialization. + void getParameters(AnyMap& node) const override; + + //! Create a rate evaluator for this reaction rate type. + unique_ptr newMultiRate() const override { + return make_unique>(); + } + + //! String identifying this reaction rate type. + const string type() const override { + return "vibrational-relaxation"; + } + + //! Set context of reaction rate evaluation. + /*! + * Vibrational relaxation rates are intended for irreversible + * non-equilibrium plasma reactions. Reversible reactions are rejected + * because the reverse rate cannot be obtained from conventional + * thermochemistry for these models. + */ + void setContext(const Reaction& rxn, const Kinetics& kin) override; + + //! Evaluate the forward rate coefficient. + double evalFromStruct(const DetailedVibData& shared_data) const { + const double invT = shared_data.recipT; + const double invT13 = std::cbrt(invT); + + return m_scaling * m_A * std::exp( + m_b * shared_data.logT + + m_B + + m_C * invT13 + + m_D * std::pow(invT, m_m) + + m_E * std::pow(invT, m_z) + ); + } + + //! Evaluate the scaled temperature derivative. + /*! + * This returns: + * + * @f[ + * \frac{1}{k_f} \frac{d k_f}{dT} + * = + * \frac{d \ln k_f}{dT} + * @f] + * + * For the internal generic expression, this is: + * + * @f[ + * \frac{b}{T} + * - \frac{C}{3} T^{-4/3} + * - m D T^{-m-1} + * - z E T^{-z-1} + * @f] + */ + double ddTScaledFromStruct(const DetailedVibData& shared_data) const; + +private: + //! Dimensionless constant in the exponential. + double m_B = 0.0; + + //! Coefficient multiplying T^(-1/3). + double m_C = 0.0; + + //! Coefficient multiplying T^(-m). + double m_D = 0.0; + + //! Dimensionless scaling factor. + double m_scaling = 1.0; + + //! Temperature exponent used by the D term. + double m_m = 2.0 / 3.0; + + //! Coefficient multiplying T^(-z). + double m_E = 0.0; + + //! Temperature exponent used by the E term. + double m_z = 1.0; + + //! Castela coefficient a. + double m_castela_a = 0.0; + + //! Castela coefficient b. + double m_castela_b = 0.0; + + //! Castela reference pressure. + double m_referencePressure = OneAtm; + + //! Selected vibrational relaxation model. + /*! + * Accepted values: + * + * - `constant` + * - `detailed-vv-vt` + * - `starikovskiy` + * - `castela` + */ + string m_vibration_model = "detailed-vv-vt"; + + //! YAML key names. + string m_B_str = "B"; + string m_C_str = "C"; + string m_D_str = "D"; + string m_scaling_str = "scaling"; + string m_m_str = "m"; + string m_E_str = "E"; + string m_z_str = "z"; + string m_reference_pressure_str = "reference-pressure"; + + //! Configure the ArrheniusBase part from an already-converted internal A value. + /*! + * This is needed for models such as Castela, where the user-facing YAML does + * not contain a standard Arrhenius A coefficient. + */ + void configureBaseFromInternalA(const AnyMap& node, + const UnitStack& rate_units, + double A, double b); + + //! Configure the ArrheniusBase part from a YAML A value and an explicit b. + /*! + * This is needed for models such as constant and Starikovskiy, where the YAML + * does not contain the standard Arrhenius pair A / b. + */ + void configureBaseFromYamlA(const AnyMap& node, + const UnitStack& rate_units, + const AnyValue& A, + double b); + + }; + +} + +#endif \ No newline at end of file diff --git a/src/kinetics/DetailedVVVTRate.cpp b/src/kinetics/DetailedVVVTRate.cpp deleted file mode 100644 index 90ea8aef863..00000000000 --- a/src/kinetics/DetailedVVVTRate.cpp +++ /dev/null @@ -1,194 +0,0 @@ -//! @file DetailedVVVTRate.cpp - -// This file is part of Cantera. See License.txt in the top-level directory or -// at https://cantera.org/license.txt for license and copyright information. - -#include "cantera/kinetics/DetailedVVVTRate.h" -#include "cantera/kinetics/Reaction.h" -#include "cantera/thermo/ThermoPhase.h" - -#include -#include - -namespace Cantera -{ - -bool DetailedVibData::update(const ThermoPhase& phase, const Kinetics& kin) -{ - double T = phase.temperature(); - - if (T == temperature) { - return false; - } - - ReactionData::update(T); - return true; -} - - -DetailedVVVTRate::DetailedVVVTRate() -{ -} - - -DetailedVVVTRate::DetailedVVVTRate(double A, double B, double C, double D, - double b, double scaling) - : ArrheniusBase(A, b, 0.0) - , m_B(B) - , m_C(C) - , m_D(D) - , m_scaling(scaling) -{ -} - - -DetailedVVVTRate::DetailedVVVTRate(const AnyMap& node, - const UnitStack& rate_units) - : DetailedVVVTRate() -{ - setParameters(node, rate_units); -} - - -void DetailedVVVTRate::setParameters(const AnyMap& node, - const UnitStack& rate_units) -{ - // Let ArrheniusBase handle: - // - // - storage of the input AnyMap - // - negative-A handling - // - conversion of A using the normal unit system - // - reading of b - // - // The Ea parameter, if present, is parsed by ArrheniusBase but is not used - // by DetailedVVVTRate. The rate expression implemented here uses B, C, D - // instead of an Arrhenius activation energy. - ArrheniusBase::setParameters(node, rate_units); - - if (!node.hasKey("rate-constant")) { - return; - } - - const auto& rate = node["rate-constant"]; - - if (!rate.is()) { - throw InputFileError("DetailedVVVTRate::setParameters", node, - "The 'rate-constant' field for a 'detailed-vv-vt' reaction " - "must be a mapping containing A, b, B, C, D, and optionally " - "scaling."); - } - - const auto& rate_map = rate.as(); - - // B is dimensionless. - if (rate_map.hasKey(m_B_str)) { - m_B = rate_map[m_B_str].asDouble(); - } - - // C is interpreted as K^(1/3), assuming T is in K. - // It is intentionally read as a raw number and is not converted by the - // Cantera unit system. - if (rate_map.hasKey(m_C_str)) { - m_C = rate_map[m_C_str].asDouble(); - } - - // D is interpreted as K^(2/3), assuming T is in K. - // It is intentionally read as a raw number and is not converted by the - // Cantera unit system. - if (rate_map.hasKey(m_D_str)) { - m_D = rate_map[m_D_str].asDouble(); - } - - // scaling is dimensionless. - if (rate_map.hasKey(m_scaling_str)) { - m_scaling = rate_map[m_scaling_str].asDouble(); - } -} - - -void DetailedVVVTRate::getParameters(AnyMap& node) const -{ - if (!valid()) { - return; - } - - if (allowNegativePreExponentialFactor()) { - node["negative-A"] = true; - } - - AnyMap rateNode; - - // Store A using the same convention as ArrheniusBase::getRateParameters. - // When the reaction has been associated with a Kinetics object, Cantera - // knows the conversion units for the leading rate coefficient. - if (conversionUnits().factor() != 0.0) { - rateNode[m_A_str].setQuantity(m_A, conversionUnits()); - } else { - // This case can occur when the rate was created outside the context of - // a Kinetics object and therefore the reaction-order-dependent units - // are not known. - rateNode[m_A_str] = m_A; - rateNode["__unconvertible__"] = true; - } - - // b is dimensionless. - rateNode[m_b_str] = m_b; - - // Custom DetailedVVVTRate parameters. - // - // B is dimensionless. - // C is interpreted as K^(1/3), assuming T is in K. - // D is interpreted as K^(2/3), assuming T is in K. - // scaling is dimensionless. - // - // These values are deliberately serialized as raw floating-point numbers. - rateNode[m_B_str] = m_B; - rateNode[m_C_str] = m_C; - rateNode[m_D_str] = m_D; - rateNode[m_scaling_str] = m_scaling; - - rateNode.setFlowStyle(); - - node["rate-constant"] = std::move(rateNode); -} - - -double DetailedVVVTRate::ddTScaledFromStruct( - const DetailedVibData& shared_data) const -{ - const double invT = shared_data.recipT; - const double invT13 = std::cbrt(invT); - - // For: - // - // k = scaling * A * exp(b*log(T) + B + C*T^(-1/3) + D*T^(-2/3)) - // - // the logarithmic derivative is: - // - // (1/k) * dk/dT - // = b/T - // - (C/3) * T^(-4/3) - // - (2D/3) * T^(-5/3) - // - // Using invT = 1/T and invT13 = T^(-1/3): - // - // T^(-4/3) = invT13 * invT - // T^(-5/3) = invT13 * invT13 * invT - return m_b * invT - - (m_C / 3.0) * invT13 * invT - - (2.0 * m_D / 3.0) * invT13 * invT13 * invT; -} - - -void DetailedVVVTRate::setContext(const Reaction& rxn, const Kinetics& kin) -{ - // DetailedVVVTRate is intended for non-equilibrium plasma kinetics. - // The reverse rate cannot be calculated from conventional thermochemistry - // without an additional non-equilibrium model. - if (rxn.reversible) { - throw InputFileError("DetailedVVVTRate::setContext", rxn.input, - "DetailedVVVTRate does not support reversible reactions."); - } -} - -} \ No newline at end of file diff --git a/src/kinetics/ReactionRateFactory.cpp b/src/kinetics/ReactionRateFactory.cpp index b8eb268bf70..cb079dc0ee1 100644 --- a/src/kinetics/ReactionRateFactory.cpp +++ b/src/kinetics/ReactionRateFactory.cpp @@ -18,7 +18,7 @@ #include "cantera/kinetics/InterfaceRate.h" #include "cantera/kinetics/PlogRate.h" #include "cantera/kinetics/TwoTempPlasmaRate.h" -#include "cantera/kinetics/DetailedVVVTRate.h" +#include "cantera/kinetics/VibrationalRelaxationRate.h" namespace Cantera { @@ -41,9 +41,9 @@ ReactionRateFactory::ReactionRateFactory() return new TwoTempPlasmaRate(node, rate_units); }); - // DetailedVVVTRate evaluator - reg("detailed-vv-vt", [](const AnyMap& node, const UnitStack& rate_units) { - return new DetailedVVVTRate(node, rate_units); + // VibrationalRelaxationRate evaluator + reg("vibrational-relaxation", [](const AnyMap& node, const UnitStack& rate_units) { + return new VibrationalRelaxationRate(node, rate_units); }); // ElectronCollisionPlasmaRate evaluator diff --git a/src/kinetics/TwoTempPlasmaRate.cpp b/src/kinetics/TwoTempPlasmaRate.cpp index bd7bbe73dbd..0b6563e98e3 100644 --- a/src/kinetics/TwoTempPlasmaRate.cpp +++ b/src/kinetics/TwoTempPlasmaRate.cpp @@ -63,9 +63,7 @@ TwoTempPlasmaRate::TwoTempPlasmaRate(double A, double b, double Ea, double EE) TwoTempPlasmaRate::TwoTempPlasmaRate(double A, double b, double Ea, double EE, double bg) : ArrheniusBase(A, b, Ea) { - m_Ea_str = "Ea-gas"; - m_E4_str = "Ea-electron"; - m_E4_R = EE / GasConstant; + TwoTempPlasmaRate(A, b, Ea, EE); m_bg = bg; } diff --git a/src/kinetics/VibrationalRelaxationRate.cpp b/src/kinetics/VibrationalRelaxationRate.cpp new file mode 100644 index 00000000000..08dec883de7 --- /dev/null +++ b/src/kinetics/VibrationalRelaxationRate.cpp @@ -0,0 +1,835 @@ +//! @file VibrationalRelaxationRate.cpp + +// This file is part of Cantera. See License.txt in the top-level directory or +// at https://cantera.org/license.txt for license and copyright information. + +#include "cantera/kinetics/VibrationalRelaxationRate.h" +#include "cantera/kinetics/Reaction.h" +#include "cantera/thermo/ThermoPhase.h" +#include "cantera/base/global.h" + +#include +#include +#include +#include +#include +#include + +namespace Cantera +{ + +namespace +{ + +void requireNoKey(const AnyMap& node, const string& key, + const string& model, const string& where) +{ + if (node.hasKey(key)) { + throw InputFileError(where, node, + "Key '{}' is not allowed for vibration_model '{}'.", key, model); + } +} + + +void requireKey(const AnyMap& node, const string& key, + const string& model, const string& where) +{ + if (!node.hasKey(key)) { + throw InputFileError(where, node, + "Missing required key '{}' for vibration_model '{}'.", key, model); + } +} + + +bool isVibrationalSpecies(const string& name) +{ + // Supported names: + // + // N2(v) + // O2(v0) + // O2(v1) + // O2(v12) + // + // The old form O2(v=1) is intentionally not required. + const auto pos = name.find("(v"); + return pos != string::npos && !name.empty() && name.back() == ')'; +} + + +string groundStateName(const string& name) +{ + const auto pos = name.find("(v"); + if (pos == string::npos) { + return name; + } + return name.substr(0, pos); +} + + +string vibrationalFamilyName(const string& name) +{ + // Collapse O2(v1), O2(v2), O2(v12), and O2(v) into O2(v). + return groundStateName(name) + "(v)"; +} + + +double compositionSum(const Composition& comp) +{ + double sum = 0.0; + for (const auto& item : comp) { + sum += item.second; + } + return sum; +} + + +Composition replaceVibrationalSpeciesByGroundState(const Composition& comp) +{ + Composition out; + for (const auto& item : comp) { + const string& name = item.first; + const double value = item.second; + + if (isVibrationalSpecies(name)) { + out[groundStateName(name)] += value; + } else { + out[name] += value; + } + } + return out; +} + + +bool sameComposition(const Composition& a, const Composition& b, + double tol = 1e-12) +{ + std::set names; + + for (const auto& item : a) { + names.insert(item.first); + } + for (const auto& item : b) { + names.insert(item.first); + } + + for (const auto& name : names) { + double av = 0.0; + double bv = 0.0; + + const auto ait = a.find(name); + if (ait != a.end()) { + av = ait->second; + } + + const auto bit = b.find(name); + if (bit != b.end()) { + bv = bit->second; + } + + if (std::abs(av - bv) > tol) { + return false; + } + } + + return true; +} + + +std::vector vibrationalSpeciesInComposition(const Composition& comp) +{ + std::vector out; + + for (const auto& item : comp) { + const string& name = item.first; + const double value = item.second; + + if (value != 0.0 && isVibrationalSpecies(name)) { + out.push_back(name); + } + } + + return out; +} + + +// Registry used to check that a given vibrational family uses exactly one +// relaxation model inside one Kinetics object. +// +// Allowed: +// +// N2(v) -> constant +// O2(v) -> detailed-vv-vt +// NH3(v) -> starikovskiy +// +// Forbidden: +// +// N2(v) + O -> castela +// N2(v) + N2 -> starikovskiy +// +std::map> s_modelByFamily; +std::set s_warnedMixedModels; + + +void registerVibrationalModelConsistency(const Kinetics& kin, + const string& family, + const string& model, + const AnyMap& input) +{ + auto& modelByFamily = s_modelByFamily[&kin]; + + const auto existing = modelByFamily.find(family); + if (existing != modelByFamily.end() && existing->second != model) { + throw InputFileError("VibrationalRelaxationRate::setContext", input, + "Inconsistent vibration_model for vibrational family '{}'. " + "This family was already registered with model '{}', but the " + "current reaction uses model '{}'. A given vibrational family " + "must use exactly one relaxation model.", + family, existing->second, model); + } + + modelByFamily[family] = model; + + std::set models; + for (const auto& item : modelByFamily) { + models.insert(item.second); + } + + if (models.size() > 1 && !s_warnedMixedModels.count(&kin)) { + std::ostringstream msg; + msg << "Multiple vibrational relaxation models were detected in the " + << "same kinetics object. This is allowed only if each " + << "vibrational family is internally consistent:\n"; + + for (const auto& item : modelByFamily) { + msg << " - " << item.first << ": " << item.second << "\n"; + } + + warn_user("VibrationalRelaxationRate::setContext", msg.str()); + s_warnedMixedModels.insert(&kin); + } +} + + +string inferRelaxingFamily(const Reaction& rxn) +{ + const auto vibReactants = vibrationalSpeciesInComposition(rxn.reactants); + + if (vibReactants.empty()) { + throw InputFileError("VibrationalRelaxationRate::setContext", rxn.input, + "A vibrational-relaxation reaction must contain at least one " + "vibrational reactant, for example O2(v1), N2(v), or NH3(v1)."); + } + + return vibrationalFamilyName(vibReactants.front()); +} + + +void validateSimpleRelaxationToGroundState(const Reaction& rxn, + const string& model) +{ + const auto vibReactants = vibrationalSpeciesInComposition(rxn.reactants); + const auto vibProducts = vibrationalSpeciesInComposition(rxn.products); + + if (vibReactants.size() != 1) { + throw InputFileError("VibrationalRelaxationRate::setContext", rxn.input, + "vibration_model '{}' expects exactly one vibrational reactant.", + model); + } + + if (!vibProducts.empty()) { + throw InputFileError("VibrationalRelaxationRate::setContext", rxn.input, + "vibration_model '{}' describes relaxation to the ground state " + "and therefore does not allow vibrational products.", model); + } + + const Composition relaxedReactants = + replaceVibrationalSpeciesByGroundState(rxn.reactants); + + if (!sameComposition(relaxedReactants, rxn.products)) { + throw InputFileError("VibrationalRelaxationRate::setContext", rxn.input, + "Invalid vibrational relaxation stoichiometry for model '{}'. " + "Expected a reaction equivalent to X(v) + M => X + M, where the " + "collider M is unchanged.", model); + } + + if (std::abs(compositionSum(rxn.reactants) - 2.0) > 1e-12 + || std::abs(compositionSum(rxn.products) - 2.0) > 1e-12) + { + throw InputFileError("VibrationalRelaxationRate::setContext", rxn.input, + "vibration_model '{}' expects a bimolecular relaxation reaction " + "of the form X(v) + M => X + M.", model); + } +} + + +void validateCastelaReaction(const Reaction& rxn) +{ + validateSimpleRelaxationToGroundState(rxn, "castela"); + + const auto vibReactants = vibrationalSpeciesInComposition(rxn.reactants); + const string family = vibrationalFamilyName(vibReactants.front()); + + if (family != "N2(v)") { + throw InputFileError("VibrationalRelaxationRate::setContext", rxn.input, + "The Castela relaxation model is only valid for N2 vibrational " + "relaxation. Found vibrational family '{}'.", family); + } + + // Castela parameters used here are only intended for the following + // colliders: N2, O2, and O. + const Composition collapsedReactants = + replaceVibrationalSpeciesByGroundState(rxn.reactants); + + string collider = ""; + for (const auto& item : collapsedReactants) { + const string& name = item.first; + const double value = item.second; + + if (name == "N2") { + if (std::abs(value - 2.0) < 1e-12) { + collider = "N2"; + } + } else if (std::abs(value - 1.0) < 1e-12) { + collider = name; + } + } + + if (collider != "N2" && collider != "O2" && collider != "O") { + throw InputFileError("VibrationalRelaxationRate::setContext", rxn.input, + "The Castela relaxation model only supports colliders 'N2', " + "'O2', and 'O'. Found collider '{}'.", collider); + } +} + + +void validateDetailedRelaxationReaction(const Reaction& rxn) +{ + const auto vibReactants = vibrationalSpeciesInComposition(rxn.reactants); + const auto vibProducts = vibrationalSpeciesInComposition(rxn.products); + + if (vibReactants.empty()) { + throw InputFileError("VibrationalRelaxationRate::setContext", rxn.input, + "A detailed vibrational relaxation reaction must contain at least " + "one vibrational reactant."); + } + + // All vibrational species in a detailed relaxation reaction are required + // to belong to the same vibrational family. + const string family = vibrationalFamilyName(vibReactants.front()); + + for (const auto& sp : vibReactants) { + const string spFamily = vibrationalFamilyName(sp); + if (spFamily != family) { + throw InputFileError("VibrationalRelaxationRate::setContext", rxn.input, + "Invalid detailed vibrational relaxation reaction: all " + "vibrational reactants must belong to the same vibrational " + "family. Found '{}' and '{}'.", family, spFamily); + } + } + + for (const auto& sp : vibProducts) { + const string spFamily = vibrationalFamilyName(sp); + if (spFamily != family) { + throw InputFileError("VibrationalRelaxationRate::setContext", rxn.input, + "Invalid detailed vibrational relaxation reaction: all " + "vibrational products must belong to the same vibrational " + "family. Found '{}' and '{}'.", family, spFamily); + } + } + + // The reaction must conserve the ground-state composition once all + // vibrational labels are removed. + const Composition collapsedReactants = + replaceVibrationalSpeciesByGroundState(rxn.reactants); + + const Composition collapsedProducts = + replaceVibrationalSpeciesByGroundState(rxn.products); + + if (!sameComposition(collapsedReactants, collapsedProducts)) { + throw InputFileError("VibrationalRelaxationRate::setContext", rxn.input, + "Invalid detailed vibrational relaxation reaction: replacing all " + "vibrational species by their ground-state species does not " + "conserve stoichiometry."); + } + + // The current detailed VV/VT formulation is bimolecular. + if (std::abs(compositionSum(rxn.reactants) - 2.0) > 1e-12 + || std::abs(compositionSum(rxn.products) - 2.0) > 1e-12) + { + throw InputFileError("VibrationalRelaxationRate::setContext", rxn.input, + "Invalid detailed vibrational relaxation reaction: expected a " + "bimolecular reaction with two reactant molecules and two product " + "molecules."); + } +} + +} // namespace + + +bool DetailedVibData::update(const ThermoPhase& phase, const Kinetics& kin) +{ + const double T = phase.temperature(); + + if (T == temperature) { + return false; + } + + ReactionData::update(T); + return true; +} + + +VibrationalRelaxationRate::VibrationalRelaxationRate() +{ +} + + +VibrationalRelaxationRate::VibrationalRelaxationRate( + double A, double B, double C, double D, + double b, double scaling, double m, double E, double z) + : ArrheniusBase(A, b, 0.0) + , m_B(B) + , m_C(C) + , m_D(D) + , m_scaling(scaling) + , m_m(m) + , m_E(E) + , m_z(z) +{ +} + + +VibrationalRelaxationRate::VibrationalRelaxationRate( + const AnyMap& node, const UnitStack& rate_units) + : VibrationalRelaxationRate() +{ + setParameters(node, rate_units); +} + + +void VibrationalRelaxationRate::configureBaseFromInternalA( + const AnyMap& node, const UnitStack& rate_units, double A, double b) +{ + // Store the original input and configure reaction-rate units. + // + // We intentionally do not call ArrheniusBase::setParameters here because + // some vibration models do not expose a standard YAML rate-constant with + // both A and b. Castela is the main example. + ReactionRate::setParameters(node, rate_units); + setRateUnits(rate_units); + + m_negativeA_ok = node.getBool("negative-A", false); + + m_A = A; + m_b = b; + + if (m_A != 0.0) { + m_logA = std::log(std::abs(m_A)); + } else { + m_logA = NAN; + } + + // Critical: ArrheniusBase::validate checks this flag. + m_valid = true; +} + + +void VibrationalRelaxationRate::configureBaseFromYamlA( + const AnyMap& node, const UnitStack& rate_units, + const AnyValue& A, double b) +{ + // Store the original input and configure reaction-rate units first, so + // conversionUnits() is available for A. + ReactionRate::setParameters(node, rate_units); + setRateUnits(rate_units); + + m_negativeA_ok = node.getBool("negative-A", false); + + // Convert the user-facing YAML A value with Cantera's standard + // rate-coefficient unit conversion. + m_A = node.units().convertRateCoeff(A, conversionUnits()); + m_b = b; + + if (m_A != 0.0) { + m_logA = std::log(std::abs(m_A)); + } else { + m_logA = NAN; + } + + // Critical: ArrheniusBase::validate checks this flag. + m_valid = true; +} + + +void VibrationalRelaxationRate::setParameters(const AnyMap& node, + const UnitStack& rate_units) +{ + // The reaction rate type is always: + // + // type: vibrational-relaxation + // + // The physical model is selected separately by: + // + // vibration_model: constant + // vibration_model: detailed-vv-vt + // vibration_model: starikovskiy + // vibration_model: castela + // + // For backward compatibility, the default model is detailed-vv-vt. + m_vibration_model = "detailed-vv-vt"; + + if (node.hasKey("vibration_model")) { + m_vibration_model = node["vibration_model"].asString(); + } else if (node.hasKey("model")) { + m_vibration_model = node["model"].asString(); + } + + if (!node.hasKey("rate-constant")) { + throw InputFileError("VibrationalRelaxationRate::setParameters", node, + "A vibrational-relaxation reaction requires a 'rate-constant' " + "mapping."); + } + + const auto& rate = node["rate-constant"]; + + if (!rate.is()) { + throw InputFileError("VibrationalRelaxationRate::setParameters", node, + "The 'rate-constant' field must be a mapping."); + } + + const auto& rateMap = rate.as(); + + if (m_vibration_model == "constant") { + // Constant model: + // + // k(T) = A + + requireKey(rateMap, m_A_str, m_vibration_model, + "VibrationalRelaxationRate::setParameters"); + + requireNoKey(rateMap, m_b_str, m_vibration_model, + "VibrationalRelaxationRate::setParameters"); + requireNoKey(rateMap, "n", m_vibration_model, + "VibrationalRelaxationRate::setParameters"); + requireNoKey(rateMap, m_B_str, m_vibration_model, + "VibrationalRelaxationRate::setParameters"); + requireNoKey(rateMap, m_C_str, m_vibration_model, + "VibrationalRelaxationRate::setParameters"); + requireNoKey(rateMap, m_D_str, m_vibration_model, + "VibrationalRelaxationRate::setParameters"); + requireNoKey(rateMap, m_m_str, m_vibration_model, + "VibrationalRelaxationRate::setParameters"); + requireNoKey(rateMap, m_E_str, m_vibration_model, + "VibrationalRelaxationRate::setParameters"); + requireNoKey(rateMap, m_z_str, m_vibration_model, + "VibrationalRelaxationRate::setParameters"); + requireNoKey(rateMap, m_scaling_str, m_vibration_model, + "VibrationalRelaxationRate::setParameters"); + + configureBaseFromYamlA(node, rate_units, rateMap[m_A_str], 0.0); + + m_B = 0.0; + m_C = 0.0; + m_D = 0.0; + m_m = 2.0 / 3.0; + m_E = 0.0; + m_z = 1.0; + m_scaling = 1.0; + } + else if (m_vibration_model == "detailed-vv-vt") { + // Detailed VV/VT model: + // + // k(T) = scaling * A * exp( + // b * log(T) + // + B + // + C * T^(-1/3) + // + D * T^(-2/3) + // ) + + requireKey(rateMap, m_A_str, m_vibration_model, + "VibrationalRelaxationRate::setParameters"); + + requireNoKey(rateMap, "n", m_vibration_model, + "VibrationalRelaxationRate::setParameters"); + requireNoKey(rateMap, m_m_str, m_vibration_model, + "VibrationalRelaxationRate::setParameters"); + requireNoKey(rateMap, m_E_str, m_vibration_model, + "VibrationalRelaxationRate::setParameters"); + requireNoKey(rateMap, m_z_str, m_vibration_model, + "VibrationalRelaxationRate::setParameters"); + + ArrheniusBase::setParameters(node, rate_units); + + m_B = rateMap.hasKey(m_B_str) ? rateMap[m_B_str].asDouble() : 0.0; + m_C = rateMap.hasKey(m_C_str) ? rateMap[m_C_str].asDouble() : 0.0; + m_D = rateMap.hasKey(m_D_str) ? rateMap[m_D_str].asDouble() : 0.0; + m_m = 2.0 / 3.0; + m_E = 0.0; + m_z = 1.0; + m_scaling = rateMap.hasKey(m_scaling_str) + ? rateMap[m_scaling_str].asDouble() + : 1.0; + } + else if (m_vibration_model == "starikovskiy") { + // User-facing formula: + // + // k(T) = A * T^n * exp( + // K + // - B * T^(-1/3) + // + C * T^(-m) + // + D * T^(-z) + // ) + // + // Internal mapping: + // + // m_b = n + // m_B = K + // m_C = -B + // m_D = C + // m_m = m + // m_E = D + // m_z = z + + requireKey(rateMap, m_A_str, m_vibration_model, + "VibrationalRelaxationRate::setParameters"); + + requireNoKey(rateMap, m_b_str, m_vibration_model, + "VibrationalRelaxationRate::setParameters"); + requireNoKey(rateMap, m_scaling_str, m_vibration_model, + "VibrationalRelaxationRate::setParameters"); + + const double n = rateMap.hasKey("n") ? rateMap["n"].asDouble() : 0.0; + configureBaseFromYamlA(node, rate_units, rateMap[m_A_str], n); + + m_B = rateMap.hasKey("K") ? rateMap["K"].asDouble() : 0.0; + m_C = rateMap.hasKey("B") ? -rateMap["B"].asDouble() : 0.0; + m_D = rateMap.hasKey("C") ? rateMap["C"].asDouble() : 0.0; + m_m = rateMap.hasKey("m") ? rateMap["m"].asDouble() : 1.0; + m_E = rateMap.hasKey("D") ? rateMap["D"].asDouble() : 0.0; + m_z = rateMap.hasKey("z") ? rateMap["z"].asDouble() : 1.0; + m_scaling = 1.0; + + if (m_m <= 0.0 || m_z <= 0.0) { + throw InputFileError("VibrationalRelaxationRate::setParameters", node, + "The Starikovskiy exponents 'm' and 'z' must be positive."); + } + } + else if (m_vibration_model == "castela") { + // Castela model: + // + // Original relaxation time: + // + // tau_k = p0 / p_k + // * exp[a_k * (T^(-1/3) - b_k) - 18.42] + // + // Equivalent bimolecular rate coefficient: + // + // k_k(T) = R T / p0 + // * exp[18.42 + a_k b_k - a_k T^(-1/3)] + // + // Internal mapping: + // + // A = R / p0 + // b = 1 + // B = 18.42 + a_k b_k + // C = -a_k + // D = 0 + // E = 0 + // scaling = 1 + + requireKey(rateMap, "a", m_vibration_model, + "VibrationalRelaxationRate::setParameters"); + requireKey(rateMap, "b", m_vibration_model, + "VibrationalRelaxationRate::setParameters"); + + requireNoKey(rateMap, m_A_str, m_vibration_model, + "VibrationalRelaxationRate::setParameters"); + requireNoKey(rateMap, "n", m_vibration_model, + "VibrationalRelaxationRate::setParameters"); + requireNoKey(rateMap, "K", m_vibration_model, + "VibrationalRelaxationRate::setParameters"); + requireNoKey(rateMap, m_B_str, m_vibration_model, + "VibrationalRelaxationRate::setParameters"); + requireNoKey(rateMap, m_C_str, m_vibration_model, + "VibrationalRelaxationRate::setParameters"); + requireNoKey(rateMap, m_D_str, m_vibration_model, + "VibrationalRelaxationRate::setParameters"); + requireNoKey(rateMap, m_m_str, m_vibration_model, + "VibrationalRelaxationRate::setParameters"); + requireNoKey(rateMap, m_E_str, m_vibration_model, + "VibrationalRelaxationRate::setParameters"); + requireNoKey(rateMap, m_z_str, m_vibration_model, + "VibrationalRelaxationRate::setParameters"); + requireNoKey(rateMap, m_scaling_str, m_vibration_model, + "VibrationalRelaxationRate::setParameters"); + + m_castela_a = rateMap["a"].asDouble(); + m_castela_b = rateMap["b"].asDouble(); + + if (rateMap.hasKey(m_reference_pressure_str)) { + m_referencePressure = rateMap.convert(m_reference_pressure_str, "Pa"); + } else { + m_referencePressure = OneAtm; + } + + if (m_referencePressure <= 0.0) { + throw InputFileError("VibrationalRelaxationRate::setParameters", node, + "Castela reference-pressure must be positive."); + } + + configureBaseFromInternalA( + node, rate_units, GasConstant / m_referencePressure, 1.0); + + m_B = 18.42 + m_castela_a * m_castela_b; + m_C = -m_castela_a; + m_D = 0.0; + m_m = 2.0 / 3.0; + m_E = 0.0; + m_z = 1.0; + m_scaling = 1.0; + } + else { + throw InputFileError("VibrationalRelaxationRate::setParameters", node, + "Unrecognized vibration_model '{}'. Expected 'detailed-vv-vt', " + "'starikovskiy', 'castela', or 'constant'.", + m_vibration_model); + } +} + + +void VibrationalRelaxationRate::getParameters(AnyMap& node) const +{ + if (!valid()) { + return; + } + + if (allowNegativePreExponentialFactor()) { + node["negative-A"] = true; + } + + node["vibration_model"] = m_vibration_model; + + AnyMap rateNode; + + auto storePreExponentialFactor = [&](AnyMap& target, double A) { + if (conversionUnits().factor() != 0.0) { + target[m_A_str].setQuantity(A, conversionUnits()); + } else { + target[m_A_str] = A; + target["__unconvertible__"] = true; + } + }; + + if (m_vibration_model == "constant") { + const double tol = 1e-12; + + if (std::abs(m_b) > tol + || std::abs(m_B) > tol + || std::abs(m_C) > tol + || std::abs(m_D) > tol + || std::abs(m_E) > tol) + { + throw InputFileError("VibrationalRelaxationRate::getParameters", node, + "Cannot serialize this rate as 'constant': the internal " + "parameters contain temperature-dependent terms."); + } + + storePreExponentialFactor(rateNode, m_scaling * m_A); + } + else if (m_vibration_model == "detailed-vv-vt") { + storePreExponentialFactor(rateNode, m_A); + + rateNode[m_b_str] = m_b; + rateNode[m_B_str] = m_B; + rateNode[m_C_str] = m_C; + rateNode[m_D_str] = m_D; + rateNode[m_scaling_str] = m_scaling; + } + else if (m_vibration_model == "starikovskiy") { + storePreExponentialFactor(rateNode, m_A); + + rateNode["n"] = m_b; + rateNode["K"] = m_B; + rateNode["B"] = -m_C; + rateNode["C"] = m_D; + rateNode["m"] = m_m; + rateNode["D"] = m_E; + rateNode["z"] = m_z; + } + else if (m_vibration_model == "castela") { + const double tol = 1e-12; + + if (std::abs(m_b - 1.0) > tol + || std::abs(m_D) > tol + || std::abs(m_E) > tol + || std::abs(m_scaling - 1.0) > tol) + { + throw InputFileError("VibrationalRelaxationRate::getParameters", node, + "Cannot serialize this rate as 'castela': the internal " + "parameters are not consistent with the Castela form. " + "Expected b = 1, D = 0, E = 0, and scaling = 1."); + } + + if (m_referencePressure <= 0.0) { + throw InputFileError("VibrationalRelaxationRate::getParameters", node, + "Cannot serialize this rate as 'castela': " + "reference-pressure must be positive."); + } + + rateNode["a"] = m_castela_a; + rateNode["b"] = m_castela_b; + rateNode[m_reference_pressure_str].setQuantity(m_referencePressure, "Pa"); + } + else { + throw InputFileError("VibrationalRelaxationRate::getParameters", node, + "Unrecognized vibration_model '{}'. Expected 'detailed-vv-vt', " + "'starikovskiy', 'castela', or 'constant'.", + m_vibration_model); + } + + rateNode.setFlowStyle(); + node["rate-constant"] = std::move(rateNode); +} + + +double VibrationalRelaxationRate::ddTScaledFromStruct( + const DetailedVibData& shared_data) const +{ + const double invT = shared_data.recipT; + const double invT13 = std::cbrt(invT); + + return m_b * invT + - (m_C / 3.0) * invT13 * invT + - m_m * m_D * std::pow(invT, m_m) * invT + - m_z * m_E * std::pow(invT, m_z) * invT; +} + + +void VibrationalRelaxationRate::setContext(const Reaction& rxn, const Kinetics& kin) +{ + if (rxn.reversible) { + throw InputFileError("VibrationalRelaxationRate::setContext", rxn.input, + "Vibrational relaxation rates do not support reversible " + "reactions."); + } + + const string family = inferRelaxingFamily(rxn); + + if (m_vibration_model == "constant") { + validateSimpleRelaxationToGroundState(rxn, "constant"); + } else if (m_vibration_model == "castela") { + validateCastelaReaction(rxn); + } else if (m_vibration_model == "starikovskiy") { + validateSimpleRelaxationToGroundState(rxn, "starikovskiy"); + } else if (m_vibration_model == "detailed-vv-vt") { + validateDetailedRelaxationReaction(rxn); + } else { + throw InputFileError("VibrationalRelaxationRate::setContext", rxn.input, + "Unrecognized vibration_model '{}'.", m_vibration_model); + } + + registerVibrationalModelConsistency( + kin, family, m_vibration_model, rxn.input); +} + +} // namespace Cantera \ No newline at end of file diff --git a/test/kinetics/kineticsFromScratch.cpp b/test/kinetics/kineticsFromScratch.cpp index bcd4e1d92f9..332ab9cb111 100644 --- a/test/kinetics/kineticsFromScratch.cpp +++ b/test/kinetics/kineticsFromScratch.cpp @@ -14,7 +14,7 @@ #include "cantera/kinetics/InterfaceRate.h" #include "cantera/kinetics/PlogRate.h" #include "cantera/kinetics/TwoTempPlasmaRate.h" -#include "cantera/kinetics/DetailedVVVTRate.h" +#include "cantera/kinetics/VibrationalRelaxationRate.h" #include "cantera/base/Array.h" #include "cantera/base/stringUtils.h" diff --git a/test/kinetics/kineticsFromYaml.cpp b/test/kinetics/kineticsFromYaml.cpp index c34a8633f8a..d03da101282 100644 --- a/test/kinetics/kineticsFromYaml.cpp +++ b/test/kinetics/kineticsFromYaml.cpp @@ -13,7 +13,7 @@ #include "cantera/kinetics/InterfaceRate.h" #include "cantera/kinetics/PlogRate.h" #include "cantera/kinetics/TwoTempPlasmaRate.h" -#include "cantera/kinetics/DetailedVVVTRate.h" +#include "cantera/kinetics/VibrationalRelaxationRate.h" #include "cantera/thermo/SurfPhase.h" #include "cantera/thermo/ThermoFactory.h" #include "cantera/base/Array.h" From 63963dcdfde4ab22b8b20efcebcfcd58ff610e41 Mon Sep 17 00:00:00 2001 From: Gaetanosaure Date: Thu, 21 May 2026 14:35:07 +0200 Subject: [PATCH 21/29] Added some details in the header of the VibrationalRelaxationRate class with some references as to where the model and rate formulas come from and where data for these can be found. --- .../kinetics/VibrationalRelaxationRate.h | 34 +++++++++++++++++++ src/thermo/PlasmaPhase.cpp | 1 + 2 files changed, 35 insertions(+) diff --git a/include/cantera/kinetics/VibrationalRelaxationRate.h b/include/cantera/kinetics/VibrationalRelaxationRate.h index 59592f6ec27..c22d21d8ed3 100644 --- a/include/cantera/kinetics/VibrationalRelaxationRate.h +++ b/include/cantera/kinetics/VibrationalRelaxationRate.h @@ -4,6 +4,40 @@ //! This file is part of Cantera. See License.txt in the top-level directory or //! at https://cantera.org/license.txt for license and copyright information. + +//! It implements the relaxation rate of vibrationally excited species in plasma kinetics. +//! Four models are currently supported by this class: 'constant', 'detailed-vv-vt', 'starikovskiy', and 'castela'. + +//! The 'constant' model relaxes the vibrational species with a constant rate coefficient. It could just as well be an +//! Arrhenius rate, but the constant model is provided for convenience and to avoid confusion with conventional Arrhenius rates. + +//! The 'detailed-vv-vt' model is meant to fully resolve vibrational relaxation by taking into account all vibrational species +//! in the phase (Ex: N2(v=1-8)) and solves for their V-T and V-V relaxation. The rates are based on the SSH theory detailed +//! in the Chapter 7 of the following book: +//! Capitelli, M., Ferreira, C. M., Gordiets, B. F., & Osipov, A. I. (2013). +//! Plasma kinetics in atmospheric gases (Vol. 31). Springer Science & Business Media. + +//! The 'castela' model is meant to be used only for N2 vibrational relaxation, by collisions with N2, O2 and O exclusively. +//! It implements the mean vibrational energy relaxation model using a fictive cantera species. +//! The Castela model is based on the following paper: +//! Castela, M., Fiorina, B., Coussement, A., Gicquel, O., Darabiha, N., & Laux, C. O. (2016). +//! Modelling the impact of non-equilibrium discharges on reactive mixtures for simulations of +//! plasma-assisted ignition in turbulent flows. Combustion and flame, 166, 133-147. + +//! The so-called 'starikovskiy' model is an extension of the Castela model to several vibrational species and additional colliders. +//! A lot of vibrational relaxation rates can be found in Table 1 of the following paper: +//! Starikovskiy, A., & Aleksandrov, N. (2013). Plasma-assisted ignition and combustion. +//! Progress in Energy and Combustion Science, 39(1), 61-110. +//! The rates for the vibrational relaxation of NH3 can be found in the reaction mechanism provided in supplementary material +//! of the following paper: +//! Zhong, H. et al. (2023) “Understanding non-equilibrium N2O/NOx chemistry in plasma-assisted low-temperature NH3 oxidation,” +//! Combustion and Flame, 256. +//! More rates for the vibrational relaxation of CH$ can be found in the following paper: +//! Popov, N.A. (2016) “Kinetics of plasma-assisted combustion: Effect of non-equilibrium excitation on the ignition +//! and oxidation of combustible mixtures,” Plasma Sources Science and Technology. Institute of Physics Publishing. + + + #ifndef CT_VIBRATIONALRELAXATIONRATE_H #define CT_VIBRATIONALRELAXATIONRATE_H diff --git a/src/thermo/PlasmaPhase.cpp b/src/thermo/PlasmaPhase.cpp index 5aa549015c6..dbddc7cdcff 100644 --- a/src/thermo/PlasmaPhase.cpp +++ b/src/thermo/PlasmaPhase.cpp @@ -1201,6 +1201,7 @@ double PlasmaPhase::jouleHeatingPower() const return sigma * E * E; // W/m^3 } +// actually be careful, this is inelastic + growth double PlasmaPhase::inelasticPower() { // Joule heating: sigma * E^2 [W/m^3] From 1097c7ca9f965ba74e4851f109200203d48831d2 Mon Sep 17 00:00:00 2001 From: Gaetanosaure Date: Fri, 22 May 2026 10:21:19 +0200 Subject: [PATCH 22/29] more comments --- include/cantera/kinetics/VibrationalRelaxationRate.h | 10 ++++++++-- 1 file changed, 8 insertions(+), 2 deletions(-) diff --git a/include/cantera/kinetics/VibrationalRelaxationRate.h b/include/cantera/kinetics/VibrationalRelaxationRate.h index c22d21d8ed3..9c4054e0c50 100644 --- a/include/cantera/kinetics/VibrationalRelaxationRate.h +++ b/include/cantera/kinetics/VibrationalRelaxationRate.h @@ -12,10 +12,16 @@ //! Arrhenius rate, but the constant model is provided for convenience and to avoid confusion with conventional Arrhenius rates. //! The 'detailed-vv-vt' model is meant to fully resolve vibrational relaxation by taking into account all vibrational species -//! in the phase (Ex: N2(v=1-8)) and solves for their V-T and V-V relaxation. The rates are based on the SSH theory detailed -//! in the Chapter 7 of the following book: +//! in the phase (Ex: N2(v=1-8)) and solves for their V-T and V-V relaxation. The scaling of the rates are based on the SSH theory +//! detailed in Chapter 7 of the following book: //! Capitelli, M., Ferreira, C. M., Gordiets, B. F., & Osipov, A. I. (2013). //! Plasma kinetics in atmospheric gases (Vol. 31). Springer Science & Business Media. +//! The simplified version of the SSH theory implemented in this model is based on the harmonic oscillator approximation and can be found +//! in the following paper (equations 18 and 19): +//! Vasco Guerra et al 2019 Plasma Sources Sci. Technol. 28 073001 +//! The k10 rates are taken from //! Zhong, H. et al. (2023) “Understanding non-equilibrium N2O/NOx chemistry in plasma-assisted +//! low-temperature NH3 oxidation,” Combustion and Flame, 256, which itself took those rates from both the Capitelli book and +//! the article from Strikovskiy and Aleksandrov (2013), see more below. //! The 'castela' model is meant to be used only for N2 vibrational relaxation, by collisions with N2, O2 and O exclusively. //! It implements the mean vibrational energy relaxation model using a fictive cantera species. From d6eec198489cf149ab2cbe59de59087fcc74042e Mon Sep 17 00:00:00 2001 From: Gaetanosaure Date: Fri, 22 May 2026 16:42:08 +0200 Subject: [PATCH 23/29] [thermo] Fixed the EEDF calculation - as shoown on the official website example, the computed EEDF was wrong. The algorithm is fixed by these changes --- src/thermo/EEDFTwoTermApproximation.cpp | 106 +++++++++++++++++++++++- 1 file changed, 103 insertions(+), 3 deletions(-) diff --git a/src/thermo/EEDFTwoTermApproximation.cpp b/src/thermo/EEDFTwoTermApproximation.cpp index 615c9f1e93c..5c793f12bdb 100644 --- a/src/thermo/EEDFTwoTermApproximation.cpp +++ b/src/thermo/EEDFTwoTermApproximation.cpp @@ -672,38 +672,138 @@ void EEDFTwoTermApproximation::updateMoleFractions() } } +// Base implementation of the function +// void EEDFTwoTermApproximation::calculateTotalCrossSection() +// { +// m_totalCrossSectionCenter.assign(m_points, 0.0); +// m_totalCrossSectionEdge.assign(m_points + 1, 0.0); +// for (size_t k = 0; k < m_phase->nCollisions(); k++) { +// auto x = m_phase->collisionRate(k)->energyLevels(); +// auto y = m_phase->collisionRate(k)->crossSections(); + +// for (size_t i = 0; i < m_points; i++) { +// m_totalCrossSectionCenter[i] += m_X_targets[m_klocTargets[k]] * +// linearInterp(m_gridCenter[i], x, y); +// } +// 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); +// } +// } +// } + +// // base implementation of the function +// void EEDFTwoTermApproximation::calculateTotalElasticCrossSection() +// { +// m_sigmaElastic.clear(); +// m_sigmaElastic.resize(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); +// 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); +// } +// } +// } + +// The old implementation counted the effective cros section as an elastic contribution +// thus double counting the inelastic contributions to the total cross section. void EEDFTwoTermApproximation::calculateTotalCrossSection() { m_totalCrossSectionCenter.assign(m_points, 0.0); m_totalCrossSectionEdge.assign(m_points + 1, 0.0); + + std::vector has_effective_for_target(m_phase->nSpecies(), false); + + // build this list first to avoid scanning through all collisions for each target twice in the loop. + for (size_t ke : m_phase->kElastic()) { + if (m_phase->collisionRate(ke)->kind() == "effective") { + has_effective_for_target[m_phase->targetIndex(ke)] = true; + } + } + for (size_t k = 0; k < m_phase->nCollisions(); k++) { + const std::string& kind = m_phase->collisionRate(k)->kind(); + const size_t target = m_phase->targetIndex(k); + + // If this target has an effective cross section, then that effective + // cross section already contains elastic + inelastic contributions. + // So only the effective cross section contributes to sigma_total. + if (has_effective_for_target[target] && kind != "effective") { + continue; + } + auto x = m_phase->collisionRate(k)->energyLevels(); auto y = m_phase->collisionRate(k)->crossSections(); + const double X_target = m_X_targets[m_klocTargets[k]]; + for (size_t i = 0; i < m_points; i++) { - m_totalCrossSectionCenter[i] += m_X_targets[m_klocTargets[k]] * + m_totalCrossSectionCenter[i] += X_target * linearInterp(m_gridCenter[i], x, y); } + for (size_t i = 0; i < m_points + 1; i++) { - m_totalCrossSectionEdge[i] += m_X_targets[m_klocTargets[k]] * + m_totalCrossSectionEdge[i] += X_target * linearInterp(m_gridEdge[i], x, y); } } } +// new implementation of the function proposed by nicolas void EEDFTwoTermApproximation::calculateTotalElasticCrossSection() { m_sigmaElastic.clear(); m_sigmaElastic.resize(m_points, 0.0); for (size_t k : m_phase->kElastic()) { + + writelog("Enter with an elastic cs\n"); auto x = m_phase->collisionRate(k)->energyLevels(); auto y = m_phase->collisionRate(k)->crossSections(); + vector y_elastic(y.data(), y.data() + y.size()); + // Note: // moleFraction(m_kTargets[k]) <=> m_X_targets[m_klocTargets[k]] double mass_ratio = ElectronMass / (m_phase->molecularWeight(m_kTargets[k]) / Avogadro); + + if (m_phase->collisionRate(k)->kind()=="effective") { + + writelog("Enter with an elastic cs that is actually an effective\n"); + + for (size_t ki : m_phase->kInelastic()) + { + if(m_phase->targetIndex(ki) == m_phase->targetIndex(k)){ + writelog("loop over inelastic processes: process {}\n", ki); + auto xi = m_phase->collisionRate(ki)->energyLevels(); + auto yi = m_phase->collisionRate(ki)->crossSections(); + for (size_t i = 0; i < x.size(); i++) + { + y_elastic[i] -= linearInterp(x[i], xi, yi); + } + } + } + // check that the reconstructed elastic cross section is non-negative. + for (size_t i = 0; i < y_elastic.size(); i++) { + if (y_elastic[i] < 0.0) { + if (y_elastic[i] > -1e-30) { + y_elastic[i] = 0.0; + } else { + writelog("Warning: reconstructed elastic cross section is negative " + "for collision {} at energy {}: {}\n", + k, x[i], y_elastic[i]); + y_elastic[i] = 0.0; + } + } + } + } + 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); + linearInterp(m_gridEdge[i], x, y_elastic); } } } From c137c1d00f8e81de6794747f11d89d1634920bc4 Mon Sep 17 00:00:00 2001 From: Gaetanosaure Date: Fri, 22 May 2026 20:56:03 +0200 Subject: [PATCH 24/29] [thermo] Added function checkSpeciesNoCrossSections for additional errors to the user during the run. --- include/cantera/thermo/EEDFTwoTermApproximation.h | 3 +++ src/thermo/EEDFTwoTermApproximation.cpp | 12 ++++++++++++ 2 files changed, 15 insertions(+) diff --git a/include/cantera/thermo/EEDFTwoTermApproximation.h b/include/cantera/thermo/EEDFTwoTermApproximation.h index 666bbfca4d6..26eb2877d98 100644 --- a/include/cantera/thermo/EEDFTwoTermApproximation.h +++ b/include/cantera/thermo/EEDFTwoTermApproximation.h @@ -283,6 +283,9 @@ class EEDFTwoTermApproximation //! Compute the total (elastic + inelastic) cross section void calculateTotalCrossSection(); + //! + void checkSpeciesNoCrossSection(); + //! Compute the L1 norm of a function f defined over a given energy grid. //! //! @param f Vector representing the function values (EEDF) diff --git a/src/thermo/EEDFTwoTermApproximation.cpp b/src/thermo/EEDFTwoTermApproximation.cpp index 5c793f12bdb..84fde1efe71 100644 --- a/src/thermo/EEDFTwoTermApproximation.cpp +++ b/src/thermo/EEDFTwoTermApproximation.cpp @@ -143,6 +143,7 @@ int EEDFTwoTermApproximation::calculateDistributionFunction() } updateMoleFractions(); + checkSpeciesNoCrossSection(); updateCrossSections(); const double EN = m_phase->reducedElectricField(); @@ -808,6 +809,17 @@ void EEDFTwoTermApproximation::calculateTotalElasticCrossSection() } } +void EEDFTwoTermApproximation::checkSpeciesNoCrossSection() +{ + // warn that a specific species needs cross-section data. + for (size_t k : m_kOthers) { + if (m_phase->moleFraction(k) > m_moleFractionThreshold) { + writelog("EEDFTwoTermApproximation:checkSpeciesNoCrossSection\n"); + writelog("Warning:The mole fraction of species {} is more than 0.01 (X = {:.3g}) but it has no cross-section data\n", m_phase->speciesName(k), m_phase->moleFraction(k)); + } + } +} + void EEDFTwoTermApproximation::setGridCache() { m_sigma.clear(); From 98e070bd4c82db20560d06e526abd0b06c6c474a Mon Sep 17 00:00:00 2001 From: Gaetanosaure Date: Fri, 22 May 2026 23:05:05 +0200 Subject: [PATCH 25/29] [thermo] added some routines to check for the potential declaration of vibrational fictive species, and in this case they send a warning if X_A(v) / ( X_A(v) + X(A)) > 1 % --- include/cantera/thermo/PlasmaPhase.h | 15 +++++ src/thermo/PlasmaPhase.cpp | 97 ++++++++++++++++++++++++++++ 2 files changed, 112 insertions(+) diff --git a/include/cantera/thermo/PlasmaPhase.h b/include/cantera/thermo/PlasmaPhase.h index 1a6f5a085c8..45feb588f79 100644 --- a/include/cantera/thermo/PlasmaPhase.h +++ b/include/cantera/thermo/PlasmaPhase.h @@ -457,6 +457,21 @@ class PlasmaPhase: public IdealGasPhase //! Calculates the power losses (W/m³) of the plasma electrons through inelastic collisions. double inelasticPower(); + struct VibrationalReservoirSpecies { + size_t reservoirIndex = npos; + size_t baseSpeciesIndex = npos; + }; + + std::vector m_vibrationalReservoirSpecies; + + bool m_vibrationalReservoirSpeciesNeedUpdate = true; + + double m_vibrationalMoleFractionThreshold = 1e-2; + double m_vibrationalAbsoluteMoleFractionThreshold = 1e-20; + + void updateVibrationalReservoirSpecies(); + void checkVibrationalReservoirMoleFractions(); + protected: void updateThermo() const override; diff --git a/src/thermo/PlasmaPhase.cpp b/src/thermo/PlasmaPhase.cpp index dbddc7cdcff..bc100e7ca45 100644 --- a/src/thermo/PlasmaPhase.cpp +++ b/src/thermo/PlasmaPhase.cpp @@ -19,6 +19,22 @@ namespace Cantera { namespace { const double gamma = sqrt(2 * ElectronCharge / ElectronMass); + + bool isVibrationalReservoirName(const string& name, string& baseName) + { + const string suffix = "(v)"; + + if (name.size() <= suffix.size()) { + return false; + } + + if (name.compare(name.size() - suffix.size(), suffix.size(), suffix) != 0) { + return false; + } + + baseName = name.substr(0, name.size() - suffix.size()); + return !baseName.empty(); + } } PlasmaPhase::PlasmaPhase(const string& inputFile, const string& id_) @@ -656,6 +672,7 @@ void PlasmaPhase::addStandaloneElectronCollision(const AnyMap& item) addCollision(R); } + bool PlasmaPhase::addSpecies(shared_ptr spec) { bool added = IdealGasPhase::addSpecies(spec); @@ -673,6 +690,11 @@ bool PlasmaPhase::addSpecies(shared_ptr spec) "Only one electron species is allowed.", spec->name); } } + + if (added) { + m_vibrationalReservoirSpeciesNeedUpdate = true; + } + return added; } @@ -1221,8 +1243,83 @@ double PlasmaPhase::intrinsicHeating() const double qJ = jouleHeatingPower(); checkFinite(qJ); + checkVibrationalReservoirMoleFractions(); + return qJ; } +void PlasmaPhase::updateVibrationalReservoirSpecies() +{ + m_vibrationalReservoirSpecies.clear(); + + for (size_t k = 0; k < nSpecies(); k++) { + string baseName; + const string& reservoirName = speciesName(k); + + if (!isVibrationalReservoirName(reservoirName, baseName)) { + continue; + } + + size_t kBase = speciesIndex(baseName, false); + + if (kBase == npos) { + warn_user("PlasmaPhase::updateVibrationalReservoirSpecies", + "Species '{}' matches the fictive vibrational-reservoir naming " + "convention, but the corresponding base species '{}' was not found. " + "This species will not be monitored as a vibrational reservoir.", + reservoirName, baseName); + continue; + } + + VibrationalReservoirSpecies reservoir; + reservoir.reservoirIndex = k; + reservoir.baseSpeciesIndex = kBase; + + m_vibrationalReservoirSpecies.push_back(reservoir); + } + + m_vibrationalReservoirSpeciesNeedUpdate = false; +} + +void PlasmaPhase::checkVibrationalReservoirMoleFractions() +{ + if (m_vibrationalReservoirSpeciesNeedUpdate) { + updateVibrationalReservoirSpecies(); + } + + for (const auto& reservoir : m_vibrationalReservoirSpecies) { + const size_t kReservoir = reservoir.reservoirIndex; + const size_t kBase = reservoir.baseSpeciesIndex; + + const double Xv = moleFraction(kReservoir); + const double Xb = moleFraction(kBase); + + const double pool = Xv + Xb; + + if (pool <= m_vibrationalAbsoluteMoleFractionThreshold) { + continue; + } + + const double reservoirFraction = Xv / pool; + + if (reservoirFraction > m_vibrationalMoleFractionThreshold) { + const string& reservoirName = speciesName(kReservoir); + const string& baseName = speciesName(kBase); + + writelog("Warning: fictive vibrational reservoir species '{}' contains " + "{:.3e} of the total '{}' pool. " + "X({}) = {:.3e}, X({}) = {:.3e}, threshold = {:.3e}. " + "Chemistry involving '{}' may be affected because part of the " + "material is stored in an inert vibrational reservoir.\n", + reservoirName, + reservoirFraction, + baseName, + reservoirName, Xv, + baseName, Xb, + m_vibrationalMoleFractionThreshold, + baseName); + } + } +} } From aea876e3a249628da4aad99f5f4f21f261819004 Mon Sep 17 00:00:00 2001 From: Gaetanosaure Date: Fri, 22 May 2026 23:25:54 +0200 Subject: [PATCH 26/29] cleaned code and added some comments --- include/cantera/thermo/PlasmaPhase.h | 9 ++++++ src/thermo/EEDFTwoTermApproximation.cpp | 41 ++----------------------- 2 files changed, 11 insertions(+), 39 deletions(-) diff --git a/include/cantera/thermo/PlasmaPhase.h b/include/cantera/thermo/PlasmaPhase.h index 45feb588f79..2017356e7b3 100644 --- a/include/cantera/thermo/PlasmaPhase.h +++ b/include/cantera/thermo/PlasmaPhase.h @@ -457,19 +457,28 @@ class PlasmaPhase: public IdealGasPhase //! Calculates the power losses (W/m³) of the plasma electrons through inelastic collisions. double inelasticPower(); + //! Allows to store info on potential vibrational reservoir / fictive species that could be declare by the user + //! These species take the form A(v) where A is a base species and v represents the general vibration. struct VibrationalReservoirSpecies { size_t reservoirIndex = npos; size_t baseSpeciesIndex = npos; }; + //! list of detected vibrational reservoir species in the phase. std::vector m_vibrationalReservoirSpecies; + //! Flag to indicate whether the vibrational reservoir species list needs to be updated bool m_vibrationalReservoirSpeciesNeedUpdate = true; + // Thresholds for identifying vibrational reservoir species based on their mole fractions double m_vibrationalMoleFractionThreshold = 1e-2; double m_vibrationalAbsoluteMoleFractionThreshold = 1e-20; + //! Update the list of vibrational reservoir / fictive species void updateVibrationalReservoirSpecies(); + + //! Checks that declared vibrational reservoir / fictive species have mole fractions + //! low enough not to impact the phase chemistry. void checkVibrationalReservoirMoleFractions(); protected: diff --git a/src/thermo/EEDFTwoTermApproximation.cpp b/src/thermo/EEDFTwoTermApproximation.cpp index 84fde1efe71..92d9f4cfd30 100644 --- a/src/thermo/EEDFTwoTermApproximation.cpp +++ b/src/thermo/EEDFTwoTermApproximation.cpp @@ -673,46 +673,9 @@ void EEDFTwoTermApproximation::updateMoleFractions() } } -// Base implementation of the function -// void EEDFTwoTermApproximation::calculateTotalCrossSection() -// { -// m_totalCrossSectionCenter.assign(m_points, 0.0); -// m_totalCrossSectionEdge.assign(m_points + 1, 0.0); -// for (size_t k = 0; k < m_phase->nCollisions(); k++) { -// auto x = m_phase->collisionRate(k)->energyLevels(); -// auto y = m_phase->collisionRate(k)->crossSections(); - -// for (size_t i = 0; i < m_points; i++) { -// m_totalCrossSectionCenter[i] += m_X_targets[m_klocTargets[k]] * -// linearInterp(m_gridCenter[i], x, y); -// } -// 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); -// } -// } -// } - -// // base implementation of the function -// void EEDFTwoTermApproximation::calculateTotalElasticCrossSection() -// { -// m_sigmaElastic.clear(); -// m_sigmaElastic.resize(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); -// 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); -// } -// } -// } - -// The old implementation counted the effective cros section as an elastic contribution +// The former implementation counted the effective cros section as an elastic contribution // thus double counting the inelastic contributions to the total cross section. +// This implemetation avoids this drawback. void EEDFTwoTermApproximation::calculateTotalCrossSection() { m_totalCrossSectionCenter.assign(m_points, 0.0); From a0318f1725cbc9e580e543f2caa0fcb449ed6a14 Mon Sep 17 00:00:00 2001 From: Gaetanosaure Date: Fri, 22 May 2026 23:38:27 +0200 Subject: [PATCH 27/29] [yaml] Added the possibility to control the geometric grid ratio from the YAML with the parameter geometric_grid_ratio --- include/cantera/thermo/EEDFTwoTermApproximation.h | 4 ++-- src/thermo/EEDFTwoTermApproximation.cpp | 4 +--- src/thermo/PlasmaPhase.cpp | 11 +++++++++++ 3 files changed, 14 insertions(+), 5 deletions(-) diff --git a/include/cantera/thermo/EEDFTwoTermApproximation.h b/include/cantera/thermo/EEDFTwoTermApproximation.h index 26eb2877d98..7a3e27eaeb4 100644 --- a/include/cantera/thermo/EEDFTwoTermApproximation.h +++ b/include/cantera/thermo/EEDFTwoTermApproximation.h @@ -78,8 +78,8 @@ class EEDFTwoTermApproximation void setQuadraticGrid(double& kTe_max, size_t& ncell); //! Sets a geometric energy grid for the EEDF solver, defined by the maximum energy and the number of grid cells. - //! @todo the geometric ratio is currently fixed to 1.05 but this parameter could be chosen by the user in the yaml file. - void setGeometricGrid(double& kTe_max, size_t& ncell); + //! the geometric ratio defaults to 1.01 if not specified. + void setGeometricGrid(double& kTe_max, size_t& ncell, double ratio = 1.01); //! Sets a custom energy grid for the EEDF solver, defined by the user-provided vector of energy levels. void setCustomGrid(span levels); diff --git a/src/thermo/EEDFTwoTermApproximation.cpp b/src/thermo/EEDFTwoTermApproximation.cpp index 92d9f4cfd30..5f71bd2f01f 100644 --- a/src/thermo/EEDFTwoTermApproximation.cpp +++ b/src/thermo/EEDFTwoTermApproximation.cpp @@ -66,7 +66,7 @@ void EEDFTwoTermApproximation::setQuadraticGrid(double& kTe_max, size_t& ncell) setGridCache(); } -void EEDFTwoTermApproximation::setGeometricGrid(double& kTe_max, size_t& ncell) +void EEDFTwoTermApproximation::setGeometricGrid(double& kTe_max, size_t& ncell, double ratio) { m_points = ncell; @@ -75,8 +75,6 @@ void EEDFTwoTermApproximation::setGeometricGrid(double& kTe_max, size_t& ncell) m_f0.resize(m_points); m_f0_edge.resize(m_points + 1); - // @todo Make the geometric-grid ratio configurable from input. - double ratio = 1.05; double denominator = std::pow(ratio, m_points) - 1.0; if (std::abs(denominator) < 1e-14) { diff --git a/src/thermo/PlasmaPhase.cpp b/src/thermo/PlasmaPhase.cpp index bc100e7ca45..28b9ee95d84 100644 --- a/src/thermo/PlasmaPhase.cpp +++ b/src/thermo/PlasmaPhase.cpp @@ -506,7 +506,18 @@ void PlasmaPhase::setParameters(const AnyMap& phaseNode, const AnyMap& rootNode) } else if (energyLevelsDistribution == "Quadratic") { m_eedfSolver->setQuadraticGrid(initialMaxEnergy, nGridCells); } else if (energyLevelsDistribution == "Geometric") { + if (eedf.hasKey("geometric_grid_ratio")) { + double ratio = eedf["geometric_grid_ratio"].asDouble(); + if (!std::isfinite(ratio) || ratio <= 1.0) { + throw CanteraError("PlasmaPhase::setParameters", + "geometric_grid_ratio must be finite and greater than 1.0."); + } + m_eedfSolver->setGeometricGrid(initialMaxEnergy, nGridCells, ratio); + } else { + writelog("No geometric_grid_ratio key found for geometric grid in the input file. " + "Defaulting to a geometric ratio of 1.01.\n"); m_eedfSolver->setGeometricGrid(initialMaxEnergy, nGridCells); + } } else { throw CanteraError("PlasmaPhase::setParameters", "energy_levels_distribution should be Linear, Quadratic or Geometric."); From 40db335f1bdd238671cb61ef57ec10bfe03de089 Mon Sep 17 00:00:00 2001 From: Gaetanosaure Date: Mon, 25 May 2026 12:06:46 +0200 Subject: [PATCH 28/29] [thermo] Corrected the extension of cross-sections using the new function linearInterpCrossSectionZeroOutside --- .../cantera/thermo/EEDFTwoTermApproximation.h | 17 +++++++++ src/thermo/EEDFTwoTermApproximation.cpp | 38 +++++++++++++++++-- 2 files changed, 52 insertions(+), 3 deletions(-) diff --git a/include/cantera/thermo/EEDFTwoTermApproximation.h b/include/cantera/thermo/EEDFTwoTermApproximation.h index 7a3e27eaeb4..b22fa500584 100644 --- a/include/cantera/thermo/EEDFTwoTermApproximation.h +++ b/include/cantera/thermo/EEDFTwoTermApproximation.h @@ -158,6 +158,23 @@ class EEDFTwoTermApproximation //! Controls the threshold below which the EEDF solver does not solve for the EEDF but imposes a Maxwellian distribution. void setReducedFieldThresholdBeforeMaxwellianTd(double threshold); + //! An extension of the linearInterp function that returns specified values when the input is + //! out of bounds instead of returning one of the extremities of the list. + double linearInterpBounded(double x, + span xpts, + span fpts, + double below_value, + double above_value); + + //! A special use case of linearInterpBounded for cross sections + //! It is not necessary as the code could just use the previous function + //! using the values 0, but the author of this code finds it clearer to + //! the reader to have an explicit name for this use case. + double linearInterpCrossSectionZeroOutside(double x, + span xpts, + span fpts); + + protected: //! Formerly options for the EEDF solver diff --git a/src/thermo/EEDFTwoTermApproximation.cpp b/src/thermo/EEDFTwoTermApproximation.cpp index 5f71bd2f01f..1961fa2a964 100644 --- a/src/thermo/EEDFTwoTermApproximation.cpp +++ b/src/thermo/EEDFTwoTermApproximation.cpp @@ -744,7 +744,7 @@ void EEDFTwoTermApproximation::calculateTotalElasticCrossSection() auto yi = m_phase->collisionRate(ki)->crossSections(); for (size_t i = 0; i < x.size(); i++) { - y_elastic[i] -= linearInterp(x[i], xi, yi); + y_elastic[i] -= linearInterpCrossSectionZeroOutside(x[i], xi, yi); } } } @@ -754,9 +754,12 @@ void EEDFTwoTermApproximation::calculateTotalElasticCrossSection() if (y_elastic[i] > -1e-30) { y_elastic[i] = 0.0; } else { + const std::string& effectiveKind = m_phase->collisionRate(k)->kind(); + std::string effectiveTarget = m_phase->speciesName(m_phase->targetIndex(k)); + writelog("Warning: reconstructed elastic cross section is negative " - "for collision {} at energy {}: {}\n", - k, x[i], y_elastic[i]); + "for collision {} kind {} target {} at energy {}: {}\n", + k, effectiveKind, effectiveTarget, x[i], y_elastic[i]); y_elastic[i] = 0.0; } } @@ -770,6 +773,35 @@ void EEDFTwoTermApproximation::calculateTotalElasticCrossSection() } } +double EEDFTwoTermApproximation::linearInterpBounded(double x, + span xpts, + span fpts, + double below_value, + double above_value) +{ + AssertThrowMsg(!xpts.empty(), "linearInterpBounded", "x data empty"); + AssertThrowMsg(!fpts.empty(), "linearInterpBounded", "f(x) data empty"); + AssertThrowMsg(xpts.size() == fpts.size(), "linearInterpBounded", + "len(xpts) = {}, len(fpts) = {}", xpts.size(), fpts.size()); + + if (x < xpts.front()) { + return below_value; + } + + if (x > xpts.back()) { + return above_value; + } + + return linearInterp(x, xpts, fpts); +} + +double EEDFTwoTermApproximation::linearInterpCrossSectionZeroOutside(double x, + span xpts, + span fpts) +{ + return linearInterpBounded(x, xpts, fpts, 0.0, 0.0); +} + void EEDFTwoTermApproximation::checkSpeciesNoCrossSection() { // warn that a specific species needs cross-section data. From 34f651648fa7cbab1babbb5b4b78551f0f14d572 Mon Sep 17 00:00:00 2001 From: Gaetanosaure Date: Tue, 26 May 2026 18:26:02 +0200 Subject: [PATCH 29/29] [kinetics] changed the subtype name in the class VibrationalRelaxationRate representing detailed vibration which was formerly known as detailed-vv-vt to multi-state-resolved for better readability --- .../kinetics/VibrationalRelaxationRate.h | 12 ++++++------ src/kinetics/VibrationalRelaxationRate.cpp | 18 +++++++++--------- 2 files changed, 15 insertions(+), 15 deletions(-) diff --git a/include/cantera/kinetics/VibrationalRelaxationRate.h b/include/cantera/kinetics/VibrationalRelaxationRate.h index 9c4054e0c50..2df4def6e8a 100644 --- a/include/cantera/kinetics/VibrationalRelaxationRate.h +++ b/include/cantera/kinetics/VibrationalRelaxationRate.h @@ -6,12 +6,12 @@ //! It implements the relaxation rate of vibrationally excited species in plasma kinetics. -//! Four models are currently supported by this class: 'constant', 'detailed-vv-vt', 'starikovskiy', and 'castela'. +//! Four models are currently supported by this class: 'constant', 'multi-state-resolved', 'starikovskiy', and 'castela'. //! The 'constant' model relaxes the vibrational species with a constant rate coefficient. It could just as well be an //! Arrhenius rate, but the constant model is provided for convenience and to avoid confusion with conventional Arrhenius rates. -//! The 'detailed-vv-vt' model is meant to fully resolve vibrational relaxation by taking into account all vibrational species +//! The 'multi-state-resolved' model is meant to fully resolve vibrational relaxation by taking into account all vibrational species //! in the phase (Ex: N2(v=1-8)) and solves for their V-T and V-V relaxation. The scaling of the rates are based on the SSH theory //! detailed in Chapter 7 of the following book: //! Capitelli, M., Ferreira, C. M., Gordiets, B. F., & Osipov, A. I. (2013). @@ -78,7 +78,7 @@ struct DetailedVibData : public ReactionData * relaxation models: * * - `constant` - * - `detailed-vv-vt` + * - `multi-state-resolved` * - `starikovskiy` * - `castela` * @@ -108,7 +108,7 @@ struct DetailedVibData : public ReactionData * * @code{.yaml} * vibration_model: constant - * vibration_model: detailed-vv-vt + * vibration_model: multi-state-resolved * vibration_model: starikovskiy * vibration_model: castela * @endcode @@ -251,11 +251,11 @@ class VibrationalRelaxationRate : public ArrheniusBase * Accepted values: * * - `constant` - * - `detailed-vv-vt` + * - `multi-state-resolved` * - `starikovskiy` * - `castela` */ - string m_vibration_model = "detailed-vv-vt"; + string m_vibration_model = "multi-state-resolved"; //! YAML key names. string m_B_str = "B"; diff --git a/src/kinetics/VibrationalRelaxationRate.cpp b/src/kinetics/VibrationalRelaxationRate.cpp index 08dec883de7..79aabf911ee 100644 --- a/src/kinetics/VibrationalRelaxationRate.cpp +++ b/src/kinetics/VibrationalRelaxationRate.cpp @@ -158,7 +158,7 @@ std::vector vibrationalSpeciesInComposition(const Composition& comp) // Allowed: // // N2(v) -> constant -// O2(v) -> detailed-vv-vt +// O2(v) -> multi-state-resolved // NH3(v) -> starikovskiy // // Forbidden: @@ -471,12 +471,12 @@ void VibrationalRelaxationRate::setParameters(const AnyMap& node, // The physical model is selected separately by: // // vibration_model: constant - // vibration_model: detailed-vv-vt + // vibration_model: multi-state-resolved // vibration_model: starikovskiy // vibration_model: castela // - // For backward compatibility, the default model is detailed-vv-vt. - m_vibration_model = "detailed-vv-vt"; + // For backward compatibility, the default model is multi-state-resolved. + m_vibration_model = "multi-state-resolved"; if (node.hasKey("vibration_model")) { m_vibration_model = node["vibration_model"].asString(); @@ -536,7 +536,7 @@ void VibrationalRelaxationRate::setParameters(const AnyMap& node, m_z = 1.0; m_scaling = 1.0; } - else if (m_vibration_model == "detailed-vv-vt") { + else if (m_vibration_model == "multi-state-resolved") { // Detailed VV/VT model: // // k(T) = scaling * A * exp( @@ -690,7 +690,7 @@ void VibrationalRelaxationRate::setParameters(const AnyMap& node, } else { throw InputFileError("VibrationalRelaxationRate::setParameters", node, - "Unrecognized vibration_model '{}'. Expected 'detailed-vv-vt', " + "Unrecognized vibration_model '{}'. Expected 'multi-state-resolved', " "'starikovskiy', 'castela', or 'constant'.", m_vibration_model); } @@ -736,7 +736,7 @@ void VibrationalRelaxationRate::getParameters(AnyMap& node) const storePreExponentialFactor(rateNode, m_scaling * m_A); } - else if (m_vibration_model == "detailed-vv-vt") { + else if (m_vibration_model == "multi-state-resolved") { storePreExponentialFactor(rateNode, m_A); rateNode[m_b_str] = m_b; @@ -782,7 +782,7 @@ void VibrationalRelaxationRate::getParameters(AnyMap& node) const } else { throw InputFileError("VibrationalRelaxationRate::getParameters", node, - "Unrecognized vibration_model '{}'. Expected 'detailed-vv-vt', " + "Unrecognized vibration_model '{}'. Expected 'multi-state-resolved', " "'starikovskiy', 'castela', or 'constant'.", m_vibration_model); } @@ -821,7 +821,7 @@ void VibrationalRelaxationRate::setContext(const Reaction& rxn, const Kinetics& validateCastelaReaction(rxn); } else if (m_vibration_model == "starikovskiy") { validateSimpleRelaxationToGroundState(rxn, "starikovskiy"); - } else if (m_vibration_model == "detailed-vv-vt") { + } else if (m_vibration_model == "multi-state-resolved") { validateDetailedRelaxationReaction(rxn); } else { throw InputFileError("VibrationalRelaxationRate::setContext", rxn.input,