More photolysis API (#3)

- Complete YAML input for photolysis x-section
- A way to modify reaction stoichiometric coefficient

TODO:
- [ ]  Add example VULCAN mechanism
- [ ]  Connecting reaction stoichiometric coefficient

include/cantera/kinetics/BulkKinetics.h

@@ -92,7 +92,14 @@ class BulkKinetics : public Kinetics
     }
 
     //! @}
+
+    //! @name Photolysis methods
+    //! @{
+
+    void updateActinicFlux(void *rt_solver) override;
 
+    //! @}
+
 protected:
     //! @name Internal service methods
     //!

include/cantera/kinetics/Kinetics.h

@@ -1398,6 +1398,62 @@ class Kinetics
         return m_root.lock();
     }
 
+    //! @name Photolysis calculation methods
+    //! @{
+
+    size_t nWavelengths() const {
+        return m_wavelength.size();
+    }
+
+    /**
+     * Set the wavelengths at which the actinic flux is calculated.
+     */
+    void setWavelength(double const* wavelength, size_t n) {
+        m_wavelength.assign(wavelength, wavelength + n);
+        m_actinicFlux.resize(n);
+        std::fill(m_actinicFlux.begin(), m_actinicFlux.end(), 0.0);
+    }
+
+    /**
+     * Get the wavelengths at which the actinic flux is calculated.
+     */
+    void getWavelength(double* wavelength) const {
+        std::copy(m_wavelength.begin(), m_wavelength.end(), wavelength);
+    }
+
+    /**
+     * Update the actinic flux for each wavelength.
+     */
+    virtual void updateActinicFlux(void *rt_solver) {
+        throw NotImplementedError("Kinetics::updateActinicFlux");
+    }
+
+    /**
+     * Check whether the actinic flux has been updated since the last call to
+     * #updateActinicFlux.
+     */
+    bool hasNewActinicFlux() const {
+        return m_hasNewActinicFlux;
+    }
+
+    /**
+     * Get the actinic flux for each wavelength. The actinic flux is the
+     * integral of the photon flux over all wavelengths, and is used to
+     * calculate the photolysis rates of reactions.
+     */
+    void getActinicFlux(double *actinic_flux) const {
+        std::copy(m_actinicFlux.begin(), m_actinicFlux.end(), actinic_flux);
+    }
+
+    /**
+     * Modify the stoichiometric coefficient of the product in a reaction.
+     * The set of products must not be changed.
+     * Only the stoichiometric coefficient is changed.
+     */
+    virtual void modifyProductStoichiometry(size_t i, Composition const& comp);
+
+    //! @}
+
 protected:
     //! Cache for saved calculations within each Kinetics object.
     ValueCache m_cache;
@@ -1518,6 +1574,15 @@ class Kinetics
 
     //! reference to Solution
     std::weak_ptr<Solution> m_root;
+
+    //! Photon wavelengths
+    vector<double> m_wavelength;
+
+    //! Photon actinic fluxes
+    vector<double> m_actinicFlux;
+
+    //! Flag indicating whether the actinic fluxes have been updated
+    bool m_hasNewActinicFlux = false;
 };
 
 }

include/cantera/kinetics/Photolysis.h

@@ -9,17 +9,13 @@
 #include "ReactionRate.h"
 #include "MultiRate.h"
 
-inline int locate(double const *xx, double x, int n) {
-  return 0;
-}
-
-inline void interpn(double *val, double const *coor, double const *data,
-             double const *axis, size_t const *len, int ndim, int nval) {
-}
-
 namespace Cantera
 {
 
+int locate(double const *xx, double x, int n);
+void interpn(double *val, double const *coor, double const *data,
+             double const *axis, size_t const *len, int ndim, int nval);
+
 class ThermoPhase;
 class Kinetics;
 
@@ -130,9 +126,25 @@ class PhotolysisRate : public PhotolysisBase {
   }
 
   double evalFromStruct(PhotolysisData const& data) {
+    if (m_crossSection.empty()) {
+      return 0.;
+    }
+
     double wmin = m_temp_wave_grid[m_ntemp];
     double wmax = m_temp_wave_grid.back();
 
+    if (wmin > data.wavelength.front()) {
+      throw CanteraError("PhotolysisRate::evalFromStruct",
+                         "Wavelength out of range: {} nm < {} nm",
+                         wmin, data.wavelength.front());
+    }
+
+    if (wmax < data.wavelength.back()) {
+      throw CanteraError("PhotolysisRate::evalFromStruct",
+                         "Wavelength out of range: {} nm > {} nm",
+                         wmax, data.wavelength.back());
+    }
+
     int iwmin = locate(data.wavelength.data(), wmin, data.wavelength.size());
     int iwmax = locate(data.wavelength.data(), wmax, data.wavelength.size());
 
@@ -150,7 +162,7 @@ class PhotolysisRate : public PhotolysisBase {
       stoich = 0.0;
 
     for (int i = iwmin; i < iwmax; i++) {
-      coord[1] = data.actinicFlux[i+1];
+      coord[1] = data.wavelength[i+1];
       interpn(cross2, coord, m_crossSection.data(), m_temp_wave_grid.data(),
           len, 2, m_branch.size());
 

include/cantera/kinetics/StoichManager.h

@@ -690,6 +690,7 @@ class StoichManagerN
         }
         if (frac || k.size() > 3) {
             m_cn_list.emplace_back(rxn, k, order, stoich);
+            m_rxn_in_cnlist[rxn] = m_cn_list.size() - 1;
         } else {
             // Try to express the reaction with unity stoichiometric
             // coefficients (by repeating species when necessary) so that the
@@ -714,6 +715,7 @@ class StoichManagerN
                 break;
             default:
                 m_cn_list.emplace_back(rxn, k, order, stoich);
+                m_rxn_in_cnlist[rxn] = m_cn_list.size() - 1;
             }
         }
         m_ready = false;
@@ -800,6 +802,25 @@ class StoichManagerN
         _scale(m_cn_list.begin(), m_cn_list.end(), in, out, factor);
     }
 
+    /** 
+     * Modify reaction stoichiometry coefficients
+     * This method is limited to reactions in #m_cn_list (fractional stoichiometry)
+     */
+    void modify(size_t rxn, const vector<size_t>& k, const vector<double>& order,
+             const vector<double>& stoich)
+    {
+        auto it = m_rxn_in_cnlist.find(rxn);
+        if (it == m_rxn_in_cnlist.end()) {
+            throw CanteraError("StoichManagerN::modify", "Reaction index {} "
+                "is not in the list of reactions with fractional stoichiometry.",
+                rxn);
+        }
+        m_cn_list[it->second] = C_AnyN(rxn, k, order, stoich);
+
+        for (size_t n = 0; n < stoich.size(); n++)
+            m_stoichCoeffs.coeffRef(k[n], rxn) = stoich[n];
+    }
+
 private:
     bool m_ready; //!< Boolean flag indicating whether object is fully configured
 
@@ -816,6 +837,9 @@ class StoichManagerN
     vector<int> m_outerIndices;
     vector<int> m_innerIndices;
     vector<double> m_values;
+
+    //! Map of reaction indices to indices in #m_cn_list
+    map<int, int> m_rxn_in_cnlist;
 };
 
 }

src/kinetics/BulkKinetics.cpp

@@ -700,4 +700,11 @@ void BulkKinetics::assertDerivativesValid(const string& name)
     }
 }
 
+void __attribute__((weak)) BulkKinetics::updateActinicFlux(void *rt_solver)
+{
+    double *flux = reinterpret_cast<double*>(rt_solver);
+    m_actinicFlux.assign(flux, flux + m_wavelength.size());
+    m_hasNewActinicFlux = true;
+}
+
 }

src/kinetics/Kinetics.cpp

@@ -693,4 +693,25 @@ shared_ptr<const Reaction> Kinetics::reaction(size_t i) const
     return m_reactions[i];
 }
 
+void Kinetics::modifyProductStoichiometry(size_t irxn, Composition const& comp)
+{
+    checkReactionIndex(irxn);
+
+    vector<size_t> pk;
+    vector<double> pstoich;
+
+    for (const auto& [name, stoich] : comp) {
+        pk.push_back(kineticsSpeciesIndex(name));
+        pstoich.push_back(stoich);
+    }
+
+    m_reactions[irxn]->products = comp;
+    m_productStoich.modify(irxn, pk, pstoich, pstoich);
+
+    m_stoichMatrix = m_productStoich.stoichCoeffs();
+    m_stoichMatrix -= m_reactantStoich.stoichCoeffs();
+
+    invalidateCache();
+}
+
 }

src/kinetics/Photolysis.cpp

@@ -1,5 +1,6 @@
 //! @file Photolysis.cpp
 
+#include "cantera/kinetics/Kinetics.h"
 #include "cantera/kinetics/Photolysis.h"
 #include "cantera/thermo/ThermoPhase.h"
 #include "cantera/base/stringUtils.h"
@@ -16,17 +17,26 @@ bool PhotolysisData::update(const ThermoPhase& thermo, const Kinetics& kin)
     changed = true;
   }
 
+  if (wavelength.empty()) {
+    size_t nwave = kin.nWavelengths();
+
+    wavelength.resize(nwave);
+    actinicFlux.resize(nwave);
+
+    kin.getWavelength(wavelength.data());
+    kin.getActinicFlux(actinicFlux.data());
+    changed = true;
+  } else if (kin.hasNewActinicFlux()) {
+    kin.getActinicFlux(actinicFlux.data());
+    changed = true;
+  }
+
   return changed;
 }
 
 bool PhotolysisData::check() const
 {
     // Check that the wavelength grid is valid
-    if (wavelength.empty()) {
-        throw CanteraError("PhotolysisData::update",
-                           "Wavelength grid is empty.");
-    }
-
     if (wavelength.size() < 2) {
         throw CanteraError("PhotolysisData::update",
                            "Wavelength grid must have at least two points.");
@@ -195,7 +205,18 @@ void PhotolysisBase::setParameters(AnyMap const& node, UnitStack const& rate_uni
   m_ntemp = temperature.size();
   m_nwave = wavelength.size();
   m_temp_wave_grid.resize(m_ntemp + m_nwave);
+
+  for (size_t i = 0; i < m_ntemp; i++) {
+    m_temp_wave_grid[i] = temperature[i];
+  }
+
+  for (size_t i = 0; i < m_nwave; i++) {
+    m_temp_wave_grid[m_ntemp + i] = wavelength[i];
+  }
+
+  // only works for one temperature range
   m_crossSection = xsection;
+  m_crossSection.insert(m_crossSection.end(), xsection.begin(), xsection.end());
 
   if (node.hasKey("rate-constant")) {
     setRateParameters(node["rate-constant"], branches);

src/kinetics/interpn.cpp

@@ -0,0 +1,91 @@
+#include <stdlib.h>
+
+namespace Cantera
+{
+
+/*! Given an array xx[0..n-1] , and given a value x ,
+ * returns a value j such that x is between xx[j] and xx[j+1].
+ * xx must be monotonic, either increasing or decreasing.
+ * j=0 or j=n is returned to indicate that x is out of range.
+ * adapted from Numerical Recipes in C, 2nd Ed., p. 117.
+ */
+int locate(double const *xx, double x, int n) {
+  xx -= 1;  // zero-offset to unit-offset
+
+  int j;
+  int ju, jm, jl;
+  int ascnd;
+
+  jl = 0;
+  ju = n + 1;
+  ascnd = (xx[n] >= xx[1]);
+  while (ju - jl > 1) {
+    jm = (ju + jl) >> 1;
+    if (x >= xx[jm] == ascnd)
+      jl = jm;
+    else
+      ju = jm;
+  }
+  if (x == xx[1])
+    j = 1;
+  else if (x == xx[n])
+    j = n - 1;
+  else
+    j = jl;
+
+  j -= 1;  // unit-offset to zero-offset
+  return j;
+}
+
+/*! Multidimensional linear interpolation
+ * val[0..nval-1]   : output values
+ * coor[0..ndim-1]  : coordinate of the interpolation point
+ * data[...]        : points to the start position of a multidimensional data
+ * table. len[0..ndim-1]   : length of each dimension axis[...]        :
+ * coordinates of each dimesnion is placed sequentially in axis
+ */
+void interpn(double *val, double const *coor, double const *data,
+             double const *axis, size_t const *len, int ndim, int nval) {
+  int i1, i2;
+  i1 = locate(axis, *coor, *len);
+
+  // if the interpolation value is out of bound
+  // use the closest value
+  if (i1 == -1) {
+    i1 = 0;
+    i2 = 0;
+  } else if (i1 == *len - 1) {
+    i1 = *len - 1;
+    i2 = *len - 1;
+  } else
+    i2 = i1 + 1;
+
+  double *v1 = (double *)malloc(nval * sizeof(double));
+  double *v2 = (double *)malloc(nval * sizeof(double));
+
+  double x1 = axis[i1];
+  double x2 = axis[i2];
+
+  if (ndim == 1) {
+    for (int j = 0; j < nval; ++j) {
+      v1[j] = data[i1 * nval + j];
+      v2[j] = data[i2 * nval + j];
+    }
+  } else {
+    int s = nval;
+    for (int j = 1; j < ndim; ++j) s *= len[j];
+    interpn(v1, coor + 1, data + i1 * s, axis + *len, len + 1, ndim - 1, nval);
+    interpn(v2, coor + 1, data + i2 * s, axis + *len, len + 1, ndim - 1, nval);
+  }
+
+  if (x2 != x1)
+    for (int j = 0; j < nval; ++j)
+      val[j] = ((*coor - x1) * v2[j] + (x2 - *coor) * v1[j]) / (x2 - x1);
+  else
+    for (int j = 0; j < nval; ++j) val[j] = (v1[j] + v2[j]) / 2.;
+
+  free(v1);
+  free(v2);
+}
+
+}  // namespace Cantera

test/data/ch4_photolysis.yaml

@@ -83,5 +83,6 @@ reactions:
     - format: YAML
       temperature-range: [0., 300.]
       data:
-        - [0., 1.e-18]
-        - [100., 2.e-18]
+        - [200., 1.e-4]
+        - [250., 2.e-4]
+        - [300., 3.e-4]