Extend the Python interface to expose Reaction & Rate Law objects

CMakeLists.txt

@@ -51,7 +51,7 @@ endif()
 #######################################################################
 
 if (SKBUILD)
-    add_subdirectory(thirdparty/pybind11)
+    find_package(pybind11 REQUIRED CONFIG)
     add_subdirectory(interface/python)
     set(CMAKE_POSITION_INDEPENDENT_CODE ON)
 endif()

interface/python/CMakeLists.txt

@@ -1,5 +1,6 @@
 pybind11_add_module(_mutationpp
         src/mutationpp_python.cpp
+	src/pyReaction.cpp
         src/pyMixtureOptions.cpp
         src/pyMixture.cpp
 )

interface/python/src/mutationpp_python.cpp

@@ -8,11 +8,15 @@
 
 namespace py = pybind11;
 
+void py_export_RateLaw(py::module &);
+void py_export_Reaction(py::module &);
 void py_export_MixtureOptions(py::module &);
 void py_export_Mixture(py::module &);
 
 PYBIND11_MODULE(_mutationpp, m) {
     m.doc() = "Mutation++ Python bindings";
+    py_export_RateLaw(m);
+    py_export_Reaction(m);
     py_export_MixtureOptions(m);
     py_export_Mixture(m);
 }

interface/python/src/pyMixture.cpp

@@ -64,36 +64,36 @@ void py_export_Mixture(py::module &m) {
           },
            "Returns the Stefan-Boltzmann constant (W/m^2-K^4)")  
 
-      .def("nElements", &Mutation::Mixture::nElements,
+      .def_property_readonly("num_elements", &Mutation::Mixture::nElements,
            "Returns the number of elements considered in the mixture.")
 
-      .def("nAtoms", &Mutation::Mixture::nAtoms,
+      .def_property_readonly("num_atoms", &Mutation::Mixture::nAtoms,
            "Returns the number of atomic species in the mixture.")
 
-      .def("nMolecules", &Mutation::Mixture::nMolecules,
+      .def_property_readonly("num_molecules", &Mutation::Mixture::nMolecules,
            "Returns the number of molecules in the mixture.")
 
-      .def("nHeavy", &Mutation::Mixture::nHeavy,
+      .def_property_readonly("num_heavy", &Mutation::Mixture::nHeavy,
            "Returns the number of heavy particles (non electrons) in the "
            "mixture.")
 
-      .def("nSpecies", &Mutation::Mixture::nSpecies,
+      .def_property_readonly("num_species", &Mutation::Mixture::nSpecies,
            "Returns the number of species considered in the mixture.")
 
-      .def("nPhases", &Mutation::Mixture::nPhases,
+      .def_property_readonly("num_phases", &Mutation::Mixture::nPhases,
            "Returns the number of phases belonging to this mixture.")
 
-      .def("nGas", &Mutation::Mixture::nGas,
+      .def_property_readonly("num_gas", &Mutation::Mixture::nGas,
            "Returns number of gas species in the mixture.")
 
-      .def("nCondensed", &Mutation::Mixture::nCondensed,
+      .def_property_readonly("numcondensed", &Mutation::Mixture::nCondensed,
            "Returns the number of condensed phase species in the mixture.")
 
-      .def("nEnergyEqns", &Mutation::Mixture::nEnergyEqns,
+      .def_property_readonly("num_energy_eqns", &Mutation::Mixture::nEnergyEqns,
            "Returns the number of energy equations associated with the mixture "
            "StateModel.")
 
-      .def("nMassEqns", &Mutation::Mixture::nMassEqns,
+      .def_property_readonly("num_mass_eqns", &Mutation::Mixture::nMassEqns,
            "Returns the number of mass equations associated with the mixture "
            "StateModel.")
 
@@ -105,6 +105,28 @@ void py_export_Mixture(py::module &m) {
       .def("hasElectrons", &Mutation::Mixture::hasElectrons,
            "Returns true if this mixture includes electrons, false otherwise.")
 
+      .def_property_readonly("reactions",
+			     &Mutation::Mixture::reactions,
+			     "Returns all reactions in this mixture.")
+
+      .def("reaction",
+	   [](const Mutation::Mixture &self, size_t idx) {
+	     return self.reactions()[idx];
+	   },
+	   "Returns reaction with index idx.")
+
+      .def_property_readonly("reactants",
+	   [](const Mutation::Mixture &self, size_t idx) {
+	     return self.reactions()[idx].reactants();
+	   },
+	   "Returns the reactants for reaction index idx.")
+
+      .def_property_readonly("products",
+	   [](const Mutation::Mixture &self, size_t idx) {
+	     return self.reactions()[idx].products();
+	   },
+	   "Returns the reactants for reaction index idx.")
+
       .def("setState", &Mutation::Mixture::setState,
            "Sets the state of the mixture using the StateModel belonging to "
            "the mixture."
@@ -590,7 +612,7 @@ void py_export_Mixture(py::module &m) {
           "model is an"
           " equilibrium one.")
 
-      .def("nReactions", &Mutation::Mixture::nReactions,
+      .def_property_readonly("num_reactions", &Mutation::Mixture::nReactions,
            "Returns the number of reactions in the mechanism.")
 
       .def(
@@ -898,4 +920,4 @@ void py_export_Mixture(py::module &m) {
           "Return an array of gibbs free energy for each species per unit mass at standard pressure")     
 
           ;
-}
+}

interface/python/src/pyReaction.cpp

@@ -0,0 +1,50 @@
+#include <Reaction.h>
+#include <RateManager.h>
+#include <pybind11/numpy.h>
+#include <pybind11/stl.h>
+
+namespace py = pybind11;
+
+/**
+ * Python wrapper for the Reaction class.
+ */
+void py_export_RateLaw(py::module &m) {
+  py::class_<Mutation::Kinetics::Arrhenius>(m, "Arrhenius")
+    .def_property_readonly("log_pre_exponential",
+			   &Mutation::Kinetics::Arrhenius::logA,
+			   "Returns the logarithm of the pre-exponential for this rate law.")
+    .def_property_readonly("exponent",
+			   &Mutation::Kinetics::Arrhenius::n,
+			   "Returns the temperature exponent for this rate law.")
+    .def_property_readonly("activation_temperature",
+			   &Mutation::Kinetics::Arrhenius::T,
+			   "Returns the activation temperature for this rate law.");
+}
+
+/**
+ * Python wrapper for the Reaction class.
+ */
+void py_export_Reaction(py::module &m) {
+  py::class_<Mutation::Kinetics::Reaction>(m, "Reaction")
+    .def_property_readonly("reactants",
+			   &Mutation::Kinetics::Reaction::reactants,
+			   "Returns the reactants for this reaction.")
+
+    .def_property_readonly("products",
+			   &Mutation::Kinetics::Reaction::products,
+			   "Returns the products for this reaction.")
+
+    .def("rate_law",
+	 &Mutation::Kinetics::Reaction::rateLaw,
+	 py::return_value_policy::reference_internal,
+	 "Returns the rate law for this reaction.")
+
+    .def_property_readonly("fwd_rate_coeff_temperature",
+			   &Mutation::Kinetics::Reaction::fwdRateTemperature,
+			   "Returns the temperature at which the fwd "
+			   "rate coeff. for this reaction is evaluated.")
+    .def_property_readonly("rev_rate_coeff_temperature",
+			   &Mutation::Kinetics::Reaction::revRateTemperature,
+			   "Returns the temperature at which the rev "
+			   "rate coeff. for this reaction is evaluated.");
+}

pyproject.toml

@@ -1,5 +1,12 @@
 [build-system]
-requires = ["setuptools", "wheel", "scikit-build", "cmake", "ninja"]
+requires = [
+	 "setuptools",
+	 "wheel",
+	 "scikit-build",
+	 "pybind11",
+	 "cmake",
+	 "ninja"
+]
 build-backend = "setuptools.build_meta"
 
 [tool.pytest.ini_options]

setup.py

@@ -1,6 +1,6 @@
 import re
 import sys
-
+import pybind11
 from pathlib import Path
 
 try:
@@ -52,5 +52,6 @@ def get_version_from_cmake(root_cmakelists):
     extras_require={
         "test": ["numpy", "pytest"],
     },
+    cmake_args=[f"-DCMAKE_PREFIX_PATH={pybind11.get_cmake_dir()}"],
     cmake_install_dir="interface/python/mutationpp",
 )

src/kinetics/RateLawGroup.h

@@ -42,6 +42,7 @@
 #include "Reaction.h"
 //#include "StateModel.h"
 #include "StoichiometryManager.h"
+#include "ReactionType.h"
 
 class StateModel;
 
@@ -172,6 +173,10 @@ class RateLawGroup1T : public RateLawGroup
         m_last_t = m_t;
     }
 
+  const RateLawTemperature evaluationTemperature(){
+    return TSelectorType().getT_id();
+  }
+
 private:
 
     /// vector of rates to evaluate

src/kinetics/RateLaws.cpp

@@ -84,14 +84,13 @@ Arrhenius::Arrhenius(const XmlElement& node, const int order)
     node.getAttribute("A", m_lnA, 
         "Arrhenius rate law must define coefficient A!");
     node.parseCheck(m_lnA > 0.0, "Pre-exponential factors must be positive > 0");
-
     // Convert to correct units based on the order of the reaction and
     // store the log value
     Units A_units = 
         (((sm_aunits[1]^3) / sm_aunits[0])^(order-1)) /
         (sm_aunits[2] * (sm_aunits[3]^m_n));
     m_lnA = std::log(A_units.convertToBase(m_lnA));
-    
+
     // Load the characteristic temperature
     if (node.hasAttribute("Ea")) {
         node.getAttribute("Ea", m_temp);
@@ -103,7 +102,7 @@ Arrhenius::Arrhenius(const XmlElement& node, const int order)
         m_temp = sm_eunits[2].convertToBase(m_temp);
     } else {
         node.parseError("Arrhenius rate law must define coefficient Ea or T!");
-    }
+    }    
 }
 
     } // namespace Kinetics

src/kinetics/RateLaws.h

@@ -33,6 +33,7 @@
 #include <cstdlib>
 
 #include "Utilities.h"
+#include "ReactionType.h"
 
 namespace Mutation {
     namespace Kinetics {
@@ -47,6 +48,7 @@ class RateLaw
 
     virtual ~RateLaw() { };
     virtual RateLaw* clone() const = 0;
+    virtual void set_evaluation_temperature(const RateLawTemperature eval_temp_id) {};
 };
 
 /**
@@ -61,7 +63,9 @@ class Arrhenius : public RateLaw
     Arrhenius(const Mutation::Utilities::IO::XmlElement& node, const int order);
 
     Arrhenius(const Arrhenius& to_copy)
-        : m_lnA(to_copy.m_lnA), m_n(to_copy.m_n), m_temp(to_copy.m_temp)
+      : m_lnA(to_copy.m_lnA),
+	m_n(to_copy.m_n),
+	m_temp(to_copy.m_temp)
     { }
 
     virtual ~Arrhenius() { };
@@ -78,6 +82,10 @@ class Arrhenius : public RateLaw
         return (k*invT*(m_n + m_temp*invT));
     }
 
+    double logA() const {
+      return m_lnA;
+    }
+
     double A() const { 
         return std::exp(m_lnA);
     }