[plasma] make changing reaction and rate update interpolated cross se…

…ctions

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
@@ -127,24 +129,22 @@ class ElectronCollisionPlasmaRate : public ReactionRate
     }
 
     //! The value of #m_energyLevels [eV]
-    const vector<double>& energyLevels() const {
+    virtual const vector<double>& energyLevels() const override {
         return m_energyLevels;
     }
 
     //! The value of #m_crossSections [m2]
-    const vector<double>& crossSections() const {
+    virtual const vector<double>& crossSections() const override {
         return m_crossSections;
     }
 
     //! The value of #m_crossSectionsInterpolated [m2]
-    const vector<double>& crossSectionInterpolated() const {
+    virtual const vector<double>& crossSectionInterpolated() const override {
         return m_crossSectionsInterpolated;
     }
 
-    //! Set the value of #m_crossSectionsInterpolated [m2]
-    void setCrossSectionInterpolated(vector<double>& cs) {
-        m_crossSectionsInterpolated = cs;
-    }
+    //! Update the value of #m_crossSectionsInterpolated [m2]
+    virtual void updateInterpolatedCrossSection(const vector<double>&) override;
 
 private:
     //! electron energy levels [eV]

include/cantera/kinetics/Kinetics.h

@@ -1407,6 +1407,26 @@ class Kinetics
         return m_root.lock();
     }
 
+    //! Register a function to be called if reaction is changed.
+    //! @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 changed.
+    //! When the callback becomes invalid (for example, the corresponding object is
+    //! being deleted, the removeReactionChangedCallback() method must be invoked.
+    //! @since New in %Cantera 3.1
+    void registerReactionChangedCallback(void* id, const function<void()>& callback)
+    {
+        m_reactionChangedCallbacks[id] = callback;
+    }
+
+    //! Remove the reaction-changed callback function associated with the specified object.
+    //! @since New in %Cantera 3.1
+    void removeReactionChangedCallback(void* id)
+    {
+        m_reactionChangedCallbacks.erase(id);
+    }
+
 protected:
     //! Cache for saved calculations within each Kinetics object.
     ValueCache m_cache;
@@ -1530,6 +1550,9 @@ class Kinetics
 
     //! reference to Solution
     std::weak_ptr<Solution> m_root;
+
+    //! Callback functions that are invoked when the reaction is changed.
+    map<void*, function<void()>> m_reactionChangedCallbacks;
 };
 
 }

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.1
+    void registerSetRateCallback(void* id, const function<void()>& callback);
+
+    //! Remove the setRate callback function associated with the specified object.
+    //! @since New in %Cantera 3.1
+    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/kinetics/ReactionRate.h

@@ -216,6 +216,35 @@ class ReactionRate
         m_composition_dependent_rate = comp_dep;
     }
 
+    //! @name Properties of the electron-collision-plasma rate
+    //! @{
+
+    //! The value of #ElectronCollisionPlasmaRate::m_energyLevels [eV]
+    virtual const vector<double>& energyLevels() const {
+        throw NotImplementedError("ReactionRate::energyLevels",
+                                  "Not implemented by '{}' object.", type());
+    }
+
+    //! The value of #ElectronCollisionPlasmaRate::m_crossSections [m2]
+    virtual const vector<double>& crossSections() const {
+        throw NotImplementedError("ReactionRate::crossSections",
+                                  "Not implemented by '{}' object.", type());
+    }
+
+    //! The value of #ElectronCollisionPlasmaRate::m_crossSectionsInterpolated
+    virtual const vector<double>& crossSectionInterpolated() const {
+        throw NotImplementedError("ReactionRate::crossSectionInterpolated",
+                                  "Not implemented by '{}' object.", type());
+    }
+
+    //! Set the value of #ElectronCollisionPlasmaRate::m_crossSectionsInterpolated
+    virtual void updateInterpolatedCrossSection(const vector<double>& sharedLevels) {
+        throw NotImplementedError("ReactionRate::updateInterpolatedCrossSection",
+                                  "Not implemented by '{}' object.", type());
+    }
+
+    //! @}
+
 protected:
     //! Get parameters
     //! @param node  AnyMap containing rate information

include/cantera/thermo/PlasmaPhase.h

@@ -94,6 +94,8 @@ class PlasmaPhase: public IdealGasPhase
      */
     explicit PlasmaPhase(const string& inputFile="", const string& id="");
 
+    ~PlasmaPhase();
+
     string type() const override {
         return "plasma";
     }
@@ -300,11 +302,6 @@ class PlasmaPhase: public IdealGasPhase
         return m_levelNum;
     }
 
-    //! Return the interpolated cross section Number #m_csNum
-    int interpCrossSectionNum() const {
-        return m_csNum;
-    }
-
     //! add Solution object
     virtual void addSolution(std::weak_ptr<Solution> soln);
 
@@ -332,10 +329,6 @@ class PlasmaPhase: public IdealGasPhase
     //! the new level.
     void electronEnergyLevelChanged();
 
-    //! When the interpolated cross sections changed, plasma properties such as
-    //! elastic power loss need to be re-evaluated.
-    void interpolatedCrossSectionsChanged();
-
     //! Check the electron energy levels
     /*!
      *  The values of electron energy levels need to be positive and
@@ -369,8 +362,8 @@ class PlasmaPhase: public IdealGasPhase
     //! Electron energy distribution norm
     void normalizeElectronEnergyDistribution();
 
-    //! Update interpolated cross sections
-    void updateInterpolatedCrossSections();
+    //! Update interpolated cross section of a collision
+    bool updateInterpolatedCrossSection(size_t k);
 
     //! Update electron energy distribution difference
     void updateElectronEnergyDistDifference();
@@ -422,6 +415,9 @@ class PlasmaPhase: public IdealGasPhase
      */
     vector<double> m_elasticElectronEnergyLossCoefficients;
 
+    //! update the elastic electron energy loss coefficient of i collision
+    void updateElasticElectronEnergyLossCoefficient(size_t i);
+
     //! update elastic electron energy loss coefficients
     void updateElasticElectronEnergyLossCoefficients();
 
@@ -434,27 +430,26 @@ class PlasmaPhase: public IdealGasPhase
     //! #m_electronEnergyLevels changes, this int is incremented.
     int m_levelNum = -1;
 
-    //! Interpolated collision cross section change variable. Whenever
-    //! #ElectronCollisionPlasmaRate::m_crossSectionsInterpolated changes,
-    //! this int is incremented.
-    int m_csNum = -1;
-
     //! The list of shared pointers of plasma collision reactions
     vector<shared_ptr<Reaction>> m_collisions;
 
-    //! The list of shared pointers of ElectronCollisionPlasmaRate
-    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;
 
-    //! Set collisions
+    //! 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);
+
 };
 
 }

src/kinetics/ElectronCollisionPlasmaRate.cpp

@@ -70,26 +70,23 @@ void ElectronCollisionPlasmaRate::getParameters(AnyMap& node) const {
     node["cross-sections"] = m_crossSections;
 }
 
+void ElectronCollisionPlasmaRate::updateInterpolatedCrossSection(
+    const vector<double>& sharedLevels) {
+    m_crossSectionsInterpolated.clear();
+    m_crossSectionsInterpolated.reserve(sharedLevels.size());
+    for (double level : sharedLevels) {
+        m_crossSectionsInterpolated.emplace_back(linearInterp(level,
+                                            m_energyLevels, m_crossSections));
+    }
+}
+
 double ElectronCollisionPlasmaRate::evalFromStruct(
     const ElectronCollisionPlasmaData& shared_data)
 {
     // Interpolate cross-sections data to the energy levels of
     // the electron energy distribution function
     if (shared_data.levelChanged) {
-        m_crossSectionsInterpolated.clear();
-        m_crossSectionsInterpolated.reserve(shared_data.energyLevels.size());
-        for (double level : shared_data.energyLevels) {
-            m_crossSectionsInterpolated.emplace_back(linearInterp(level,
-                                                m_energyLevels, m_crossSections));
-        }
-    } else {
-        // The interpolated cross section is set by the PlasmaPhase
-        // check the size of the vector
-        if (m_crossSectionsInterpolated.size() != shared_data.energyLevels.size()) {
-            throw CanteraError("ElectronCollisionPlasmaRate::evalFromStruct",
-                               "Energy levels and interpolated cross section ",
-                               "must have the same length.");
-        }
+        updateInterpolatedCrossSection(shared_data.energyLevels);
     }
 
     // Map cross sections to Eigen::ArrayXd

src/kinetics/Kinetics.cpp

@@ -703,13 +703,24 @@ bool Kinetics::addReaction(shared_ptr<Reaction> r, bool resize)
         m_ready = false;
     }
 
+    for (const auto& [id, callback] : m_reactionChangedCallbacks) {
+        callback();
+    }
+
     return true;
 }
 
 void Kinetics::modifyReaction(size_t i, shared_ptr<Reaction> rNew)
 {
     checkReactionIndex(i);
     shared_ptr<Reaction>& rOld = m_reactions[i];
+
+    if (rNew->rate()->type() == "electron-collision-plasma") {
+        throw CanteraError("Kinetics::modifyReaction",
+            "Type electron-collision-plasma is not supported. "
+            "Use the rate object of the reaction to modify the data.");
+    }
+
     if (rNew->type() != rOld->type()) {
         throw CanteraError("Kinetics::modifyReaction",
             "Reaction types are different: {} != {}.",

src/kinetics/Reaction.cpp

@@ -316,6 +316,19 @@ void Reaction::setRate(shared_ptr<ReactionRate> rate)
             "Reaction rate for reaction '{}' must not be empty.", equation());
     }
     m_rate = rate;
+    for (const auto& [id, callback] : m_setRateCallbacks) {
+        callback();
+    }
+}
+
+void Reaction::registerSetRateCallback(void* id, const function<void()>& callback)
+{
+    m_setRateCallbacks[id] = callback;
+}
+
+void Reaction::removeSetRateCallback(void* id)
+{
+    m_setRateCallbacks.erase(id);
 }
 
 string Reaction::reactantString() const