diff --git a/CMakeLists.txt b/CMakeLists.txt index 1b67d871..e8fc2fd9 100644 --- a/CMakeLists.txt +++ b/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() diff --git a/interface/python/CMakeLists.txt b/interface/python/CMakeLists.txt index 074ce029..bac2711f 100644 --- a/interface/python/CMakeLists.txt +++ b/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 ) diff --git a/interface/python/src/mutationpp_python.cpp b/interface/python/src/mutationpp_python.cpp index 05fba317..3e7285e2 100644 --- a/interface/python/src/mutationpp_python.cpp +++ b/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); } diff --git a/interface/python/src/pyMixture.cpp b/interface/python/src/pyMixture.cpp index d3e83824..2ec2b79c 100644 --- a/interface/python/src/pyMixture.cpp +++ b/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") ; -} \ No newline at end of file +} diff --git a/interface/python/src/pyReaction.cpp b/interface/python/src/pyReaction.cpp new file mode 100644 index 00000000..56bad6da --- /dev/null +++ b/interface/python/src/pyReaction.cpp @@ -0,0 +1,50 @@ +#include +#include +#include +#include + +namespace py = pybind11; + +/** + * Python wrapper for the Reaction class. + */ +void py_export_RateLaw(py::module &m) { + py::class_(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_(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."); +} diff --git a/pyproject.toml b/pyproject.toml index 7bd830f8..6a4606b8 100644 --- a/pyproject.toml +++ b/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] diff --git a/setup.py b/setup.py index 3197aef9..ae6e9fe6 100644 --- a/setup.py +++ b/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", ) diff --git a/src/kinetics/RateLawGroup.h b/src/kinetics/RateLawGroup.h index 78b7b679..796f3e15 100644 --- a/src/kinetics/RateLawGroup.h +++ b/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 diff --git a/src/kinetics/RateLaws.cpp b/src/kinetics/RateLaws.cpp index 012eaea2..d471e46b 100644 --- a/src/kinetics/RateLaws.cpp +++ b/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 diff --git a/src/kinetics/RateLaws.h b/src/kinetics/RateLaws.h index 8f9ac6b8..de80308b 100644 --- a/src/kinetics/RateLaws.h +++ b/src/kinetics/RateLaws.h @@ -33,6 +33,7 @@ #include #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); } diff --git a/src/kinetics/RateManager.cpp b/src/kinetics/RateManager.cpp index 40d4460b..ef49362e 100644 --- a/src/kinetics/RateManager.cpp +++ b/src/kinetics/RateManager.cpp @@ -35,82 +35,44 @@ namespace Mutation { namespace Kinetics { -//============================================================================== - -// Simple macro to create a temperature selector type -#define TEMPERATURE_SELECTOR(__NAME__,__T__)\ -class __NAME__\ -{\ -public:\ - inline double getT(const Thermodynamics::StateModel* const state) const {\ - return ( __T__ );\ - }\ +// // Simple macro to create a temperature selector type +// #define TEMPERATURE_SELECTOR(__NAME__,__T__,__Tname__) \ +// class __NAME__\ +// {\ +// public:\ +// inline double getT(const Thermodynamics::StateModel* const state) const {\ +// return ( __T__ );\ +// }\ +// inline std::string getT_name() const {\ +// return ( __Tname__ );\ +// }\ +// }; + +// /// Temperature selector which returns the current translational temperature +// TEMPERATURE_SELECTOR(TSelector, state->T(), "translational") + +// /// Temperature selector which returns the current electron temperature +// //TEMPERATURE_SELECTOR(TeSelector, std::min(state->Te(), 10000.0)) +// TEMPERATURE_SELECTOR(TeSelector, state->Te(), "electron") + +// /// Temperature selector which returns the current value of sqrt(T*Tv) +// TEMPERATURE_SELECTOR(ParkSelector, std::sqrt(state->T()*state->Tv()), "geometric") + +// #undef TEMPERATURE_SELECTOR + +inline double TSelector::getT(const Thermodynamics::StateModel* const state) const +{ + return state->T(); }; - -/// Temperature selector which returns the current translational temperature -TEMPERATURE_SELECTOR(TSelector, state->T()) - -/// Temperature selector which returns the current electron temperature -//TEMPERATURE_SELECTOR(TeSelector, std::min(state->Te(), 10000.0)) -TEMPERATURE_SELECTOR(TeSelector, state->Te()) - -/// Temperature selector which returns the current value of sqrt(T*Tv) -TEMPERATURE_SELECTOR(ParkSelector, std::sqrt(state->T()*state->Tv())) - -#undef TEMPERATURE_SELECTOR - -/// Arrhenius group evaluated at T -typedef RateLawGroup1T ArrheniusT; - -/// Arrhenius group evaluated at Te -typedef RateLawGroup1T ArrheniusTe; - -/// Arrhenius group evaluated at sqrt(T*Tv) -typedef RateLawGroup1T ArrheniusPark; - -//============================================================================== - -/** - * Used to define which rate law groups the forward and reverse rate laws are - * evaluated in for a given ReactionType value. The default is an Arrhenius - * rate law with Tf = Tb = T. - */ -template -struct RateSelector { - typedef ArrheniusT ForwardGroup; - typedef ArrheniusT ReverseGroup; +inline double TeSelector::getT(const Thermodynamics::StateModel* const state) const +{ + return state->Te(); }; - -#define SELECT_RATE_LAWS(__TYPE__,__FORWARD__,__REVERSE__)\ -template <> struct RateSelector<__TYPE__> {\ - typedef __FORWARD__ ForwardGroup;\ - typedef __REVERSE__ ReverseGroup;\ +inline double ParkSelector::getT(const Thermodynamics::StateModel* const state) const +{ + return std::sqrt(state->T()*state->Tv()); }; -// Default rate law groups for non (kf(T), kb(T)) reaction types -SELECT_RATE_LAWS(ASSOCIATIVE_IONIZATION, ArrheniusT, ArrheniusTe) -SELECT_RATE_LAWS(DISSOCIATIVE_RECOMBINATION, ArrheniusTe, ArrheniusT) -SELECT_RATE_LAWS(ASSOCIATIVE_DETACHMENT, ArrheniusT, ArrheniusTe) -SELECT_RATE_LAWS(DISSOCIATIVE_ATTACHMENT, ArrheniusTe, ArrheniusT) -SELECT_RATE_LAWS(DISSOCIATION_E, ArrheniusTe, ArrheniusTe) -SELECT_RATE_LAWS(RECOMBINATION_E, ArrheniusTe, ArrheniusTe) -SELECT_RATE_LAWS(DISSOCIATION_M, ArrheniusPark, ArrheniusT) -SELECT_RATE_LAWS(RECOMBINATION_M, ArrheniusT, ArrheniusPark) -SELECT_RATE_LAWS(IONIZATION_E, ArrheniusTe, ArrheniusT) -SELECT_RATE_LAWS(ION_RECOMBINATION_E, ArrheniusT, ArrheniusTe) -SELECT_RATE_LAWS(IONIZATION_M, ArrheniusT, ArrheniusT) -SELECT_RATE_LAWS(ION_RECOMBINATION_M, ArrheniusT, ArrheniusT) -SELECT_RATE_LAWS(ELECTRONIC_ATTACHMENT_M, ArrheniusTe, ArrheniusT) -SELECT_RATE_LAWS(ELECTRONIC_DETACHMENT_M, ArrheniusT, ArrheniusTe) -SELECT_RATE_LAWS(ELECTRONIC_ATTACHMENT_E, ArrheniusTe, ArrheniusTe) -SELECT_RATE_LAWS(ELECTRONIC_DETACHMENT_E, ArrheniusTe, ArrheniusTe) -SELECT_RATE_LAWS(EXCHANGE, ArrheniusT, ArrheniusT) -SELECT_RATE_LAWS(EXCITATION_M, ArrheniusT, ArrheniusT) -SELECT_RATE_LAWS(EXCITATION_E, ArrheniusTe, ArrheniusTe) - -#undef SELECT_RATE_LAWS - -//============================================================================== RateManager::RateManager(size_t ns, const std::vector& reactions) : m_ns(ns), m_nr(reactions.size()), mp_lnkf(NULL), mp_lnkb(NULL), @@ -197,11 +159,12 @@ struct is_same { template void RateManager::addRate(const size_t rxn, const Reaction& reaction) -{ +{ + m_rate_groups.addRateCoefficient(rxn, reaction.rateLaw()); if (reaction.isReversible()) { - + // Make use of forward computation when possible if (is_same::value) m_to_copy.push_back(rxn); @@ -212,7 +175,7 @@ void RateManager::addRate(const size_t rxn, const Reaction& reaction) rxn+m_nr, reaction.rateLaw()); m_rate_groups.addReaction(rxn, reaction); - + } else { m_irr.push_back(rxn); } diff --git a/src/kinetics/RateManager.h b/src/kinetics/RateManager.h index 5d93148c..235f33f0 100644 --- a/src/kinetics/RateManager.h +++ b/src/kinetics/RateManager.h @@ -31,9 +31,66 @@ #include "RateLawGroup.h" #include "Reaction.h" + namespace Mutation { - namespace Kinetics { + namespace Kinetics { + + class TSelector + { + public: + inline double getT(const Thermodynamics::StateModel* const state) const; + }; + + class TeSelector + { + public: + inline double getT(const Thermodynamics::StateModel* const state) const; + }; + + class ParkSelector + { + public: + inline double getT(const Thermodynamics::StateModel* const state) const; + }; + +/** + * Used to define which rate law groups the forward and reverse rate laws are + * evaluated in for a given ReactionType value. The default is an Arrhenius + * rate law with Tf = Tb = T. + */ +/// Arrhenius group evaluated at T +typedef RateLawGroup1T ArrheniusT; + +/// Arrhenius group evaluated at Te +typedef RateLawGroup1T ArrheniusTe; +/// Arrhenius group evaluated at sqrt(T*Tv) +typedef RateLawGroup1T ArrheniusPark; + +// Map RateLawTemperature to RateLawGroup1T +template struct RateLawGroupForTemp; +template<> struct RateLawGroupForTemp { using type = ArrheniusT; }; +template<> struct RateLawGroupForTemp { using type = ArrheniusTe; }; +template<> struct RateLawGroupForTemp { using type = ArrheniusPark; }; + +// template +// struct RateSelector { +// typedef ArrheniusT ForwardGroup; +// typedef ArrheniusT ReverseGroup; +// }; + +// Rate selector +template struct RateSelector { +private: + static constexpr auto temps = getRateLawT(T); + +public: + using ForwardGroup = typename RateLawGroupForTemp::type; + using ReverseGroup = typename RateLawGroupForTemp::type; +}; + + + class Reaction; /** @@ -78,7 +135,7 @@ class RateManager const std::vector& irrReactions() const { return m_irr; } - + private: /** @@ -128,6 +185,7 @@ class RateManager /// Stores the indices of non-reversible reactions std::vector m_irr; + }; diff --git a/src/kinetics/Reaction.cpp b/src/kinetics/Reaction.cpp index 5cdc20c1..5a8df24d 100644 --- a/src/kinetics/Reaction.cpp +++ b/src/kinetics/Reaction.cpp @@ -325,62 +325,70 @@ void Reaction::determineType(const class Thermodynamics& thermo) if (ecount > 0) { if (chgcount > 0) { - if (m_inert_e) - m_type = ION_RECOMBINATION_E; - else if (m_inert) - m_type = ION_RECOMBINATION_M; - else - m_type = DISSOCIATIVE_RECOMBINATION; + if (m_inert_e){ + m_type = ION_RECOMBINATION_E; + }else if (m_inert){ + m_type = ION_RECOMBINATION_M; + }else { + m_type = DISSOCIATIVE_RECOMBINATION; + } } else { - if (m_inert_e) - m_type = ELECTRONIC_ATTACHMENT_E; - else if (m_inert) - m_type = ELECTRONIC_ATTACHMENT_M; - else - m_type = DISSOCIATIVE_ATTACHMENT; + if (m_inert_e){ + m_type = ELECTRONIC_ATTACHMENT_E; + }else if (m_inert){ + m_type = ELECTRONIC_ATTACHMENT_M; + }else { + m_type = DISSOCIATIVE_ATTACHMENT; + } } } else if (ecount < 0) { if (chgcount > 0) { - if (m_inert_e) - m_type = ELECTRONIC_DETACHMENT_E; - else if (m_inert) - m_type = ELECTRONIC_DETACHMENT_M; - else - m_type = ASSOCIATIVE_DETACHMENT; + if (m_inert_e){ + m_type = ELECTRONIC_DETACHMENT_E; + }else if (m_inert){ + m_type = ELECTRONIC_DETACHMENT_M; + }else{ + m_type = ASSOCIATIVE_DETACHMENT; + } + } else { - if (m_inert_e) - m_type = IONIZATION_E; - else if (m_inert) - m_type = IONIZATION_M; - else - m_type = ASSOCIATIVE_IONIZATION; + if (m_inert_e){ + m_type = IONIZATION_E; + }else if (m_inert){ + m_type = IONIZATION_M; + }else{ + m_type = ASSOCIATIVE_IONIZATION; + } } } else { if (nReactants() > nProducts()) { - if (m_inert_e) - m_type = RECOMBINATION_E; - else - m_type = RECOMBINATION_M; // including double recombination + if (m_inert_e){ + m_type = RECOMBINATION_E; + }else{ + m_type = RECOMBINATION_M; // including double recombination + } } else if (nReactants() < nProducts()) { - if (m_inert_e) - m_type = DISSOCIATION_E; - else - m_type = DISSOCIATION_M; // including double dissociation + if (m_inert_e){ + m_type = DISSOCIATION_E; + }else{ + m_type = DISSOCIATION_M; // including double dissociation + } } else { - if (m_inert_e) - m_type = EXCITATION_E; - else if (m_inert) - m_type = EXCITATION_M; - else - m_type = EXCHANGE; // exchange of charge, atom, or both simultaneously + if (m_inert_e){ + m_type = EXCITATION_E; + }else if (m_inert){ + m_type = EXCITATION_M; + }else{ + m_type = EXCHANGE; // exchange of charge, atom, or both + } } } } diff --git a/src/kinetics/Reaction.h b/src/kinetics/Reaction.h index 6f671a48..7c6359b2 100644 --- a/src/kinetics/Reaction.h +++ b/src/kinetics/Reaction.h @@ -203,7 +203,15 @@ class Reaction const std::vector >& efficiencies() const { return m_thirdbodies; } - + + const std::string fwdRateTemperature() { + return rateLawTemperatureString(getRateLawT(m_type).first); + } + + const std::string revRateTemperature() { + return rateLawTemperatureString(getRateLawT(m_type).second); + } + friend void swap(Reaction& left, Reaction& right); private: @@ -254,6 +262,8 @@ class Reaction */ void determineType(const Mutation::Thermodynamics::Thermodynamics& thermo); + + private: std::string m_formula; diff --git a/src/kinetics/ReactionType.cpp b/src/kinetics/ReactionType.cpp index ef28eb4c..af44d266 100644 --- a/src/kinetics/ReactionType.cpp +++ b/src/kinetics/ReactionType.cpp @@ -78,6 +78,14 @@ const char* const reactionTypeString(const ReactionType type) } } +const std::string rateLawTemperatureString(const RateLawTemperature eval_temp){ + switch (eval_temp){ + case TTRANSLATIONAL: return "translational"; + case TELECTRON: return "electron"; + case TGEOMETRIC_TTv: return "geometric_ttv"; + } +} + //============================================================================== } // namespace Kinetics diff --git a/src/kinetics/ReactionType.h b/src/kinetics/ReactionType.h index 6261f640..2dc60c1d 100644 --- a/src/kinetics/ReactionType.h +++ b/src/kinetics/ReactionType.h @@ -28,6 +28,7 @@ #ifndef KINETICS_REACTION_TYPE_H #define KINETICS_REACTION_TYPE_H +#include namespace Mutation { namespace Kinetics { @@ -59,13 +60,50 @@ enum ReactionType EXCITATION_M, MAX_REACTION_TYPES /// Number of ReactionType values (leave at end of list) }; - + /** * Simpy returns a string which represents the reaction type given by the * ReactionType argument. */ const char* const reactionTypeString(const ReactionType type); +/** + * Enumerate possible temperatures at which + * rate laws are evaluated. + */ +enum RateLawTemperature{ + TTRANSLATIONAL, + TELECTRON, + TGEOMETRIC_TTv +}; + +const std::string rateLawTemperatureString(const RateLawTemperature eval_temp); + +constexpr std::pair getRateLawT(int rxn_type) { + switch (rxn_type) { + case ASSOCIATIVE_IONIZATION: return {TTRANSLATIONAL, TELECTRON}; + case DISSOCIATIVE_RECOMBINATION: return {TELECTRON, TTRANSLATIONAL}; + case ASSOCIATIVE_DETACHMENT: return {TTRANSLATIONAL, TELECTRON}; + case DISSOCIATIVE_ATTACHMENT: return {TELECTRON, TTRANSLATIONAL}; + case DISSOCIATION_E: return {TELECTRON, TELECTRON}; + case RECOMBINATION_E: return {TELECTRON, TELECTRON}; + case DISSOCIATION_M: return {TGEOMETRIC_TTv, TTRANSLATIONAL}; + case RECOMBINATION_M: return {TTRANSLATIONAL, TGEOMETRIC_TTv}; + case IONIZATION_E: return {TELECTRON, TTRANSLATIONAL}; + case ION_RECOMBINATION_E: return {TTRANSLATIONAL, TELECTRON}; + case IONIZATION_M: return {TTRANSLATIONAL, TTRANSLATIONAL}; + case ION_RECOMBINATION_M: return {TTRANSLATIONAL, TTRANSLATIONAL}; + case ELECTRONIC_ATTACHMENT_M: return {TELECTRON, TTRANSLATIONAL}; + case ELECTRONIC_DETACHMENT_M: return {TTRANSLATIONAL, TELECTRON}; + case ELECTRONIC_ATTACHMENT_E: return {TELECTRON, TELECTRON}; + case ELECTRONIC_DETACHMENT_E: return {TELECTRON, TELECTRON}; + case EXCHANGE: return {TTRANSLATIONAL, TTRANSLATIONAL}; + case EXCITATION_M: return {TTRANSLATIONAL, TTRANSLATIONAL}; + case EXCITATION_E: return {TELECTRON, TELECTRON}; + default: return {TTRANSLATIONAL, TTRANSLATIONAL}; + } +}; + } // namespace Kinetics } // namespace Mutation