Enhance the PlasmaPhase to use the cross section data with an example of calculating the elastic collision energy loss rate

include/cantera/kinetics/ElectronCollisionPlasmaRate.h

@@ -31,6 +31,8 @@ struct ElectronCollisionPlasmaData : public ReactionData
         ReactionData::invalidateCache();
         energyLevels.resize(0);
         distribution.resize(0);
+        m_dist_number = -1;
+        m_level_number = -1;
     }
 
     vector<double> energyLevels; //!< electron energy levels
@@ -89,6 +91,7 @@ struct ElectronCollisionPlasmaData : public ReactionData
 class ElectronCollisionPlasmaRate : public ReactionRate
 {
 public:
+
     ElectronCollisionPlasmaRate() = default;
 
     ElectronCollisionPlasmaRate(const AnyMap& node,
@@ -135,6 +138,14 @@ class ElectronCollisionPlasmaRate : public ReactionRate
         return m_crossSections;
     }
 
+    //! The value of #m_crossSectionsInterpolated [m2]
+    const vector<double>& crossSectionInterpolated() const {
+        return m_crossSectionsInterpolated;
+    }
+
+    //! Update the value of #m_crossSectionsInterpolated [m2]
+    void updateInterpolatedCrossSection(const vector<double>&);
+
 private:
     //! electron energy levels [eV]
     vector<double> m_energyLevels;

include/cantera/kinetics/Kinetics.h

@@ -1407,6 +1407,26 @@
         return m_root.lock();
     }
 
+    //! Register a function to be called if reaction is added.
+    //! @param id  A unique ID corresponding to the object affected by the callback.
+    //!   Typically, this is a pointer to an object that also holds a reference to the
+    //!   Kinetics object.
+    //! @param callback  The callback function to be called after any reaction is added.
+    //! When the callback becomes invalid (for example, the corresponding object is
+    //! being deleted, the removeReactionAddedCallback() method must be invoked.
+    //! @since New in %Cantera 3.1
+    void registerReactionAddedCallback(void* id, const function<void()>& callback)
+    {
+        m_reactionAddedCallbacks[id] = callback;
+    }
+
+    //! Remove the reaction-changed callback function associated with the specified object.
+    //! @since New in %Cantera 3.1
+    void removeReactionAddedCallback(void* id)
+    {
+        m_reactionAddedCallbacks.erase(id);
+    }
+
 protected:
     //! Cache for saved calculations within each Kinetics object.
     ValueCache m_cache;
@@ -1530,6 +1550,9 @@
 
     //! reference to Solution
     std::weak_ptr<Solution> m_root;
+
+    //! Callback functions that are invoked when the reaction is added.
+    map<void*, function<void()>> m_reactionAddedCallbacks;
 };
 
 }

include/cantera/kinetics/Reaction.h

@@ -162,6 +162,18 @@ class Reaction
         return bool(m_third_body);
     }
 
+    //! Register a function to be called if setRate is used.
+    //! @param id  A unique ID corresponding to the object affected by the callback.
+    //! @param callback  The callback function to be called after setRate is called.
+    //! When the callback becomes invalid (for example, the corresponding object is
+    //! being deleted), the removeSetRateCallback() method must be invoked.
+    //! @since New in %Cantera 3.2
+    void registerSetRateCallback(void* id, const function<void()>& callback);
+
+    //! Remove the setRate callback function associated with the specified object.
+    //! @since New in %Cantera 3.2
+    void removeSetRateCallback(void* id);
+
 protected:
     //! Store the parameters of a Reaction needed to reconstruct an identical
     //! object using the newReaction(AnyMap&, Kinetics&) function. Does not
@@ -183,6 +195,9 @@ class Reaction
     //! Relative efficiencies of third-body species in enhancing the reaction rate
     //! (if applicable)
     shared_ptr<ThirdBody> m_third_body;
+
+    //! Callback functions that are invoked when the rate is replaced.
+    map<void*, function<void()>> m_setRateCallbacks;
 };
 
 

include/cantera/thermo/Phase.h

@@ -18,6 +18,7 @@ namespace Cantera
 
 class Solution;
 class Species;
+class Kinetics;
 
 //! Class Phase is the base class for phases of matter, managing the species and
 //! elements in a phase, as well as the independent variables of temperature,

include/cantera/thermo/PlasmaPhase.h

@@ -15,38 +15,63 @@
 namespace Cantera
 {
 
-/**
- * Base class for a phase with plasma properties. This class manages the
- * plasma properties such as electron energy distribution function (EEDF).
- * There are two ways to define the electron distribution and electron
- * temperature. The first method uses setElectronTemperature() to set
- * the electron temperature which is used to calculate the electron energy
- * distribution with isotropic-velocity model. The generalized electron
- * energy distribution for isotropic-velocity distribution can be
- * expressed as [1,2],
+class Reaction;
+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
+ * using two primary methods for defining the electron distribution and electron
+ * temperature.
+ *
+ * The first method utilizes setElectronTemperature(), which sets the electron
+ * temperature and calculates the electron energy distribution assuming an
+ * isotropic-velocity model. Note that all units in PlasmaPhase are in SI, except
+ * for electron energy, which is measured in volts.
+ *
+ * The generalized electron energy distribution for an isotropic-velocity
+ * distribution (as described by Gudmundsson @cite gudmundsson2001 and Khalilpour
+ * and Foroutan @cite khalilpour2020)
+ * is given by:
  *   @f[
  *          f(\epsilon) = c_1 \frac{\sqrt{\epsilon}}{\epsilon_m^{3/2}}
- *          \exp(-c_2 (\frac{\epsilon}{\epsilon_m})^x),
+ *          \exp \left(-c_2 \left(\frac{\epsilon}{\epsilon_m}\right)^x \right),
  *   @f]
- * where @f$ x = 1 @f$ and @f$ x = 2 @f$ correspond to the Maxwellian and
- * Druyvesteyn (default) electron energy distribution, respectively.
- * @f$ \epsilon_m = 3/2 T_e @f$ [eV] (mean electron energy). The second
- * method uses setDiscretizedElectronEnergyDist() to manually set electron
- * energy distribution and calculate electron temperature from mean electron
- * energy, which is calculated as [3],
+ * where @f$ x = 1 @f$ corresponds to a Maxwellian distribution and
+ * @f$ x = 2 @f$ corresponds to a Druyvesteyn distribution (which is the
+ * default). Here, @f$ \epsilon_m = \frac{3}{2} T_e @f$ [V] represents the
+ * mean electron energy.
+ *
+ * The total probability distribution integrates to one:
  *   @f[
- *          \epsilon_m = \int_0^{\infty} \epsilon^{3/2} f(\epsilon) d\epsilon,
+ *           \int_0^{\infty} f(\epsilon) d\epsilon = 1.
  *   @f]
- * which can be calculated using trapezoidal rule,
+ * According to Hagelaar and Pitchford @cite hagelaar2005, the electron energy
+ * probability function can be defined as
+ * @f$ F(\epsilon) = \frac{f(\epsilon)}{\sqrt{\epsilon}} @f$ with units of
+ * [V@f$^{-3/2}@f$]. The generalized form of the electron energy probability
+ * function for isotropic-velocity distributions is:
  *   @f[
- *          \epsilon_m = \sum_i (\epsilon^{5/2}_{i+1} - \epsilon^{5/2}_i)
- *                       (f(\epsilon_{i+1}) + f(\epsilon_i)) / 2,
+ *          F(\epsilon) = c_1 \frac{1}{\epsilon_m^{3/2}}
+ *          \exp\left(-c_2 \left(\frac{\epsilon}{\epsilon_m}\right)^x\right),
  *   @f]
- * where @f$ i @f$ is the index of energy levels.
+ * and this form is used to model the isotropic electron energy distribution
+ * in PlasmaPhase.
  *
- * For references, see Gudmundsson @cite gudmundsson2001; Khalilpour and Foroutan
- * @cite khalilpour2020; Hagelaar and Pitchford @cite hagelaar2005, and BOLOS
- * @cite BOLOS.
+ * The second method allows for manual definition of the electron energy
+ * distribution using setDiscretizedElectronEnergyDist(). In this approach,
+ * the electron temperature is derived from the mean electron energy,
+ * @f$ \epsilon_m @f$, which can be calculated as follows @cite hagelaar2005 :
+ *   @f[
+ *          \epsilon_m = \int_0^{\infty} \epsilon^{3/2} F(\epsilon) d\epsilon.
+ *   @f]
+ * This integral can be approximated using the trapezoidal rule,
+ *   @f[
+ *          \epsilon_m = \sum_i \left(\epsilon_{i+1}^{5/2} - \epsilon_i^{5/2}\right)
+ *                       \frac{F(\epsilon_{i+1}) + F(\epsilon_i)}{2},
+ *   @f]
+ * where @f$ i @f$ is the index of discrete energy levels, or Simpson's rule.
  *
  * @warning  This class is an experimental part of %Cantera and may be
  *           changed or removed without notice.
@@ -71,6 +96,8 @@ class PlasmaPhase: public IdealGasPhase
      */
     explicit PlasmaPhase(const string& inputFile="", const string& id="");
 
+    ~PlasmaPhase();
+
     string type() const override {
         return "plasma";
     }
@@ -81,7 +108,9 @@ class PlasmaPhase: public IdealGasPhase
     //! @param  levels The vector of electron energy levels (eV).
     //!                Length: #m_nPoints.
     //! @param  length The length of the @c levels.
-    void setElectronEnergyLevels(const double* levels, size_t length);
+    //! @param  updateEnergyDist update electron energy distribution
+    void setElectronEnergyLevels(const double* levels, size_t length,
+                                 bool updateEnergyDist=true);
 
     //! Get electron energy levels.
     //! @param  levels The vector of electron energy levels (eV). Length: #m_nPoints
@@ -99,13 +128,6 @@ class PlasmaPhase: public IdealGasPhase
                                           const double* distrb,
                                           size_t length);
 
-    //! Set discretized electron energy distribution.
-    //! @param  distrb The vector of electron energy distribution.
-    //!                Length: #m_nPoints.
-    //! @param  length The length of the vectors, which equals #m_nPoints.
-    void setDiscretizedElectronEnergyDist(const double* distrb,
-                                          size_t length);
-
     //! Get electron energy distribution.
     //! @param  distrb The vector of electron energy distribution.
     //!                Length: #m_nPoints.
@@ -198,6 +220,11 @@ class PlasmaPhase: public IdealGasPhase
         return m_nPoints;
     }
 
+    //! Number of collisions
+    size_t nCollisions() const {
+        return m_collisions.size();
+    }
+
     //! Electron Species Index
     size_t electronSpeciesIndex() const {
         return m_electronSpeciesIndex;
@@ -271,6 +298,19 @@ class PlasmaPhase: public IdealGasPhase
         return m_levelNum;
     }
 
+    virtual void setSolution(std::weak_ptr<Solution> soln) override;
+
+    /**
+     * The elastic power loss (J/s/m³)
+     *   @f[
+     *     P_k = N_A N_A C_e e \sum_k C_k K_k,
+     *   @f]
+     * where @f$ C_k @f$ and @f$ C_e @f$ are the concentration (kmol/m³) of the
+     * target species and electrons, respectively. @f$ K_k @f$ is the elastic
+     * electron energy loss coefficient (eV-m³/s).
+     */
+    double elasticPowerLoss();
+
 protected:
     void updateThermo() const override;
 
@@ -317,6 +357,12 @@ class PlasmaPhase: public IdealGasPhase
     //! Electron energy distribution norm
     void normalizeElectronEnergyDistribution();
 
+    //! Update interpolated cross section of a collision
+    bool updateInterpolatedCrossSection(size_t k);
+
+    //! Update electron energy distribution difference
+    void updateElectronEnergyDistDifference();
+
     // Electron energy order in the exponential term
     double m_isotropicShapeFactor = 2.0;
 
@@ -345,6 +391,39 @@ class PlasmaPhase: public IdealGasPhase
     //! Flag of normalizing electron energy distribution
     bool m_do_normalizeElectronEnergyDist = true;
 
+    //! Data for initiate reaction
+    AnyMap m_root;
+
+    //! Electron energy distribution Difference dF/dε (V^-5/2)
+    Eigen::ArrayXd m_electronEnergyDistDiff;
+
+    //! Elastic electron energy loss coefficients (eV m3/s)
+    /*! The elastic electron energy loss coefficient for species k is,
+     *   @f[
+     *     K_k = \frac{2 m_e}{m_k} \sqrt{\frac{2 e}{m_e}} \int_0^{\infty} \sigma_k
+     *           \epsilon^2 \left( F_0 + \frac{k_B T}{e}
+     *           \frac{\partial F_0}{\partial \epsilon} \right) d \epsilon,
+     *   @f]
+     * where @f$ m_e @f$ [kg] is the electron mass, @f$ \epsilon @f$ [V] is the
+     * electron energy, @f$ \sigma_k @f$ [m2] is the reaction collision cross section,
+     * @f$ F_0 @f$ [V^(-3/2)] is the normalized electron energy distribution function.
+     */
+    vector<double> m_elasticElectronEnergyLossCoefficients;
+
+    //! Updates the elastic electron energy loss coefficient for collision index i
+    /*! Calculates the elastic energy loss coefficient using the current electron
+        energy distribution and cross sections.
+    */
+    void updateElasticElectronEnergyLossCoefficient(size_t i);
+
+    //! Update elastic electron energy loss coefficients
+    /*! Used by elasticPowerLoss() and other plasma property calculations that
+        depends on #m_elasticElectronEnergyLossCoefficients. This function calls
+        updateInterpolatedCrossSection() before calling
+        updateElasticElectronEnergyLossCoefficient()
+    */
+    void updateElasticElectronEnergyLossCoefficients();
+
 private:
     //! Electron energy distribution change variable. Whenever
     //! #m_electronEnergyDist changes, this int is incremented.
@@ -353,6 +432,30 @@ class PlasmaPhase: public IdealGasPhase
     //! Electron energy level change variable. Whenever
     //! #m_electronEnergyLevels changes, this int is incremented.
     int m_levelNum = -1;
+
+    //! The list of shared pointers of plasma collision reactions
+    vector<shared_ptr<Reaction>> m_collisions;
+
+    //! The list of shared pointers of collision rates
+    vector<shared_ptr<ElectronCollisionPlasmaRate>> m_collisionRates;
+
+    //! The collision-target species indices of #m_collisions
+    vector<size_t> m_targetSpeciesIndices;
+
+    //! Interpolated cross sections. This is used for storing
+    //! interpolated cross sections temporarily.
+    vector<double> m_interp_cs;
+
+    //! The list of whether the interpolated cross sections is ready
+    vector<bool> m_interp_cs_ready;
+
+    //! Set collisions. This function sets the list of collisions and
+    //! the list of target species using #addCollision.
+    void setCollisions();
+
+    //! Add a collision and record the target species
+    void addCollision(std::shared_ptr<Reaction> collision);
+
 };
 
 }

include/cantera/thermo/ThermoPhase.h

@@ -15,6 +15,7 @@
 #include "MultiSpeciesThermo.h"
 #include "cantera/base/Units.h"
 #include "cantera/base/AnyMap.h"
+#include "cantera/base/Solution.h"
 
 namespace Cantera
 {
@@ -1957,6 +1958,12 @@ class ThermoPhase : public Phase
 
     //! @}
 
+    //! Set the link to the Solution object that owns this ThermoPhase
+    //! @param soln Weak pointer to the parent Solution object
+    virtual void setSolution(std::weak_ptr<Solution> soln) {
+        m_soln = soln;
+    }
+
 protected:
     //! Store the parameters of a ThermoPhase object such that an identical
     //! one could be reconstructed using the newThermo(AnyMap&) function. This
@@ -1993,6 +2000,9 @@ class ThermoPhase : public Phase
 
     //! last value of the temperature processed by reference state
     mutable double m_tlast = 0.0;
+
+    //! reference to Solution
+    std::weak_ptr<Solution> m_soln;
 };
 
 }

interfaces/cython/cantera/thermo.pxd

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

interfaces/cython/cantera/thermo.pyx

@@ -1904,6 +1904,15 @@ cdef class ThermoPhase(_SolutionBase):
                 raise ThermoModelMethodError(self.thermo_model)
             return pystr(self.plasma.electronSpeciesName())
 
+    property elastic_power_loss:
+        """
+        Elastic power loss (J/s/m3)
+        .. versionadded:: 3.2
+        """
+        def __get__(self):
+            if not self._enable_plasma:
+                raise ThermoModelMethodError(self.thermo_model)
+            return self.plasma.elasticPowerLoss()
 
 cdef class InterfacePhase(ThermoPhase):
     """ A class representing a surface, edge phase """