From 2e03c918458972daf367f6804f3e5259f77f2f12 Mon Sep 17 00:00:00 2001 From: Cheng Li <69489965+chengcli@users.noreply.github.com> Date: Mon, 8 Jul 2024 16:45:08 -0400 Subject: [PATCH] Add condensation to Kinetics (#16) - Add condensation class to Kinetics - Add vapor pressures and Nucleation --- include/cantera/kinetics/Condensation.h | 111 ++ include/cantera/kinetics/Freezing.h | 55 + include/cantera/kinetics/Kinetics.h | 5 + include/cantera/kinetics/Nucleation.h | 69 ++ include/cantera/kinetics/svp_funcs.h | 8 + include/cantera/kinetics/vapors/ammonia.dat | 59 + .../kinetics/vapors/ammonia_vapors.hpp | 78 ++ .../vapors/ammonium_hydrosulfide_vapors.hpp | 22 + .../kinetics/vapors/carbon_dioxide_vapors.hpp | 26 + .../cantera/kinetics/vapors/fit_ammonia.py | 18 + .../vapors/hydrogen_sulfide_vapors.hpp | 75 ++ .../kinetics/vapors/methane_vapors.hpp | 29 + .../kinetics/vapors/potassium_vapors.hpp | 9 + .../cantera/kinetics/vapors/sodium_vapors.hpp | 10 + .../kinetics/vapors/sulfur_dioxide_vapors.hpp | 16 + .../cantera/kinetics/vapors/water_vapors.hpp | 110 ++ include/cantera/thermo/IdealGasSurfPhase.h | 47 + include/cantera/thermo/IdealMoistPhase.h | 100 ++ include/cantera/thermo/LatticeSolidPhase.h | 4 +- include/cantera/thermo/Phase.h | 12 +- include/cantera/thermo/SingleSpeciesTP.h | 2 +- include/cantera/thermo/ThermoPhase.h | 11 + install_cantera.sh | 23 + src/kinetics/Condensation.cpp | 455 +++++++ src/kinetics/Freezing.cpp | 62 + src/kinetics/KineticsFactory.cpp | 2 + src/kinetics/Nucleation.cpp | 157 +++ src/kinetics/ReactionRateFactory.cpp | 14 + src/thermo/IdealGasSurfPhase.cpp | 21 + src/thermo/IdealMoistPhase.cpp | 87 ++ src/thermo/Phase.cpp | 32 +- src/thermo/ThermoFactory.cpp | 4 + z.ab_patch | 1043 +++++++++++++++++ 33 files changed, 2767 insertions(+), 9 deletions(-) create mode 100644 include/cantera/kinetics/Condensation.h create mode 100644 include/cantera/kinetics/Freezing.h create mode 100644 include/cantera/kinetics/Nucleation.h create mode 100644 include/cantera/kinetics/svp_funcs.h create mode 100644 include/cantera/kinetics/vapors/ammonia.dat create mode 100644 include/cantera/kinetics/vapors/ammonia_vapors.hpp create mode 100644 include/cantera/kinetics/vapors/ammonium_hydrosulfide_vapors.hpp create mode 100644 include/cantera/kinetics/vapors/carbon_dioxide_vapors.hpp create mode 100755 include/cantera/kinetics/vapors/fit_ammonia.py create mode 100644 include/cantera/kinetics/vapors/hydrogen_sulfide_vapors.hpp create mode 100644 include/cantera/kinetics/vapors/methane_vapors.hpp create mode 100644 include/cantera/kinetics/vapors/potassium_vapors.hpp create mode 100644 include/cantera/kinetics/vapors/sodium_vapors.hpp create mode 100644 include/cantera/kinetics/vapors/sulfur_dioxide_vapors.hpp create mode 100644 include/cantera/kinetics/vapors/water_vapors.hpp create mode 100644 include/cantera/thermo/IdealGasSurfPhase.h create mode 100644 include/cantera/thermo/IdealMoistPhase.h create mode 100755 install_cantera.sh create mode 100644 src/kinetics/Condensation.cpp create mode 100644 src/kinetics/Freezing.cpp create mode 100644 src/kinetics/Nucleation.cpp create mode 100644 src/thermo/IdealGasSurfPhase.cpp create mode 100644 src/thermo/IdealMoistPhase.cpp create mode 100644 z.ab_patch diff --git a/include/cantera/kinetics/Condensation.h b/include/cantera/kinetics/Condensation.h new file mode 100644 index 0000000000..2a348f71f4 --- /dev/null +++ b/include/cantera/kinetics/Condensation.h @@ -0,0 +1,111 @@ +#ifndef CT_CONDENSATION_H +#define CT_CONDENSATION_H + +#include "Kinetics.h" + +namespace Cantera +{ + +class Condensation : public Kinetics { + public: + Condensation() = default; + + ~Condensation() override {} + + void setQuantityMoleFraction() { + if (!m_use_mole_fraction) { + m_ROP_ok = false; + } + m_use_mole_fraction = true; + } + + void setQuantityConcentration() { + if (m_use_mole_fraction) { + m_ROP_ok = false; + } + m_use_mole_fraction = false; + } + + void resizeReactions() override; + + string kineticsType() const override { + return "nucleation"; + } + + void getActivityConcentrations(double* const pdata) override; + + bool isReversible(size_t i) override { + return true; + } + + void getFwdRateConstants(double* kfwd) override; + + void getFwdRateConstants_ddT(double* dkfwd) override; + + void resizeSpecies() override; + + bool addReaction(shared_ptr r, bool resize=true) override; + + void updateROP() override; + + Eigen::SparseMatrix netRatesOfProgress_ddCi() override; + Eigen::SparseMatrix netRatesOfProgress_ddX() override; + + protected: + void _update_rates_T(double *pdata, double *pdata_ddT); + void _update_rates_C(double *pdata); + void _update_rates_X(double *pdata); + + //! Process derivatives + //! @param stoich stoichiometry manager + //! @param in rate expression used for the derivative calculation + //! @param ddX true: w.r.t mole fractions false: w.r.t species concentrations + //! @return a sparse matrix of derivative contributions for each reaction of + //! dimensions nTotalReactions by nTotalSpecies + Eigen::SparseMatrix calculateCompositionDerivatives( + StoichManagerN& stoich, const vector& in, bool ddX=true); + + //! This variable has two interpretations. + //! If m_use_mole_fraction is true, then it is the vector of mole fractions. + //! If m_use_mole_fraction is false, then it is the vector of concentrations. + Eigen::VectorXd m_conc; + Eigen::VectorXd m_intEng; + Eigen::VectorXd m_cv; + + //! Number of dimensions of reacting phase (2 for InterfaceKinetics, 1 for + //! EdgeKinetics) + size_t m_nDim = 2; + + //! Vector of rate handlers for interface reactions + vector> m_interfaceRates; + map m_interfaceTypes; //!< Rate handler mapping + + //! reaction index for x <=> y + vector m_jxy; + + //! reaction index for x1 + x2 <=> y + vector m_jxxy; + + //! reaction index for freezing reaction x(l) <=> x(s) + vector m_jyy; + + //! rate jacobian matrix + Eigen::SparseMatrix m_jac; + + //! rate jacobian with respect to temperature + vector m_rfn_ddT; + + bool m_use_mole_fraction = false; + bool m_ROP_ok = false; + + //! Current temperature of the data + double m_temp = 0.0; + + //! Buffers for partial rop results with length nReactions() + vector m_rbuf0; + vector m_rbuf1; +}; + +} + +#endif diff --git a/include/cantera/kinetics/Freezing.h b/include/cantera/kinetics/Freezing.h new file mode 100644 index 0000000000..ab6e25df35 --- /dev/null +++ b/include/cantera/kinetics/Freezing.h @@ -0,0 +1,55 @@ +#ifndef CT_FREEZING_H +#define CT_FREEZING_H + +#include + +#include "cantera/base/ct_defs.h" +#include "cantera/base/Units.h" +#include "cantera/kinetics/Arrhenius.h" +#include "ReactionRate.h" +#include "MultiRate.h" + +namespace Cantera +{ + +class AnyValue; +class AnyMap; + +class FreezingRate : public ReactionRate { + public: + FreezingRate() = default; + FreezingRate(const AnyMap& node, const UnitStack& rate_units); + + unique_ptr newMultiRate() const override { + return make_unique>(); + } + + //! Set the rate parameters for this reaction. + void setRateParameters(const AnyValue& equation, + const AnyValue& rate, + const AnyMap& node); + + //! return the rate coefficient type + const string type() const override { return "freezing"; } + + void getParameters(AnyMap& rateNode, const Units& rate_units=Units(0.)) const; + using ReactionRate::getParameters; + + void validate(const string& equation, const Kinetics& kin) override; + + double evalFromStruct(const ArrheniusData& shared_data) const; + + double ddTScaledFromStruct(const ArrheniusData& shared_data) const { + return 0.; + } + + protected: + string m_formula_str = "formula"; + std::function m_meltfunc; + + double m_t3 = 0.0; +}; + +} + +#endif diff --git a/include/cantera/kinetics/Kinetics.h b/include/cantera/kinetics/Kinetics.h index 4a4f91c14b..f8f36882cd 100644 --- a/include/cantera/kinetics/Kinetics.h +++ b/include/cantera/kinetics/Kinetics.h @@ -303,6 +303,11 @@ class Kinetics */ size_t kineticsSpeciesIndex(const string& nm) const; + //! shorthand for kineticsSpeciesIndex + size_t speciesIndex(const string& nm) const { + return kineticsSpeciesIndex(nm); + } + /** * This function looks up the name of a species and returns a * reference to the ThermoPhase object of the phase where the species diff --git a/include/cantera/kinetics/Nucleation.h b/include/cantera/kinetics/Nucleation.h new file mode 100644 index 0000000000..0a93428f0f --- /dev/null +++ b/include/cantera/kinetics/Nucleation.h @@ -0,0 +1,69 @@ +#ifndef CT_NUCLEATION_H +#define CT_NUCLEATION_H + +#include + +#include "cantera/base/ct_defs.h" +#include "cantera/base/Units.h" +#include "cantera/kinetics/Arrhenius.h" +#include "ReactionRate.h" +#include "MultiRate.h" + +namespace Cantera +{ + +class AnyValue; +class AnyMap; + +class NucleationRate : public ReactionRate { + public: + NucleationRate() = default; + NucleationRate(const AnyMap& node, const UnitStack& rate_units); + + unique_ptr newMultiRate() const override { + return make_unique>(); + } + + //! Set the rate parameters for this reaction. + void setRateParameters(const AnyValue& equation, + const AnyValue& rate, + const AnyMap& node); + + //! return the rate coefficient type + const string type() const override { return "nucleation"; } + + void getParameters(AnyMap& rateNode, const Units& rate_units=Units(0.)) const; + using ReactionRate::getParameters; + + void validate(const string& equation, const Kinetics& kin) override; + + double evalFromStruct(const ArrheniusData& shared_data) const; + + double ddTScaledFromStruct(const ArrheniusData& shared_data) const; + + protected: + string m_formula_str = "formula"; + + //! returns s = svp(T)/RT or svp(T)/(RT)^2 + std::function m_svp; + + //! returns d(log(s))/dT + std::function m_logsvp_ddT; + + size_t m_order = 1; + + double m_t3 = 0.0; + double m_p3 = 0.0; + double m_beta = 0.0; + double m_delta = 0.0; + + double m_min_temp = 0.0; + double m_max_temp = 1.0e30; +}; + +std::function find_svp_function(const Composition& reactants, + const string& svp_name); + +} + +#endif diff --git a/include/cantera/kinetics/svp_funcs.h b/include/cantera/kinetics/svp_funcs.h new file mode 100644 index 0000000000..fe5e854f5c --- /dev/null +++ b/include/cantera/kinetics/svp_funcs.h @@ -0,0 +1,8 @@ +#ifndef SVP_FUNCS_INC_H +#define SVP_FUNCS_INC_H + +#include "vapors/water_vapors.hpp" +#include "vapors/ammonia_vapors.hpp" +#include "vapors/ammonium_hydrosulfide_vapors.hpp" + +#endif // SVP_FUNCS_INC_H diff --git a/include/cantera/kinetics/vapors/ammonia.dat b/include/cantera/kinetics/vapors/ammonia.dat new file mode 100644 index 0000000000..951eb7b381 --- /dev/null +++ b/include/cantera/kinetics/vapors/ammonia.dat @@ -0,0 +1,59 @@ +# Ammonia Thermodynamic Properties from Engineering Toolbox +# https://www.engineeringtoolbox.com/ammonia-d_971.html +# +# T - temperature (K) +# P - pressure (bar) +# alpha - specific volume (m^3/kg) +# Hf(l) - liquid specific enthalpy (kJ/kg) +# Hf(v) - vapor specific enthalpy (kJ/kg) +# Sf(l) - liquid specific entropy (kJ/kg) +# Sf(v) - vapor specific entropy (kJ/kg) +# +# Saturated Properties Superheated Properties +# ---------------------------- ------------------------- +# T P alpha Hf(l) Hf(v) Sf(l) Sf(v) H(50K) S(50K) H(100K) S(100K) +-50 0.4089 2.625 -44.4 1373.3 -0.194 6.159 1479.8 6.592 1585.9 6.948 +-45 0.5454 2.005 -22.3 1381.6 -0.096 6.057 1489.3 6.486 1596.1 6.839 +-40 0.7177 1.552 0 1390.0 0 5.962 1498.6 6.387 1606.3 6.736 +-35 0.9322 1.216 22.3 1397.9 0.095 5.872 1507.9 6.293 1616.3 6.639 +-30 1.196 0.9633 44.7 1405.6 0.188 5.785 1517.0 6.203 1626.3 6.547 +-28 1.317 0.8809 53.6 1408.5 0.224 5.751 1520.7 6.169 1630.3 6.512 +-26 1.447 0.8058 62.6 1411.4 0.261 5.718 1524.3 6.135 1634.2 6.477 +-24 1.588 0.7389 71.7 1414.3 0.297 5.686 1527.9 6.103 1638.2 6.444 +-22 1.740 0.6783 80.8 1417.3 0.333 5.655 1531.4 6.071 1642.2 6.411 +-20 1.902 0.6237 89.8 1420.0 0.368 5.623 1534.8 6.039 1646.0 6.379 +-18 2.077 0.5743 98.8 1422.7 0.404 5.593 1538.2 6.008 1650.0 6.347 +-16 2.265 0.5296 107.9 1425.3 0.440 5.563 1541.7 5.978 1653.8 6.316 +-14 2.465 0.4890 117.0 1427.9 0.475 5.533 1545.1 5.948 1657.7 6.286 +-12 2.680 0.4521 126.2 1430.5 0.510 5.504 1548.5 5.919 1661.5 6.256 +-10 2.908 0.4185 135.4 1433.0 0.544 5.475 1551.7 5.891 1665.3 6.227 +-8 3.153 0.3879 144.5 1435.3 0.579 5.447 1554.9 5.863 1669.0 6.199 +-6 3.413 0.3599 153.6 1437.6 0.613 5.419 1558.2 5.836 1672.8 6.171 +-4 3.691 0.3344 162.8 1439.9 0.647 5.392 1561.4 5.808 1676.4 6.143 +-2 3.983 0.3110 172.0 1442.2 0.681 5.365 1564.6 5.782 1680.1 6.116 +0 4.295 0.2895 181.2 1444.4 0.715 5.340 1567.8 5.756 1683.9 6.090 +2 4.625 0.2699 190.4 1446.5 0.749 5.314 1570.9 5.731 1687.5 6.065 +4 4.975 0.2517 199.7 1448.5 0.782 5.288 1574.0 5.706 1691.2 6.040 +6 5.346 0.2351 209.1 1450.6 0.816 5.263 1577.0 5.682 1694.9 6.015 +8 5.736 0.2198 218.5 1452.5 0.849 5.238 1580.1 5.658 1698.4 5.991 +10 6.149 0.2056 227.8 1454.3 0.881 5.213 1583.1 5.634 1702.2 5.967 +12 6.585 0.1926 237.2 1456.1 0.914 5.189 1586.0 5.611 1705.7 5.943 +14 7.045 0.1805 246.6 1457.8 0.947 5.165 1588.9 5.588 1709.1 5.920 +16 7.529 0.1693 256.0 1459.5 0.979 5.141 1591.7 5.565 1712.5 5.898 +18 8.035 0.1590 265.5 1461.1 1.012 5.118 1594.4 5.543 1715.9 5.876 +20 8.570 0.1494 275.1 1462.6 1.044 5.095 1597.2 5.521 1719.3 5.854 +22 9.134 0.1405 284.6 1463.9 1.076 5.072 1600.0 5.499 1722.8 5.832 +24 9.722 0.1322 294.1 1465.2 1.108 5.049 1602.7 5.478 1726.3 5.811 +26 10.34 0.1245 303.7 1466.5 1.140 5.027 1605.3 5.458 1729.6 5.790 +28 10.99 0.1173 313.4 1467.8 1.172 5.005 1608.0 5.437 1732.7 5.770 +30 11.67 0.1106 323.1 1468.9 1.204 4.984 1610.5 5.417 1735.9 5.750 +32 12.37 0.1044 332.8 1469.9 1.235 4.962 1613.0 5.397 1739.3 5.731 +34 13.11 0.0986 342.5 1470.8 1.267 4.940 1615.4 5.378 1742.6 5.711 +36 13.89 0.0931 352.3 1471.8 1.298 4.919 1617.8 5.358 1745.7 5.692 +38 14.70 0.0880 362.1 1472.6 1.329 4.898 1620.1 5.340 1748.7 5.674 +40 15.54 0.0833 371.9 1473.3 1.360 4.877 1622.4 5.321 1751.9 5.655 +42 16.42 0.0788 381.8 1473.8 1.391 4.856 1624.6 5.302 1755.0 5.637 +44 17.34 0.0746 391.8 1474.2 1.422 4.835 1626.8 5.284 1758.0 5.619 +46 18.30 0.0706 401.8 1474.5 1.453 4.814 1629.0 5.266 1761.0 5.602 +48 19.29 0.0670 411.9 1474.7 1.484 4.793 1631.1 5.248 1764.0 5.584 +50 20.33 0.0635 421.9 1474.7 1.515 4.773 1633.1 5.230 1766.8 5.567 diff --git a/include/cantera/kinetics/vapors/ammonia_vapors.hpp b/include/cantera/kinetics/vapors/ammonia_vapors.hpp new file mode 100644 index 0000000000..eb7f4a1910 --- /dev/null +++ b/include/cantera/kinetics/vapors/ammonia_vapors.hpp @@ -0,0 +1,78 @@ +#ifndef SRC_SNAP_THERMODYNAMICS_VAPORS_AMMONIA_VAPORS_HPP_ +#define SRC_SNAP_THERMODYNAMICS_VAPORS_AMMONIA_VAPORS_HPP_ + +// C/C++ +#include + +inline double svpnh3(double t, double p, double beta, double gamma) { + return p * exp((1. - 1. / t) * beta - gamma * log(t)); +} + +inline double svp_nh3_UMich(double T) { + double x = -1790.00 / T - 1.81630 * log10(T) + 14.97593; + return pow(10.0, x) * 1.E5 / 760.; +} + +// (164, 371.5) +inline double svp_nh3_Antoine(double T) { + double result; + if (T < 239.6) + result = pow(10., 3.18757 - (506.713 / (T - 80.78))); + else + result = pow(10., 4.86886 - (1113.928 / (T - 10.409))); + return 1.E5 * result; +} + +// (130, 200) +inline double sat_vapor_p_NH3_Hubner(double T) { + double A = 24.3037, B = -1766.28, C = -5.64472, D = 0.00740241; + + double x = A + B / T + C * log10(T) + D * T; + return pow(10.0, x); +} + +inline double sat_vapor_p_NH3_BriggsS(double T) { + double a[6], x; + if (T < 195) { + a[1] = -4122.; + a[2] = 41.67871; + a[3] = -1.81630; + a[4] = 0.; + a[5] = 0.; + } else { + a[1] = -4409.3512; + a[2] = 76.864252; + a[3] = -8.4598340; + a[4] = 5.51029e-03; + a[5] = 6.80463e-06; + } + x = a[1] / T + a[2] + a[3] * log(T) + a[4] * T + a[5] * pow(T, 2); + return exp(x) / 10.; +} + +// (15, 195.4) +inline double sat_vapor_p_NH3_Fray(double T) { + double a[7], x = 0; + a[0] = 1.596e+01; + a[1] = -3.537e+03; + a[2] = -3.310e+04; + a[3] = 1.742e+06; + a[4] = -2.995e+07; + a[5] = 0.; + a[6] = 0.; + + for (int i = 1; i < 7; i++) + x = x + a[i] / pow(T, i); // best fit in [15K; 195.41K] + + return 1.E5 * exp(x + a[0]); +} + +inline double sat_vapor_p_NH3_Ideal(double T) { + double betal = 20.08, gammal = 5.62, betas = 20.64, gammas = 1.43, tr = 195.4, + pr = 6060.; + + return T > tr ? svpnh3(T / tr, pr, betal, gammal) + : svpnh3(T / tr, pr, betas, gammas); +} + +#endif // SRC_SNAP_THERMODYNAMICS_VAPORS_AMMONIA_VAPORS_HPP_ diff --git a/include/cantera/kinetics/vapors/ammonium_hydrosulfide_vapors.hpp b/include/cantera/kinetics/vapors/ammonium_hydrosulfide_vapors.hpp new file mode 100644 index 0000000000..ca89556652 --- /dev/null +++ b/include/cantera/kinetics/vapors/ammonium_hydrosulfide_vapors.hpp @@ -0,0 +1,22 @@ +#ifndef VAPORS_AMMONIUM_HYDROSULFIDE_VAPORS_HPP_ +#define VAPORS_AMMONIUM_HYDROSULFIDE_VAPORS_HPP_ + +// C/C++ +#include + +const double Pcgs_of_atm = 1013250.0; // atmospheres to dynes/cm**2 + +inline double svp_nh3_h2s_Umich(double T) { + double const GOLB2 = (14.83 - (4715.0 / T)); + return (pow(10.0, GOLB2)) * Pcgs_of_atm * Pcgs_of_atm; +} + +inline double svp_nh3_h2s_Lewis(double T) { + return pow(10., 14.82 - 4705. / T) * 101325. * 101325.; +} + +inline double logsvp_ddT_nh3_h2s_Lewis(double T) { + return 4705. * log(10) / (T * T); +} + +#endif // VAPORS_AMMONIUM_HYDROSULFIDE_VAPORS_HPP_ diff --git a/include/cantera/kinetics/vapors/carbon_dioxide_vapors.hpp b/include/cantera/kinetics/vapors/carbon_dioxide_vapors.hpp new file mode 100644 index 0000000000..c8a8257d70 --- /dev/null +++ b/include/cantera/kinetics/vapors/carbon_dioxide_vapors.hpp @@ -0,0 +1,26 @@ +#ifndef SRC_SNAP_THERMODYNAMICS_VAPORS_CARBON_DIOXIDE_VAPORS_HPP_ +#define SRC_SNAP_THERMODYNAMICS_VAPORS_CARBON_DIOXIDE_VAPORS_HPP_ + +// C/C++ +#include + +// best fit in [154.26, 195.89] K +inline double sat_vapor_p_CO2_Antoine(double T) { + double A = 6.81228; + double B = 1301.679; + double C = -34.94; + + double result = pow(10, A - B / (T + C)); + return 1.E5 * result; + + // best fit in [154.2, 204] K + + // double A = 9.8106; + // double B = 1347.79; + // double C = 273; + + // double result = pow(10, A - B / ((T-273.2) + C)); + // return 133.3 * result; +} + +#endif // SRC_SNAP_THERMODYNAMICS_VAPORS_CARBON_DIOXIDE_VAPORS_HPP_ diff --git a/include/cantera/kinetics/vapors/fit_ammonia.py b/include/cantera/kinetics/vapors/fit_ammonia.py new file mode 100755 index 0000000000..1f33f24a01 --- /dev/null +++ b/include/cantera/kinetics/vapors/fit_ammonia.py @@ -0,0 +1,18 @@ +#! /usr/bin/env python2.7 +from scipy.stats import linregress +from pylab import * + +data = genfromtxt('ammonia.dat') +temp = data[:,0] + 273.15 +hfv = data[:,4] +a, b, r, p, e = linregress(temp, hfv) +print 'cpv = %.2f kJ/kg' % a + +hfl = data[:,3] +a, b, r, p, e = linregress(temp, hfl) +print 'cpl = %.2f kJ/kg' % a + +hfl = data[:,3] +a, b, r, p, e = linregress(temp, hfv - hfl) +print 'delta_cp = %.2f kJ/kg' % -a +print 'mu = %.2f kJ/kg' % b diff --git a/include/cantera/kinetics/vapors/hydrogen_sulfide_vapors.hpp b/include/cantera/kinetics/vapors/hydrogen_sulfide_vapors.hpp new file mode 100644 index 0000000000..21e894bf70 --- /dev/null +++ b/include/cantera/kinetics/vapors/hydrogen_sulfide_vapors.hpp @@ -0,0 +1,75 @@ +#ifndef SRC_SNAP_THERMODYNAMICS_VAPORS_HYDROGEN_SULFIDE_VAPORS_HPP_ +#define SRC_SNAP_THERMODYNAMICS_VAPORS_HYDROGEN_SULFIDE_VAPORS_HPP_ + +// C/C++ +#include + +inline double Pcgs_from_torr(double P) { return (1013250.0 / 760.0) * P; } + +inline double sat_vapor_p_H2S_solid_UMich(double T) { + // Vapor of H2S over H2S liquid, T > 188 K, over H2S ice for T < 188 K + double x; + double log10T = log10(T); + if (T > 212.75) + x = 7.7547 - 976.0 / T - 0.12058 * log10T; + else if (T > 188.0) + x = 15.4859 - 1264.3574 / T - 2.86206 * log10T; + else + x = 57.19 - 2461.84 / T - 18.443 * log10T; + + return Pcgs_from_torr(pow(10.0, x)) * 0.1; +} + +// best fit in [138.8K; 349.5K] (webbook.nist.gov) +inline double sat_vapor_p_H2S_Antoine(double T) { + double result; + if (T < 212.8) + result = pow(10., 4.43681 - (829.439 / (T - 25.412))); + else + result = pow(10., 4.52887 - (958.587 / (T - 0.539))); + + return 1.E5 * result; +} + +// best fit in [120K; 190K] +inline double sat_vapor_p_H2S_solid_Hubner(double T) { + double A, B, C, D; + double log10T = log10(T); + A = 6.96156; + B = -903.815; + C = 0.258812; + D = 0.00873804; + + double x = A + B / T + C * log10(T) + D * T; + return 1.E5 * pow(10.0, x); +} + +// best fit in [110K; 187.57K] +inline double sat_vapor_p_H2S_solid_Fray(double T) { + double a[7], x; + x = 0.; + if (T > 127) { + a[0] = 8.933; + a[1] = -7.260e+02; + a[2] = -3.504e+05; + a[3] = 2.724e+07; + a[4] = -8.582e+08; + a[5] = 0.; + a[6] = 0.; + } else { + a[0] = 1.298e+01; + a[1] = -2.707e+03; + a[2] = 0.; + a[3] = 0.; + a[4] = 0.; + a[5] = 0.; + a[6] = 0.; + } + + for (int i = 1; i < 7; i++) { + x = x + a[i] / pow(T, i); + } + return 1.E5 * exp(x + a[0]); +} + +#endif // SRC_SNAP_THERMODYNAMICS_VAPORS_HYDROGEN_SULFIDE_VAPORS_HPP_ diff --git a/include/cantera/kinetics/vapors/methane_vapors.hpp b/include/cantera/kinetics/vapors/methane_vapors.hpp new file mode 100644 index 0000000000..157c2caec0 --- /dev/null +++ b/include/cantera/kinetics/vapors/methane_vapors.hpp @@ -0,0 +1,29 @@ +#ifndef SRC_SNAP_THERMODYNAMICS_VAPORS_METHANE_VAPORS_HPP_ +#define SRC_SNAP_THERMODYNAMICS_VAPORS_METHANE_VAPORS_HPP_ + +// C/C++ +#include + +inline double svpch4(double t, double p, double beta, double gamma) { + return p * exp((1. - 1. / t) * beta - gamma * log(t)); +} + +inline double sat_vapor_p_CH4_Ideal(double T) { + double betal = 10.15, gammal = 2.1; + double betas = 10.41, gammas = 0.9; + double tr = 90.67, pr = 11690.; + + return T > tr ? svpch4(T / tr, pr, betal, gammal) + : svpch4(T / tr, pr, betas, gammas); +} + +// best fit in [90.99, 189.99] K +inline double sat_vapor_p_CH4_Antoine(double T) { + double A = 3.9895; + double B = 443.028; + double C = -0.49; + double result = pow(10, A - B / (T + C)); + return 1.E5 * result; +} + +#endif // SRC_SNAP_THERMODYNAMICS_VAPORS_METHANE_VAPORS_HPP_ diff --git a/include/cantera/kinetics/vapors/potassium_vapors.hpp b/include/cantera/kinetics/vapors/potassium_vapors.hpp new file mode 100644 index 0000000000..90d7816efc --- /dev/null +++ b/include/cantera/kinetics/vapors/potassium_vapors.hpp @@ -0,0 +1,9 @@ +#ifndef SRC_SNAP_THERMODYNAMICS_VAPORS_POTASSIUM_VAPORS_HPP_ +#define SRC_SNAP_THERMODYNAMICS_VAPORS_POTASSIUM_VAPORS_HPP_ + +inline double sat_vapor_p_KCl_Lodders(double T) { + double logp = 7.611 - 11382. / T; + return 1.E5 * exp(logp); +} + +#endif // SRC_SNAP_THERMODYNAMICS_VAPORS_POTASSIUM_VAPORS_HPP_ diff --git a/include/cantera/kinetics/vapors/sodium_vapors.hpp b/include/cantera/kinetics/vapors/sodium_vapors.hpp new file mode 100644 index 0000000000..8b1ed62a64 --- /dev/null +++ b/include/cantera/kinetics/vapors/sodium_vapors.hpp @@ -0,0 +1,10 @@ +#ifndef SRC_SNAP_THERMODYNAMICS_VAPORS_SODIUM_VAPORS_HPP_ +#define SRC_SNAP_THERMODYNAMICS_VAPORS_SODIUM_VAPORS_HPP_ + +// TODO(cli) check this equation +inline double sat_vapor_p_Na_H2S_Visscher(double T, double pH2S) { + double log10p = 8.55 - 13889. / T - 0.5 * log10(pH2S / 1E5); + return 1.E5 * pow(10., log10p); +} + +#endif // SRC_SNAP_THERMODYNAMICS_VAPORS_SODIUM_VAPORS_HPP_ diff --git a/include/cantera/kinetics/vapors/sulfur_dioxide_vapors.hpp b/include/cantera/kinetics/vapors/sulfur_dioxide_vapors.hpp new file mode 100644 index 0000000000..aba9663336 --- /dev/null +++ b/include/cantera/kinetics/vapors/sulfur_dioxide_vapors.hpp @@ -0,0 +1,16 @@ +#ifndef SRC_SNAP_THERMODYNAMICS_VAPORS_SULFUR_DIOXIDE_VAPORS_HPP_ +#define SRC_SNAP_THERMODYNAMICS_VAPORS_SULFUR_DIOXIDE_VAPORS_HPP_ + +// C/C++ +#include + +// best fit in [177.7, 263] K +inline double sat_vapor_p_SO2_Antoine(double T) { + double A = 3.48586; + double B = 668.225; + double C = -72.252; + double result = pow(10, A - B / (T + C)); + return 1.E5 * result; +} + +#endif // SRC_SNAP_THERMODYNAMICS_VAPORS_SULFUR_DIOXIDE_VAPORS_HPP_ diff --git a/include/cantera/kinetics/vapors/water_vapors.hpp b/include/cantera/kinetics/vapors/water_vapors.hpp new file mode 100644 index 0000000000..34ebc68b10 --- /dev/null +++ b/include/cantera/kinetics/vapors/water_vapors.hpp @@ -0,0 +1,110 @@ +#ifndef SRC_SNAP_THERMODYNAMICS_VAPORS_WATER_VAPORS_HPP_ +#define SRC_SNAP_THERMODYNAMICS_VAPORS_WATER_VAPORS_HPP_ + +#include + +inline double svph2o(double t, double p, double beta, double gamma) { + return p * exp((1. - 1. / t) * beta - gamma * log(t)); +} + +inline double svp_h2o_Umich(double T) { + double x = -2445.5646 / T + 8.2312 * log10(T) - 0.01677006 * T + + 1.20514e-05 * T * T - 6.757169; + return pow(10.0, x) * 1.E5 / 760.; +} + +// (273, 373) +inline double svp_h2o_Antoine(double T) { + double result; + + if (T < 303.) + result = pow(10., 5.40221 - (1838.675 / (T - 31.737))); + else if (T < 333.) + result = pow(10., 5.20389 - (1733.926 / (T - 39.485))); + else if (T < 363.) + result = pow(10., 5.0768 - (1659.793 / (T - 45.854))); + else + result = pow(10., 5.08354 - (1663.125 / (T - 45.622))); + + return 1.E5 * result; +} + +// (100, 373.16) +inline double sat_vapor_p_H2O_Hubner(double T) { + double A, B, C, D; + A = 4.07023; + B = -2484.98; + C = 3.56654; + D = -0.00320981; + double x = A + B / T + C * log10(T) + D * T; // best fit in [100K; 273.16K] + return pow(10.0, x); +} + +// (0, 273.16) +inline double sat_vapor_p_H2O_Fray(double T) { + double a[7], x, pt, Tt, tr; + pt = 6.11657e-03; + Tt = 273.16; + tr = T / Tt; + x = 0.; + a[0] = 20.9969665107897; + a[1] = 3.72437478271362; + a[2] = -13.9205483215524; + a[3] = 29.6988765013566; + a[4] = -40.1972392635944; + a[5] = 29.7880481050215; + a[6] = -9.13050963547721; + + for (int i = 0; i < 7; i++) x = x + a[i] * pow(tr, i); + return 1.E5 * pt * exp(1.5 * log(tr) + (1 - 1 / tr) * x); +} + +inline double sat_vapor_p_H2O_BriggsS(double T) { + double a[6], x; + if (T < 273.16) { + a[1] = -5631.1206; + a[2] = -8.363602; + a[3] = 8.2312; + a[4] = -3.861449e-02; + a[5] = 2.77494e-05; + } else { + a[1] = -2313.0338; + a[2] = -164.03307; + a[3] = 38.053682; + a[4] = -1.3844344e-01; + a[5] = 7.4465367e-05; + } + x = a[1] / T + a[2] + a[3] * log(T) + a[4] * T + a[5] * pow(T, 2); + return exp(x) / 10.; +} + +inline double sat_vapor_p_H2O_Bolton(double T) { + double result; + result = 612.2 * exp(17.67 * (T - 273.15) / (T - 29.65)); + return result; +} + +inline double sat_vapor_p_H2O_Smithsonian(double T) { + double result; + result = 100. * exp(23.33086 - 6111.72784 / T + 0.15215 * log(T)); + return result; +} + +inline double sat_vapor_p_H2O_Ideal(double T) { + double betal = 24.845, gammal = 4.986009, betas = 22.98, gammas = 0.52, + tr = 273.16, pr = 611.7; + return T > tr ? svph2o(T / tr, pr, betal, gammal) + : svph2o(T / tr, pr, betas, gammas); +} + +inline double sat_vapor_p_H2O_liquid_Ideal(double T) { + double betal = 24.845, gammal = 4.986009, tr = 273.16, pr = 611.7; + return svph2o(T / tr, pr, betal, gammal); +} + +inline double sat_vapor_p_H2O_solid_Ideal(double T) { + double betas = 22.98, gammas = 0.52, tr = 273.16, pr = 611.7; + return svph2o(T / tr, pr, betas, gammas); +} + +#endif // SRC_SNAP_THERMODYNAMICS_VAPORS_WATER_VAPORS_HPP_ diff --git a/include/cantera/thermo/IdealGasSurfPhase.h b/include/cantera/thermo/IdealGasSurfPhase.h new file mode 100644 index 0000000000..edd207e78e --- /dev/null +++ b/include/cantera/thermo/IdealGasSurfPhase.h @@ -0,0 +1,47 @@ +#ifndef CT_IDEALGASSURF_H +#define CT_IDEALGASSURF_H + +#include "IdealGasPhase.h" +#include "SurfPhase.h" + +namespace Cantera +{ + +class IdealGasSurfPhase: public SurfPhase +{ + public: + explicit IdealGasSurfPhase(const string& inputFile="", const string& id=""): + SurfPhase(inputFile, id) + {} + + bool isCompressible() const override { + return true; + } + + string type() const override { + return "ideal-gas-surface"; + } + + string phaseOfMatter() const override { + return "gas"; + } + + double pressure() const override { + return GasConstant * molarDensity() * temperature(); + } + + void setPressure(double p) override { + setDensity(p * meanMolecularWeight() / RT()); + } + + protected: + // dimensionless parameters + double m_h0_RT; + double m_cp0_R; + double m_g0_RT; + double m_s0_R; +}; + +} + +#endif diff --git a/include/cantera/thermo/IdealMoistPhase.h b/include/cantera/thermo/IdealMoistPhase.h new file mode 100644 index 0000000000..2109c9d207 --- /dev/null +++ b/include/cantera/thermo/IdealMoistPhase.h @@ -0,0 +1,100 @@ +#ifndef CT_IDEALMOISTPHASE_H +#define CT_IDEALMOISTPHASE_H + +#include "IdealGasPhase.h" + +namespace Cantera +{ + +class IdealMoistPhase : public IdealGasPhase { + public: + explicit IdealMoistPhase(const std::string& infile = "", const std::string& id=""): + IdealGasPhase(infile, id) {} + + string type() const override { + return "ideal-moist"; + } + + void setTemperature(double temp) override { + IdealGasPhase::setTemperature(temp); + m_pressure_sets_temperature = false; + } + + void setDensity(double dens) override { + IdealGasPhase::setDensity(dens); + m_pressure_sets_temperature = true; + } + + /*! + * @f[ + * \hat s(T, P) = \sum_k X_k \hat s^0_k(T) - \hat R \ln \frac{P}{P^0}. + * - \sum_{k \elem G} X_k \ln X_k + * @f] + */ + double entropy_mole() const override { + return GasConstant * (mean_X(entropy_R_ref()) - _sum_xlogx_g() + - std::log(pressure() / refPressure())); + } + + /*! + * @f[ + * \hat c_v = \hat c_p - g * \hat R. + * @f] + */ + double cv_mole() const override { + return cp_mole() - _g_ov_mu() * meanMolecularWeight() * GasConstant; + } + + double pressure() const override; + + void setPressure(double p) override; + + bool addSpecies(shared_ptr spec) override; + + void setState_DP(double rho, double p) override { + setDensity(rho); + setPressure(p); + } + + void getChemPotentials(double* mu) const override { + getStandardChemPotentials(mu); + } + + void getCv_R(double* cvr) const override; + + void getIntEnergy_RT_ref(double* urt) const override; + + size_t nGas() const { + return nSpecies() - m_ncloud; + } + + protected: + /*! + * @f[ + * \frac{\sum_{i \elem G} X_i}{\bar{\mu}} + * @f] + */ + double _g_ov_mu() const; + + /*! + * @f[ + * \frac{\sum_{i \elem G} X_i \ln X_i}{\sum_{i \elem G} X_i} + * @f] + */ + double _sum_xlogx_g() const; + + bool m_pressure_sets_temperature = false; + + //! Number of cloud species + size_t m_ncloud = 0; + + //! Vector of molar volumes for each species in the solution + /** + * Species molar volumes (@f$ m^3 kmol^-1 @f$) at the current mixture state. + */ + mutable vector m_speciesMolarVolume; +}; + +} + +#endif diff --git a/include/cantera/thermo/LatticeSolidPhase.h b/include/cantera/thermo/LatticeSolidPhase.h index 040d7c8708..94d6a2df37 100644 --- a/include/cantera/thermo/LatticeSolidPhase.h +++ b/include/cantera/thermo/LatticeSolidPhase.h @@ -292,11 +292,11 @@ class LatticeSolidPhase : public ThermoPhase */ void getMoleFractions(double* const x) const; - void setMassFractions(const double* const y) override { + void setMassFractions(const double* const y, size_t stride = 1) override { throw NotImplementedError("LatticeSolidPhase::setMassFractions"); } - void setMassFractions_NoNorm(const double* const y) override { + void setMassFractions_NoNorm(const double* const y, size_t stride = 1) override { throw NotImplementedError("LatticeSolidPhase::setMassFractions_NoNorm"); } diff --git a/include/cantera/thermo/Phase.h b/include/cantera/thermo/Phase.h index 116edb4be4..50b1d92b8a 100644 --- a/include/cantera/thermo/Phase.h +++ b/include/cantera/thermo/Phase.h @@ -448,14 +448,18 @@ class Phase //! must be greater than or equal to the number of //! species. The Phase object will normalize this vector //! before storing its contents. - virtual void setMassFractions(const double* const y); + virtual void setMassFractions(const double* const y, size_t stride = 1); //! Set the mass fractions to the specified values without normalizing. //! This is useful when the normalization condition is being handled by //! some other means, for example by a constraint equation as part of a //! larger set of equations. //! @param y Input vector of mass fractions. Length is m_kk. - virtual void setMassFractions_NoNorm(const double* const y); + virtual void setMassFractions_NoNorm(const double* const y, size_t stride = 1); + + + //! Set the mass fractions to the specified values without normalizing. + virtual void setMassFractionsPartial(const double* const y, size_t stride = 1); //! Get the species concentrations (kmol/m^3). /*! @@ -874,6 +878,10 @@ class Phase //! Flag determining whether case sensitive species names are enforced bool m_caseSensitiveSpecies = false; + vector const& ym() const { + return m_ym; + } + private: //! Find lowercase species name in m_speciesIndices when case sensitive //! species names are not enforced and a user specifies a non-canonical diff --git a/include/cantera/thermo/SingleSpeciesTP.h b/include/cantera/thermo/SingleSpeciesTP.h index 24647ddbcc..6be89d2df9 100644 --- a/include/cantera/thermo/SingleSpeciesTP.h +++ b/include/cantera/thermo/SingleSpeciesTP.h @@ -215,7 +215,7 @@ class SingleSpeciesTP : public ThermoPhase //! @{ //! Mass fractions are fixed, with Y[0] = 1.0. - void setMassFractions(const double* const y) override {}; + void setMassFractions(const double* const y, size_t stride = 1) override {}; //! Mole fractions are fixed, with x[0] = 1.0. void setMoleFractions(const double* const x) override {}; diff --git a/include/cantera/thermo/ThermoPhase.h b/include/cantera/thermo/ThermoPhase.h index 635015b57a..a3daf73188 100644 --- a/include/cantera/thermo/ThermoPhase.h +++ b/include/cantera/thermo/ThermoPhase.h @@ -925,6 +925,17 @@ class ThermoPhase : public Phase throw NotImplementedError("ThermoPhase::getCp_R"); } + //! Get the nondimensional Heat Capacities at constant volume for the + //! species standard states at the current *T* and *P* of the + //! solution + /*! + * @param cpr Output vector of nondimensional standard state heat + * capacities. Length: m_kk. + */ + virtual void getCv_R(double* cvr) const { + throw NotImplementedError("ThermoPhase::getCv_R"); + } + //! Get the molar volumes of the species standard states at the current //! *T* and *P* of the solution. /*! diff --git a/install_cantera.sh b/install_cantera.sh new file mode 100755 index 0000000000..7f03c6facc --- /dev/null +++ b/install_cantera.sh @@ -0,0 +1,23 @@ +#! /bin/bash + +CXX=g++ +CC=gcc +cxx_flags=-std=c++17 +prefix=${HOME}/opt/ +python_package=y +f90_interface=n +system_eigen=n +system_blas_lapack=n +boost_inc_dir=`pwd`/cantera/ext/cliboost/ + +cd cantera +git submodule init ext/cliboost +git submodule update ext/cliboost + +scons build CXX=${CXX} CC=${CC} cxx_flags=${cxx_flags} prefix=${prefix} \ + python_package=${python_package} f90_interface=${f90_interface} \ + system_eigen=${system_eigen} \ + system_blas_lapack=${system_blas_lapack} \ + boost_inc_dir=${boost_inc_dir} -j8 +scons install +cd .. diff --git a/src/kinetics/Condensation.cpp b/src/kinetics/Condensation.cpp new file mode 100644 index 0000000000..de96835529 --- /dev/null +++ b/src/kinetics/Condensation.cpp @@ -0,0 +1,455 @@ +// C/C++ +#include +#include + +#include "cantera/kinetics/Reaction.h" +#include "cantera/kinetics/Condensation.h" +#include "cantera/thermo/ThermoPhase.h" +#include "cantera/thermo/IdealMoistPhase.h" + +namespace Cantera +{ + +// x -> y at constant volume (mole concentration) +inline pair satfunc1v(double s, double x, double y, + double logs_ddT = 0.) +{ + if (y < 0.) return {-y, 0.}; + + double rate = x - s; + if (rate > 0. || (rate < 0. && y > - rate)) { + return {rate, - s * logs_ddT}; + } + return {-y, 0.}; +} + +inline void set_jac1v( + Eigen::SparseMatrix &jac, double s, double const *conc, + int j, int ix, int iy) +{ + double const& x = conc[ix]; + double const& y = conc[iy]; + double rate = x - s; + if (rate > 0. || (rate < 0. && y > - rate)) { + jac.coeffRef(j, ix) = 1.; + } else { + jac.coeffRef(j, iy) = -1.; + } +} + +// x1 + x2 -> y at constant volume (mole concentration) +inline pair satfunc2v(double s, double x1, double x2, double y, + double logs_ddT = 0.) +{ + if (y < 0.) return {-y, 0.}; + + double delta = (x1 - x2) * (x1 - x2) + 4 * s; + double rate = (x1 + x2 - sqrt(delta)) / 2.; + + if (rate > 0. || (rate < 0. && y > - rate)) { + return {rate, - s * logs_ddT / sqrt(delta)}; + } + return {-y, 0.}; +} + +inline void set_jac2v( + Eigen::SparseMatrix &jac, double s, double const *conc, + int j, int ix1, int ix2, int iy) +{ + double const& x1 = conc[ix1]; + double const& x2 = conc[ix2]; + double const& y = conc[iy]; + double delta = (x1 - x2) * (x1 - x2) + 4 * s; + double rate = (x1 + x2 - sqrt(delta)) / 2.; + + if (rate > 0. || (rate < 0. && y > - rate)) { + jac.coeffRef(j, ix1) = (1. - (x1 - x2) / sqrt(delta)) / 2.; + jac.coeffRef(j, ix2) = (1. - (x2 - x1) / sqrt(delta)) / 2.; + } else { + jac.coeffRef(j, iy) = -1.; + } +} + +// x -> y at constant pressure (mole fraction) +inline double satfunc1p(double s, double x, double y, double g) +{ + // boil all condensates + if (s > 1.) { + return -y; + } + + //std::cout << "s = " << s << ", x = " << x << ", y = " << y << ", g = " << g << std::endl; + + double rate = (x - s * g) / (1. - s); + + if (rate > 0. || (rate < 0. && y > - rate)) { + return rate; + } + return -y; +} + +inline void set_jac1p(Eigen::SparseMatrix &jac, + double s, double const* frac, double g, + int j, int ix, int iy) +{ + // boil all condensates + if (s > 1.) { + jac.coeffRef(j, iy) = -1.; + return; + } + + double const& x = frac[ix]; + double const& y = frac[iy]; + + double rate = (x - s * g) / (1. - s); + + if (rate > 0. || (rate < 0. && y > - rate)) { + jac.coeffRef(j, ix) = 1.; + } else { + jac.coeffRef(j, iy) = -1.; + } +} + +// x1 + x2 -> y at constant pressure (mole fraction) +inline double satfunc2p(double s, double x1, double x2, double y, double g) +{ + double delta = (x1 - x2) * (x1 - x2) + 4 * s * (g - 2. * x1) * (g - 2. * x2); + double rate = (x1 + x2 - 4 * g * s - sqrt(delta)) / (2. * (1. - 4. * s)); + + if (rate > 0. || (rate < 0. && y > - rate)) { + return rate; + } + return -y; +} + +inline void set_jac2p(Eigen::SparseMatrix &jac, + double s, double const* frac, double g, + int j, int ix1, int ix2, int iy) +{ + double const& x1 = frac[ix1]; + double const& x2 = frac[ix2]; + double const& y = frac[iy]; + + double delta = (x1 - x2) * (x1 - x2) + 4 * s * (g - 2. * x1) * (g - 2. * x2); + double rate = (x1 + x2 - 4 * g * s - sqrt(delta)) / (2. * (1. - 4. * s)); + + if (rate > 0. || (rate < 0. && y > - rate)) { + jac.coeffRef(j, ix1) = (1. - (x1 - x2) / sqrt(delta)) / 2.; + jac.coeffRef(j, ix2) = (1. - (x2 - x1) / sqrt(delta)) / 2.; + } else { + jac.coeffRef(j, iy) = -1.; + } +} + +inline Eigen::VectorXd linear_solve_rop( + Eigen::SparseMatrix const& jac, + Eigen::SparseMatrix const& stoich, + Eigen::VectorXd const& b) +{ + auto A = jac * stoich; + + Eigen::SparseQR, Eigen::COLAMDOrdering> solver; + solver.compute(A); + Eigen::VectorXd r = solver.solve(b); + + if (solver.info() != Eigen::Success) { + throw CanteraError("Condensation::updateROP", + "Failed to solve for net rates of progress."); + } + + /*std::cout << jac << std::endl; + std::cout << A << std::endl; + std::cout << A.transpose() * A << std::endl; + std::cout << "b = " << b << std::endl; + std::cout << "r = " << -r << std::endl;*/ + + return r; +} + +void Condensation::resizeReactions() +{ + Kinetics::resizeReactions(); + + for (auto& rates : m_interfaceRates) { + rates->resize(nTotalSpecies(), nReactions(), nPhases()); + } + m_jac.resize(nReactions(), nTotalSpecies()); + m_rfn_ddT.resize(nReactions()); +} + +void Condensation::resizeSpecies() +{ + size_t kOld = m_kk; + Kinetics::resizeSpecies(); + if (m_kk != kOld && nReactions()) { + throw CanteraError("InterfaceKinetics::resizeSpecies", "Cannot add" + " species to InterfaceKinetics after reactions have been added."); + } + + m_conc.resize(m_kk); + m_intEng.resize(m_kk); + m_cv.resize(m_kk); +} + +void Condensation::getActivityConcentrations(double* const pdata) +{ + if (m_use_mole_fraction) { + _update_rates_X(pdata); + } else { + _update_rates_C(pdata); + } +} + +void Condensation::getFwdRateConstants(double* kfwd) +{ + _update_rates_T(kfwd, nullptr); +} + +void Condensation::getFwdRateConstants_ddT(double* kfwd) +{ + _update_rates_T(nullptr, kfwd); +} + +bool Condensation::addReaction(shared_ptr r_base, bool resize) +{ + size_t i = nReactions(); + bool added = Kinetics::addReaction(r_base, resize); + if (!added) { + return false; + } + + // Set index of rate to number of reaction within kinetics + shared_ptr rate = r_base->rate(); + rate->setRateIndex(nReactions() - 1); + rate->setContext(*r_base, *this); + + string rtype = rate->subType(); + if (rtype == "") { + rtype = rate->type(); + } + + if (rtype == "nucleation") { + if (r_base->reactants.size() == 1) { + m_jxy.push_back(i); + } else if (r_base->reactants.size() == 2) { + m_jxxy.push_back(i); + } + } else if (rtype == "freezing") { + m_jyy.push_back(i); + } else { + throw CanteraError("Condensation::addReaction", + "Unknown reaction type '{}'", rtype); + } + + // If necessary, add new interface MultiRate evaluator + if (m_interfaceTypes.find(rtype) == m_interfaceTypes.end()) { + m_interfaceTypes[rtype] = m_interfaceRates.size(); + m_interfaceRates.push_back(rate->newMultiRate()); + m_interfaceRates.back()->resize(m_kk, nReactions(), nPhases()); + } + + // Add reaction rate to evaluator + size_t index = m_interfaceTypes[rtype]; + m_interfaceRates[index]->add(nReactions() - 1, *rate); + + return true; +} + +void Condensation::updateROP() { + _update_rates_T(m_rfn.data(), m_rfn_ddT.data()); + if (m_use_mole_fraction) { + _update_rates_X(m_conc.data()); + } else { + _update_rates_C(m_conc.data()); + } + + if (m_ROP_ok) { + return; + } + + m_jac.setZero(); + + double pres = thermo().pressure(); + double temp = thermo().temperature(); + double dens = pres / (GasConstant * temp); + double xgas = 0.; + + if (m_use_mole_fraction) { + size_t ngas = static_cast(thermo()).nGas(); + for (size_t i = 0; i < ngas; i++) + xgas += m_conc[i]; + //std::cout << "xgas = " << xgas << std::endl; + } + + Eigen::VectorXd b(nReactions()); + Eigen::VectorXd b_ddT(nReactions()); + Eigen::SparseMatrix stoich(m_stoichMatrix); + Eigen::SparseMatrix rate_ddT(nReactions(), nTotalSpecies()); + + b.setZero(); + b_ddT.setZero(); + rate_ddT.setZero(); + + // nucleation: x <=> y + for (auto j : m_jxy) { + // std::cout << "jxy = " << j << std::endl; + // inactive reactions + if (m_rfn[j] < 0.0) { + for (int i = 0; i < nTotalSpecies(); i++) + stoich.coeffRef(i,j) = 0.0; + continue; + } + + auto& R = m_reactions[j]; + size_t ix = kineticsSpeciesIndex(R->reactants.begin()->first); + size_t iy = kineticsSpeciesIndex(R->products.begin()->first); + + if (m_use_mole_fraction) { + b(j) = satfunc1p(m_rfn[j] / dens, m_conc[ix], m_conc[iy], xgas); + set_jac1p(m_jac, m_rfn[j] / dens, m_conc.data(), xgas, j, ix, iy); + } else { + auto result = satfunc1v(m_rfn[j], m_conc[ix], m_conc[iy], m_rfn_ddT[j]); + b(j) = result.first; + b_ddT(j) = result.second; + set_jac1v(m_jac, m_rfn[j], m_conc.data(), j, ix, iy); + } + } + + // nucleation: x1 + x2 <=> y + for (auto j : m_jxxy) { + //std::cout << "jxxy = " << j << std::endl; + // inactive reactions + if (m_rfn[j] < 0.0) { + for (int i = 0; i < nTotalSpecies(); i++) + stoich.coeffRef(i,j) = 0.0; + continue; + } + + auto& R = m_reactions[j]; + size_t ix1 = kineticsSpeciesIndex(R->reactants.begin()->first); + size_t ix2 = kineticsSpeciesIndex(next(R->reactants.begin())->first); + size_t iy = kineticsSpeciesIndex(R->products.begin()->first); + + if (m_use_mole_fraction) { + b(j) = satfunc2p(m_rfn[j] / (dens * dens), m_conc[ix1], m_conc[ix2], m_conc[iy], xgas); + set_jac2p(m_jac, m_rfn[j] / (dens * dens), m_conc.data(), xgas, j, ix1, ix2, iy); + } else { + auto result = satfunc2v(m_rfn[j], m_conc[ix1], m_conc[ix2], m_conc[iy], m_rfn_ddT[j]); + b(j) = result.first; + b_ddT(j) = result.second; + set_jac2v(m_jac, m_rfn[j], m_conc.data(), j, ix1, ix2, iy); + } + } + + // freezing: y1 <=> y2 + for (auto j : m_jyy) { + auto& R = m_reactions[j]; + size_t iy1 = kineticsSpeciesIndex(R->reactants.begin()->first); + size_t iy2 = kineticsSpeciesIndex(R->products.begin()->first); + + if (temp > m_rfn[j]) { // higher than freezing temperature + if (m_conc[iy2] > 0.) { + b(j) = - m_conc[iy2]; + m_jac.coeffRef(j, iy2) = -1.; + } + } else { // lower than freezing temperature + if (m_conc[iy1] > 0.) { + b(j) = m_conc[iy1]; + m_jac.coeffRef(j, iy1) = 1.; + } + } + } + + //std::cout << "b_ddt = " << b_ddT << std::endl; + + // set up temperature gradient + if (!m_use_mole_fraction) { + for (size_t j = 0; j != nReactions(); ++j) { + if (b_ddT[j] != 0.) { + for (size_t i = 0; i != nTotalSpecies(); ++i) + rate_ddT.coeffRef(j, i) = b_ddT[j] * m_intEng[i]; + } + } + + /*std::cout << "u = " << m_intEng << std::endl; + std::cout << "cc = " << m_conc.dot(m_cv) << std::endl; + std::cout << "cv = " << m_cv << std::endl;*/ + m_jac -= rate_ddT / m_conc.dot(m_cv); + } + + // solve the optimal net rates + auto r = linear_solve_rop(m_jac, stoich, b); + + // scale rate down if some species becomes negative + Eigen::VectorXd rates = - stoich * r; + + for (size_t j = 0; j != nReactions(); ++j) { + m_ropf[j] = std::max(0., -r(j)); + m_ropr[j] = std::max(0., r(j)); + m_ropnet[j] = m_ropf[j] - m_ropr[j]; + } + + m_ROP_ok = true; +} + +void Condensation::_update_rates_T(double *pdata, double *pdata_ddT) +{ + // Go find the temperature from the surface + double T = thermo().temperature(); + + if (T != m_temp) { + m_temp = T; + m_ROP_ok = false; + } + + if (pdata_ddT != nullptr) { + std::fill(pdata_ddT, pdata_ddT + nReactions(), 1.); + } + + // loop over interface MultiRate evaluators for each reaction type + for (auto& rates : m_interfaceRates) { + bool changed = rates->update(thermo(), *this); + if (changed) { + if (pdata != nullptr) { + rates->getRateConstants(pdata); + } + if (pdata_ddT != nullptr) { + rates->processRateConstants_ddT(pdata_ddT, nullptr, 0.); + } + m_ROP_ok = false; + } + } +} + +void Condensation::_update_rates_C(double *pdata) +{ + thermo().getActivityConcentrations(pdata); + thermo().getIntEnergy_RT(m_intEng.data()); + thermo().getCv_R(m_cv.data()); + + for (size_t i = 0; i < m_kk; i++) { + m_intEng[i] *= thermo().temperature(); + } + + m_ROP_ok = false; +} + +void Condensation::_update_rates_X(double *pdata) +{ + thermo().getMoleFractions(pdata); + m_ROP_ok = false; +} + +Eigen::SparseMatrix Condensation::netRatesOfProgress_ddX() +{ + updateROP(); + return m_jac; +} + +Eigen::SparseMatrix Condensation::netRatesOfProgress_ddCi() +{ + updateROP(); + return m_jac; +} + +} diff --git a/src/kinetics/Freezing.cpp b/src/kinetics/Freezing.cpp new file mode 100644 index 0000000000..3102f25963 --- /dev/null +++ b/src/kinetics/Freezing.cpp @@ -0,0 +1,62 @@ +#include "cantera/numerics/Func1.h" +#include "cantera/base/stringUtils.h" +#include "cantera/kinetics/Freezing.h" +#include "cantera/kinetics/Reaction.h" + +namespace Cantera +{ + +FreezingRate::FreezingRate(const AnyMap& node, const UnitStack& rate_units) + : FreezingRate() +{ + setParameters(node, rate_units); + + if (!node.hasKey("rate-constant")) { + throw CanteraError("FreezingRate::FreezingRate", + "Missing 'rate-constant' key in rate node."); + } + + setRateParameters(node["equation"], node["rate-constant"], node); +} + +void FreezingRate::setRateParameters( + const AnyValue& equation, + const AnyValue& rate, + const AnyMap& node) +{ + if (rate.empty()) { + throw InputFileError("FreezingRate::setRateParameters", rate, + "Missing rate constant data."); + } + + if (!rate.is()) { + throw InputFileError("FreezingRate::setRateParameters", rate, + "Expected a parameter map."); + } + + auto& rate_map = rate.as(); + + m_t3 = rate_map["T3"].asDouble(); + m_valid = true; +} + +void FreezingRate::validate(const string& equation, const Kinetics& kin) +{ + if (!m_valid) { + throw InputFileError("FreezingRate::validate", m_input, + "Rate object for reaction '{}' is not configured.", equation); + } +} + +double FreezingRate::evalFromStruct(const ArrheniusData& shared_data) const +{ + return m_t3; +} + +void FreezingRate::getParameters(AnyMap& rateNode, const Units& rate_units) const +{ + throw NotImplementedError("FreezingRate::getParameters", + "Not implemented by '{}' object.", type()); +} + +} diff --git a/src/kinetics/KineticsFactory.cpp b/src/kinetics/KineticsFactory.cpp index db4d1b6a38..ae91489247 100644 --- a/src/kinetics/KineticsFactory.cpp +++ b/src/kinetics/KineticsFactory.cpp @@ -9,6 +9,7 @@ #include "cantera/kinetics/BulkKinetics.h" #include "cantera/kinetics/InterfaceKinetics.h" #include "cantera/kinetics/EdgeKinetics.h" +#include "cantera/kinetics/Condensation.h" #include "cantera/kinetics/Reaction.h" #include "cantera/thermo/ThermoPhase.h" #include "cantera/base/stringUtils.h" @@ -35,6 +36,7 @@ KineticsFactory::KineticsFactory() { addDeprecatedAlias("surface", "surf"); reg("edge", []() { return new EdgeKinetics(); }); addDeprecatedAlias("edge", "Edge"); + reg("condensation", []() { return new Condensation(); }); } KineticsFactory* KineticsFactory::factory() { diff --git a/src/kinetics/Nucleation.cpp b/src/kinetics/Nucleation.cpp new file mode 100644 index 0000000000..c25bc8eec9 --- /dev/null +++ b/src/kinetics/Nucleation.cpp @@ -0,0 +1,157 @@ +#include "cantera/numerics/Func1.h" +#include "cantera/base/stringUtils.h" +#include "cantera/kinetics/Nucleation.h" +#include "cantera/kinetics/svp_funcs.h" +#include "cantera/kinetics/Reaction.h" + +namespace Cantera +{ + +string concatenate(const Composition& comp, char sep) { + if (comp.empty()) return ""; + if (comp.size() == 1) return comp.begin()->first; + + std::string result = comp.begin()->first; + + for (auto it = std::next(comp.begin()); it != comp.end(); ++it) { + result += sep + it->first; + } + + return result; +} + +std::function find_svp(const Composition& reactants, + const string& svp_name) +{ + string name = concatenate(reactants, '-') + '-' + svp_name; + + if (name == "NH3-H2S-lewis" || name == "H2S-NH3-lewis") { + return svp_nh3_h2s_Lewis; + } else if (name == "H2O-antoine") { + return svp_h2o_Antoine; + } else if (name == "NH3-antoine") { + return svp_nh3_Antoine; + } + + throw CanteraError("find_svp", + "No SVP function found for reaction '{}'.", name); +} + +std::function find_logsvp_ddT(const Composition& reactants, + const string& svp_name) +{ + string name = concatenate(reactants, '-') + '-' + svp_name; + + if (name == "NH3-H2S-lewis" || name == "H2S-NH3-lewis") { + return logsvp_ddT_nh3_h2s_Lewis; + } else if (name == "H2O-antoine") { + return nullptr; // logsvp_ddT_h2o_Antoine; + } else if (name == "NH3-antoine") { + return nullptr; // logsvp_ddT_nh3_Antoine; + } + + throw CanteraError("find_logsvp_ddT", + "No SVP_DDT function found for reaction '{}'.", name); +} + +NucleationRate::NucleationRate(const AnyMap& node, const UnitStack& rate_units) + : NucleationRate() +{ + setParameters(node, rate_units); + + if (!node.hasKey("rate-constant")) { + throw CanteraError("NucleationRate::NucleationRate", + "Missing 'rate-constant' key in rate node."); + } + + setRateParameters(node["equation"], node["rate-constant"], node); +} + +void NucleationRate::setRateParameters( + const AnyValue& equation, + const AnyValue& rate, + const AnyMap& node) +{ + if (rate.empty()) { + throw InputFileError("NucleationRate::setRateParameters", rate, + "Missing rate constant data."); + } + + if (!rate.is()) { + throw InputFileError("NucleationRate::setRateParameters", rate, + "Expected a parameter map."); + } + + auto& rate_map = rate.as(); + string svp_name = rate_map[m_formula_str].asString(); + + if (rate.hasKey("minT")) { + m_min_temp = rate_map["minT"].asDouble(); + }; + + if (rate.hasKey("maxT")) { + m_max_temp = rate_map["maxT"].asDouble(); + }; + + Reaction rtmp; + parseReactionEquation(rtmp, equation.asString(), node, nullptr); + m_order = rtmp.reactants.size(); + + if (svp_name == "ideal") { + m_t3 = rate_map["T3"].asDouble(); + m_p3 = rate_map["P3"].asDouble(); + m_beta = rate_map["beta"].asDouble(); + m_delta = rate_map["delta"].asDouble(); + m_svp = [this](double T) { + return m_p3 * exp((1. - m_t3 / T) * m_beta - m_delta * log(T / m_t3)); + }; + m_logsvp_ddT = [this](double T) { + return m_beta * m_t3 / (T * T) - m_delta / T; + }; + if (m_delta > 0.) { + m_max_temp = m_beta * m_t3 / m_delta; + } + } else { + m_svp = find_svp(rtmp.reactants, svp_name); + m_logsvp_ddT = find_logsvp_ddT(rtmp.reactants, svp_name); + } + + m_valid = true; +} + +void NucleationRate::validate(const string& equation, const Kinetics& kin) +{ + if (!m_valid) { + throw InputFileError("NucleationRate::validate", m_input, + "Rate object for reaction '{}' is not configured.", equation); + } +} + +double NucleationRate::evalFromStruct(const ArrheniusData& shared_data) const +{ + double T = shared_data.temperature; + if (T < m_min_temp || T > m_max_temp) { + return -1; + } + + double RT = GasConstant * T; + return m_svp(T) / pow(RT, m_order); +} + +double NucleationRate::ddTScaledFromStruct(const ArrheniusData& shared_data) const +{ + double T = shared_data.temperature; + if (T < m_min_temp || T > m_max_temp) { + return 0.; + } + + return m_logsvp_ddT(T) - m_order / T; +} + +void NucleationRate::getParameters(AnyMap& rateNode, const Units& rate_units) const +{ + throw NotImplementedError("NucleationRate::getParameters", + "Not implemented by '{}' object.", type()); +} + +} diff --git a/src/kinetics/ReactionRateFactory.cpp b/src/kinetics/ReactionRateFactory.cpp index 605447768d..3653a02a3a 100644 --- a/src/kinetics/ReactionRateFactory.cpp +++ b/src/kinetics/ReactionRateFactory.cpp @@ -18,6 +18,8 @@ #include "cantera/kinetics/PlogRate.h" #include "cantera/kinetics/TwoTempPlasmaRate.h" #include "cantera/kinetics/Photolysis.h" +#include "cantera/kinetics/Nucleation.h" +#include "cantera/kinetics/Freezing.h" namespace Cantera { @@ -110,6 +112,18 @@ ReactionRateFactory::ReactionRateFactory() reg("photolysis", [](const AnyMap& node, const UnitStack& rate_units) { return new PhotolysisRate(node, rate_units); }); + + // Nucleation evaluator + reg("nucleation", [](const AnyMap& node, const UnitStack& rate_units) { + return new NucleationRate(node, rate_units); + }); + addAlias("nucleation", "interface-nucleation"); + + // Freezing evaluator + reg("freezing", [](const AnyMap& node, const UnitStack& rate_units) { + return new FreezingRate(node, rate_units); + }); + addAlias("freezing", "interface-freezing"); } ReactionRateFactory* ReactionRateFactory::factory() { diff --git a/src/thermo/IdealGasSurfPhase.cpp b/src/thermo/IdealGasSurfPhase.cpp new file mode 100644 index 0000000000..518da38f5f --- /dev/null +++ b/src/thermo/IdealGasSurfPhase.cpp @@ -0,0 +1,21 @@ +#include "cantera/thermo/IdealGasSurfPhase.h" + +namespace Cantera +{ + +/*double IdealGasSurfPhase::entropy_mole() const +{ + return GasConstant * (mean_X(entropy_R_ref()) - sum_xlogx() - std::log(pressure() / refPressure())); +} + +double IdealGasSurfPhase::cp_mole() const +{ + return GasConstant * mean_X(cp_R_ref()); +} + +double IdealGasSurfPhase::cv_mole() const +{ + return cp_mole() - GasConstant; +}*/ + +} diff --git a/src/thermo/IdealMoistPhase.cpp b/src/thermo/IdealMoistPhase.cpp new file mode 100644 index 0000000000..3de2e93996 --- /dev/null +++ b/src/thermo/IdealMoistPhase.cpp @@ -0,0 +1,87 @@ +#include "cantera/thermo/IdealMoistPhase.h" +#include "cantera/thermo/Species.h" +#include "cantera/base/utilities.h" + +namespace Cantera +{ + +double IdealMoistPhase::_g_ov_mu() const { + return std::accumulate(ym().begin(), ym().begin() + nGas(), 0.0); +} + +double IdealMoistPhase::_sum_xlogx_g() const { + double sumxlogx = 0; + double sumx = 0; + for (size_t k = 0; k < nGas(); k++) { + sumxlogx += ym()[k] * std::log(std::max(ym()[k], SmallNumber)); + sumx += ym()[k]; + } + return sumxlogx / sumx + std::log(meanMolecularWeight()); +} + +void IdealMoistPhase::setPressure(double p) { + if (p <= 0) { + throw CanteraError("IdealGasPhase::setState_DP", + "pressure must be positive"); + } + + if (m_pressure_sets_temperature) { + setTemperature(p / (_g_ov_mu() * GasConstant * density())); + } else { + setDensity(p / (GasConstant * temperature() * _g_ov_mu())); + } +} + +double IdealMoistPhase::pressure() const { + return GasConstant * temperature() * density() * _g_ov_mu(); +} + +bool IdealMoistPhase::addSpecies(shared_ptr spec) { + bool added = IdealGasPhase::addSpecies(spec); + if (spec->input.hasKey("equation-of-state")) { + auto& eos = spec->input["equation-of-state"].getMapWhere("model", "constant-volume"); + double mv; + if (eos.hasKey("density")) { + mv = molecularWeight(m_kk-1) / eos.convert("density", "kg/m^3"); + } else if (eos.hasKey("molar-density")) { + mv = 1.0 / eos.convert("molar-density", "kmol/m^3"); + } else if (eos.hasKey("molar-volume")) { + mv = eos.convert("molar-volume", "m^3/kmol"); + } else { + throw CanteraError("IdealSolidSolnPhase::addSpecies", + "equation-of-state entry for species '{}' is missing " + "'density', 'molar-volume', or 'molar-density' " + "specification", spec->name); + } + m_speciesMolarVolume.push_back(mv); + m_ncloud++; + } else { + m_speciesMolarVolume.push_back(0.); + } +} + +void IdealMoistPhase::getCv_R(double* cvr) const +{ + const vector& _cpr = cp_R_ref(); + // gas + for (size_t k = 0; k < nGas(); k++) { + cvr[k] = _cpr[k] - 1.0; + } + + // clouds + copy(_cpr.begin() + nGas(), _cpr.end(), cvr + nGas()); +} + +void IdealMoistPhase::getIntEnergy_RT_ref(double* urt) const { + const vector& _h = enthalpy_RT_ref(); + + // gas + for (size_t k = 0; k < nGas(); k++) { + urt[k] = _h[k] - 1.0; + } + + // clouds + copy(_h.begin() + nGas(), _h.end(), urt + nGas()); +} + +} diff --git a/src/thermo/Phase.cpp b/src/thermo/Phase.cpp index dc0f6a4fd3..f8d9f16a7f 100644 --- a/src/thermo/Phase.cpp +++ b/src/thermo/Phase.cpp @@ -338,10 +338,10 @@ void Phase::setMoleFractionsByName(const string& x) setMoleFractionsByName(parseCompString(x)); } -void Phase::setMassFractions(const double* const y) +void Phase::setMassFractions(const double* const y, size_t stride) { for (size_t k = 0; k < m_kk; k++) { - m_y[k] = std::max(y[k], 0.0); // Ignore negative mass fractions + m_y[k] = std::max(y[k*stride], 0.0); // Ignore negative mass fractions } double norm = accumulate(m_y.begin(), m_y.end(), 0.0); scale(m_y.begin(), m_y.end(), m_y.begin(), 1.0/norm); @@ -352,10 +352,17 @@ void Phase::setMassFractions(const double* const y) compositionChanged(); } -void Phase::setMassFractions_NoNorm(const double* const y) +void Phase::setMassFractions_NoNorm(const double* const y, size_t stride) { double sum = 0.0; - copy(y, y + m_kk, m_y.begin()); + if (stride == 1) { + copy(y, y + m_kk, m_y.begin()); + } else { + for (size_t k = 0; k < m_kk; k++) { + m_y[k] = std::max(y[k*stride], 0.0); // Ignore negative mass fractions + } + } + transform(m_y.begin(), m_y.end(), m_rmolwts.begin(), m_ym.begin(), multiplies()); sum = accumulate(m_ym.begin(), m_ym.end(), 0.0); @@ -363,6 +370,23 @@ void Phase::setMassFractions_NoNorm(const double* const y) compositionChanged(); } +void Phase::setMassFractionsPartial(const double* const y, size_t stride) +{ + double sum = 0.0; + + // fix negative concentration + for (size_t k = 1; k < m_kk; k++) { + m_y[k] = std::max(y[(k - 1)*stride], 0.0); + sum += m_y[k]; + } + m_y[0] = 1.0 - sum; + + transform(m_y.begin(), m_y.end(), m_rmolwts.begin(), + m_ym.begin(), multiplies()); + m_mmw = 1.0 / accumulate(m_ym.begin(), m_ym.end(), 0.0); + compositionChanged(); +} + void Phase::setMassFractionsByName(const Composition& yMap) { vector mf = getCompositionFromMap(yMap); diff --git a/src/thermo/ThermoFactory.cpp b/src/thermo/ThermoFactory.cpp index c3d76bc4e6..06aec34915 100644 --- a/src/thermo/ThermoFactory.cpp +++ b/src/thermo/ThermoFactory.cpp @@ -15,6 +15,8 @@ #include "cantera/thermo/PDSSFactory.h" #include "cantera/thermo/MultiSpeciesThermo.h" #include "cantera/thermo/IdealGasPhase.h" +#include "cantera/thermo/IdealGasSurfPhase.h" +#include "cantera/thermo/IdealMoistPhase.h" #include "cantera/thermo/PlasmaPhase.h" #include "cantera/thermo/IdealSolidSolnPhase.h" @@ -53,6 +55,8 @@ ThermoFactory::ThermoFactory() addDeprecatedAlias("none", "None"); reg("ideal-gas", []() { return new IdealGasPhase(); }); addDeprecatedAlias("ideal-gas", "IdealGas"); + reg("ideal-gas-surface", []() { return new IdealGasSurfPhase(); }); + reg("ideal-moist", []() { return new IdealMoistPhase(); }); reg("plasma", []() { return new PlasmaPhase(); }); reg("ideal-surface", []() { return new SurfPhase(); }); addDeprecatedAlias("ideal-surface", "Surface"); diff --git a/z.ab_patch b/z.ab_patch new file mode 100644 index 0000000000..3d8999773e --- /dev/null +++ b/z.ab_patch @@ -0,0 +1,1043 @@ +diff --git a/INSTALL.sh b/INSTALL.sh +new file mode 100644 +index 000000000..0e158dff6 +--- /dev/null ++++ b/INSTALL.sh +@@ -0,0 +1,20 @@ ++CXX=g++ ++CC=gcc ++cxx_flags="-std=c++17" ++prefix=${HOME}/opt/ ++python_package=n ++f90_interface=n ++system_eigen=n ++extra_inc_dirs="/opt/include /usr/include/yaml-cpp" # Include both directories ++system_blas_lapack=n ++boost_inc_dir=$(pwd)/ext/cliboost ++system_sundials=n ++system_yamlcpp=n ++ ++scons build CXX=${CXX} CC=${CC} cxx_flags="${cxx_flags}" prefix=${prefix} \ ++ python_package=${python_package} f90_interface=${f90_interface} \ ++ system_eigen=${system_eigen} extra_inc_dirs="${extra_inc_dirs}" \ ++ system_blas_lapack=${system_blas_lapack} boost_inc_dir=${boost_inc_dir} system_sundials=${system_sundials} \ ++ system_yamlcpp=${system_yamlcpp} -j8 ++ ++scons install +diff --git a/include/cantera/kinetics/Arrhenius.h b/include/cantera/kinetics/Arrhenius.h +index 2791faca8..473cfb262 100644 +--- a/include/cantera/kinetics/Arrhenius.h ++++ b/include/cantera/kinetics/Arrhenius.h +@@ -1,212 +1,215 @@ +-/** +- * @file Arrhenius.h +- * Header for reaction rates that involve Arrhenius-type kinetics. +- */ +- +-// This file is part of Cantera. See License.txt in the top-level directory or +-// at https://cantera.org/license.txt for license and copyright information. +- +-#ifndef CT_ARRHENIUS_H +-#define CT_ARRHENIUS_H +- +-#include "cantera/base/ct_defs.h" +-#include "cantera/base/Units.h" +-#include "cantera/kinetics/ReactionData.h" +-#include "ReactionRate.h" +-#include "MultiRate.h" +- +-namespace Cantera +-{ +- +-class AnyValue; +-class AnyMap; +- +-//! Data container holding shared data specific to ArrheniusRate +-/** +- * The data container `ArrheniusData` holds precalculated data common to +- * all `ArrheniusRate` objects. +- */ +-struct ArrheniusData : public ReactionData +-{ +- bool update(const ThermoPhase& phase, const Kinetics& kin) override; +- using ReactionData::update; +-}; +- +- +-//! Base class for Arrhenius-type Parameterizations +-/*! +- * This base class provides a minimally functional interface that allows for parameter +- * access from derived classes as well as classes that use Arrhenius-type expressions +- * internally, for example FalloffRate and PlogRate. +- * @ingroup arrheniusGroup +- */ +-class ArrheniusBase : public ReactionRate +-{ +-public: +- //! Default constructor. +- ArrheniusBase() {} +- +- //! Constructor. +- /*! +- * @param A Pre-exponential factor. The unit system is (kmol, m, s); actual units +- * depend on the reaction order and the dimensionality (surface or bulk). +- * @param b Temperature exponent (non-dimensional) +- * @param Ea Activation energy in energy units [J/kmol] +- */ +- ArrheniusBase(double A, double b, double Ea); +- +- //! Constructor based on AnyValue content +- ArrheniusBase(const AnyValue& rate, const UnitSystem& units, +- const UnitStack& rate_units); +- +- explicit ArrheniusBase(const AnyMap& node, const UnitStack& rate_units={}); +- +- //! Perform object setup based on AnyValue node information +- /*! +- * Used to set parameters from a child of the reaction node, which may have +- * different names for different rate parameterizations, such as falloff rates. +- * +- * @param rate Child of the reaction node containing Arrhenius rate parameters. +- * For example, the `rate-coefficient` node for a standard Arrhenius reaction. +- * @param units Unit system +- * @param rate_units Unit definitions specific to rate information +- */ +- void setRateParameters(const AnyValue& rate, +- const UnitSystem& units, +- const UnitStack& rate_units); +- +- //! Get Arrhenius parameters used to populate the `rate-coefficient` or +- //! equivalent field +- void getRateParameters(AnyMap& node) const; +- +- void setParameters(const AnyMap& node, const UnitStack& rate_units) override; +- +- void getParameters(AnyMap& node) const override; +- +- //! Check rate expression +- void check(const string& equation) override; +- +- void validate(const string& equation, const Kinetics& kin) override; +- +- //! Return the pre-exponential factor *A* (in m, kmol, s to powers depending +- //! on the reaction order) +- /*! +- * Class specializations may provide alternate definitions that describe +- * an effective pre-exponential factor that depends on the thermodynamic state. +- */ +- virtual double preExponentialFactor() const { +- return m_A; +- } +- +- //! Return the temperature exponent *b* +- /*! +- * Class specializations may provide alternate definitions that describe +- * an effective temperature exponent that depends on the thermodynamic state. +- */ +- virtual double temperatureExponent() const { +- return m_b; +- } +- +- //! Return the activation energy *Ea* [J/kmol] +- //! The value corresponds to the constant specified by input parameters; +- /*! +- * Class specializations may provide alternate definitions that describe +- * an effective activation energy that depends on the thermodynamic state. +- */ +- virtual double activationEnergy() const { +- return m_Ea_R * GasConstant; +- } +- +- //! Return reaction order associated with the reaction rate +- double order() const { +- return m_order; +- } +- +- //! Set units of the reaction rate expression +- void setRateUnits(const UnitStack& rate_units) override { +- ReactionRate::setRateUnits(rate_units); +- if (rate_units.size() > 1) { +- m_order = 1 - rate_units.product().dimension("quantity"); +- } else { +- m_order = NAN; +- } +- } +- +- //! Get flag indicating whether negative A values are permitted +- bool allowNegativePreExponentialFactor() const { +- return m_negativeA_ok; +- } +- +- //! Set flag indicating whether negative A values are permitted +- void setAllowNegativePreExponentialFactor(bool value) { +- m_negativeA_ok = value; +- } +- +-protected: +- bool m_negativeA_ok = false; //!< Permissible negative A values +- double m_A = NAN; //!< Pre-exponential factor +- double m_b = NAN; //!< Temperature exponent +- double m_Ea_R = 0.; //!< Activation energy (in temperature units) +- double m_E4_R = 0.; //!< Optional 4th energy parameter (in temperature units) +- double m_logA = NAN; //!< Logarithm of pre-exponential factor +- double m_order = NAN; //!< Reaction order +- string m_A_str = "A"; //!< The string for the pre-exponential factor +- string m_b_str = "b"; //!< The string for temperature exponent +- string m_Ea_str = "Ea"; //!< The string for activation energy +- string m_E4_str = ""; //!< The string for an optional 4th parameter +-}; +- +-//! Arrhenius reaction rate type depends only on temperature +-/*! +- * A reaction rate coefficient of the following form. +- * +- * @f[ +- * k_f = A T^b \exp (-Ea/RT) +- * @f] +- * +- * @ingroup arrheniusGroup +- */ +-class ArrheniusRate : public ArrheniusBase +-{ +-public: +- using ArrheniusBase::ArrheniusBase; // inherit constructors +- +- unique_ptr newMultiRate() const override { +- return make_unique>(); +- } +- +- const string type() const override { +- return "Arrhenius"; +- } +- +- //! Evaluate reaction rate +- double evalRate(double logT, double recipT) const { +- return m_A * std::exp(m_b * logT - m_Ea_R * recipT); +- } +- +- //! Evaluate natural logarithm of the rate constant. +- double evalLog(double logT, double recipT) const { +- return m_logA + m_b * logT - m_Ea_R * recipT; +- } +- +- //! Evaluate reaction rate +- /*! +- * @param shared_data data shared by all reactions of a given type +- */ +- double evalFromStruct(const ArrheniusData& shared_data) const { +- return m_A * std::exp(m_b * shared_data.logT - m_Ea_R * shared_data.recipT); +- } +- +- //! Evaluate derivative of reaction rate with respect to temperature +- //! divided by reaction rate +- /*! +- * @param shared_data data shared by all reactions of a given type +- */ +- double ddTScaledFromStruct(const ArrheniusData& shared_data) const { +- return (m_Ea_R * shared_data.recipT + m_b) * shared_data.recipT; +- } +-}; +- +-} +- +-#endif ++/** ++ * @file Arrhenius.h ++ * Header for reaction rates that involve Arrhenius-type kinetics. ++ */ ++ ++// This file is part of Cantera. See License.txt in the top-level directory or ++// at https://cantera.org/license.txt for license and copyright information. ++ ++#ifndef CT_ARRHENIUS_H ++#define CT_ARRHENIUS_H ++ ++#include "cantera/base/ct_defs.h" ++#include "cantera/base/Units.h" ++#include "cantera/kinetics/ReactionData.h" ++#include "ReactionRate.h" ++#include "MultiRate.h" ++ ++namespace Cantera ++{ ++ ++class AnyValue; ++class AnyMap; ++ ++//! Data container holding shared data specific to ArrheniusRate ++/** ++ * The data container `ArrheniusData` holds precalculated data common to ++ * all `ArrheniusRate` objects. ++ */ ++struct ArrheniusData : public ReactionData ++{ ++ bool update(const ThermoPhase& phase, const Kinetics& kin) override; ++ using ReactionData::update; ++}; ++ ++ ++//! Base class for Arrhenius-type Parameterizations ++/*! ++ * This base class provides a minimally functional interface that allows for parameter ++ * access from derived classes as well as classes that use Arrhenius-type expressions ++ * internally, for example FalloffRate and PlogRate. ++ * @ingroup arrheniusGroup ++ */ ++class ArrheniusBase : public ReactionRate ++{ ++public: ++ //! Default constructor. ++ ArrheniusBase() {} ++ ++ //! Constructor. ++ /*! ++ * @param A Pre-exponential factor. The unit system is (kmol, m, s); actual units ++ * depend on the reaction order and the dimensionality (surface or bulk). ++ * @param T0 Reference temperature (K) ++ * @param b Temperature exponent (non-dimensional) ++ * @param Ea Activation energy in energy units [J/kmol] ++ */ ++ ArrheniusBase(double A, double b,double T0,double Ea); ++ ++ //! Constructor based on AnyValue content ++ ArrheniusBase(const AnyValue& rate, const UnitSystem& units, ++ const UnitStack& rate_units); ++ ++ explicit ArrheniusBase(const AnyMap& node, const UnitStack& rate_units={}); ++ ++ //! Perform object setup based on AnyValue node information ++ /*! ++ * Used to set parameters from a child of the reaction node, which may have ++ * different names for different rate parameterizations, such as falloff rates. ++ * ++ * @param rate Child of the reaction node containing Arrhenius rate parameters. ++ * For example, the `rate-coefficient` node for a standard Arrhenius reaction. ++ * @param units Unit system ++ * @param rate_units Unit definitions specific to rate information ++ */ ++ void setRateParameters(const AnyValue& rate, ++ const UnitSystem& units, ++ const UnitStack& rate_units); ++ ++ //! Get Arrhenius parameters used to populate the `rate-coefficient` or ++ //! equivalent field ++ void getRateParameters(AnyMap& node) const; ++ ++ void setParameters(const AnyMap& node, const UnitStack& rate_units) override; ++ ++ void getParameters(AnyMap& node) const override; ++ ++ //! Check rate expression ++ void check(const string& equation) override; ++ ++ void validate(const string& equation, const Kinetics& kin) override; ++ ++ //! Return the pre-exponential factor *A* (in m, kmol, s to powers depending ++ //! on the reaction order) ++ /*! ++ * Class specializations may provide alternate definitions that describe ++ * an effective pre-exponential factor that depends on the thermodynamic state. ++ */ ++ virtual double preExponentialFactor() const { ++ return m_A; ++ } ++ ++ //! Return the temperature exponent *b* ++ /*! ++ * Class specializations may provide alternate definitions that describe ++ * an effective temperature exponent that depends on the thermodynamic state. ++ */ ++ virtual double temperatureExponent() const { ++ return m_b; ++ } ++ ++ //! Return the activation energy *Ea* [J/kmol] ++ //! The value corresponds to the constant specified by input parameters; ++ /*! ++ * Class specializations may provide alternate definitions that describe ++ * an effective activation energy that depends on the thermodynamic state. ++ */ ++ virtual double activationEnergy() const { ++ return m_Ea_R * GasConstant; ++ } ++ ++ //! Return reaction order associated with the reaction rate ++ double order() const { ++ return m_order; ++ } ++ ++ //! Set units of the reaction rate expression ++ void setRateUnits(const UnitStack& rate_units) override { ++ ReactionRate::setRateUnits(rate_units); ++ if (rate_units.size() > 1) { ++ m_order = 1 - rate_units.product().dimension("quantity"); ++ } else { ++ m_order = NAN; ++ } ++ } ++ ++ //! Get flag indicating whether negative A values are permitted ++ bool allowNegativePreExponentialFactor() const { ++ return m_negativeA_ok; ++ } ++ ++ //! Set flag indicating whether negative A values are permitted ++ void setAllowNegativePreExponentialFactor(bool value) { ++ m_negativeA_ok = value; ++ } ++ ++protected: ++ bool m_negativeA_ok = false; //!< Permissible negative A values ++ double m_A = NAN; //!< Pre-exponential factor ++ double m_b = NAN; //!< Temperature exponent ++ double m_T0 = 1.0; //!< Reference temperature ++ double m_Ea_R = 0.; //!< Activation energy (in temperature units) ++ double m_E4_R = 0.; //!< Optional 4th energy parameter (in temperature units) ++ double m_logA = NAN; //!< Logarithm of pre-exponential factor ++ double m_order = NAN; //!< Reaction order ++ string m_A_str = "A"; //!< The string for the pre-exponential factor ++ string m_b_str = "b"; //!< The string for temperature exponent ++ string m_T0_str = "T0"; //!< The string for reference temperature ++ string m_Ea_str = "Ea"; //!< The string for activation energy ++ string m_E4_str = ""; //!< The string for an optional 4th parameter ++}; ++ ++//! Arrhenius reaction rate type depends only on temperature ++/*! ++ * A reaction rate coefficient of the following form. ++ * ++ * @f[ ++ * k_f = A (T/T0)^b \exp (-Ea/RT) ++ * @f] ++ * ++ * @ingroup arrheniusGroup ++ */ ++class ArrheniusRate : public ArrheniusBase ++{ ++public: ++ using ArrheniusBase::ArrheniusBase; // inherit constructors ++ ++ unique_ptr newMultiRate() const override { ++ return make_unique>(); ++ } ++ ++ const string type() const override { ++ return "T-Arrhenius"; ++ } ++ ++ //! Evaluate reaction rate ++ double evalRate(double logT, double recipT) const { ++ return m_A * std::exp((m_b * (logT - std::log(m_T0))) - (m_Ea_R * recipT)); ++ } ++ ++ //! Evaluate natural logarithm of the rate constant. ++ double evalLog(double logT, double recipT) const { ++ return m_logA + (m_b * (logT - std::log(m_T0))) - (m_Ea_R * recipT); ++ } ++ ++ //! Evaluate reaction rate ++ /*! ++ * @param shared_data data shared by all reactions of a given type ++ */ ++ double evalFromStruct(const ArrheniusData& shared_data) const { ++ return m_A * std::exp((m_b * (shared_data.logT - std::log(m_T0))) - (m_Ea_R * shared_data.recipT)); ++ } ++ ++ //! Evaluate derivative of reaction rate with respect to temperature ++ //! divided by reaction rate ++ /*! ++ * @param shared_data data shared by all reactions of a given type ++ */ ++ double ddTScaledFromStruct(const ArrheniusData& shared_data) const { ++ return (m_Ea_R * shared_data.recipT + m_b) * shared_data.recipT; ++ } ++}; ++ ++} ++ ++#endif +diff --git a/include/cantera/kinetics/BlowersMaselRate.h b/include/cantera/kinetics/BlowersMaselRate.h +index d10db0b72..cbef45a7a 100644 +--- a/include/cantera/kinetics/BlowersMaselRate.h ++++ b/include/cantera/kinetics/BlowersMaselRate.h +@@ -81,7 +81,7 @@ public: + * @param w Average bond dissociation energy of the bond being formed and + * broken in the reaction, in energy units [J/kmol] + */ +- BlowersMaselRate(double A, double b, double Ea0, double w); ++ BlowersMaselRate(double A, double b, double T0, double Ea0, double w); + + explicit BlowersMaselRate(const AnyMap& node, + const UnitStack& rate_units={}); +diff --git a/include/cantera/kinetics/Photolysis.h b/include/cantera/kinetics/Photolysis.h +index ee9b4c0cf..358e9f5b3 100644 +--- a/include/cantera/kinetics/Photolysis.h ++++ b/include/cantera/kinetics/Photolysis.h +@@ -1,4 +1,6 @@ +-//! @file Photolysis.h ++/** @file Photolysis.h ++ * Header for reaction rates that involve Photochemical reactions ++ */ + + #ifndef CT_PHOTOLYSIS_H + #define CT_PHOTOLYSIS_H +@@ -69,12 +71,16 @@ class PhotolysisBase : public ReactionRate { + //! @param branch_map Map of branch names to branch indices + void setRateParameters(const AnyValue& rate, map const& branch_map); + ++ //! Get the parameters corresponding to node rate-constant + void getParameters(AnyMap& node) const override; + ++ //! Get the parameters for a given node with flow style output + void getRateParameters(AnyMap& node) const; + ++ //! Checks for temperature range, and wavelength data + void check(string const& equation) override; + ++ //! Checks for valid species, stoichiometric balance, and consistency with photolysis branches + void validate(const string& equation, const Kinetics& kin) override; + + vector getCrossSection(double temp, double wavelength) const; +@@ -122,14 +128,21 @@ class PhotolysisRate : public PhotolysisBase { + return make_unique>(); + } + ++ //! reaction string type for photolysis reactions + const string type() const override { + return "Photolysis"; + } +- ++ ++ //! net stoichiometric coefficients of photolysis products + Composition const& photoProducts() const override { + return m_net_products; + } + ++/** ++ * @brief Calculates the photolysis reaction rate and updates stoichiometric concentration of products ++ * ++ * @return total photolysis rate from all the branches ++ */ + double evalFromStruct(PhotolysisData const& data); + + protected: +diff --git a/include/cantera/kinetics/TwoTempPlasmaRate.h b/include/cantera/kinetics/TwoTempPlasmaRate.h +index 7f17e9de7..f09f11b55 100644 +--- a/include/cantera/kinetics/TwoTempPlasmaRate.h ++++ b/include/cantera/kinetics/TwoTempPlasmaRate.h +@@ -67,10 +67,11 @@ public: + * @param A Pre-exponential factor. The unit system is (kmol, m, s); actual units + * depend on the reaction order and the dimensionality (surface or bulk). + * @param b Temperature exponent (non-dimensional) ++ * @param T0 Reference temperature (K) + * @param Ea Activation energy in energy units [J/kmol] + * @param EE Activation electron energy in energy units [J/kmol] + */ +- TwoTempPlasmaRate(double A, double b, double Ea=0.0, double EE=0.0); ++ TwoTempPlasmaRate(double A, double b, double T0, double Ea=0.0, double EE=0.0); + + TwoTempPlasmaRate(const AnyMap& node, const UnitStack& rate_units={}); + +diff --git a/src/kinetics/Arrhenius.cpp b/src/kinetics/Arrhenius.cpp +index ecd89b203..42f82d618 100644 +--- a/src/kinetics/Arrhenius.cpp ++++ b/src/kinetics/Arrhenius.cpp +@@ -1,156 +1,160 @@ +-//! @file Arrhenius.cpp +- +-// This file is part of Cantera. See License.txt in the top-level directory or +-// at https://cantera.org/license.txt for license and copyright information. +- +-#include "cantera/kinetics/Arrhenius.h" +-#include "cantera/thermo/ThermoPhase.h" +- +-namespace Cantera +-{ +- +-ArrheniusBase::ArrheniusBase(double A, double b, double Ea) +- : m_A(A) +- , m_b(b) +- , m_Ea_R(Ea / GasConstant) +-{ +- if (m_A > 0.0) { +- m_logA = std::log(m_A); +- } +- m_valid = true; +-} +- +-ArrheniusBase::ArrheniusBase(const AnyValue& rate, const UnitSystem& units, +- const UnitStack& rate_units) +-{ +- setRateUnits(rate_units); +- setRateParameters(rate, units, rate_units); +-} +- +-ArrheniusBase::ArrheniusBase(const AnyMap& node, const UnitStack& rate_units) +-{ +- setParameters(node, rate_units); +-} +- +-void ArrheniusBase::setRateParameters( +- const AnyValue& rate, const UnitSystem& units, const UnitStack& rate_units) +-{ +- m_Ea_R = 0.; // assume zero if not provided +- m_E4_R = 0.; // assume zero if not provided +- if (rate.empty()) { +- m_A = NAN; +- m_b = NAN; +- m_logA = NAN; +- setRateUnits(Units(0.)); +- return; +- } +- +- if (rate.is()) { +- +- auto& rate_map = rate.as(); +- m_A = units.convertRateCoeff(rate_map[m_A_str], conversionUnits()); +- m_b = rate_map[m_b_str].asDouble(); +- if (rate_map.hasKey(m_Ea_str)) { +- m_Ea_R = units.convertActivationEnergy(rate_map[m_Ea_str], "K"); +- } +- if (rate_map.hasKey(m_E4_str)) { +- m_E4_R = units.convertActivationEnergy(rate_map[m_E4_str], "K"); +- } +- } else { +- auto& rate_vec = rate.asVector(2, 4); +- m_A = units.convertRateCoeff(rate_vec[0], conversionUnits()); +- m_b = rate_vec[1].asDouble(); +- if (rate_vec.size() > 2) { +- m_Ea_R = units.convertActivationEnergy(rate_vec[2], "K"); +- } +- if (rate_vec.size() > 3) { +- m_E4_R = units.convertActivationEnergy(rate_vec[3], "K"); +- } +- } +- if (m_A > 0.0) { +- m_logA = std::log(m_A); +- } +- m_valid = true; +-} +- +-void ArrheniusBase::getRateParameters(AnyMap& node) const +-{ +- if (!valid()) { +- // Return empty/unmodified AnyMap +- return; +- } +- +- if (conversionUnits().factor() != 0.0) { +- node[m_A_str].setQuantity(m_A, conversionUnits()); +- } else { +- node[m_A_str] = m_A; +- // This can't be converted to a different unit system because the dimensions of +- // the rate constant were not set. Can occur if the reaction was created outside +- // the context of a Kinetics object and never added to a Kinetics object. +- node["__unconvertible__"] = true; +- } +- node[m_b_str] = m_b; +- node[m_Ea_str].setQuantity(m_Ea_R, "K", true); +- if (m_E4_str != "") { +- node[m_E4_str].setQuantity(m_E4_R, "K", true); +- } +- node.setFlowStyle(); +-} +- +-void ArrheniusBase::setParameters(const AnyMap& node, const UnitStack& rate_units) +-{ +- ReactionRate::setParameters(node, rate_units); +- m_negativeA_ok = node.getBool("negative-A", false); +- if (!node.hasKey("rate-constant")) { +- setRateParameters(AnyValue(), node.units(), rate_units); +- return; +- } +- setRateParameters(node["rate-constant"], node.units(), rate_units); +-} +- +-void ArrheniusBase::getParameters(AnyMap& node) const { +- if (m_negativeA_ok) { +- node["negative-A"] = true; +- } +- AnyMap rateNode; +- getRateParameters(rateNode); +- if (!rateNode.empty()) { +- // RateType object is configured +- node["rate-constant"] = std::move(rateNode); +- } +-} +- +-void ArrheniusBase::check(const string& equation) +-{ +- if (!m_negativeA_ok && m_A < 0) { +- if (equation == "") { +- throw CanteraError("ArrheniusBase::check", +- "Detected negative pre-exponential factor (A={}).\n" +- "Enable 'allowNegativePreExponentialFactor' to suppress " +- "this message.", m_A); +- } +- throw InputFileError("ArrheniusBase::check", m_input, +- "Undeclared negative pre-exponential factor found in reaction '{}'", +- equation); +- } +-} +- +-void ArrheniusBase::validate(const string& equation, const Kinetics& kin) +-{ +- if (!valid()) { +- throw InputFileError("ArrheniusBase::validate", m_input, +- "Rate object for reaction '{}' is not configured.", equation); +- } +-} +- +-bool ArrheniusData::update(const ThermoPhase& phase, const Kinetics& kin) +-{ +- double T = phase.temperature(); +- if (T == temperature) { +- return false; +- } +- update(T); +- return true; +-} +- +-} ++//! @file Arrhenius.cpp ++ ++// This file is part of Cantera. See License.txt in the top-level directory or ++// at https://cantera.org/license.txt for license and copyright information. ++ ++#include "cantera/kinetics/Arrhenius.h" ++#include "cantera/thermo/ThermoPhase.h" ++ ++namespace Cantera ++{ ++ ++ArrheniusBase::ArrheniusBase(double A, double b, double T0, double Ea) ++ : m_A(A) ++ , m_b(b) ++ , m_T0(T0) ++ , m_Ea_R(Ea / GasConstant) ++{ ++ if (m_A > 0.0) { ++ m_logA = std::log(m_A); ++ } ++ m_valid = true; ++} ++ ++ArrheniusBase::ArrheniusBase(const AnyValue& rate, const UnitSystem& units, ++ const UnitStack& rate_units) ++{ ++ setRateUnits(rate_units); ++ setRateParameters(rate, units, rate_units); ++} ++ ++ArrheniusBase::ArrheniusBase(const AnyMap& node, const UnitStack& rate_units) ++{ ++ setParameters(node, rate_units); ++} ++ ++void ArrheniusBase::setRateParameters( ++ const AnyValue& rate, const UnitSystem& units, const UnitStack& rate_units) ++{ ++ m_Ea_R = 0.; // assume zero if not provided ++ m_E4_R = 0.; // assume zero if not provided ++ if (rate.empty()) { ++ m_A = NAN; ++ m_b = NAN; ++ m_T0 = 1.0; ++ m_logA = NAN; ++ setRateUnits(Units(0.)); ++ return; ++ } ++ ++ if (rate.is()) { ++ ++ auto& rate_map = rate.as(); ++ m_A = units.convertRateCoeff(rate_map[m_A_str], conversionUnits()); ++ m_b = rate_map[m_b_str].asDouble(); ++ m_T0 = rate_map[m_T0_str].asDouble(); ++ if (rate_map.hasKey(m_Ea_str)) { ++ m_Ea_R = units.convertActivationEnergy(rate_map[m_Ea_str], "K"); ++ } ++ if (rate_map.hasKey(m_E4_str)) { ++ m_E4_R = units.convertActivationEnergy(rate_map[m_E4_str], "K"); ++ } ++ } else { ++ auto& rate_vec = rate.asVector(2, 4); ++ m_A = units.convertRateCoeff(rate_vec[0], conversionUnits()); ++ m_b = rate_vec[1].asDouble(); ++ if (rate_vec.size() > 2) { ++ m_Ea_R = units.convertActivationEnergy(rate_vec[2], "K"); ++ } ++ if (rate_vec.size() > 3) { ++ m_E4_R = units.convertActivationEnergy(rate_vec[3], "K"); ++ } ++ } ++ if (m_A > 0.0) { ++ m_logA = std::log(m_A); ++ } ++ m_valid = true; ++} ++ ++void ArrheniusBase::getRateParameters(AnyMap& node) const ++{ ++ if (!valid()) { ++ // Return empty/unmodified AnyMap ++ return; ++ } ++ ++ if (conversionUnits().factor() != 0.0) { ++ node[m_A_str].setQuantity(m_A, conversionUnits()); ++ } else { ++ node[m_A_str] = m_A; ++ // This can't be converted to a different unit system because the dimensions of ++ // the rate constant were not set. Can occur if the reaction was created outside ++ // the context of a Kinetics object and never added to a Kinetics object. ++ node["__unconvertible__"] = true; ++ } ++ node[m_b_str] = m_b; ++ node[m_T0_str] = m_T0; ++ node[m_Ea_str].setQuantity(m_Ea_R, "K", true); ++ if (m_E4_str != "") { ++ node[m_E4_str].setQuantity(m_E4_R, "K", true); ++ } ++ node.setFlowStyle(); ++} ++ ++void ArrheniusBase::setParameters(const AnyMap& node, const UnitStack& rate_units) ++{ ++ ReactionRate::setParameters(node, rate_units); ++ m_negativeA_ok = node.getBool("negative-A", false); ++ if (!node.hasKey("rate-constant")) { ++ setRateParameters(AnyValue(), node.units(), rate_units); ++ return; ++ } ++ setRateParameters(node["rate-constant"], node.units(), rate_units); ++} ++ ++void ArrheniusBase::getParameters(AnyMap& node) const { ++ if (m_negativeA_ok) { ++ node["negative-A"] = true; ++ } ++ AnyMap rateNode; ++ getRateParameters(rateNode); ++ if (!rateNode.empty()) { ++ // RateType object is configured ++ node["rate-constant"] = std::move(rateNode); ++ } ++} ++ ++void ArrheniusBase::check(const string& equation) ++{ ++ if (!m_negativeA_ok && m_A < 0) { ++ if (equation == "") { ++ throw CanteraError("ArrheniusBase::check", ++ "Detected negative pre-exponential factor (A={}).\n" ++ "Enable 'allowNegativePreExponentialFactor' to suppress " ++ "this message.", m_A); ++ } ++ throw InputFileError("ArrheniusBase::check", m_input, ++ "Undeclared negative pre-exponential factor found in reaction '{}'", ++ equation); ++ } ++} ++ ++void ArrheniusBase::validate(const string& equation, const Kinetics& kin) ++{ ++ if (!valid()) { ++ throw InputFileError("ArrheniusBase::validate", m_input, ++ "Rate object for reaction '{}' is not configured.", equation); ++ } ++} ++ ++bool ArrheniusData::update(const ThermoPhase& phase, const Kinetics& kin) ++{ ++ double T = phase.temperature(); ++ if (T == temperature) { ++ return false; ++ } ++ update(T); ++ return true; ++} ++ ++} +diff --git a/src/kinetics/BlowersMaselRate.cpp b/src/kinetics/BlowersMaselRate.cpp +index c51d8d718..528d17a13 100644 +--- a/src/kinetics/BlowersMaselRate.cpp ++++ b/src/kinetics/BlowersMaselRate.cpp +@@ -42,8 +42,8 @@ BlowersMaselRate::BlowersMaselRate() + m_E4_str = "w"; + } + +-BlowersMaselRate::BlowersMaselRate(double A, double b, double Ea0, double w) +- : ArrheniusBase(A, b, Ea0) ++BlowersMaselRate::BlowersMaselRate(double A, double b, double T0, double Ea0, double w) ++ : ArrheniusBase(A, b, T0, Ea0) + { + m_Ea_str = "Ea0"; + m_E4_str = "w"; +diff --git a/src/kinetics/Photolysis.cpp b/src/kinetics/Photolysis.cpp +index 67d2ccb87..a87a8802d 100644 +--- a/src/kinetics/Photolysis.cpp ++++ b/src/kinetics/Photolysis.cpp +@@ -44,11 +44,13 @@ bool PhotolysisData::check() const + "Wavelength grid must have at least two points."); + } + ++ // Check that wavelength grid values are positive + if (wavelength[0] <= 0.0) { + throw CanteraError("PhotolysisData::update", + "Wavelength grid must be positive."); + } + ++ // Check that wavelength grid values are monotonic and increasing + for (size_t i = 1; i < wavelength.size(); i++) { + if (wavelength[i] <= wavelength[i-1]) { + throw CanteraError("PhotolysisData::update", +@@ -61,12 +63,14 @@ bool PhotolysisData::check() const + throw CanteraError("PhotolysisData::update", + "Actinic flux is empty."); + } +- ++ ++ // Check that actinic flux grid should have same size as wavelength grid + if (actinicFlux.size() != wavelength.size()) { + throw CanteraError("PhotolysisData::update", + "Actinic flux must have the same size as the wavelength grid."); + } + ++ // Check that actinic flux values are positive + for (size_t i = 0; i < actinicFlux.size(); i++) { + if (actinicFlux[i] < 0.0) { + throw CanteraError("PhotolysisData::update", +@@ -87,6 +91,7 @@ PhotolysisBase::PhotolysisBase( + m_ntemp = temp.size(); + m_nwave = wavelength.size(); + ++ // Grid for temperature and wavelength + m_temp_wave_grid.resize(m_ntemp + m_nwave); + for (size_t i = 0; i < m_ntemp; i++) { + m_temp_wave_grid[i] = temp[i]; +@@ -100,6 +105,7 @@ PhotolysisBase::PhotolysisBase( + m_branch.push_back(parseCompString(branch)); + } + ++ // Check if cross-section data size + if (m_ntemp * m_nwave * branches.size() != m_crossSection.size()) { + throw CanteraError("PhotolysisBase::PhotolysisBase", + "Cross-section data size does not match the temperature, " +@@ -182,16 +188,19 @@ void PhotolysisBase::setParameters(AnyMap const& node, UnitStack const& rate_uni + } else if (rtmp.products != rtmp.reactants) { // this is not photoabsorption + m_branch.push_back(rtmp.products); + } +- ++ + if (node.hasKey("cross-section")) { + for (auto const& data: node["cross-section"].asVector()) { + auto format = data["format"].asString(); + auto temp = data["temperature-range"].asVector(2, 2); ++ ++ //Check temperature range to be monotonically increasing + if (temp[0] >= temp[1]) { + throw CanteraError("PhotolysisBase::setParameters", + "Temperature range must be strictly increasing."); + } +- ++ ++ //Check for gaps in temperature range + if (temperature.empty()) { + temperature = temp; + } else { +@@ -209,9 +218,12 @@ void PhotolysisBase::setParameters(AnyMap const& node, UnitStack const& rate_uni + result.first.push_back(entry[0]); + result.second.push_back(entry[1]); + } ++ ++ //Read file names for VULCAN photochemistry database + } else if (format == "VULCAN") { + auto files = data["filenames"].asVector(); + result = load_xsection_vulcan(files, m_branch); ++ //Read file name for KINETICS photochemistry database + } else if (format == "KINETICS7") { + auto files = data["filenames"].asVector(); + result = load_xsection_kinetics7(files, m_branch); +@@ -267,11 +279,13 @@ void PhotolysisBase::setParameters(AnyMap const& node, UnitStack const& rate_uni + m_valid = true; + } + ++ //Set the rate parameters to flow style ?? (TBD) + void PhotolysisBase::getRateParameters(AnyMap& node) const + { + node.setFlowStyle(); + } + ++ //Get rate parameters for photolysis reaction + void PhotolysisBase::getParameters(AnyMap& node) const + { + AnyMap rateNode; +@@ -282,6 +296,7 @@ void PhotolysisBase::getParameters(AnyMap& node) const + } + } + ++//Check temperature range, and wavelength data for photolysis reactions + void PhotolysisBase::check(string const& equation) + { + if (m_ntemp < 2) { +@@ -296,6 +311,7 @@ void PhotolysisBase::check(string const& equation) + } + } + ++//Check rate coefficient + void PhotolysisBase::validate(string const& equation, Kinetics const& kin) + { + if (!valid()) { +@@ -330,6 +346,7 @@ void PhotolysisBase::validate(string const& equation, Kinetics const& kin) + } + } + ++ // Check for consistency in species in reaction string, and photolysis branches + if (species_from_equation != species_from_branches) { + throw InputFileError("PhotolysisBase::validate", m_input, + "Reaction '{}' has different products than the photolysis branches.", equation); +@@ -353,6 +370,7 @@ vector PhotolysisBase::getCrossSection(double temp, double wavelength) c + return cross; + } + ++// Evaluate the photolysis rate, and update the stoichiometric coefficient for different branches + double PhotolysisRate::evalFromStruct(PhotolysisData const& data) { + double wmin = m_temp_wave_grid[m_ntemp]; + double wmax = m_temp_wave_grid.back(); +@@ -382,10 +400,9 @@ double PhotolysisRate::evalFromStruct(PhotolysisData const& data) { + + double coord[2] = {data.temperature, data.wavelength[0]}; + size_t len[2] = {m_ntemp, m_nwave}; ++ + +- // debug +- //std::cout << "coord = " << coord[0] << " " << coord[1] << std::endl; +- ++ // N-space interpolation to determine photolysis cross section + interpn(cross1, coord, m_crossSection.data(), m_temp_wave_grid.data(), + len, 2, m_branch.size()); + +@@ -417,6 +434,7 @@ double PhotolysisRate::evalFromStruct(PhotolysisData const& data) { + + // photodissociation only + for (size_t n = 1; n < m_branch.size(); n++) { ++ //Photochemical rate constant + double rate = 0.5 * (data.wavelength[i+1] - data.wavelength[i]) + * (cross1[n] * data.actinicFlux[i] + cross2[n] * data.actinicFlux[i+1]); + +diff --git a/src/kinetics/TwoTempPlasmaRate.cpp b/src/kinetics/TwoTempPlasmaRate.cpp +index ce3f46949..1a930ada5 100644 +--- a/src/kinetics/TwoTempPlasmaRate.cpp ++++ b/src/kinetics/TwoTempPlasmaRate.cpp +@@ -51,8 +51,8 @@ TwoTempPlasmaRate::TwoTempPlasmaRate() + m_E4_str = "Ea-electron"; + } + +-TwoTempPlasmaRate::TwoTempPlasmaRate(double A, double b, double Ea, double EE) +- : ArrheniusBase(A, b, Ea) ++TwoTempPlasmaRate::TwoTempPlasmaRate(double A, double b, double T0, double Ea, double EE) ++ : ArrheniusBase(A, b, T0, Ea) + { + m_Ea_str = "Ea-gas"; + m_E4_str = "Ea-electron";