diff --git a/AUTHORS.md b/AUTHORS.md index 1461b0421b1..610a6edde5a 100644 --- a/AUTHORS.md +++ b/AUTHORS.md @@ -16,6 +16,7 @@ update, please report on Cantera's - **Philip Berndt** - **Guus Bertens** [@guusbertens](https://github.com/guusbertens) - Eindhoven University of Technology - **Wolfgang Bessler** [@wbessler](https://github.com/wbessler) - Offenburg University of Applied Science +- **Gaétan Bizot** [@Gaetanosaure] (https://github.com/Gaetanosaure) - University of Florence - **Paul Blum** [@paulblum](https://github.com/paulblum) - **Tilman Bremer** - **Victor Brunini** [@vbrunini](https://github.com/vbrunini) - Sandia National Laboratory diff --git a/include/cantera/kinetics/ElectronCollisionPlasmaRate.h b/include/cantera/kinetics/ElectronCollisionPlasmaRate.h index 1cddd48df41..7b8508f6c00 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; + + //! 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 (still accepted but deprecated) 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; @@ -240,11 +252,19 @@ 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 + 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/kinetics/TwoTempPlasmaRate.h b/include/cantera/kinetics/TwoTempPlasmaRate.h index 7f17e9de723..441da3bf33e 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 */ @@ -66,10 +68,13 @@ class TwoTempPlasmaRate : public ArrheniusBase /*! * @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 Temperature exponent (non-dimensional) + * @param b Electron 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, double EE, double bg); + 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/include/cantera/kinetics/VibrationalRelaxationRate.h b/include/cantera/kinetics/VibrationalRelaxationRate.h new file mode 100644 index 00000000000..2df4def6e8a --- /dev/null +++ b/include/cantera/kinetics/VibrationalRelaxationRate.h @@ -0,0 +1,293 @@ +//! @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. + + +//! It implements the relaxation rate of vibrationally excited species in plasma kinetics. +//! 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 '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). +//! 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. +//! 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 + +#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` + * - `multi-state-resolved` + * - `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: multi-state-resolved + * 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` + * - `multi-state-resolved` + * - `starikovskiy` + * - `castela` + */ + string m_vibration_model = "multi-state-resolved"; + + //! 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/include/cantera/thermo/EEDFTwoTermApproximation.h b/include/cantera/thermo/EEDFTwoTermApproximation.h index fddc47d5448..b22fa500584 100644 --- a/include/cantera/thermo/EEDFTwoTermApproximation.h +++ b/include/cantera/thermo/EEDFTwoTermApproximation.h @@ -61,28 +61,119 @@ 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. + //! 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); + + //! 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; } - double getElectronMobility() const { - return m_electronMobility; - } + //! 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); + + //! 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: @@ -192,7 +283,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(); @@ -209,13 +300,16 @@ 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) //! @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] @@ -236,7 +330,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) @@ -249,13 +343,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 @@ -270,10 +364,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 @@ -281,6 +379,43 @@ 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; + + // 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; + + //! 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); + + //! 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 } // end of namespace Cantera diff --git a/include/cantera/thermo/PlasmaPhase.h b/include/cantera/thermo/PlasmaPhase.h index 61bb85c9338..2017356e7b3 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. * @@ -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. @@ -186,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) @@ -337,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; @@ -368,10 +375,15 @@ 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; } - //! 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³] @@ -385,7 +397,18 @@ class PlasmaPhase: public IdealGasPhase //! Set reduced electric field given in [V·m²] void setReducedElectricField(double EN) { - m_electricField = EN * molarDensity() * Avogadro; // [V/m] + 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; } virtual void setSolution(std::weak_ptr soln) override; @@ -428,8 +451,36 @@ 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(); + + //! 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: void updateThermo() const override; @@ -438,7 +489,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(); @@ -452,7 +503,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. @@ -527,6 +578,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; @@ -560,6 +612,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 @@ -600,6 +657,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/interfaces/cython/cantera/lxcat2yaml.py b/interfaces/cython/cantera/lxcat2yaml.py index eda3b873da0..1f17a27cb34 100644 --- a/interfaces/cython/cantera/lxcat2yaml.py +++ b/interfaces/cython/cantera/lxcat2yaml.py @@ -1,380 +1,1157 @@ #!/usr/bin/env python3 +""" +lxcat2yaml.py -# 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. +Converts the latest version to date of LXCat XML files (08/06/2026) to a CanteraPlasma 4.0 YAML file. + +By default, the script writes only the top-level `electron-collisions` block. +If `--lxcat2phase converter.yaml` is provided, it also writes a `reactions` +block. The converter file maps LXCat product/state names to the species names +actually present in the target Cantera phase. Mapping values may also be full +product-side expressions, for example ``O + O`` when a complex product appears or +``0.876 N2 + 0.124 N2(v)`` should the user want to use the mean vibrational energy framework. + +If `--mech mechanism.yaml` is also provided, generated reactions are appended +to the existing top-level `reactions` section of the mechanism, and generated +`electron-collisions` are appended as the final top-level section. The merged +mechanism is then validated by loading it with Cantera. +Warning: the provided mechanism is expected to have been already pre-processed by the user +that should already have included the block related to the EEDF calculation in the YAML +and declared the phase as a PlasmaPhase. -""" -lxcat2yaml.py: Convert the LXCat integral cross-section data in XML format to YAML format. -The cross-section data is used to calculate the reaction rate of a electron-collision -process in a plasma. The data can be downloaded at https://nl.lxcat.net/data/download.php. Usage: - lxcat2yaml [--input=] - [--database=] - [--mech=] - [--phase=] - [--insert] - [--output=] - -Example: - lxcat2yaml --input=mycs.xml --database=itikawa --mech=oxygen-plasma.yaml - --phase=isotropic-electron-energy-plasma --insert - --output=oxygen-itikawa-plasma.yaml + lxcat2yaml --input test.xml --output collisions.yaml + lxcat2yaml --input test.xml --database Phelps \ + --lxcat2phase converter.yaml --output plasma_collisions.yaml + lxcat2yaml --input test.xml --database Phelps \ + --lxcat2phase converter.yaml --mech base_mechanism.yaml \ + --phase plasma --output base_mechanism_lxcat.yaml + +Example converter.yaml: + species: + N2(A3,v0-4): N2(A) + N2(A3,v5-9): N2(A) + N2(A3,v10-): N2(A) + N2(B3): N2(B) + N2(C3): N2(C) + O2(6.0eV): O + O + N2(v1): 0.8765957447 N2 + 0.1234042553 N2(v) (For mean vibrational energy frameworks). + +A flat mapping is also accepted: + N2(A3,v0-4): N2(A) + N2(A3,v5-9): N2(A) """ + from __future__ import annotations import argparse -import sys -import textwrap -import xml.etree.ElementTree as etree -from collections.abc import Iterable, Sequence +import re +import xml.etree.ElementTree as ET +from collections import defaultdict from pathlib import Path -from typing import TypeAlias, TypeVar, cast - -from ruamel import yaml -from ruamel.yaml.comments import CommentedMap, CommentedSeq -from ruamel.yaml.nodes import MappingNode -from ruamel.yaml.representer import SafeRepresenter +from typing import Any try: - import cantera as ct - OptionalSolutionType: TypeAlias = ct.Solution | None - Solution: type[ct.Solution] | None = ct.Solution -except ImportError: - print("The Cantera Python module was not found" - ", so the mechanism file cannot be used.") - Solution = None - -BlockMap: type[CommentedMap] = CommentedMap - -class Process: - """A class of YAML data for collision of a target species""" - def __init__(self, equation: str, energy_levels: list[float], cross_sections: list[float]) -> None: - self.equation = equation - self.energy_levels = energy_levels - self.cross_sections = cross_sections - - @classmethod - def to_yaml(cls, representer: SafeRepresenter, node: Process) -> MappingNode: - out = BlockMap([('equation', node.equation), - ('type', 'electron-collision-plasma'), - ('energy-levels', node.energy_levels), - ('cross-sections', node.cross_sections), - ]) - return representer.represent_dict(out) - -# Define YAML emitter -emitter: yaml.YAML = yaml.YAML() -emitter.register_class(Process) - -# Return indices of a child name -def get_children(parent: etree.Element[str], child_name: str) -> list[etree.Element[str]]: - return [child for child in parent if child.tag.find(child_name) != -1] - -_VT = TypeVar("_VT") # Value type. - -def Flowlist(*args: Iterable[_VT], **kwargs: _VT) -> list[_VT]: - """A YAML sequence that flows onto one line.""" - lst: CommentedSeq = CommentedSeq(*args, **kwargs) - lst.fa.set_flow_style() - return cast(list[_VT], lst) - -class IncorrectXMLNode(LookupError): - def __init__(self, message: str = "", node: etree.Element | None = None) -> None: - """Error raised when a required node is incorrect in the XML tree. - - :param message: - The error message to be displayed to the user. - :param node: - The XML node from which the requested node is incorrect. - """ - if node is not None: - # node str - node_str = etree.tostring(node, encoding="unicode") - - # Print the XML node - if message: - message += "\n" + node_str + from ruamel import yaml + from ruamel.yaml.comments import CommentedSeq +except ImportError as exc: + raise SystemExit( + "This script requires ruamel.yaml. Install it with: pip install ruamel.yaml" + ) from exc + + +# --------------------------------------------------------------------------- +# Small XML utilities +# --------------------------------------------------------------------------- + +def local_name(tag: str) -> str: + """Return the local XML tag name, stripping any namespace if present.""" + if "}" in tag: + return tag.split("}", 1)[1] + return tag + + +def child(parent: ET.Element, name: str) -> ET.Element | None: + """Return the first direct child whose local tag name is `name`.""" + for elem in parent: + if local_name(elem.tag) == name: + return elem + return None + + +def children(parent: ET.Element, name: str) -> list[ET.Element]: + """Return all direct children whose local tag name is `name`.""" + return [elem for elem in parent if local_name(elem.tag) == name] + + +def text_of(elem: ET.Element | None, default: str = "") -> str: + """Return stripped text from an XML element, or `default` if missing.""" + if elem is None or elem.text is None: + return default + return elem.text.strip() + + +def parse_float_list(text: str) -> list[float]: + """Parse a whitespace-separated list of floats.""" + return [float(x) for x in text.split()] + + +# --------------------------------------------------------------------------- +# YAML utilities +# --------------------------------------------------------------------------- + +def flow_list(values: list[float]) -> CommentedSeq: + """Force a YAML sequence to be emitted in compact flow style.""" + seq = CommentedSeq(values) + seq.fa.set_flow_style() + return seq + + +def make_yaml_emitter() -> yaml.YAML: + """Create the YAML emitter used for the output file.""" + emitter = yaml.YAML() + emitter.default_flow_style = False + emitter.width = 1000 + return emitter + + +# --------------------------------------------------------------------------- +# LXCat -> CanteraPlasma normalization +# --------------------------------------------------------------------------- + +ELECTRON_ALIASES = {"e", "E", "electron", "Electron"} +KINDS_WITHOUT_REACTIONS = {"effective", "elastic"} + + +def normalize_species_name(name: str) -> str: + """ + Normalize an LXCat species name toward a Cantera-compatible name. + + Examples: + e -> Electron + E -> Electron + O2^+ -> O2+ + O^- -> O- + N2^+(B2SIGMA) -> N2^+(B2SIGMA) + " O2 " -> O2 + + State-resolved charged species are intentionally kept in LXCat notation. + For example, ``N2^+(B2SIGMA)`` is not converted to ``N2+(B2SIGMA)`` so + that users can map it explicitly in the LXCat-to-phase converter file. + Simple ions without state labels are still converted to the usual + Cantera-like charge notation. + """ + name = name.strip() + + if name in ELECTRON_ALIASES: + return "Electron" + + # Remove stray internal whitespace, especially around charges. + name = re.sub(r"\s+", "", name) + + # Keep state-resolved charged species in LXCat notation, for example + # N2^+(B2SIGMA). The user can then provide an explicit mapping such as: + # N2^+(B2SIGMA): N2+ + # or any other phase-species name in the converter YAML file. + if re.search(r"\^[+-]\(", name): + return name + + # LXCat charge notation for simple ions O2^+ / O^- -> O2+ / O-. + name = name.replace("^+", "+") + name = name.replace("^-", "-") + + return name + + +def is_product_separator_plus(text: str, index: int) -> bool: + """ + Return True if `text[index]` is a plus sign separating two products. + + LXCat product strings may contain plus signs both as species separators and + as ionic charge markers. For example: + O^-+H2 -> separator between O^- and H2 + O+O -> separator between O and O + N2^+(B2SIGMA) -> charge marker, not a separator + N2+(B2SIGMA) -> charge marker, not a separator + O2+ -> final charge marker, not a separator + """ + if text[index] != "+": + return False + + previous_char = text[index - 1] if index > 0 else "" + next_char = text[index + 1] if index + 1 < len(text) else "" + + # Caret charge notation: N2^+(B2SIGMA), O^+, etc. + if previous_char == "^": + return False + + # Final plus sign: O2+, Ar+, etc. + if not next_char: + return False + + # Charge followed by an excited-state/state label: N2+(B2SIGMA). + if next_char == "(": + return False + + return True + + +def split_compact_products(product_text: str) -> list[str]: + """ + Split products that LXCat sometimes stores as a single compact string. + + Example: + "O^-+H2" -> ["O^-", "H2"] + "O+O" -> ["O", "O"] + + Note: + Positive ions with state labels such as "N2^+(B2SIGMA)" must remain a + single species and must not be split into "N2^" and "(B2SIGMA)". + """ + product_text = product_text.strip() + + if not product_text: + return [] + + parts: list[str] = [] + start = 0 + + for index, char in enumerate(product_text): + if char == "+" and is_product_separator_plus(product_text, index): + part = product_text[start:index].strip() + if part: + parts.append(part) + start = index + 1 + + final_part = product_text[start:].strip() + if final_part: + parts.append(final_part) + + return parts + + +def normalize_kind(process_type: str) -> str: + """Normalize LXCat process type names to lower-case CanteraPlasma kinds.""" + mapping = { + "Effective": "effective", + "Excitation": "excitation", + "Ionization": "ionization", + "Attachment": "attachment", + "Elastic": "elastic", + } + return mapping.get(process_type, process_type.strip().lower()) + + +def format_threshold(value: float) -> str: + """Return a stable, readable threshold representation for collision names.""" + return f"{value:.12g}" + + +def sanitize_for_name(text: str) -> str: + """ + Clean a string so that it can be used in a collision `name` field. + + Parentheses and + / - signs are intentionally preserved because they keep + names readable, for example: + N2(v1) + O+O + O2+ + """ + text = text.strip() + text = text.replace(" ", "") + text = text.replace("=>", "_") + text = text.replace("->", "_") + text = text.replace("/", "-") + text = text.replace(",", "") + text = text.replace(";", "") + text = text.replace(":", "") + return text + + +def make_collision_name( + database_id: str, + target: str, + kind: str, + product: str, + threshold: float, +) -> str: + """Create a deterministic collision name from LXCat process metadata.""" + db = sanitize_for_name(database_id.lower()) + target_s = sanitize_for_name(target) + product_s = sanitize_for_name(product) + thr_s = sanitize_for_name(format_threshold(threshold)) + + return f"{db}_{target_s}_{kind}_{product_s}_{thr_s}" + + +# --------------------------------------------------------------------------- +# Parsing one LXCat Process node +# --------------------------------------------------------------------------- + +def parse_threshold( + process: ET.Element, + kind: str, + energy_levels: list[float], + cross_sections: list[float], +) -> float: + """ + Determine the threshold / energy loss of an LXCat process. + + Rules: + 1. effective / elastic -> threshold = 0.0 + 2. if is strictly positive -> use it + 3. otherwise, if sigma starts at zero, use the last energy whose + sigma is zero before the first positive sigma value + 4. otherwise -> 0.0 + """ + if kind in KINDS_WITHOUT_REACTIONS: + return 0.0 + + parameters = child(process, "Parameters") + if parameters is not None: + e_node = child(parameters, "E") + if e_node is not None and e_node.text: + explicit_threshold = float(e_node.text.strip()) + + if explicit_threshold > 0.0: + return explicit_threshold + + if energy_levels and cross_sections and cross_sections[0] == 0.0: + last_zero_energy = energy_levels[0] + + for energy, sigma in zip(energy_levels, cross_sections): + if sigma == 0.0: + last_zero_energy = energy else: - message = "\n" + node_str + return last_zero_energy - super().__init__(message) + return 0.0 -def convert( - inpfile: str | Path | None = None, - database: str | None = None, - mechfile: str | None = None, - phase: str | None = None, - insert: bool | None = True, - outfile: str | Path | None = None, -) -> None: - """Convert an LXCat XML file to a YAML file. - - :param inpfile: - The input LXCat file name. - :param database: - The name of the database. For example, "itikawa". - :param mechfile: - The reaction mechanism file. This option requires using the Cantera library. - :param phase: - The phase name of the mechanism file. This option requires a ``mechfile`` to - also be specified. - :param insert: - The flag of whether to insert the collision reactions or not. - :param outfile: - The output YAML file name. - - All files are assumed to be relative to the current working directory of the Python - process running this script. - """ - if inpfile is not None: - inpfile = Path(inpfile) - lxcat_text = inpfile.read_text().lstrip() - if outfile is None: - outfile = inpfile.with_suffix(".yaml") - else: - raise ValueError("'inpfile' must be specified") - - if insert and mechfile is None: - raise ValueError("'mech' must be specified if 'insert' is used") - - gas: OptionalSolutionType = None - if mechfile is not None: - if Solution is None: - print("Cantera is not used, so the mechanism file cannot be used.") - sys.exit(1) - elif phase is not None: - gas = Solution(mechfile, phase, transport_model=None) + return 0.0 + + +def parse_species_block(process: ET.Element) -> tuple[list[str], list[str]]: + """ + Extract Reactant/Product entries from a block such as: + + + e + N2 + E + N2(v1) + + """ + species = child(process, "Species") + if species is None: + return [], [] + + reactants: list[str] = [] + products: list[str] = [] + + for node in children(species, "Reactant"): + value = text_of(node) + if value: + reactants.append(normalize_species_name(value)) + + for node in children(species, "Product"): + raw = text_of(node) + for part in split_compact_products(raw): + products.append(normalize_species_name(part)) + + return reactants, products + + +def choose_target(reactants: list[str]) -> str | None: + """Return the first non-electron reactant.""" + for species in reactants: + if species != "Electron": + return species + return None + + +def choose_product(target: str, products: list[str]) -> str: + """ + Return the collision product, excluding electrons. + + Examples: + [Electron, N2(v1)] -> N2(v1) + [Electron, Electron, O2+] -> O2+ + [CO, O-] -> CO + O- + [] -> target + """ + non_electron_products = [p for p in products if p != "Electron"] + + if not non_electron_products: + return target + + return " + ".join(non_electron_products) + + +def parse_process( + process: ET.Element, + database_id: str, + used_names: defaultdict[str, int], +) -> tuple[dict[str, Any], dict[str, Any]] | None: + """ + Convert one LXCat node into internal collision/reaction data. + + Return None if the process cannot be converted. + """ + data_x = child(process, "DataX") + data_y = child(process, "DataY") + + if data_x is None or data_y is None: + return None + + energy_levels = parse_float_list(text_of(data_x)) + cross_sections = parse_float_list(text_of(data_y)) + + if len(energy_levels) != len(cross_sections): + raise ValueError( + "DataX and DataY do not have the same length for process:\n" + + ET.tostring(process, encoding="unicode") + ) + + process_type = process.attrib.get("type", "unknown") + kind = normalize_kind(process_type) + + reactants, products = parse_species_block(process) + target = choose_target(reactants) + + if target is None: + return None + + product = choose_product(target, products) + threshold = parse_threshold( + process=process, + kind=kind, + energy_levels=energy_levels, + cross_sections=cross_sections, + ) + + # Same convention as the original Cantera script: if the first energy point + # is above the threshold, explicitly insert a zero-cross-section point at + # the threshold. + if energy_levels and energy_levels[0] > threshold: + energy_levels = [threshold, *energy_levels] + cross_sections = [0.0, *cross_sections] + elif cross_sections: + cross_sections[0] = 0.0 + + name = make_collision_name( + database_id=database_id, + target=target, + kind=kind, + product=product, + threshold=threshold, + ) + + # Ensure that the collision name is unique if duplicates occur. + used_names[name] += 1 + if used_names[name] > 1: + name = f"{name}_{used_names[name]}" + + collision: dict[str, Any] = { + "name": name, + "target": target, + "product": product, + "kind": kind, + "threshold": threshold, + "energy-levels": flow_list(energy_levels), + "cross-sections": flow_list(cross_sections), + } + + comment = text_of(child(process, "Comment")) + updated = text_of(child(process, "Updated")) + reaction = text_of(child(process, "Reaction")) + + # Useful metadata for checking the conversion. These fields can be removed + # later if a more minimal YAML output is desired. + if reaction: + collision["lxcat-reaction"] = reaction + if updated: + collision["updated"] = updated + if comment: + collision["comment"] = comment + + reaction_candidate: dict[str, Any] = { + "collision": name, + "kind": kind, + "reactants": reactants, + "products": products, + "target": target, + "product": product, + "lxcat-reaction": reaction, + } + + return collision, reaction_candidate + + +# --------------------------------------------------------------------------- +# LXCat-to-phase species mapping and reaction generation +# --------------------------------------------------------------------------- + +def split_reaction_expression(expression: str) -> list[str]: + """ + Split a product-side mapping expression into additive terms. + + The same plus-sign handling as LXCat product parsing is used, so ionic + charge markers are not treated as species separators. Examples: + O + O -> ["O", "O"] + O^-+H2 -> ["O^-", "H2"] + 0.876 N2 + 0.124 N2(v) -> ["0.876 N2", "0.124 N2(v)"] + N2^+(B2SIGMA) -> ["N2^+(B2SIGMA)"] + """ + return split_compact_products(expression) + + +def normalize_phase_expression(expression: str) -> str: + """ + Normalize a converter mapping value while preserving product expressions. + + Mapping values may be simple species names: + N2(A3,v0-4): N2(A3) + + or complete product-side expressions: + O2(6.0eV): O + O + N2(v1): 0.8765957447 N2 + 0.1234042553 N2(v) + + Each species term is normalized independently, but stoichiometric + coefficients and explicit product separators are preserved. + """ + expression = str(expression).strip() + if not expression: + return expression + + normalized_terms: list[str] = [] + + for term in split_reaction_expression(expression): + term = term.strip() + if not term: + continue + + coefficient_match = re.fullmatch( + r"([+-]?(?:\d+(?:\.\d*)?|\.\d+)(?:[eE][+-]?\d+)?)\s+(.+)", + term, + ) + + if coefficient_match: + coefficient = coefficient_match.group(1) + species = normalize_species_name(coefficient_match.group(2)) + normalized_terms.append(f"{coefficient} {species}") else: - gas = Solution(mechfile, transport_model=None) + normalized_terms.append(normalize_species_name(term)) + + return " + ".join(normalized_terms) + + +def load_lxcat2phase_mapping(path: Path) -> dict[str, str]: + """ + Load a mapping from LXCat species names to phase species names or product + expressions. + + Accepted formats: + species: + N2(A3,v0-4): N2(A) + O2(6.0eV): O + O + N2(v1): 0.8765957447 N2 + 0.1234042553 N2(v) + + or directly: + N2(A3,v0-4): N2(A) + O2(6.0eV): O + O + """ + loader = yaml.YAML(typ="safe") + with path.open("r", encoding="utf-8") as stream: + data = loader.load(stream) or {} + + if not isinstance(data, dict): + raise ValueError("The lxcat2phase converter YAML must contain a mapping.") + + if "species" in data: + data = data["species"] or {} + + if not isinstance(data, dict): + raise ValueError("The `species` entry in the converter YAML must be a mapping.") + + mapping: dict[str, str] = {} + for lxcat_name, phase_expression in data.items(): + if phase_expression is None: + continue + lxcat_key = normalize_species_name(str(lxcat_name)) + mapping[lxcat_key] = normalize_phase_expression(str(phase_expression)) + + return mapping + + +def is_simple_phase_species(species: str) -> bool: + """ + Return True for species names that can reasonably map to themselves. + + This covers ground-state atoms/molecules and simple ions such as O, O2, + CO2, H2O, O-, O2+, Ar+, etc. Excited-state labels like N2(A3,v0-4), + pseudo-products such as 6.3eVloss, or star labels require an explicit + lxcat2phase mapping. + """ + if species == "Electron": + return True + + return bool(re.fullmatch(r"[A-Z][A-Za-z0-9]*[+-]?", species)) + + +def project_species_to_phase( + species: str, + lxcat2phase: dict[str, str], +) -> tuple[str | None, bool]: + """ + Project one LXCat species name to a phase species name or product-side + expression. + + Return `(projected_expression, missing_mapping)`. If `missing_mapping` is + True, no usable phase species name could be determined. + """ + if species in lxcat2phase: + return lxcat2phase[species], False + + if is_simple_phase_species(species): + return species, False + + return None, True + + +def format_species_side(species: list[str]) -> str: + """ + Format a reaction side from a list of species names or expressions. + + Expressions returned by the converter mapping may already contain product + separators, for example ``O + O`` or ``0.876 N2 + 0.124 N2(v)``. Joining + them with `` + `` preserves the intended full reaction side. + """ + return " + ".join(species) + + +def mark_duplicate_reactions(reactions: list[dict[str, Any]]) -> int: + """ + Add ``duplicate: true`` to reactions that share the same equation. + + Cantera requires duplicate reactions to be explicitly marked. This pass is + applied after all LXCat species have been projected to phase species, since + two distinct LXCat processes may collapse to the same phase-level equation. + + Existing reactions from a mechanism file are also supported. Entries that + do not have an ``equation`` field are ignored. + + Return the number of reactions marked as duplicates. + """ + equation_counts: defaultdict[str, int] = defaultdict(int) + for reaction in reactions: + equation = reaction.get("equation") + if equation: + equation_counts[equation] += 1 + + duplicate_count = 0 + for reaction in reactions: + equation = reaction.get("equation") + if equation and equation_counts[equation] > 1: + reaction["duplicate"] = True + duplicate_count += 1 + + return duplicate_count + + +def build_reactions( + reaction_candidates: list[dict[str, Any]], + lxcat2phase: dict[str, str], +) -> tuple[list[dict[str, Any]], dict[str, list[dict[str, str]]], int, int]: + """ + Build CanteraPlasma reaction entries from parsed LXCat process metadata. + + Reactions are generated only for processes whose kind is not effective or + elastic. If a non-trivial LXCat species cannot be mapped to the phase, the + reaction is skipped and the missing species is recorded. Reactions that end + up with identical phase-level equations are marked with ``duplicate: true``. + """ + reactions: list[dict[str, Any]] = [] + missing: dict[str, list[dict[str, str]]] = defaultdict(list) + skipped = 0 + + for candidate in reaction_candidates: + kind = candidate["kind"] + if kind in KINDS_WITHOUT_REACTIONS: + continue + + projected_reactants: list[str] = [] + projected_products: list[str] = [] + candidate_missing = False + + for species in candidate["reactants"]: + projected, is_missing = project_species_to_phase(species, lxcat2phase) + if is_missing or projected is None: + candidate_missing = True + missing[species].append( + { + "collision": candidate["collision"], + "lxcat-reaction": candidate.get("lxcat-reaction", ""), + } + ) + else: + projected_reactants.append(projected) + + for species in candidate["products"]: + projected, is_missing = project_species_to_phase(species, lxcat2phase) + if is_missing or projected is None: + candidate_missing = True + missing[species].append( + { + "collision": candidate["collision"], + "lxcat-reaction": candidate.get("lxcat-reaction", ""), + } + ) + else: + projected_products.append(projected) + + if candidate_missing: + skipped += 1 + continue + + # Some LXCat processes may omit products. In that case, keep the target + # unchanged while preserving the electron, matching the original intent + # of a non-reactive energy-loss channel. + if not projected_products: + projected_products = list(projected_reactants) + + equation = f"{format_species_side(projected_reactants)} => {format_species_side(projected_products)}" + + reactions.append( + { + "equation": equation, + "type": "electron-collision-plasma", + "collision": candidate["collision"], + } + ) + + duplicate_count = mark_duplicate_reactions(reactions) + + return reactions, dict(missing), skipped, duplicate_count + + +def print_missing_mapping_report( + missing: dict[str, list[dict[str, str]]], + skipped: int, +) -> None: + """Print a detailed terminal report for missing lxcat2phase mappings.""" + if not missing: + print("All non-trivial LXCat species used in reactions were mapped successfully.") + return + + total_occurrences = sum(len(items) for items in missing.values()) + print() + print("Missing LXCat-to-phase species mappings") + print("---------------------------------------") + print(f"Skipped reaction candidates: {skipped}") + print(f"Missing species: {len(missing)}") + print(f"Missing occurrences: {total_occurrences}") + print() + print("Add entries such as the following to your converter YAML:") + print("species:") + for species in sorted(missing): + print(f" {species}: ") + print() + print("Detailed occurrences:") + for species in sorted(missing): + print(f"- {species}") + for item in missing[species]: + reaction = item.get("lxcat-reaction") or "" + print(f" collision: {item['collision']}") + print(f" lxcat-reaction: {reaction}") + + +# --------------------------------------------------------------------------- +# Full-file conversion +# --------------------------------------------------------------------------- + +def database_matches(database_node: ET.Element, requested: str | None) -> bool: + """Return True if the XML database node matches the requested database.""" + if requested is None: + return True + + db_id = database_node.attrib.get("id", "") + db_name = database_node.attrib.get("name", "") + + return requested in {db_id, db_name} + + +def parse_lxcat_xml( + input_path: Path, + requested_database: str | None = None, +) -> tuple[list[dict[str, Any]], list[dict[str, Any]]]: + """Parse the LXCat XML file and return collisions plus reaction candidates.""" + xml_text = input_path.read_text(encoding="utf-8").lstrip() + root = ET.fromstring(xml_text) - elif phase is not None: - raise ValueError("'mech' must be specified if 'phase' is used.") + collisions: list[dict[str, Any]] = [] + reaction_candidates: list[dict[str, Any]] = [] + used_names: defaultdict[str, int] = defaultdict(int) - xml_tree = etree.fromstring(lxcat_text) + databases = [node for node in root if local_name(node.tag) == "Database"] - # If insert key word is used, create a process list, - # and append all processes together - process_list: list[Process] | None = None - if not insert: - process_list = [] + if not databases: + raise ValueError("No node found in the XML file.") - for database_node in xml_tree: - if database is not None: - if database_node.attrib["id"] != database: + for database_node in databases: + if not database_matches(database_node, requested_database): + continue + + database_id = ( + database_node.attrib.get("id") + or database_node.attrib.get("name") + or "lxcat" + ) + + groups_node = child(database_node, "Groups") + if groups_node is None: + continue + + for group_node in children(groups_node, "Group"): + processes_node = child(group_node, "Processes") + if processes_node is None: continue - # Get groups node - groups_node = get_children(database_node, "groups")[0] + for process_node in children(processes_node, "Process"): + parsed = parse_process( + process=process_node, + database_id=database_id, + used_names=used_names, + ) + if parsed is not None: + collision, reaction_candidate = parsed + collisions.append(collision) + reaction_candidates.append(reaction_candidate) - for group in groups_node: - for process in get_children(group, "processes")[0]: - registerProcess(process, process_list, gas) + return collisions, reaction_candidates - if not insert: - # Put process list in collision node - collision_node = {"collisions": process_list} - with Path(outfile).open("w") as output_file: - emitter.dump(collision_node, output_file) - else: - # Get mechanism file unit system - units = None - assert mechfile is not None - with open(mechfile, "r") as mech: - data = yaml.YAML(typ="rt").load(mech) - if "units" in data: - units = data["units"] - assert gas is not None - gas.write_yaml(outfile, units=units) - -def registerProcess(process: etree.Element, - process_list: list[Process] | None, - gas: OptionalSolutionType) -> None: - """ - Add a collision process (electron collision reaction) to process_list - and gas object if it exists. - - :param process: - The collision process (electron collision reaction) - :param process_list: - The list of collision processes - :param gas: - The Cantera Solution object - """ - # Get electron specie name - electron_name = gas.electron_species_name if gas is not None else "e" - - # Parse the threshold - threshold = 0.0 - parameters_node = get_children(process, "parameters")[0] - if len(get_children(parameters_node, "parameter")) == 1: - parameter = get_children(parameters_node, "parameter")[0] - if parameter.attrib["name"] == 'E': - assert parameter.text is not None - threshold = float(parameter.text) - - # Parse the equation - product_array: list[str] = [] - - products: list[etree.Element[str]] = get_children(process, "products") - if products: - for product_node in products[0]: - if product_node.tag.find("electron") != -1: - product_array.append(electron_name) - - if product_node.tag.find("molecule") != -1: - assert product_node.text is not None - product_name = product_node.text - if "state" in product_node.attrib: - state = product_node.attrib["state"].replace(" ","-") - # State is appended in a parenthesis - product_name += f"({state})" - if "charge" in product_node.attrib: - charge = int(product_node.attrib["charge"]) - if charge > 0: - product_name += charge*"+" - else: - product_name += -charge*"-" - - # Filter the collision based on the existed species in the mechanism file - if gas is not None and not product_name in gas.species_names: - return - - product_array.append(product_name) - - for reactant_node in get_children(process, "reactants")[0]: - if reactant_node.tag.find("molecule") != -1: - reactant = reactant_node.text - # Filter the collision based on the existed species in the mechanism file - if gas is not None and not reactant in gas.species_names: - return - - if product_array: # not empty - products_string = " + ".join(product_array) + +def ensure_top_level_sequence(data: dict[str, Any], key: str) -> list[Any]: + """Return a top-level YAML sequence, creating it if necessary.""" + if key not in data or data[key] is None: + data[key] = CommentedSeq() + + if not isinstance(data[key], list): + raise ValueError(f"Top-level `{key}` entry must be a YAML sequence.") + + return data[key] + + +def move_key_to_end(data: dict[str, Any], key: str) -> None: + """Move a top-level key to the end of the YAML mapping if it exists.""" + if key not in data: + return + + value = data.pop(key) + data[key] = value + + +def load_mechanism_yaml(mech_path: Path) -> dict[str, Any]: + """Load an existing Cantera YAML mechanism in round-trip mode.""" + loader = yaml.YAML(typ="rt") + with mech_path.open("r", encoding="utf-8") as stream: + data = loader.load(stream) + + if data is None: + data = {} + + if not isinstance(data, dict): + raise ValueError("The mechanism YAML must contain a top-level mapping.") + + return data + + +def write_merged_mechanism_yaml( + mech_path: Path, + output_path: Path, + reactions: list[dict[str, Any]], + collisions: list[dict[str, Any]], +) -> int: + """ + Append generated reactions/collisions to an existing mechanism YAML file. + + The existing top-level ``reactions`` section is preserved and extended. The + top-level ``electron-collisions`` section is created or extended, then moved + to the end of the file. Duplicate reaction equations are marked after the + merge, so duplicates between existing and generated reactions are handled. + + Return the number of reactions marked with ``duplicate: true`` in the final + merged mechanism. + """ + data = load_mechanism_yaml(mech_path) + + mechanism_reactions = ensure_top_level_sequence(data, "reactions") + mechanism_reactions.extend(reactions) + + mechanism_collisions = ensure_top_level_sequence(data, "electron-collisions") + mechanism_collisions.extend(collisions) + move_key_to_end(data, "electron-collisions") + + duplicate_count = mark_duplicate_reactions(mechanism_reactions) + + emitter = make_yaml_emitter() + with output_path.open("w", encoding="utf-8") as stream: + emitter.dump(data, stream) + + return duplicate_count + + +def validate_cantera_mechanism(output_path: Path, phase: str | None = None) -> None: + """Validate a generated mechanism by loading it with Cantera.""" + try: + import cantera as ct + except ImportError as exc: + raise SystemExit( + "Cantera validation was requested because --mech was provided, " + "but the Cantera Python module could not be imported." + ) from exc + + try: + if phase: + ct.Solution(str(output_path), phase, transport_model=None) + else: + ct.Solution(str(output_path), transport_model=None) + except Exception as exc: + raise SystemExit( + f"Cantera validation FAILED for {output_path}.\n{exc}" + ) from exc + + print(f"Cantera validation PASSED for: {output_path}") + + +def write_yaml( + collisions: list[dict[str, Any]], + output_path: Path, + reactions: list[dict[str, Any]] | None = None, +) -> None: + """Write the converted data to a YAML file.""" + emitter = make_yaml_emitter() + + if reactions is None: + data = { + "electron-collisions": collisions, + } else: - # No product is identified. Use the reactant as the product. - products_string = f"{reactant} + {electron_name}" + data = { + "reactions": reactions, + "electron-collisions": collisions, + } - equation = f"{reactant} + {electron_name} => {products_string}" + with output_path.open("w", encoding="utf-8") as stream: + emitter.dump(data, stream) - # Parse the cross-section data - data_x_node = get_children(process, "data_x")[0] - if data_x_node is None: - raise IncorrectXMLNode("The 'process' node requires the 'data_x' node.", process) - data_y_node = get_children(process, "data_y")[0] - if data_y_node is None: - raise IncorrectXMLNode("The 'process' node requires the 'data_y' node.", process) - assert data_x_node.text is not None - assert data_y_node.text is not None - energy_levels = Flowlist(map(float, data_x_node.text.split())) - cross_sections = Flowlist(map(float, data_y_node.text.split())) +# --------------------------------------------------------------------------- +# Public conversion API +# --------------------------------------------------------------------------- - # Edit energy levels and cross section - if len(energy_levels) != len(cross_sections): - raise IncorrectXMLNode("Energy levels (data_x) and cross section " - "(data_y) must have the same length.", process) +def convert( + inpfile: str | Path | None = None, + database: str | None = None, + lxcat2phase: str | Path | None = None, + mechfile: str | Path | None = None, + phase: str | None = None, + outfile: str | Path | None = None, + *, + report: bool = False, + validate: bool = True, +) -> Path: + """ + Convert a modern LXCat XML file to CanteraPlasma 4.0 YAML. + + Parameters + ---------- + inpfile + Input LXCat XML file. + database + Optional LXCat database id or name used as a filter. + lxcat2phase + Optional YAML mapping from LXCat species names to phase species names + or product-side expressions. If provided, a ``reactions`` section is + generated. If omitted, only ``electron-collisions`` are written. + mechfile + Optional existing Cantera YAML mechanism. If provided, generated + reactions are appended to its top-level ``reactions`` section and + generated ``electron-collisions`` are appended at the end of the file. + This option requires ``lxcat2phase``. + phase + Optional phase name used when validating a merged mechanism with + Cantera. This option is only meaningful with ``mechfile``. + outfile + Output YAML file. If omitted, the input XML suffix is replaced by + ``.yaml``. With ``mechfile``, the default is + ``_lxcat.yaml`` in the same directory as the mechanism. + report + If True, print the same status and missing-mapping messages as the + command-line interface. + validate + If True, validate merged mechanisms by loading the output with Cantera. + + Returns + ------- + pathlib.Path + Path to the generated YAML file. + """ + if inpfile is None: + raise ValueError("'inpfile' must be specified.") - if energy_levels[0] > threshold: - # Use Flowlist again to ensure correct YAML format - energy_levels = Flowlist([threshold, *energy_levels]) - cross_sections = Flowlist([0.0, *cross_sections]) + input_path = Path(inpfile) + mech_path = Path(mechfile) if mechfile is not None else None + + if mech_path is not None and lxcat2phase is None: + raise ValueError("'lxcat2phase' must be specified when 'mechfile' is used.") + + if phase is not None and mech_path is None: + raise ValueError("'mechfile' must be specified when 'phase' is used.") + + if outfile is not None: + output_path = Path(outfile) + elif mech_path is not None: + output_path = mech_path.with_name(f"{mech_path.stem}_lxcat.yaml") else: - cross_sections[0] = 0.0 + output_path = input_path.with_suffix(".yaml") + + collisions, reaction_candidates = parse_lxcat_xml( + input_path=input_path, + requested_database=database, + ) + + if lxcat2phase is None: + write_yaml(collisions, output_path) + if report: + print(f"Wrote {len(collisions)} electron-collisions to: {output_path}") + print("No reactions were generated because --lxcat2phase was not provided.") + return output_path + + converter_path = Path(lxcat2phase) + lxcat2phase_mapping = load_lxcat2phase_mapping(converter_path) + reactions, missing, skipped, duplicate_count = build_reactions( + reaction_candidates, + lxcat2phase_mapping, + ) - # If insert mode is on, add the process as a reaction to the gas object. - if gas is not None: - R = ct.Reaction( - equation=equation, - rate=ct.ElectronCollisionPlasmaRate(energy_levels=energy_levels, - cross_sections=cross_sections)) - gas.add_reaction(R) + if mech_path is not None: + final_duplicate_count = write_merged_mechanism_yaml( + mech_path=mech_path, + output_path=output_path, + reactions=reactions, + collisions=collisions, + ) + if report: + print(f"Appended {len(reactions)} generated reactions to mechanism: {output_path}") + print(f"Appended {len(collisions)} electron-collisions to mechanism: {output_path}") + print( + f"Marked {final_duplicate_count} duplicate reactions with " + "duplicate: true in the final merged mechanism." + ) + print_missing_mapping_report(missing, skipped) + if validate: + validate_cantera_mechanism(output_path, phase) + return output_path - # If insert mode is off, process_list is used to store the data. - if process_list is not None: - process_list.append(Process(equation=equation, - energy_levels=energy_levels, - cross_sections=cross_sections)) + write_yaml(collisions, output_path, reactions=reactions) + if report: + print(f"Wrote {len(reactions)} reactions to: {output_path}") + print(f"Wrote {len(collisions)} electron-collisions to: {output_path}") + print(f"Marked {duplicate_count} duplicate reactions with duplicate: true.") + print_missing_mapping_report(missing, skipped) + return output_path + +# --------------------------------------------------------------------------- +# Command-line interface +# --------------------------------------------------------------------------- def create_argparser() -> argparse.ArgumentParser: + """Create the command-line argument parser.""" parser = argparse.ArgumentParser( description=( - "Convert the LXCat integral cross-section data in XML format (LXCATML) to " - "YAML format"), - epilog=textwrap.dedent( - """ - Example:: - - lxcat2yaml --input=mycs.xml --database=itikawa --mech=oxygen-plasma.yaml - --phase=isotropic-electron-energy-plasma --insert - --output=oxygen-itikawa-plasma.yaml - - If the **lxcat2yaml** script is not on your path but the Cantera Python - module is, **lxcat2yaml** can also be invoked by running:: - - python -m cantera.lxcat2yaml --input=mycs.xml - - In both cases, the equal signs in the options are optional. In the - second case, the xml file is converted to yaml without inserting the - collision reactions into the mechanism file. - """), - formatter_class=argparse.RawDescriptionHelpFormatter + "Convert a modern LXCat XML file to a CanteraPlasma 4.0 YAML file. " + "By default only `electron-collisions` are written. If --lxcat2phase " + "is provided, a `reactions` section is generated as well. If --mech " + "is provided, the generated entries are merged into that mechanism and " + "the final YAML is validated with Cantera." + ) ) + parser.add_argument( - "--input", default=None, - help=("LXCat electron-collision cross sections input file, containing " - "a list of the electron-collision plasma reactions with the electron " - "energy levels and corresponding cross sections. Must be specified")) + "--input", + required=True, + help="Input LXCat XML file.", + ) + parser.add_argument( - "--database", default=None, - help=("The name of the database. Optional. Use it when multiple databases " - "exist in the input file.")) + "--output", + default=None, + help=( + "Output YAML file. Default: input filename with a .yaml extension, " + "or _lxcat.yaml when --mech is provided." + ), + ) + parser.add_argument( - "--mech", default=None, - help=("Cantera yaml-format reaction mechanism file. The list of the species is " - "used as the filter to determine which electron-collision reactions in " - "the input file are parsed. In addition, the electron-collision reactions " - "can be inserted automatically into the mechanism file with the argument " - "--insert and become the output file.")) + "--database", + default=None, + help="Optional filter on the LXCat database id or name, for example Phelps.", + ) + parser.add_argument( - "--phase", default=None, - help=("This specifies the name of the phase in the mechanism file. Optional.")) + "--lxcat2phase", + default=None, + help=( + "Optional YAML mapping from LXCat species names to phase species names or product-side expressions. " + "When provided, the script also writes a `reactions` section." + ), + ) + parser.add_argument( - "--insert", action="store_true", default=False, - help=("Enable inserting the electron-collision reactions into the mechanism file." - "Need to use with the argument --mech to provide the mechanism file" - "Optional.")) + "--mech", + default=None, + help=( + "Optional existing Cantera YAML mechanism. When provided, generated " + "reactions are appended to its top-level `reactions` section and " + "generated `electron-collisions` are appended at the end. This option " + "requires --lxcat2phase and triggers Cantera validation of the output." + ), + ) + parser.add_argument( - "--output", default=None, - help=("Specifies the OUTPUT file name. By default, the output file name is the " - "input file name with the extension changed to **.yaml**.")) + "--phase", + default=None, + help=( + "Optional phase name used when validating the merged mechanism with " + "Cantera. If omitted, Cantera's default phase loading behavior is used." + ), + ) return parser -def main(argv: Sequence[str] | None = None) -> None: - """Parse command line arguments and pass them to `convert`.""" + +def main() -> None: + """Run the command-line converter.""" parser = create_argparser() - if argv is None and len(sys.argv) < 2: - parser.print_help(sys.stderr) - sys.exit(1) - args = parser.parse_args(argv) - - input_file = Path(args.input) - - output_file: Path | str = args.output or input_file.with_suffix(".yaml") - convert(input_file, args.database, args.mech, args.phase, args.insert, output_file) - - if args.insert and Solution is not None: - # Test mechanism can be loaded back into Cantera - try: - print("Validating mechanism...", end="") - Solution(output_file, args.phase, transport_model=None) - print("PASSED.") - except RuntimeError as e: - print("FAILED.") - print(e) - sys.exit(1) + args = parser.parse_args() + + try: + convert( + inpfile=args.input, + database=args.database, + lxcat2phase=args.lxcat2phase, + mechfile=args.mech, + phase=args.phase, + outfile=args.output, + report=True, + validate=True, + ) + except ValueError as exc: + parser.error(str(exc)) + + if __name__ == "__main__": main() diff --git a/interfaces/cython/cantera/thermo.pxd b/interfaces/cython/cantera/thermo.pxd index ef8be3b98e1..3346f934524 100644 --- a/interfaces/cython/cantera/thermo.pxd +++ b/interfaces/cython/cantera/thermo.pxd @@ -215,6 +215,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 b13a608f544..eba60dc141c 100644 --- a/interfaces/cython/cantera/thermo.pyx +++ b/interfaces/cython/cantera/thermo.pyx @@ -2007,6 +2007,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): diff --git a/src/kinetics/ElectronCollisionPlasmaRate.cpp b/src/kinetics/ElectronCollisionPlasmaRate.cpp index ce4a1ecbcd6..a21f4fd5e0f 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()); @@ -45,35 +45,56 @@ 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) { + 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("target")) { - m_target = node["target"].asString(); + if (node.hasKey("collision")) { + m_collisionName = node["collision"].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); } -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( @@ -142,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 @@ -159,19 +180,20 @@ 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) { + const ThermoPhase& thermo = kin.thermo(); // get electron species name string electronName; @@ -215,16 +237,7 @@ 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; - } - } - } + setDefaultThreshold(); if (!rxn.reversible) { return; // end checking of forward reaction @@ -245,4 +258,118 @@ 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; + + 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(); + } + + 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); + + setDefaultThreshold(); + + 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/kinetics/ReactionRateFactory.cpp b/src/kinetics/ReactionRateFactory.cpp index 1019fe4f4ed..cb079dc0ee1 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/VibrationalRelaxationRate.h" namespace Cantera { @@ -40,6 +41,11 @@ ReactionRateFactory::ReactionRateFactory() return new TwoTempPlasmaRate(node, rate_units); }); + // VibrationalRelaxationRate evaluator + reg("vibrational-relaxation", [](const AnyMap& node, const UnitStack& rate_units) { + return new VibrationalRelaxationRate(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/src/kinetics/TwoTempPlasmaRate.cpp b/src/kinetics/TwoTempPlasmaRate.cpp index ce3f4694969..0b6563e98e3 100644 --- a/src/kinetics/TwoTempPlasmaRate.cpp +++ b/src/kinetics/TwoTempPlasmaRate.cpp @@ -57,19 +57,34 @@ 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) +{ + TwoTempPlasmaRate(A, b, Ea, EE); + 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/kinetics/VibrationalRelaxationRate.cpp b/src/kinetics/VibrationalRelaxationRate.cpp new file mode 100644 index 00000000000..79aabf911ee --- /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) -> multi-state-resolved +// 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: multi-state-resolved + // vibration_model: starikovskiy + // vibration_model: castela + // + // 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(); + } 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 == "multi-state-resolved") { + // 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 'multi-state-resolved', " + "'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 == "multi-state-resolved") { + 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 'multi-state-resolved', " + "'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 == "multi-state-resolved") { + 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/src/thermo/EEDFTwoTermApproximation.cpp b/src/thermo/EEDFTwoTermApproximation.cpp index 7fa38c916e5..1961fa2a964 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; @@ -43,6 +43,96 @@ 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, double ratio) +{ + 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 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) { @@ -51,33 +141,133 @@ int EEDFTwoTermApproximation::calculateDistributionFunction() } updateMoleFractions(); + checkSpeciesNoCrossSection(); 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); + 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", + "Invalid kTe value for Maxwellian first guess."); + } + + const double fFloor = 1e-300; + + 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)); } - } else { + } + + double fnorm = norm(m_f0, m_gridCenter); + + if (!std::isfinite(fnorm) || fnorm <= 0.0) { throw CanteraError("EEDFTwoTermApproximation::calculateDistributionFunction", - " unknown EEDF first guess"); + "Invalid norm for Maxwellian initialization."); } - } - converge(m_f0); + m_f0 /= fnorm; + }; + + // At very low reduced electric field, force a Maxwellian at the gas + // 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 { + // 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", + "Unknown EEDF first guess '{}'.", m_firstguess); + } + } + + 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; + + 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") { + setMaxwellian(m_init_kTe); + } else { + throw CanteraError("EEDFTwoTermApproximation::calculateDistributionFunction", + "Unknown EEDF first guess '{}'.", m_firstguess); + } + + updateCrossSections(); + converge(m_f0); - // write the EEDF at grid edges + } 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") { + setMaxwellian(m_init_kTe); + } else { + throw CanteraError("EEDFTwoTermApproximation::calculateDistributionFunction", + "Unknown EEDF first guess '{}'.", m_firstguess); + } + + 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()); + 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 - m_electronMobility = electronMobility(m_f0); + // Update electron mobility. + m_electronMobility = computeElectronMobility(m_f0); + return 0; } @@ -115,7 +305,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"); } } } @@ -132,7 +322,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); @@ -278,7 +468,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 { @@ -342,7 +532,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++) { @@ -394,7 +584,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); @@ -407,7 +597,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() @@ -464,6 +656,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; @@ -478,38 +671,144 @@ void EEDFTwoTermApproximation::updateMoleFractions() } } +// 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); 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] -= linearInterpCrossSectionZeroOutside(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 { + 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 {} kind {} target {} at energy {}: {}\n", + k, effectiveKind, effectiveTarget, 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); + } + } +} + +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. + 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)); } } } @@ -596,4 +895,111 @@ 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); + } + + // put back the flag at false because since the grid has changed the EEDF must be computed again. + 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; +} + +// 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 26f83aeb469..3c72d5bcd00 100644 --- a/src/thermo/PlasmaPhase.cpp +++ b/src/thermo/PlasmaPhase.cpp @@ -19,25 +19,56 @@ 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_) { - 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() @@ -60,10 +91,28 @@ 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(); 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."); @@ -78,6 +127,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 { @@ -88,16 +140,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) @@ -127,7 +185,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(); } @@ -155,7 +219,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(); } @@ -179,49 +249,102 @@ 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); } } } 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(); @@ -245,15 +368,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(); } @@ -318,32 +448,242 @@ 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 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 initial grid. + 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 initialMaxEnergy = eedf["initial_max_energy_level"].asDouble(); + size_t nGridCells = static_cast( + eedf["initial_number_of_energy_grid_cells"].asInt()); + + if (!std::isfinite(initialMaxEnergy) || initialMaxEnergy <= 0.0) { + throw CanteraError("PlasmaPhase::setParameters", + "initial_max_energy_level must be finite and 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"); + } + + m_eedfSolver->setGridType(energyLevelsDistribution); + m_eedfSolver->setInitialGridParameters(initialMaxEnergy, nGridCells); + + if (energyLevelsDistribution == "Linear") { + m_eedfSolver->setLinearGrid(initialMaxEnergy, nGridCells); + } 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."); + } + + 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(); + + 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. + m_nPoints = nGridCells + 1; + } + + m_electronEnergyLevels = Eigen::Map( + m_eedfSolver->getGridEdge().data(), m_eedfSolver->getGridEdge().size()); + + m_nPoints = m_eedfSolver->getGridEdge().size(); + + m_electronEnergyDist.resize(m_nPoints); + m_electronEnergyDist.setZero(); + + checkElectronEnergyLevels(); + electronEnergyLevelChanged(); } + + } 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; - } else { - products[item["target"].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; } - products[electronSpeciesName()] = 1; - if (rate->kind() == "ionization") { - products[electronSpeciesName()] += 1; - } else if (rate->kind() == "attachment") { - products[electronSpeciesName()] -= 1; + } + } + + // 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(); + + if (!m_electronCollisionDefinitions.count(name)) { + throw InputFileError("setParameters", rootNode, + "Reaction references unknown electron collision '{}'.", name); + } + m_referencedElectronCollisions.insert(name); } - auto R = make_shared(reactants, products, rate); - addCollision(R); } } + + 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)) { + continue; + } + + } + bool hasName = item.hasKey("name"); + + if (hasName) { + string name = item["name"].asString(); + + // New YAML format: if the collision is referenced by a reaction, do not create a synthetic reaction here. + if (m_referencedElectronCollisions.count(name)) { + continue; + } + } + + // Old YAML format for anonymous collisions or unreferenced named collisions: they remain as standalone EEDF collisions. + 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); @@ -361,6 +701,11 @@ bool PlasmaPhase::addSpecies(shared_ptr spec) "Only one electron species is allowed.", spec->name); } } + + if (added) { + m_vibrationalReservoirSpeciesNeedUpdate = true; + } + return added; } @@ -422,6 +767,7 @@ void PlasmaPhase::setCollisions() void PlasmaPhase::addCollision(shared_ptr collision) { + size_t i = nCollisions(); // setup callback to signal updating the cross-section-related @@ -447,6 +793,77 @@ void PlasmaPhase::addCollision(shared_ptr collision) " collision with equation '{}'", collision->equation()); } + // management of the new data file format + + auto ratePtr = std::dynamic_pointer_cast(collision->rate()); + + 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()); + } + + ratePtr->applyCollisionData(it->second); + } + + 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); + } + + string kindFromReaction = inferElectronCollisionKind(collision); + string kindFromCollision = ratePtr->kind(); + + bool compatibleKind = kindFromReaction == kindFromCollision; + + if ((kindFromReaction == "elastic" || kindFromReaction == "effective") && + (kindFromCollision == "elastic" || kindFromCollision == "effective")) { + compatibleKind = true; + } + + // Allow collapsed inelastic channels. + // Example: + // reactions: + // equation: Electron + N2 => Electron + N2 + // collision: phelps_N2_excitation_N2(rot)_0.02 + // + // electron-collisions: + // name: phelps_N2_excitation_N2(rot)_0.02 + // target: N2 + // product: N2(rot) + // energy-levels: [0.02, 0.03, 0.4, 0.8, 1.2, 1.6, 1.7, 1.8, 1.9, 2.0, 2.1, 2.2, 2.3, + // 2.4, 2.5, 2.6, 2.7, 2.8, 2.9, 3.0, 3.1, 3.2, 3.3, 3.6, 5.0] + // cross-sections: [0.0, 2.5e-22, 2.5e-22, 2.5e-22, 4.7e-22, 8.6e-22, 1.5e-21, 2.35e-21, + // 1.08e-20, 1.9e-20, 2.03e-20, 2.77e-20, 2.5e-20, 2.19e-20, 2.4e-20, 2.17e-20, + // 1.62e-20, 1.38e-20, 1.18e-20, 1.03e-20, 8.4e-21, 6.9e-21, 5e-21, 1.7e-21, 0.0] + // kind: excitation, + // threshold: 0.02 + + 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 '{}'.", + ratePtr->collisionName(), kindFromCollision, + collision->equation(), kindFromReaction); + } + m_collisions.emplace_back(collision); m_collisionRates.emplace_back( std::dynamic_pointer_cast(collision->rate())); @@ -457,7 +874,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")) { @@ -475,11 +891,50 @@ void PlasmaPhase::addCollision(shared_ptr collision) m_kInelastic.push_back(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(); + +} + +// 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 +{ + 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) @@ -526,7 +981,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()); @@ -553,7 +1008,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 @@ -587,6 +1052,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]; } @@ -777,7 +1248,8 @@ double PlasmaPhase::jouleHeatingPower() const return sigma * E * E; // W/m^3 } -double PlasmaPhase::intrinsicHeating() +// actually be careful, this is inelastic + growth +double PlasmaPhase::inelasticPower() { // Joule heating: sigma * E^2 [W/m^3] const double qJ = jouleHeatingPower(); @@ -787,8 +1259,92 @@ 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); + + 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); + } + } +} } diff --git a/test/data/air-plasma.yaml b/test/data/air-plasma.yaml index 4ff6bb09284..4059237a5ec 100644 --- a/test/data/air-plasma.yaml +++ b/test/data/air-plasma.yaml @@ -1,16 +1,113 @@ description: |- - Incomplete, simplified subset of the Phelps cross sections for O2/N2. - https://www.lxcat.net/Phelps. For testing purposes only. + A compilation of atomic and molecular cross-section data assembled by A. V. Phelps in + the 1970s and 1980s for gases such as O₂, N₂, He, Ar, etc. The compilation itself + is unpublished; data used from it are cited as: A. V. Phelps, private communication + (compilation of electron cross-sections), retrieved [8/1/25], from Phelps collection. + (see http://www.lxcat.net/contributors/#d19). + + This dataset is provided in a Cantera-compatible format only for demonstration + purposes. Any scientific use of this dataset should note the recommended citation + approach for data obtained from the LXCat project as described at + https://nl.lxcat.net/instructions/how_reference.php. units: {length: cm, quantity: molec, activation-energy: K} +date: Fri, 01 Aug 2025 00:00:00 -0400 + +species: +- name: Electron + composition: {E: 1} + thermo: + model: NASA7 + temperature-ranges: [200.0, 6000.0] + data: + - [2.5, 0.0, 0.0, 0.0, 0.0, -745.375, -11.7246902] + note: gas L10/92 +- name: N2 + composition: {N: 2} + thermo: + model: NASA7 + temperature-ranges: [200.0, 1000.0, 6000.0] + data: + - [3.53603521, -1.58270944e-04, -4.26984251e-07, 2.3754259e-09, -1.39708206e-12, + -1047.49645, 2.94603724] + - [2.9380297, 1.4183803e-03, -5.03281045e-07, 8.07555464e-11, -4.76064275e-15, + -917.18099, 5.95521985] + note: ATcT3E +- name: O2 + composition: {O: 2} + thermo: + model: NASA7 + temperature-ranges: [200.0, 1000.0, 6000.0] + data: + - [3.78498258, -3.02002233e-03, 9.92029171e-06, -9.77840434e-09, 3.28877702e-12, + -1064.13589, 3.64780709] + - [3.65980488, 6.59877372e-04, -1.44158172e-07, 2.14656037e-11, -1.36503784e-15, + -1216.03048, 3.42074148] + note: ATcT3E +- name: O2+ + composition: {O: 2, E: -1} + thermo: + model: NASA7 + temperature-ranges: [298.15, 1000.0, 6000.0] + data: + - [4.61017167, -6.35951952e-03, 1.42425624e-05, -1.20997923e-08, 3.70956878e-12, + 1.39742229e+05, -0.201326941] + - [3.31675922, 1.11522244e-03, -3.83492556e-07, 5.72784687e-11, -2.77648381e-15, + 1.39876823e+05, 5.44726469] + note: RUS 89 +- name: N2+ + composition: {N: 2, E: -1} + thermo: + model: NASA7 + temperature-ranges: [298.15, 1000.0, 6000.0] + data: + - [3.77540711, -2.06459157e-03, 4.75752301e-06, -3.15664228e-09, 6.70509973e-13, + 1.80481115e+05, 2.69322178] + - [3.58661363, 2.53071949e-04, 1.84778214e-07, -4.55257223e-11, 3.26818029e-15, + 1.80390994e+05, 3.09584143] + note: tpis89 +- name: O + composition: {O: 1} + thermo: + model: NASA7 + temperature-ranges: [200.0, 1000.0, 6000.0] + data: + - [3.1682671, -3.27931884e-03, 6.64306396e-06, -6.12806624e-09, 2.11265971e-12, + 2.91222592e+04, 2.05193346] + - [2.54363697, -2.73162486e-05, -4.1902952e-09, 4.95481845e-12, -4.79553694e-16, + 2.9226012e+04, 4.92229457] + note: L 1/90 +- name: O2- + composition: {O: 2, E: 1} + thermo: + model: NASA7 + temperature-ranges: [298.15, 1000.0, 6000.0] + data: + - [3.66442522, -9.28741138e-04, 6.45477082e-06, -7.7470338e-09, 2.93332662e-12, + -6870.76983, 4.35140674] + - [3.95666294, 5.98141823e-04, -2.12133905e-07, 3.63267581e-11, -2.24989228e-15, + -7062.87229, 2.2787101] + note: L 4/89 + phases: - name: air-plasma-Phelps thermo: plasma + elements: [N, O, E] kinetics: gas - species: - - nasa_gas.yaml/species: [Electron, O2, N2, O2+, N2+, O, O2-] + species: [N2, N2+, O2, O, O2-, O2+, Electron] state: {T: 300.0, P: 1 atm, X: {O2: 0.21, N2: 0.79}} + # electron-energy-distribution: + # type: Boltzmann-two-term + # initial_max_energy_level: 40.0 + # initial_number_of_energy_grid_cells: 23 + # energy_levels_distribution: Linear + # energy_grid_adaptation: + # enabled: false + # min_decay_decades: 10.0 + # max_decay_decades: 11.0 + # update_factor: 0.1 + # max_iterations: 1000 electron-energy-distribution: type: Boltzmann-two-term energy-levels: [0.0, 2.0, 4.0, 6.0, 8.0, 10.0, 12.0, 14.0, 16.0, 18.0, 20.0, 22.0, @@ -19,42 +116,27 @@ phases: reactions: - equation: O2 + Electron => Electron + Electron + O2+ type: electron-collision-plasma - energy-levels: [12.06, 13.0, 18.0, 28.0, 38.0, 48.0, 58.0, 68.0, 78.0, 88.0, 100.0, - 150.0, 200.0, 300.0, 500.0, 700.0, 1000.0, 1500.0, 2000.0, 3000.0, 5000.0, 7000.0, - 10000.0] - cross-sections: [0.0, 2.3e-22, 2e-21, 7.4e-21, 1.32e-20, 1.8e-20, 2.1e-20, 2.33e-20, - 2.5e-20, 2.6e-20, 2.7e-20, 2.7e-20, 2.5e-20, 2.17e-20, 1.66e-20, 1.35e-20, 1.04e-20, - 7.6e-21, 6e-21, 4.2e-21, 2.7e-21, 2e-21, 1.4e-21] - threshold: 12.06 -- equation: N2 + Electron => N2+ + 2 Electron - type: electron-collision-plasma - energy-levels: [0.0, 15.6, 16.0, 16.5, 17.0, 17.5, 18.0, 18.5, 19.0, 19.5, 20.0, 21.0, - 22.0, 23.0, 25.0, 30.0, 34.0, 45.0, 60.0, 75.0, 100.0, 150.0, 200.0] - cross-sections: [0.0, 0.0, 1.95e-22, 4.28e-22, 6.6e-22, 9.11e-22, 1.2e-21, 1.516e-21, - 1.841e-21, 2.13e-21, 2.502e-21, 3.181e-21, 3.869e-21, 4.557e-21, 5.924e-21, - 9.579e-21, 1.1718e-20, 1.6461e-20, 2.0181e-20, 2.2134e-20, 2.3436e-20, 2.2692e-20, - 2.1018e-20] - threshold: 15.6 + collision: phelps_O2_ionization_O2^+_12.06 + - equation: O2 + Electron => O2- type: electron-collision-plasma - energy-levels: [0.0, 0.058, 0.073, 0.083, 0.089, 0.095, 0.103, 0.109, 0.15, 0.17, 0.2, - 0.21, 0.23, 0.32, 0.33, 0.35, 0.44, 0.45, 0.47, 0.56, 0.57, 0.59, 0.68, 0.69, 0.71, - 0.79, 0.8, 0.82, 0.9, 0.91, 0.93, 1.02, 1.03, 1.05, 1.5, 100.0] - cross-sections: [0.0, 0.0, 5.60795e-41, 1.80256e-40, 4.20596e-41, 8.41193e-41, - 1.80256e-40, 0.0, 0.0, 0.0, 0.0, 3.56506e-41, 0.0, 0.0, 2.30327e-41, 0.0, 0.0, - 1.45206e-41, 0.0, 0.0, 1.10156e-41, 0.0, 0.0, 8.01136e-42, 0.0, 0.0, 7.00994e-42, - 0.0, 0.0, 5.50781e-42, 0.0, 0.0, 4.20596e-42, 0.0, 0.0, 0.0] - threshold: 0.0 + collision: phelps_O2_attachment_O2^-_0.0 + +- equation: N2 + Electron => N2+ + 2 Electron + type: electron-collision-plasma + collision: phelps_N2_ionization_N2^+_15.6 electron-collisions: -- target: N2 +- name: phelps_N2_effective_N2_0.0 + target: N2 energy-levels: [0.0, 0.015, 0.03, 0.05, 0.1, 0.15, 0.2, 0.3, 0.4, 0.7, 1.2, 1.5, 1.9, 2.2, 2.8, 3.3, 4.0, 5.0, 7.0, 10.0, 15.0, 20.0, 30.0, 75.0, 150.0] cross-sections: [1.1e-20, 2.55e-20, 3.4e-20, 4.33e-20, 5.95e-20, 7.1e-20, 7.9e-20, 9e-20, 9.7e-20, 1e-19, 1.04e-19, 1.2e-19, 1.96e-19, 2.85e-19, 2.8e-19, 1.72e-19, 1.26e-19, 1.09e-19, 1.01e-19, 1.04e-19, 1.1e-19, 1.02e-19, 9e-20, 6.6e-20, 4.9e-20] kind: effective -- target: N2 +- name: phelps_N2_excitation_N2(rot)_0.02 + target: N2 product: N2(rot) energy-levels: [0.02, 0.03, 0.4, 0.8, 1.2, 1.6, 1.7, 1.8, 1.9, 2.0, 2.1, 2.2, 2.3, 2.4, 2.5, 2.6, 2.7, 2.8, 2.9, 3.0, 3.1, 3.2, 3.3, 3.6, 5.0] @@ -63,7 +145,8 @@ electron-collisions: 1.62e-20, 1.38e-20, 1.18e-20, 1.03e-20, 8.4e-21, 6.9e-21, 5e-21, 1.7e-21, 0.0] threshold: 0.02 kind: excitation -- target: N2 +- name: phelps_N2_excitation_N2(v1)_0.29 + target: N2 product: N2(v1) energy-levels: [0.29, 0.3, 0.33, 0.4, 0.75, 0.9, 1.0, 1.1, 1.16, 1.2, 1.22, 1.4, 1.5, 1.6, 1.65, 3.6, 4.0, 5.0, 15.0, 18.0, 20.0, 22.0, 23.0, 25.0, 29.0, 32.0, 50.0, @@ -73,7 +156,8 @@ electron-collisions: 3.5e-22, 4e-22, 6.5e-22, 8.5e-22, 8.5e-22, 6e-22, 3e-22, 1.5e-22, 1.2e-22, 0.0] threshold: 0.29 kind: excitation -- target: N2 +- name: phelps_N2_excitation_N2(v1res)_0.291 + target: N2 product: N2(v1res) energy-levels: [0.0, 0.291, 1.6, 1.65, 1.7, 1.8, 1.9, 2.0, 2.1, 2.2, 2.3, 2.4, 2.5, 2.6, 2.7, 2.75, 2.8, 2.9, 3.0, 3.1, 3.2, 3.3, 3.4, 3.5, 3.6, 4.0, 100.0] @@ -83,7 +167,8 @@ electron-collisions: 0.0, 0.0] threshold: 0.291 kind: excitation -- target: N2 +- name: phelps_N2_excitation_N2(v2)_0.59 + target: N2 product: N2(v2) energy-levels: [0.0, 0.59, 1.7, 1.8, 1.9, 2.0, 2.1, 2.2, 2.3, 2.4, 2.5, 2.6, 2.7, 2.75, 2.8, 2.9, 3.0, 3.1, 3.2, 3.3, 3.4, 3.5, 3.6, 100.0] @@ -92,7 +177,8 @@ electron-collisions: 6.45e-21, 3.75e-21, 3.45e-21, 3e-21, 2.13e-21, 0.0, 0.0] threshold: 0.59 kind: excitation -- target: N2 +- name: phelps_N2_excitation_N2(v3)_0.88 + target: N2 product: N2(v3) energy-levels: [0.0, 0.88, 1.9, 2.0, 2.1, 2.2, 2.3, 2.4, 2.5, 2.6, 2.7, 2.75, 2.8, 2.9, 3.0, 3.1, 3.2, 3.3, 3.4, 100.0] @@ -101,7 +187,8 @@ electron-collisions: 2.91e-21, 6.15e-22, 0.0, 0.0] threshold: 0.88 kind: excitation -- target: N2 +- name: phelps_N2_excitation_N2(C3)_11.03 + target: N2 product: N2(C3) energy-levels: [11.03, 11.5, 12.0, 12.5, 13.0, 13.5, 13.8, 14.0, 14.2, 14.5, 15.0, 16.0, 17.0, 18.0, 19.0, 20.0, 22.0, 24.0, 26.0, 28.0, 30.0, 36.0, 40.0, 50.0, 70.0, @@ -112,7 +199,8 @@ electron-collisions: 0.0] threshold: 11.03 kind: excitation -- target: N2 +- name: phelps_N2_excitation_N2(B3)_7.35 + target: N2 product: N2(B3) energy-levels: [0.0, 7.35, 8.0, 9.0, 10.0, 11.0, 12.0, 13.0, 14.0, 15.0, 16.0, 17.0, 18.0, 20.0, 22.0, 26.0, 30.0, 34.0, 40.0, 50.0, 70.0, 150.0] @@ -121,7 +209,8 @@ electron-collisions: 9.51e-22, 8.04e-22, 6.77e-22, 5.63e-22, 4.29e-22, 2.68e-22, 6.7e-23, 0.0] threshold: 7.35 kind: excitation -- target: N2 +- name: phelps_N2_excitation_N2(a1)_8.55 + target: N2 product: N2(a1) energy-levels: [0.0, 8.55, 9.0, 14.0, 15.0, 16.0, 17.0, 18.0, 19.0, 20.0, 24.0, 26.0, 30.0, 40.0, 50.0, 70.0, 100.0, 150.0, 200.0] @@ -130,7 +219,8 @@ electron-collisions: 8.51e-22, 6.03e-22, 4.02e-22, 2.68e-22, 2.01e-22] threshold: 8.55 kind: excitation -- target: N2 +- name: phelps_N2_excitation_N2(w1)_8.89 + target: N2 product: N2(w1) energy-levels: [0.0, 8.89, 9.0, 10.0, 11.0, 12.0, 13.0, 14.0, 15.0, 16.0, 17.0, 18.0, 20.0, 22.0, 30.0, 38.0, 50.0, 150.0] @@ -139,23 +229,37 @@ electron-collisions: 4.7e-23, 0.0] threshold: 8.89 kind: excitation -- target: O2 +- name: phelps_N2_ionization_N2^+_15.6 + target: N2 + product: N2^+ + energy-levels: [0.0, 15.6, 16.0, 16.5, 17.0, 17.5, 18.0, 18.5, 19.0, 19.5, 20.0, 21.0, + 22.0, 23.0, 25.0, 30.0, 34.0, 45.0, 60.0, 75.0, 100.0, 150.0, 200.0] + cross-sections: [0.0, 0.0, 1.95e-22, 4.28e-22, 6.6e-22, 9.11e-22, 1.2e-21, 1.516e-21, + 1.841e-21, 2.13e-21, 2.502e-21, 3.181e-21, 3.869e-21, 4.557e-21, 5.924e-21, + 9.579e-21, 1.1718e-20, 1.6461e-20, 2.0181e-20, 2.2134e-20, 2.3436e-20, 2.2692e-20, + 2.1018e-20] + threshold: 15.6 + kind: ionization +- name: phelps_O2_attachment_O^-+O_0.0 + target: O2 + product: O^-+O energy-levels: [4.4, 4.9, 5.38, 5.86, 6.1, 6.48, 6.77, 7.05, 7.3, 7.53, 7.77, 8.0, 8.25, 8.73, 9.2, 9.68, 10.15, 11.35, 100.0] cross-sections: [0.0, 0.0, 2.3e-23, 7.2e-23, 1.08e-22, 1.38e-22, 1.52e-22, 1.56e-22, 1.48e-22, 1.31e-22, 1.1e-22, 8.4e-23, 5.4e-23, 2.8e-23, 1.4e-23, 8e-24, 8e-24, 8e-24, 0.0] threshold: 0.0 - product: O^-+O kind: attachment -- target: O2 +- name: phelps_O2_effective_O2_0.0 + target: O2 energy-levels: [0.0, 0.015, 0.03, 0.05, 0.1, 0.15, 0.2, 0.3, 0.4, 0.7, 1.2, 1.5, 1.9, 2.2, 2.8, 3.3, 4.0, 5.0, 7.0, 10.0, 15.0, 20.0, 30.0, 75.0, 150.0, 200.0] cross-sections: [3.5e-21, 8.7e-21, 1.24e-20, 1.6e-20, 2.5e-20, 3.1e-20, 3.6e-20, 4.5e-20, 5.2e-20, 6.1e-20, 7.9e-20, 7.6e-20, 6.9e-20, 6.5e-20, 5.8e-20, 5.5e-20, 5.5e-20, 5.6e-20, 6.6e-20, 8e-20, 8.8e-20, 8.6e-20, 8e-20, 6.8e-20, 6.7e-20, 6e-20] kind: effective -- target: O2 +- name: phelps_O2_excitation_O2(rot)_0.02 + target: O2 product: O2(rot) energy-levels: [0.07, 0.08, 0.1, 0.2, 0.21, 0.22, 0.32, 0.33, 0.35, 0.44, 0.45, 0.47, 0.56, 0.57, 0.59, 0.68, 0.69, 0.71, 0.79, 0.8, 0.81, 0.9, 0.91, 0.93, 1.02, 1.03, @@ -167,17 +271,42 @@ electron-collisions: 0.0, 0.0, 2.4e-22, 0.0, 0.0, 1.2e-22, 0.0, 0.0, 4.8e-23, 0.0] threshold: 0.02 kind: excitation -- target: O2 +- name: phelps_O2_excitation_O2(a1)_0.977 + target: O2 product: O2(a1) energy-levels: [0.977, 1.5, 3.5, 5.62, 6.53, 7.89, 13.0, 20.5, 41.0, 100.0] cross-sections: [0.0, 5.8e-23, 4.9e-22, 8.25e-22, 9.08e-22, 8.63e-22, 5.27e-22, 3.24e-22, 1.37e-22, 0.0] threshold: 0.977 kind: excitation -- target: O2 +- name: phelps_O2_excitation_O2(b1)_1.627 + target: O2 product: O2(b1) - energy-levels: [1.627, 3.0, 4.0, 7.34, 9.26, 13.0, 17.0, 20.7, 24.0, 35.1, 45.1, 100] + energy-levels: [1.627, 3.0, 4.0, 7.34, 9.26, 13.0, 17.0, 20.7, 24.0, 35.1, 45.1, 100.0] cross-sections: [0.0, 9.7e-23, 1.49e-22, 1.91e-22, 1.74e-22, 1.3e-22, 1.3e-22, 1.25e-22, 1e-22, 6.3e-23, 5e-24, 0.0] threshold: 1.627 kind: excitation +- name: phelps_O2_ionization_O2^+_12.06 + target: O2 + product: O2^+ + energy-levels: [12.06, 13.0, 18.0, 28.0, 38.0, 48.0, 58.0, 68.0, 78.0, 88.0, 100.0, + 150.0, 200.0, 300.0, 500.0, 700.0, 1000.0, 1500.0, 2000.0, 3000.0, 5000.0, 7000.0, + 10000.0] + cross-sections: [0.0, 2.3e-22, 2e-21, 7.4e-21, 1.32e-20, 1.8e-20, 2.1e-20, 2.33e-20, + 2.5e-20, 2.6e-20, 2.7e-20, 2.7e-20, 2.5e-20, 2.17e-20, 1.66e-20, 1.35e-20, 1.04e-20, + 7.6e-21, 6e-21, 4.2e-21, 2.7e-21, 2e-21, 1.4e-21] + threshold: 12.06 + kind: ionization +- name: phelps_O2_attachment_O2^-_0.0 + target: O2 + product: O2^- + energy-levels: [0.0, 0.058, 0.073, 0.083, 0.089, 0.095, 0.103, 0.109, 0.15, 0.17, 0.2, + 0.21, 0.23, 0.32, 0.33, 0.35, 0.44, 0.45, 0.47, 0.56, 0.57, 0.59, 0.68, 0.69, 0.71, + 0.79, 0.8, 0.82, 0.9, 0.91, 0.93, 1.02, 1.03, 1.05, 1.5, 100.0] + cross-sections: [0.0, 0.0, 5.60795e-41, 1.80256e-40, 4.20596e-41, 8.41193e-41, + 1.80256e-40, 0.0, 0.0, 0.0, 0.0, 3.56506e-41, 0.0, 0.0, 2.30327e-41, 0.0, 0.0, + 1.45206e-41, 0.0, 0.0, 1.10156e-41, 0.0, 0.0, 8.01136e-42, 0.0, 0.0, 7.00994e-42, + 0.0, 0.0, 5.50781e-42, 0.0, 0.0, 4.20596e-42, 0.0, 0.0, 0.0] + threshold: 0.0 + kind: attachment \ No newline at end of file diff --git a/test/data/lxcat40-test-converter.yaml b/test/data/lxcat40-test-converter.yaml new file mode 100644 index 00000000000..8784f7480ad --- /dev/null +++ b/test/data/lxcat40-test-converter.yaml @@ -0,0 +1,3 @@ +species: + O2(6.0eV): O + O + O2(8.4eV): O + O \ No newline at end of file diff --git a/test/data/lxcat40-test-mech.yaml b/test/data/lxcat40-test-mech.yaml new file mode 100644 index 00000000000..7b2d32eefd1 --- /dev/null +++ b/test/data/lxcat40-test-mech.yaml @@ -0,0 +1,83 @@ +description: |- + Minimal CanteraPlasma 4.0 mechanism used to test lxcat2yaml merge mode. + +units: + length: cm + time: s + quantity: mol + activation-energy: cal/mol + +phases: +- name: isotropic-electron-energy-plasma + thermo: plasma + elements: [E, O, C] + species: [Electron, O2, O2+, O, O-, CO, CO2] + kinetics: gas + reactions: all + electron-energy-distribution: + type: isotropic + shape-factor: 2.0 + mean-electron-energy: 5.0 eV + state: {T: 300.0, P: 1 atm} + +species: +- name: Electron + composition: {E: 1} + thermo: + model: constant-cp + T0: 300.0 + h0: 0.0 + s0: 0.0 + cp0: 20.786 +- name: O2 + composition: {O: 2} + thermo: + model: constant-cp + T0: 300.0 + h0: 0.0 + s0: 0.0 + cp0: 29.355 +- name: O2+ + composition: {O: 2, E: -1} + thermo: + model: constant-cp + T0: 300.0 + h0: 0.0 + s0: 0.0 + cp0: 29.355 +- name: O + composition: {O: 1} + thermo: + model: constant-cp + T0: 300.0 + h0: 0.0 + s0: 0.0 + cp0: 20.786 +- name: O- + composition: {O: 1, E: 1} + thermo: + model: constant-cp + T0: 300.0 + h0: 0.0 + s0: 0.0 + cp0: 20.786 +- name: CO + composition: {C: 1, O: 1} + thermo: + model: constant-cp + T0: 300.0 + h0: 0.0 + s0: 0.0 + cp0: 29.100 +- name: CO2 + composition: {C: 1, O: 2} + thermo: + model: constant-cp + T0: 300.0 + h0: 0.0 + s0: 0.0 + cp0: 37.135 + +reactions: +- equation: O + CO => CO2 + rate-constant: {A: 1.0e12, b: 0.0, Ea: 0.0} diff --git a/test/data/lxcat40-test.xml b/test/data/lxcat40-test.xml new file mode 100644 index 00000000000..cd0d3fb8945 --- /dev/null +++ b/test/data/lxcat40-test.xml @@ -0,0 +1,93 @@ + + + Minimal LXCat test file. + + + + + + + + e + CO2 + E + CO2 + + E + CO2 -> E + CO2 + + 1.0e-5 + + Effective momentum-transfer-like process. + 2026-06-08 00:00:00 + 0.0 1.0 + 1.0e-20 2.0e-20 + + + + e + O2 + E + E + O2^+ + + E + O2 -> E + E + O2^+ + + 12.0 + + Simple ionization process. + 2026-06-08 00:00:00 + 12.0 20.0 + 0.0 5.5e-22 + + + + e + O2 + E + O2(6.0eV) + + E + O2 -> E + O2(6.0eV) + + 6.0 + + Dissociation represented through converter mapping. + 2026-06-08 00:00:00 + 6.0 10.0 + 0.0 1.0e-21 + + + + e + O2 + E + O2(8.4eV) + + E + O2 -> E + O2(8.4eV) + + 8.4 + + Second dissociation channel mapped to the same phase-level equation. + 2026-06-08 00:00:00 + 8.4 12.0 + 0.0 2.0e-21 + + + + e + CO2 + CO+O^- + + E + CO2 -> CO + O^- + + 0.0 + + Dissociative attachment, threshold 3.85eV. + 2026-06-08 00:00:00 + 3.85 4.30 9.70 + 0.0 1.4e-23 0.0 + + + + + + \ No newline at end of file diff --git a/test/data/oxygen-plasma.yaml b/test/data/oxygen-plasma.yaml index 9c172c208bb..b2c4c905a8d 100644 --- a/test/data/oxygen-plasma.yaml +++ b/test/data/oxygen-plasma.yaml @@ -50,10 +50,17 @@ reactions: - equation: E + O2 + O2 => O2- + O2 type: two-temperature-plasma rate-constant: {A: 4.2e-27, b: -1.0, Ea-gas: 600, Ea-electron: 700} - -- equation: O2 + E <=> E + O2 +- equation: O2 + E <=> O2 + E type: electron-collision-plasma - note: This is a electron collision process of plasma + collision: test_O2_effective_O2_0 + +electron-collisions: +- name: test_O2_effective_O2_0 + target: O2 + product: O2 energy-levels: [0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0] cross-sections: [0.0, 5.97e-20, 6.45e-20, 6.74e-20, 6.93e-20, 7.2e-20, 7.52e-20, 7.86e-20, 8.21e-20, 8.49e-20, 8.8e-20] + kind: effective + threshold: 0 + diff --git a/test/kinetics/kineticsFromScratch.cpp b/test/kinetics/kineticsFromScratch.cpp index b4d120099db..332ab9cb111 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/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 fbe699b3889..d03da101282 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/VibrationalRelaxationRate.h" #include "cantera/thermo/SurfPhase.h" #include "cantera/thermo/ThermoFactory.h" #include "cantera/base/Array.h" diff --git a/test/python/test_convert.py b/test/python/test_convert.py index ac9d5ceada1..1b3392080ad 100644 --- a/test/python/test_convert.py +++ b/test/python/test_convert.py @@ -1678,84 +1678,105 @@ class Testlxcat2yaml: def inject_fixtures(self, test_data_path): self.test_data_path = test_data_path - def convert(self, inputFile=None, database=None, mechFile=None, phase=None, - insert=True, output=None): + def convert(self, inputFile=None, database=None, lxcat2phase=None, + mechFile=None, phase=None, output=None): if inputFile is not None: inputFile = self.test_data_path / inputFile + if lxcat2phase is not None: + lxcat2phase = self.test_data_path / lxcat2phase if mechFile is not None: mechFile = self.test_data_path / mechFile if output is None: - output = Path(inputFile).stem # strip '.xml' - # output to work dir + output = Path(inputFile).with_suffix(".yaml").name output = self.test_work_path / output + output.unlink(missing_ok=True) - lxcat2yaml.convert(inputFile, database, mechFile, phase, insert, output) + lxcat2yaml.convert( + inpfile=inputFile, + database=database, + lxcat2phase=lxcat2phase, + mechfile=mechFile, + phase=phase, + outfile=output, + ) return output - def test_mechanism_with_lxcat(self): - # get Solution from the mechanism file - phase = "isotropic-electron-energy-plasma" - mechFile = "lxcat-test-convert.yaml" - gas1 = ct.Solution(self.test_data_path / mechFile, - phase=phase, transport_model=None) - - # get a stand-alone collisions - standAloneFile = "stand-alone-lxcat.yaml" - self.convert(inputFile='lxcat-test-convert.xml', database="test", - mechFile=mechFile, insert=False, - output=standAloneFile) - - # add collisions to the reaction list - rxn_list = ct.Reaction.list_from_file(self.test_work_path / standAloneFile, - gas1, section="collisions") - for R in rxn_list: - gas1.add_reaction(R) - - # get Solution from the output file - output = "output-lxcat.yaml" - self.convert(inputFile="lxcat-test-convert.xml", database="test", - mechFile=mechFile, insert=True, - output=output) - gas2 = ct.Solution(self.test_work_path / output, - phase=phase, transport_model=None) - - # check number of reactions - assert gas1.n_reactions == gas2.n_reactions == 4 - for i in range(1, gas1.n_reactions): - assert (gas1.reaction(i).rate.energy_levels - == approx(gas2.reaction(i).rate.energy_levels)) - assert (gas1.reaction(i).rate.cross_sections - == approx(gas2.reaction(i).rate.cross_sections)) - - def test_stand_alone_lxcat(self): - outfile = "stand-alone-lxcat-without-mech.yaml" - self.convert(inputFile='lxcat-test-convert.xml', - database="test", insert=False, - output=outfile) - - # get Solution from the input file + def test_stand_alone_lxcat40(self): + output = self.convert( + inputFile="lxcat40-test.xml", + database="test", + output="stand-alone-lxcat40.yaml", + ) + + data = load_yaml(output) + assert "reactions" not in data + assert "electron-collisions" in data + + collisions = data["electron-collisions"] + assert len(collisions) == 5 + + by_name = {collision["name"]: collision for collision in collisions} + + effective = by_name["test_CO2_effective_CO2_0"] + assert effective["kind"] == "effective" + assert effective["target"] == "CO2" + assert effective["product"] == "CO2" + assert effective["threshold"] == approx(0.0) + assert effective["energy-levels"] == approx([0.0, 1.0]) + assert effective["cross-sections"] == approx([0.0, 2.0e-20]) + + ionization = by_name["test_O2_ionization_O2+_12"] + assert ionization["kind"] == "ionization" + assert ionization["target"] == "O2" + assert ionization["product"] == "O2+" + assert ionization["threshold"] == approx(12.0) + assert ionization["energy-levels"] == approx([12.0, 20.0]) + assert ionization["cross-sections"] == approx([0.0, 5.5e-22]) + + attachment = by_name["test_CO2_attachment_CO+O-_3.85"] + assert attachment["kind"] == "attachment" + assert attachment["target"] == "CO2" + assert attachment["product"] == "CO + O-" + assert attachment["threshold"] == approx(3.85) + assert attachment["energy-levels"] == approx([3.85, 4.30, 9.70]) + assert attachment["cross-sections"] == approx([0.0, 1.4e-23, 0.0]) + + def test_lxcat40_merge_into_mechanism(self): phase = "isotropic-electron-energy-plasma" - inputFile = "lxcat-test-convert-species.yaml" - gas = ct.Solution(self.test_data_path / inputFile, - phase=phase, transport_model=None) - - # add collisions to the reaction list - rxn_list = ct.Reaction.list_from_file(self.test_work_path / outfile, - gas, section="collisions") - - # verify the data - assert len(rxn_list) == 3 - assert rxn_list[0].equation == "CO2 + e => CO2 + e" - assert rxn_list[0].reaction_type == "electron-collision-plasma" - assert rxn_list[0].rate.energy_levels == approx([0.0, 1.0]) - assert rxn_list[0].rate.cross_sections == approx([0.0, 1.0e-22]) - - assert rxn_list[1].equation == "O2 + e => O2(Total-Ionization)+ + 2 e" - assert rxn_list[1].reaction_type == "electron-collision-plasma" - assert rxn_list[1].rate.energy_levels == approx([15., 20.]) - assert rxn_list[1].rate.cross_sections == approx([0.0, 5.5e-22]) - - assert rxn_list[2].equation == "O2 + e => O2-" - assert rxn_list[2].reaction_type == "electron-collision-plasma" - assert rxn_list[2].rate.energy_levels == approx([0.0, 1.0]) - assert rxn_list[2].rate.cross_sections == approx([0.0, 1.0e-22]) + output = self.convert( + inputFile="lxcat40-test.xml", + database="test", + lxcat2phase="lxcat40-test-converter.yaml", + mechFile="lxcat40-test-mech.yaml", + phase=phase, + output="lxcat40-test-merged.yaml", + ) + + data = load_yaml(output) + assert "reactions" in data + assert "electron-collisions" in data + assert len(data["electron-collisions"]) == 5 + + plasma_reactions = [ + reaction for reaction in data["reactions"] + if reaction.get("type") == "electron-collision-plasma" + ] + assert len(plasma_reactions) == 4 + + equations = [reaction["equation"] for reaction in plasma_reactions] + assert "Electron + O2 => Electron + Electron + O2+" in equations + assert "Electron + O2 => Electron + O + O" in equations + assert "Electron + CO2 => CO + O-" in equations + + duplicate_reactions = [ + reaction for reaction in plasma_reactions + if reaction.get("duplicate") is True + ] + assert len(duplicate_reactions) == 2 + assert all( + reaction["equation"] == "Electron + O2 => Electron + O + O" + for reaction in duplicate_reactions + ) + + gas = ct.Solution(output, phase=phase, transport_model=None) + assert gas.n_reactions == 5 \ No newline at end of file diff --git a/test/python/test_thermo.py b/test/python/test_thermo.py index 8e71ce310ef..1a68417cff1 100644 --- a/test/python/test_thermo.py +++ b/test/python/test_thermo.py @@ -1496,20 +1496,11 @@ def test_eedf_solver(self): grid = phase.electron_energy_levels eedf = phase.electron_energy_distribution - reference_grid = np.logspace(-1, np.log10(60)) - - reference_eedf = np.array([ - 9.1027381e-02, 9.1026393e-02, 9.1025267e-02, 9.1023985e-02, 9.1022523e-02, - 9.1020858e-02, 9.1015025e-02, 9.1006713e-02, 9.0997242e-02, 9.0986450e-02, - 9.0974154e-02, 9.0954654e-02, 9.0923885e-02, 9.0888824e-02, 9.0842837e-02, - 9.0775447e-02, 9.0695937e-02, 9.0578309e-02, 9.0398980e-02, 9.0118320e-02, - 8.9293838e-02, 8.7498617e-02, 8.3767419e-02, 7.5765714e-02, 6.4856820e-02, - 5.5592157e-02, 4.9309310e-02, 4.5268611e-02, 4.2261381e-02, 3.9440745e-02, - 3.6437762e-02, 3.3181527e-02, 2.9616717e-02, 2.5795007e-02, 2.1676205e-02, - 1.7347058e-02, 1.3022044e-02, 8.9705614e-03, 5.5251937e-03, 3.1894295e-03, - 1.7301525e-03, 8.4647152e-04, 3.6030983e-04, 1.2894755e-04, 3.7416645e-05, - 8.4693678e-06, 1.4299900e-06, 1.7026957e-07, 1.3992350e-08, 1.5340110e-09 - ]) + reference_grid = np.array([ 0., 2., 4., 6., 8., 10., 12., 14., 16., 18., 20., 22., 24., 26., 28., 30., 32., 34., 36., 38., 40.]) + + reference_eedf = np.array([9.8075881952895655e-02, 7.5100960725059007e-02, 4.6382335178763986e-02, 3.5561250551010629e-02, 2.5998954792519895e-02, 1.7784714129416464e-02, 1.1249679269582400e-02, 6.6780916799679298e-03, + 3.9075622604817523e-03, 2.3070689055146782e-03, 1.3538410735400440e-03, 7.8369063569500214e-04, 4.4631163146510136e-04, 2.4997717421735031e-04, 1.3779808478281695e-04, 7.5027829771846555e-05, + 4.0530766221145227e-05, 2.1899425412039534e-05, 1.2236213446791840e-05, 7.6958974488623866e-06, 6.3734125240941400e-06]) interp = np.interp(reference_grid, grid, eedf, left=0.0, right=0.0) diff --git a/test/zeroD/test_zeroD.cpp b/test/zeroD/test_zeroD.cpp index abf513e1fba..f56f050950b 100644 --- a/test/zeroD/test_zeroD.cpp +++ b/test/zeroD/test_zeroD.cpp @@ -212,7 +212,7 @@ TEST(zerodim, plasma_reactor_energy) const double t_end = 1e-3; ASSERT_NO_THROW(net.advance(t_end)); const double T_final = reactor->temperature(); - const double T_expected = 300.04674410693019; + const double T_expected = 300.07632822936; const double rtol = 1e-6; // Simple regression test EXPECT_NEAR(T_final, T_expected, rtol * T_expected);