[reaction, interface] add ElectronCollisionPlasmaRate

include/cantera/kinetics/ElectronCollisionPlasmaRate.h

@@ -0,0 +1,121 @@
+//! @file ElectronCollisionPlasmaRate.h   Header for plasma reaction rates parameterized
+//!     by electron collision cross section and electron energy distribution.
+
+// This file is part of Cantera. See License.txt in the top-level directory or
+// at https://cantera.org/license.txt for license and copyright information.
+
+#ifndef CT_ELECTRONCOLLISIONPLASMARATE_H
+#define CT_ELECTRONCOLLISIONPLASMARATE_H
+
+#include "cantera/thermo/PlasmaPhase.h"
+#include "cantera/kinetics/ReactionData.h"
+#include "ReactionRate.h"
+#include "MultiRate.h"
+
+namespace Cantera
+{
+
+//! Data container holding shared data specific to ElectronCollisionPlasmaRate
+/**
+ * The data container `ElectronCollisionPlasmaData` holds precalculated data common to
+ * all `ElectronCollisionPlasmaRate` objects.
+ */
+struct ElectronCollisionPlasmaData : public ReactionData
+{
+    ElectronCollisionPlasmaData();
+
+    virtual bool update(const ThermoPhase& phase, const Kinetics& kin) override;
+    using ReactionData::update;
+
+    virtual void invalidateCache() override {
+        ReactionData::invalidateCache();
+        energyLevels.resize(0);
+        distribution.resize(0);
+    }
+
+    vector<double> energyLevels; //!< electron energy levels
+    vector<double> distribution; //!< electron energy distribution
+};
+
+
+//! Electron collision plasma reaction rate type
+/*!
+ * The electron collision plasma reaction rate uses the electron collision
+ * data and the electron energy distribution to calculate the reaction rate.
+ *
+ *   @f[
+ *        k_f =  \gamma \int_0^{\infty} \epsilon \sigma F_0 d\epsilon
+ *   @f]
+ *
+ * where @f$ \gamma = (2e/m)^{1/2} @f$ is a constant, @f$ \epsilon @f$ is
+ * the electron energy, @f$ \sigma @f$ is the reaction collision cross
+ * section, and @f$ F_0 @f$ is the normalized electron energy distribution
+ * function.
+ *
+ * References:
+ *
+ * [1] Hagelaar and Pitchford @cite hagelaar2005.
+ *
+ * @ingroup plasma
+ */
+class ElectronCollisionPlasmaRate : public ReactionRate
+{
+public:
+    ElectronCollisionPlasmaRate() = default;
+
+    ElectronCollisionPlasmaRate(const AnyMap& node,
+                                const UnitStack& rate_units={})
+        : ElectronCollisionPlasmaRate()
+    {
+        setParameters(node, rate_units);
+    }
+
+    virtual void setParameters(const AnyMap& node, const UnitStack& units) override;
+
+    virtual void getParameters(AnyMap& node) const override;
+
+    unique_ptr<MultiRateBase> newMultiRate() const override {
+        return unique_ptr<MultiRateBase>(
+            new MultiRate<ElectronCollisionPlasmaRate, ElectronCollisionPlasmaData>);
+    }
+
+    virtual const std::string type() const override {
+        return "electron-collision-plasma";
+    }
+
+    virtual void setContext(const Reaction& rxn, const Kinetics& kin) override;
+
+    //! Evaluate reaction rate
+    /*!
+     *  @param shared_data  data shared by all reactions of a given type
+     */
+    double evalFromStruct(const ElectronCollisionPlasmaData& shared_data) const;
+
+    //! Evaluate derivative of reaction rate with respect to temperature
+    //! divided by reaction rate
+    //! @param shared_data  data shared by all reactions of a given type
+    double ddTScaledFromStruct(const ElectronCollisionPlasmaData& shared_data) const {
+        throw NotImplementedError("ElectronCollisionPlasmaRate::ddTScaledFromStruct");
+    }
+
+    //! The value of #m_energyLevels [eV]
+    const vector<double>& energyLevels() const {
+        return m_energyLevels;
+    }
+
+    //! The value of #m_crossSections [m2]
+    const vector<double>& crossSections() const {
+        return m_crossSections;
+    }
+
+private:
+    //! electron energy levels [eV]
+    vector<double> m_energyLevels;
+
+    //! collision cross sections [m2] at #m_energyLevels
+    vector<double> m_crossSections;
+};
+
+}
+
+#endif

include/cantera/thermo/PlasmaPhase.h

@@ -249,6 +249,11 @@ public:
     void setParameters(const AnyMap& phaseNode,
                        const AnyMap& rootNode=AnyMap()) override;
 
+    //! Electron species name
+    string electronSpeciesName() const {
+        return speciesName(m_electronSpeciesIndex);
+    }
+
 protected:
     void updateThermo() const override;
 

interfaces/cython/cantera/reaction.pxd

@@ -73,6 +73,11 @@ cdef extern from "cantera/kinetics/TwoTempPlasmaRate.h" namespace "Cantera":
         CxxTwoTempPlasmaRate(double, double, double, double)
         double activationElectronEnergy()
 
+cdef extern from "cantera/kinetics/ElectronCollisionPlasmaRate.h" namespace "Cantera":
+    cdef cppclass CxxElectronCollisionPlasmaRate "Cantera::ElectronCollisionPlasmaRate" (CxxReactionRate):
+        CxxElectronCollisionPlasmaRate(CxxAnyMap) except +translate_exception
+        vector[double]& energyLevels()
+        vector[double]& crossSections()
 
 cdef extern from "cantera/base/Array.h" namespace "Cantera":
     cdef cppclass CxxArray2D "Cantera::Array2D":
@@ -251,6 +256,10 @@ cdef class ReactionRate:
 cdef class ArrheniusRateBase(ReactionRate):
     cdef CxxArrheniusBase* base
 
+cdef class ElectronCollisionPlasmaRate(ReactionRate):
+    cdef CxxElectronCollisionPlasmaRate* base
+    cdef set_cxx_object(self)
+
 cdef class FalloffRate(ReactionRate):
     cdef CxxFalloffRate* falloff
     cdef set_cxx_object(self)

interfaces/cython/cantera/reaction.pyx

@@ -370,6 +370,67 @@ cdef class TwoTempPlasmaRate(ArrheniusRateBase):
             return self.cxx_object().activationElectronEnergy()
 
 
+cdef class ElectronCollisionPlasmaRate(ReactionRate):
+    r"""
+    A reaction rate coefficient which depends on electron collision cross section
+    and electron energy distribution
+
+    """
+    _reaction_rate_type = "electron-collision-plasma"
+
+    def __cinit__(self, energy_levels=None, cross_sections=None,
+                  input_data=None, init=True):
+        if init:
+            if isinstance(input_data, dict):
+                self._from_dict(input_data)
+            elif energy_levels is not None and cross_sections is not None:
+                self._from_parameters(energy_levels, cross_sections)
+            elif input_data is None:
+                self._from_dict({})
+            else:
+                raise TypeError("Invalid input parameters")
+            self.set_cxx_object()
+
+    def _from_dict(self, dict input_data):
+        self._rate.reset(
+            new CxxElectronCollisionPlasmaRate(py_to_anymap(input_data, hyphenize=True))
+        )
+
+    def _from_parameters(self, energy_levels, cross_sections):
+        # check length
+        if len(energy_levels) != len(cross_sections):
+            raise ValueError('Length of energy levels and '
+                             'cross sections are different')
+        cdef np.ndarray[np.double_t, ndim=1] data_energy_levels = \
+            np.ascontiguousarray(energy_levels, dtype=np.double)
+        cdef np.ndarray[np.double_t, ndim=1] data_cross_sections = \
+            np.ascontiguousarray(cross_sections, dtype=np.double)
+        input_data = {'type': 'electron-collision-plasma',
+                      'energy-levels': data_energy_levels,
+                      'cross-sections': data_cross_sections}
+        self._from_dict(input_data)
+
+    cdef set_cxx_object(self):
+        self.rate = self._rate.get()
+        self.base = <CxxElectronCollisionPlasmaRate*>self.rate
+
+    property energy_levels:
+        """
+        The energy levels [eV]. Each level corresponds to a cross section
+        of `cross_sections`.
+        """
+        def __get__(self):
+            return np.fromiter(self.base.energyLevels(), np.double)
+
+    property cross_sections:
+        """
+        The cross sections [m2]. Each cross section corresponds to a energy
+        level of `energy_levels`.
+        """
+        def __get__(self):
+            return np.fromiter(self.base.crossSections(), np.double)
+
+
 cdef class FalloffRate(ReactionRate):
     """
     Base class for parameterizations used to describe the fall-off in reaction rates
@@ -1463,7 +1524,7 @@ cdef class Reaction:
         Directories on Cantera's input file path will be searched for the
         specified file.
         """
-        root = AnyMapFromYamlFile(stringify(filename))
+        root = AnyMapFromYamlFile(stringify(str(filename)))
         cxx_reactions = CxxGetReactions(root[stringify(section)],
                                         deref(kinetics.kinetics))
         return [Reaction.wrap(r) for r in cxx_reactions]

interfaces/cython/cantera/thermo.pxd

@@ -209,6 +209,7 @@ cdef extern from "cantera/thermo/PlasmaPhase.h":
         double meanElectronEnergy()
         size_t nElectronEnergyLevels()
         double electronPressure()
+        string electronSpeciesName()
 
 
 cdef extern from "cantera/cython/thermo_utils.h":

interfaces/cython/cantera/thermo.pyx

@@ -1865,7 +1865,7 @@ cdef class ThermoPhase(_SolutionBase):
         def __set__(self, distribution_type):
             if not self._enable_plasma:
                 raise ThermoModelMethodError(self.thermo_model)
-            self.plasma.setElectronEnergyDistributionType(distribution_type)
+            self.plasma.setElectronEnergyDistributionType(stringify(distribution_type))
 
     property mean_electron_energy:
         """ Mean electron energy [eV] """
@@ -1900,6 +1900,13 @@ cdef class ThermoPhase(_SolutionBase):
                 raise ThermoModelMethodError(self.thermo_model)
             self.plasma.enableNormalizeElectronEnergyDist(enable)
 
+    property electron_species_name:
+        """ Electron species name """
+        def __get__(self):
+            if not self._enable_plasma:
+                raise ThermoModelMethodError(self.thermo_model)
+            return pystr(self.plasma.electronSpeciesName())
+
 
 cdef class InterfacePhase(ThermoPhase):
     """ A class representing a surface, edge phase """

src/kinetics/ElectronCollisionPlasmaRate.cpp

@@ -0,0 +1,155 @@
+//! @file ElectronCollisionPlasmaRate.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/ElectronCollisionPlasmaRate.h"
+#include "cantera/kinetics/Reaction.h"
+#include "cantera/kinetics/Kinetics.h"
+#include "cantera/numerics/funcs.h"
+
+namespace Cantera
+{
+
+ElectronCollisionPlasmaData::ElectronCollisionPlasmaData()
+{
+    energyLevels.assign(1, 0.0);
+    distribution.assign(1, 0.0);
+}
+
+bool ElectronCollisionPlasmaData::update(const ThermoPhase& phase,
+                                         const Kinetics& kin)
+{
+    const PlasmaPhase& pp = dynamic_cast<const PlasmaPhase&>(phase);
+
+    vector<double> levels(pp.nElectronEnergyLevels());
+    pp.getElectronEnergyLevels(levels.data());
+
+    vector<double> dist(pp.nElectronEnergyLevels());
+    pp.getElectronEnergyDistribution(dist.data());
+
+    bool changed = false;
+
+    if (levels != energyLevels) {
+        energyLevels = std::move(levels);
+        changed = true;
+    }
+
+    if (dist != distribution) {
+        distribution = std::move(dist);
+        changed = true;
+    }
+
+    return changed;
+}
+
+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("energy-levels")) {
+        m_energyLevels = node["energy-levels"].asVector<double>();
+    }
+    if (node.hasKey("cross-sections")) {
+        m_crossSections = node["cross-sections"].asVector<double>();
+    }
+    if (m_energyLevels.size() != m_crossSections.size()) {
+        throw CanteraError("ElectronCollisionPlasmaRate::setParameters",
+            "Energy levels and cross section must have the same length.");
+    }
+}
+
+void ElectronCollisionPlasmaRate::getParameters(AnyMap& node) const {
+    node["type"] = type();
+    node["energy-levels"] = m_energyLevels;
+    node["cross-sections"] = m_crossSections;
+}
+
+double ElectronCollisionPlasmaRate::evalFromStruct(const ElectronCollisionPlasmaData& shared_data) const
+{
+    // Interpolate cross-sections data to the energy levels of
+    // the electron energy distribution function
+    vector<double> cross_sections;
+    for (auto i = 0; i < shared_data.energyLevels.size(); i ++ ) {
+        cross_sections.push_back(linearInterp(shared_data.energyLevels[i],
+                                        m_energyLevels, m_crossSections));
+    }
+
+    // Map cross sections to Eigen::ArrayXd
+    auto cs_array = Eigen::Map<const Eigen::ArrayXd>(
+        cross_sections.data(), cross_sections.size()
+    );
+
+    // Map energyLevels in Eigen::ArrayXd
+    auto eps = Eigen::Map<const Eigen::ArrayXd>(
+        shared_data.energyLevels.data(), shared_data.energyLevels.size()
+    );
+
+    // Map energyLevels in Eigen::ArrayXd
+    auto distribution = Eigen::Map<const Eigen::ArrayXd>(
+        shared_data.distribution.data(), shared_data.distribution.size()
+    );
+
+    // unit in kmol/m3/s
+    return 0.5 * pow(2.0 * ElectronCharge / ElectronMass, 0.5) * Avogadro *
+           simpson(distribution.cwiseProduct(cs_array), eps.pow(2.0));
+}
+
+void ElectronCollisionPlasmaRate::setContext(const Reaction& rxn, const Kinetics& kin)
+{
+    // ElectronCollisionPlasmaReaction is for a non-equilibrium plasma, and the reverse rate
+    // cannot be calculated from the conventional thermochemistry.
+    // @todo implement the reversible rate for non-equilibrium plasma
+    if (rxn.reversible) {
+        throw InputFileError("ElectronCollisionPlasmaRate::setContext", rxn.input,
+            "ElectronCollisionPlasmaRate does not support reversible reactions");
+    }
+
+    // get electron species name
+    string electronName;
+    if (kin.thermo().type() == "plasma") {
+        electronName = dynamic_cast<const PlasmaPhase&>(kin.thermo()).electronSpeciesName();
+    } else {
+        throw CanteraError("ElectronCollisionPlasmaRate::setContext",
+                           "ElectronCollisionPlasmaRate requires plasma phase");
+    }
+
+    // Check rxn.reactantString (rxn.reactants gives reduced reactants)
+    std::istringstream iss(rxn.reactantString());
+    string token;
+    int nElectron = 0;
+    int nReactants = 0;
+
+    // Count number of electron and reactants
+    // Since the reactants are one electron and one molecule,
+    // no token can be a number.
+    while (iss >> token) {
+        if (isdigit(token[0])) {
+            throw InputFileError("ElectronCollisionPlasmaRate::setContext", rxn.input,
+                "ElectronCollisionPlasmaRate requires one electron and one molecule as reactants");
+        }
+        if (token == electronName) {
+            nElectron++;
+        }
+        if (token != "+") {
+            nReactants++;
+        }
+    }
+
+    // Number of reactants needs to be two
+    if (nReactants != 2) {
+        throw InputFileError("ElectronCollisionPlasmaRate::setContext", rxn.input,
+            "ElectronCollisionPlasmaRate requires exactly two reactants");
+    }
+
+    // Must have only one electron
+    // @todo add electron-electron collision rate
+    if (nElectron != 1) {
+        throw InputFileError("ElectronCollisionPlasmaRate::setContext", rxn.input,
+            "ElectronCollisionPlasmaRate requires one and only one electron");
+    }
+}
+
+}

src/kinetics/ReactionRateFactory.cpp

@@ -12,6 +12,7 @@
 #include "cantera/kinetics/Arrhenius.h"
 #include "cantera/kinetics/ChebyshevRate.h"
 #include "cantera/kinetics/Custom.h"
+#include "cantera/kinetics/ElectronCollisionPlasmaRate.h"
 #include "cantera/kinetics/Falloff.h"
 #include "cantera/kinetics/InterfaceRate.h"
 #include "cantera/kinetics/PlogRate.h"
@@ -38,6 +39,11 @@ ReactionRateFactory::ReactionRateFactory()
         return new TwoTempPlasmaRate(node, rate_units);
     });
 
+    // ElectronCollisionPlasmaRate evaluator
+    reg("electron-collision-plasma", [](const AnyMap& node, const UnitStack& rate_units) {
+        return new ElectronCollisionPlasmaRate(node, rate_units);
+    });
+
     // BlowersMaselRate evaluator
     reg("Blowers-Masel", [](const AnyMap& node, const UnitStack& rate_units) {
         return new BlowersMaselRate(node, rate_units);