From c0040ef1b194d8c0b8e28639bba27944fcdcaf7e Mon Sep 17 00:00:00 2001
From: Cheng Li <69489965+chengcli@users.noreply.github.com>
Date: Thu, 4 Jan 2024 14:59:29 -0500
Subject: [PATCH] Disable balance check for photolysis reactions (#6)

- Disable balance check for photolysis reactions when there are more than one products
- Any photolysis product must have a unit stoichiometric coefficient
---
 include/cantera/kinetics/Kinetics.h      |  2 +-
 include/cantera/kinetics/StoichManager.h |  4 ++--
 src/kinetics/Kinetics.cpp                |  3 ++-
 src/kinetics/Reaction.cpp                | 12 ++++++++++++
 test/data/ch4_photolysis.yaml            | 12 ++++++------
 5 files changed, 23 insertions(+), 10 deletions(-)

diff --git a/include/cantera/kinetics/Kinetics.h b/include/cantera/kinetics/Kinetics.h
index 759b1b080d..2f139d7992 100644
--- a/include/cantera/kinetics/Kinetics.h
+++ b/include/cantera/kinetics/Kinetics.h
@@ -1453,7 +1453,7 @@ class Kinetics
     virtual void modifyProductStoichCoeff(size_t i, Composition const& comp);
 
     virtual bool isPhotolysis(size_t i) const {
-        throw NotImplementedError("Kinetics::isPhotolysis");
+        return false;
     }
 
     //! @}
diff --git a/include/cantera/kinetics/StoichManager.h b/include/cantera/kinetics/StoichManager.h
index c80899bd56..1330091998 100644
--- a/include/cantera/kinetics/StoichManager.h
+++ b/include/cantera/kinetics/StoichManager.h
@@ -669,9 +669,10 @@ class StoichManagerN
      *     expression involving the species vector.
      *  @param stoich  This is used to handle fractional stoichiometric
      *     coefficients on the product side of irreversible reactions.
+     *  @param frac Specify that the reaction is a fractional reaction (photolysis).
      */
     void add(size_t rxn, const vector<size_t>& k, const vector<double>& order,
-             const vector<double>& stoich) {
+             const vector<double>& stoich, bool frac = false) {
         if (order.size() != k.size()) {
             throw CanteraError(
                 "StoichManagerN::add()", "size of order and species arrays differ");
@@ -680,7 +681,6 @@ class StoichManagerN
             throw CanteraError(
                 "StoichManagerN::add()", "size of stoich and species arrays differ");
         }
-        bool frac = false;
         for (size_t n = 0; n < stoich.size(); n++) {
             m_coeffList.emplace_back(
                 static_cast<int>(k[n]), static_cast<int>(rxn), stoich[n]);
diff --git a/src/kinetics/Kinetics.cpp b/src/kinetics/Kinetics.cpp
index d02b6e5d5a..e2a3692673 100644
--- a/src/kinetics/Kinetics.cpp
+++ b/src/kinetics/Kinetics.cpp
@@ -626,7 +626,8 @@ bool Kinetics::addReaction(shared_ptr<Reaction> r, bool resize)
 
     m_reactantStoich.add(irxn, rk, rorder, rstoich);
     // product orders = product stoichiometric coefficients
-    m_productStoich.add(irxn, pk, pstoich, pstoich);
+    bool frac = isPhotolysis(irxn) && r->products.size() > 1;
+    m_productStoich.add(irxn, pk, pstoich, pstoich, frac);
     if (r->reversible) {
         m_revProductStoich.add(irxn, pk, pstoich, pstoich);
     }
diff --git a/src/kinetics/Reaction.cpp b/src/kinetics/Reaction.cpp
index 6fe5ad3c0f..8df890836f 100644
--- a/src/kinetics/Reaction.cpp
+++ b/src/kinetics/Reaction.cpp
@@ -619,6 +619,18 @@ void Reaction::checkBalance(const Kinetics& kin) const
 {
     Composition balr, balp;
 
+    // disable balance check for photolysis reactions with more than one product
+    if (products.size() > 1 && type() == "Photolysis") {
+      for (const auto& [name, stoich] : products) {
+        if (stoich != 1.0) {
+          throw InputFileError("Reaction::checkBalance", input,
+              "Every product of a photolysis reaction '{}' must have "
+              "a unit stoichiometric coefficient", equation());
+        }
+      }
+      return;
+    }
+
     // iterate over products and reactants
     for (const auto& [name, stoich] : products) {
         const ThermoPhase& ph = kin.speciesPhase(name);
diff --git a/test/data/ch4_photolysis.yaml b/test/data/ch4_photolysis.yaml
index 47198cbd28..0ef5310f3e 100644
--- a/test/data/ch4_photolysis.yaml
+++ b/test/data/ch4_photolysis.yaml
@@ -55,18 +55,18 @@ species:
     model: constant-cp
 
 reactions:
-- equation: CH4 => 0.2 CH4 + 0.2 CH3 + 0.2 (1)CH2 + 0.2 (3)CH2 + 0.2 CH + 0.4 H2 + 0.8 H 
+- equation: CH4 => CH4 + CH3 + (1)CH2 + (3)CH2 + CH + H2 + H
   type: photolysis
   branches:
-    - name: b1  # 0.2
+    - name: b1
       product: "CH4:1"
-    - name: b2  # 0.2
+    - name: b2
       product: "CH3:1 H:1"
-    - name: b3  # 0.2
+    - name: b3
       product: "(1)CH2:1 H2:1"
-    - name: b4  # 0.2
+    - name: b4
       product: "(3)CH2:1 H:2"
-    - name: b5  # 0.2
+    - name: b5
       product: "CH:1 H2:1 H:1"
   cross-section:
     - format: KINETICS7