Revise photolysis on cantera (#7)

- Remove the photoabsorption from the photolysis equation
- More checks regarding the equation

TODO:
- [ ] Clean up code (comments, debug statements, unused functions)
- [ ] Think about how to obtain the photoabsorption rate

include/cantera/kinetics/Photolysis.h

@@ -52,7 +52,7 @@ class PhotolysisBase : public ReactionRate {
   /*!
    * @param temp Temperature grid
    * @param wavelength Wavelength grid
-   * @param branches Branch strings of the photolysis reaction
+   * @param branches Branch strings of the photolysis products
    * @param xsection Cross-section data
    */
   PhotolysisBase(vector<double> const& temp, vector<double> const& wavelength,
@@ -78,7 +78,7 @@ class PhotolysisBase : public ReactionRate {
   void validate(const string& equation, const Kinetics& kin) override;
 
  protected:
-  //! composition of branches
+  //! composition of photolysis branch products
   vector<Composition> m_branch;
 
   //! number of temperature grid points
@@ -96,7 +96,7 @@ class PhotolysisBase : public ReactionRate {
   //! The first dimension is the number of temperature grid points, the second dimension
   //! is the number of wavelength grid points, and the third dimension is the number of
   //! branches of the photolysis reaction.
-  //! Default units are nanometers, cm^2, cm^2, and cm^2, respectively.
+  //! Default units are SI units such as m, m^2, and m^2/m.
   vector<double> m_crossSection;
 };
 
@@ -128,70 +128,27 @@ class PhotolysisRate : public PhotolysisBase {
     return m_net_products;
   }
 
-  double evalFromStruct(PhotolysisData const& data) {
-    double wmin = m_temp_wave_grid[m_ntemp];
-    double wmax = m_temp_wave_grid.back();
-
-    if (m_crossSection.empty() ||
-        wmin > data.wavelength.back() || 
-        wmax < data.wavelength.front()) 
-    {
-      return 0.;
-    }
-
-    int iwmin = locate(data.wavelength.data(), wmin, data.wavelength.size());
-    int iwmax = locate(data.wavelength.data(), wmax, data.wavelength.size());
-
-    double* cross1 = new double [m_branch.size()];
-    double* cross2 = new double [m_branch.size()];
-
-    double coord[2] = {data.temperature, data.wavelength[iwmin]};
-    size_t len[2] = {m_ntemp, m_nwave};
-
-    interpn(cross1, coord, m_crossSection.data(), m_temp_wave_grid.data(),
-        len, 2, m_branch.size());
-
-    double total_rate = 0.0;
-    for (auto const& branch : m_branch)
-        for (auto const& [name, stoich] : branch)
-            m_net_products[name] = 0.;
-
-    for (int i = iwmin; i < iwmax; i++) {
-      coord[1] = data.wavelength[i+1];
-      interpn(cross2, coord, m_crossSection.data(), m_temp_wave_grid.data(),
-          len, 2, m_branch.size());
-
-      for (size_t n = 0; n < m_branch.size(); n++) {
-        double rate = 0.5 * (data.wavelength[i+1] - data.wavelength[i])
-          * (cross1[n] * data.actinicFlux[i] + cross2[n] * data.actinicFlux[i+1]);
-        for (auto const& [name, stoich] : m_branch[n]) {
-          m_net_products.at(name) += rate * stoich;
-        }
-        total_rate += rate;
-        cross1[n] = cross2[n];
-      }
-    }
-
-    for (auto& [name, stoich] : m_net_products)
-      stoich /= total_rate;
-
-    delete [] cross1;
-    delete [] cross2;
-
-    return total_rate;
-  }
+  double evalFromStruct(PhotolysisData const& data);
 
  protected:
   //! net stoichiometric coefficients of products
   Composition m_net_products;
+
+  //! photoabsorption rate coefficient
+  double m_photoabsorption_rate;
 };
 
 /**
- * Read the cross-section data from VULCAN format files
+ * @breif Read the cross-section data from VULCAN format files
+ *
+ * @param branches Composition of the photodissociation products (no
+ * photoabsorption branch). 
  *
  * @param files Vector of filenames.
  * There are two files for each photolysis reaction. The first one is for
  * cross-section data and the second one for the branch ratios.
+ *
+ * @return a pair of vectors containing the wavelength (m) and cross section data (m^2)
  */
 pair<vector<double>, vector<double>> 
 load_xsection_vulcan(vector<string> const& files, vector<Composition> const& branches);

src/kinetics/Kinetics.cpp

@@ -626,7 +626,7 @@ bool Kinetics::addReaction(shared_ptr<Reaction> r, bool resize)
 
     m_reactantStoich.add(irxn, rk, rorder, rstoich);
     // product orders = product stoichiometric coefficients
-    bool frac = isPhotolysis(irxn) && r->products.size() > 1;
+    bool frac = r->photolysis && (r->products.size() > 1);
     m_productStoich.add(irxn, pk, pstoich, pstoich, frac);
     if (r->reversible) {
         m_revProductStoich.add(irxn, pk, pstoich, pstoich);

src/kinetics/Photolysis.cpp

@@ -2,6 +2,8 @@
 
 #include "cantera/kinetics/Kinetics.h"
 #include "cantera/kinetics/Photolysis.h"
+#include "cantera/kinetics/Reaction.h"
+#include "cantera/kinetics/ReactionRateFactory.h"
 #include "cantera/thermo/ThermoPhase.h"
 #include "cantera/base/stringUtils.h"
 
@@ -151,15 +153,34 @@ void PhotolysisBase::setParameters(AnyMap const& node, UnitStack const& rate_uni
 
   ReactionRate::setParameters(node, rate_units);
 
+  // set up a dummy reaction to parse the reaction equation
+  Reaction rtmp;
+  parseReactionEquation(rtmp, node["equation"].asString(), node, nullptr);
+
+  if (rtmp.reactants.size() != 1 || rtmp.reactants.begin()->second != 1) {
+    throw CanteraError("PhotolysisBase::setParameters",
+                       "Photolysis reaction must have one reactant with stoichiometry 1.");
+  }
+
+  // b0 is reserved for the photoabsorption cross section
+  branch_map["b0"] = 0;
+  m_branch.push_back(rtmp.reactants);
+
   if (node.hasKey("branches")) {
     for (auto const& branch : node["branches"].asVector<AnyMap>()) {
-      branch_map[branch["name"].asString()] = m_branch.size();
+      std::string branch_name = branch["name"].asString();
+
+      // check duplicated branch name
+      if (branch_map.find(branch_name) != branch_map.end()) {
+        throw CanteraError("PhotolysisBase::setParameters",
+                           "Duplicated branch name '{}'.", branch_name);
+      }
+
+      branch_map[branch_name] = m_branch.size();
       m_branch.push_back(parseCompString(branch["product"].asString()));
-    } 
-  } else {
-    vector<string> tokens;
-    tokenizeString(node["equation"].asString(), tokens);
-    m_branch.push_back(parseCompString(tokens[0] + ":1"));
+    }
+  } else if (rtmp.products != rtmp.reactants) { // this is not photoabsorption
+    m_branch.push_back(rtmp.products);
   }
 
   if (node.hasKey("cross-section")) {
@@ -216,14 +237,33 @@ void PhotolysisBase::setParameters(AnyMap const& node, UnitStack const& rate_uni
     m_temp_wave_grid[m_ntemp + i] = wavelength[i];
   }
 
-  // only works for one temperature range
+  // TODO(cli): 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"], branch_map);
   }
 
+  /* debug
+  std::cout << "number of temperature: " << m_ntemp << std::endl;
+  std::cout << "number of wavelength: " << m_nwave << std::endl;
+  std::cout << "number of branches: " << m_branch.size() << std::endl;
+  for (auto const& branch : branch_map) {
+    std::cout << "branch: " << branch.first << std::endl;
+    for (auto const& [name, stoich] : m_branch[branch.second]) {
+      std::cout << name << " " << stoich << std::endl;
+    }
+  }
+  std::cout << "number of cross-section: " << m_crossSection.size() << std::endl;
+  */
+
+  if (m_ntemp * m_nwave * m_branch.size() != m_crossSection.size()) {
+    throw CanteraError("PhotolysisBase::PhotolysisBase",
+                       "Cross-section data size does not match the temperature, "
+                       "wavelength, and branch grid sizes.");
+  }
+
   m_valid = true;
 }
 
@@ -262,6 +302,123 @@ void PhotolysisBase::validate(string const& equation, Kinetics const& kin)
     throw InputFileError("PhotolysisBase::validate", m_input,
                        "Rate object for reaction '{}' is not configured.", equation);
   }
+
+  std::vector<std::string> tokens;
+  tokenizeString(equation, tokens);
+  auto reactor_comp = parseCompString(tokens[0] + ":1");
+
+  // set up a dummy reaction to parse the reaction equation
+  Reaction rtmp;
+  parseReactionEquation(rtmp, equation, m_input, nullptr);
+
+  std::set<std::string> species_from_equation, species_from_branches;
+
+  species_from_equation.insert(rtmp.reactants.begin()->first);
+  for (auto const& [name, stoich] : rtmp.products) {
+    species_from_equation.insert(name);
+  }
+
+  species_from_branches.insert(rtmp.reactants.begin()->first);
+  for (auto const& branch : m_branch) {
+    // create a Arrhenius reaction placeholder to check balance
+    Reaction rtmp(reactor_comp, branch, newReactionRate("Arrhenius"));
+    rtmp.reversible = false;
+    rtmp.checkSpecies(kin);
+    rtmp.checkBalance(kin);
+    for (auto const& [name, stoich] : branch) {
+      species_from_branches.insert(name);
+    }
+  }
+
+  if (species_from_equation != species_from_branches) {
+    throw InputFileError("PhotolysisBase::validate", m_input,
+                       "Reaction '{}' has different products than the photolysis branches.", equation);
+  }
+}
+
+double PhotolysisRate::evalFromStruct(PhotolysisData const& data) {
+    double wmin = m_temp_wave_grid[m_ntemp];
+    double wmax = m_temp_wave_grid.back();
+
+    if (m_crossSection.empty() ||
+        wmin > data.wavelength.back() || 
+        wmax < data.wavelength.front()) 
+    {
+      m_photoabsorption_rate = 0.;
+      return 0.;
+    }
+
+    /* debug
+    std::cout << "wavelength data range: " << wmin << " " << wmax << std::endl;
+    std::cout << "temperature = " << data.temperature << std::endl;
+    std::cout << "wavelength = " << std::endl;
+    for (auto w : data.wavelength) {
+      std::cout << w << std::endl;
+    }*/
+
+    double* cross1 = new double [m_branch.size()];
+    double* cross2 = new double [m_branch.size()];
+
+    double coord[2] = {data.temperature, data.wavelength[0]};
+    size_t len[2] = {m_ntemp, m_nwave};
+
+    interpn(cross1, coord, m_crossSection.data(), m_temp_wave_grid.data(),
+        len, 2, m_branch.size());
+
+    double total_rate = 0.0;
+    // prevent division by zero
+    double eps = 1.e-30;
+    double total_rate_eps = 0.;
+
+    // first branch is photoabsorption
+    m_photoabsorption_rate = 0.;
+    for (size_t n = 1; n < m_branch.size(); n++) {
+      for (auto const& [name, stoich] : m_branch[n])
+        m_net_products[name] = 0.;
+    }
+
+    for (size_t i = 0; i < data.wavelength.size() - 1; ++i) {
+      // debug
+      //std::cout << "wavelength = " << data.wavelength[i] << " " << data.wavelength[i+1] << std::endl;
+      coord[1] = data.wavelength[i+1];
+      interpn(cross2, coord, m_crossSection.data(), m_temp_wave_grid.data(),
+          len, 2, m_branch.size());
+
+      for (size_t n = 0; n < m_branch.size(); n++) {
+        double rate = 0.5 * (data.wavelength[i+1] - data.wavelength[i])
+          * (cross1[n] * data.actinicFlux[i] + cross2[n] * data.actinicFlux[i+1]);
+
+        // debug
+        //std::cout << "cross section [ " << n << "] = " << cross1[n] << " " << cross2[n] << std::endl;
+
+        if (n == 0) { // photoabsorption
+          m_photoabsorption_rate += rate;
+        } else {  // photodissociation
+          for (auto const& [name, stoich] : m_branch[n]) {
+            m_net_products.at(name) += (rate + eps) * stoich;
+          }
+          total_rate += rate;
+          total_rate_eps += rate + eps;
+        }
+
+        cross1[n] = cross2[n];
+      }
+    }
+
+    for (auto& [name, stoich] : m_net_products)
+      stoich /= total_rate_eps;
+
+    /* debug
+    for (auto const& [name, stoich] : m_net_products)
+      std::cout << name << " " << stoich << std::endl;
+    std::cout << "photodissociation rate: " << total_rate << std::endl;
+    std::cout << "photoabsorption rate: " << m_photoabsorption_rate << std::endl;
+    */
+
+    delete [] cross1;
+    delete [] cross2;
+
+    return total_rate;
 }
 
 }

src/kinetics/load_xsection_kinetics7.cpp

@@ -1,5 +1,5 @@
 #include <cmath>
-#include <stdio.h>
+#include <cstdio>
 
 #include "cantera/base/stringUtils.h"
 #include "cantera/kinetics/Reaction.h"
@@ -28,6 +28,7 @@ load_xsection_kinetics7(vector<string> const& files, vector<Composition> const&
   vector<double> wavelength;
   vector<double> xsection;
 
+  // first cross section data is always the photoabsorption cross section (no dissociation)
   int nbranch = branches.size();
   int min_is = 9999, max_ie = 0;
 
@@ -64,6 +65,7 @@ load_xsection_kinetics7(vector<string> const& files, vector<Composition> const&
     int ncols = 7;
     int nrows = ceil(1. * nwave / ncols);
 
+    equation[60] = '\0';
     auto product = parseCompString(equation);
 
     auto it = std::find(branches.begin(), branches.end(), product);
@@ -87,8 +89,10 @@ load_xsection_kinetics7(vector<string> const& files, vector<Composition> const&
           int k = i * ncols + j;
 
           if (k >= nwave) break;
-          wavelength[k] = wave / 10.; // Angstrom to nm
-          xsection[k * nbranch + b] = cross * 10.;  // cm^2 / Angstrom to cm^2 / nm
+          // Angstrom -> m
+          wavelength[k] = wave * 1.e-10;
+          // cm^2 -> m^2
+          xsection[k * nbranch + b] = cross * 1.e-4;
         }
       }
     }
@@ -114,6 +118,15 @@ load_xsection_kinetics7(vector<string> const& files, vector<Composition> const&
   free(line);
   fclose(file);
 
+  /* debug
+  for (size_t i = 0; i < wavelength.size(); i++) {
+    printf("%g ", wavelength[i]);
+    for (int j = 0; j < nbranch; j++) {
+      printf("%g ",xsection[i * nbranch + j]);
+    }
+    printf("\n");
+  }*/
+
   return {std::move(wavelength), std::move(xsection)};
 }