add ternary rate constant class

include/micm/process/arrhenius_rate_constant.hpp

@@ -1,7 +1,6 @@
-/* Copyright (C) 2023 National Center for Atmospheric Research,
- *
- * SPDX-License-Identifier: Apache-2.0
- */
+// Copyright (C) 2023 National Center for Atmospheric Research,
+//
+// SPDX-License-Identifier: Apache-2.0
 #pragma once
 
 #include <cmath>
@@ -11,24 +10,22 @@ namespace micm
 {
   struct ArrheniusRateConstantParameters
   {
-    /// @brief Pre-exponential factor, (cm−3)^(−(𝑛−1))s−1
+    /// @brief Pre-exponential factor [(mol m−3)^(−(𝑛−1)) s−1]
     double A_{ 1 };
     /// @brief Unitless exponential factor
     double B_{ 0 };
-    /// @brief Activation threshold, expected to be the negative activation energy divided by the boltzman constant (-E_a /
-    /// k_b), K
+    /// @brief Activation threshold, expected to be the negative activation energy divided by the boltzman constant
+    ///        [-E_a / k_b), K]
     double C_{ 0 };
-    /// @brief A factor that determines temperature dependence, (K)
+    /// @brief A factor that determines temperature dependence [K]
     double D_{ 300 };
-    /// @brief A factor that determines pressure dependence (Pa-1)
+    /// @brief A factor that determines pressure dependence [Pa-1]
     double E_{ 0 };
   };
 
-  /**
-   * @brief An arrhenius rate constant dependent on temperature and pressure
-   *
-   * More information can be found here: https://open-atmos.github.io/camp/html/camp_rxn_arrhenius.html
-   */
+  /// @brief An arrhenius rate constant dependent on temperature and pressure
+  ///
+  /// More information can be found here: https://open-atmos.github.io/camp/html/camp_rxn_arrhenius.html
   class ArrheniusRateConstant : public RateConstant
   {
    public:

include/micm/process/ternary_chemical_activation_rate_constant.hpp

@@ -0,0 +1,98 @@
+// Copyright (C) 2023 National Center for Atmospheric Research,
+//
+// SPDX-License-Identifier: Apache-2.0
+#pragma once
+
+#include <cmath>
+#include <micm/process/rate_constant.hpp>
+
+namespace micm
+{
+
+  struct TernaryChemicalActivationRateConstantParameters
+  {
+    /// @brief low-pressure pre-exponential factor
+    double k0_A_;
+    /// @brief low-pressure temperature-scaling parameter
+    double k0_B_ = 0.0;
+    /// @brief low-pressure exponential factor
+    double k0_C_ = 0.0;
+    /// @brief high-pressure pre-exponential factor
+    double kinf_A_;
+    /// @brief high-pressure temperature-scaling parameter
+    double kinf_B_ = 0.0;
+    /// @brief high-pressure exponential factor
+    double kinf_C_ = 0.0;
+    /// @brief TernaryChemicalActivation F_c parameter
+    double Fc_ = 0.6;
+    /// @brief TernaryChemicalActivation N parameter
+    double N_ = 1.0;
+  };
+
+  /**
+   * @brief A TernaryChemicalActivation rate constant
+   *
+   */
+  class TernaryChemicalActivationRateConstant : public RateConstant
+  {
+   public:
+    const TernaryChemicalActivationRateConstantParameters parameters_;
+
+   public:
+    /// @brief Default constructor
+    TernaryChemicalActivationRateConstant();
+
+    /// @brief An explicit constructor
+    /// @param parameters A set of troe rate constants
+    TernaryChemicalActivationRateConstant(const TernaryChemicalActivationRateConstantParameters& parameters);
+
+    /// @brief Deep copy
+    std::unique_ptr<RateConstant> clone() const override;
+
+    /// @brief Calculate the rate constant
+    /// @param conditions The current environmental conditions of the chemical system
+    /// @param custom_parameters User-defined rate constant parameters
+    /// @return A rate constant based off of the conditions in the system
+    double calculate(const Conditions& conditions, const std::vector<double>::const_iterator& custom_parameters)
+        const override;
+
+    /// @brief Calculate the rate constant
+    /// @param temperature Temperature in [K]
+    /// @param air_number_density Number density in [mol m-3]
+    /// @return
+    double calculate(const double& temperature, const double& air_number_density) const;
+  };
+
+  inline TernaryChemicalActivationRateConstant::TernaryChemicalActivationRateConstant()
+      : parameters_()
+  {
+  }
+
+  inline TernaryChemicalActivationRateConstant::TernaryChemicalActivationRateConstant(const TernaryChemicalActivationRateConstantParameters& parameters)
+      : parameters_(parameters)
+  {
+  }
+
+  inline std::unique_ptr<RateConstant> TernaryChemicalActivationRateConstant::clone() const
+  {
+    return std::unique_ptr<RateConstant>{ new TernaryChemicalActivationRateConstant{ *this } };
+  }
+
+  inline double TernaryChemicalActivationRateConstant::calculate(
+      const Conditions& conditions,
+      const std::vector<double>::const_iterator& custom_parameters) const
+  {
+    return calculate(conditions.temperature_, conditions.air_density_);
+  }
+
+  inline double TernaryChemicalActivationRateConstant::calculate(const double& temperature, const double& air_number_density) const
+  {
+    double k0 = parameters_.k0_A_ * std::exp(parameters_.k0_C_ / temperature) * pow(temperature / 300.0, parameters_.k0_B_);
+    double kinf =
+        parameters_.kinf_A_ * std::exp(parameters_.kinf_C_ / temperature) * pow(temperature / 300.0, parameters_.kinf_B_);
+
+    return k0 / (1.0 + k0 * air_number_density / kinf) *
+           pow(parameters_.Fc_, 1.0 / (1.0 + 1.0 / parameters_.N_ * pow(log10(k0 * air_number_density / kinf), 2)));
+  }
+
+}  // namespace micm

include/micm/process/troe_rate_constant.hpp

@@ -1,7 +1,6 @@
-/* Copyright (C) 2023 National Center for Atmospheric Research,
- *
- * SPDX-License-Identifier: Apache-2.0
- */
+// Copyright (C) 2023 National Center for Atmospheric Research,
+//
+// SPDX-License-Identifier: Apache-2.0
 #pragma once
 
 #include <cmath>
@@ -59,7 +58,7 @@ namespace micm
 
     /// @brief Calculate the rate constant
     /// @param temperature Temperature in [K]
-    /// @param air_number_density Number density in [# cm-3]
+    /// @param air_number_density Number density in [mol m-3]
     /// @return
     double calculate(const double& temperature, const double& air_number_density) const;
   };

test/unit/process/CMakeLists.txt

@@ -8,6 +8,7 @@ include(test_util)
 
 create_standard_test(NAME arrhenius_rate_constant SOURCES test_arrhenius_rate_constant.cpp)
 create_standard_test(NAME photolysis_rate_constant SOURCES test_photolysis_rate_constant.cpp)
+create_standard_test(NAME ternary_chemical_activation_rate_constant SOURCES test_ternary_chemical_activation_rate_constant.cpp)
 create_standard_test(NAME troe_rate_constant SOURCES test_troe_rate_constant.cpp)
 create_standard_test(NAME process_set SOURCES test_process_set.cpp)
 

test/unit/process/test_ternary_chemical_activation_rate_constant.cpp

@@ -0,0 +1,42 @@
+#include <gtest/gtest.h>
+
+#include <micm/process/ternary_chemical_activation_rate_constant.hpp>
+#include <micm/solver/state.hpp>
+#include <micm/system/system.hpp>
+
+TEST(TernaryChemicalActivationRateConstant, CalculateWithMinimalArugments)
+{
+  micm::State<micm::Matrix> state{ 0, 0, 1 };
+  state.conditions_[0].temperature_ = 301.24;  // [K]
+  state.conditions_[0].air_density_ = 42.2;    // [mol mol-1]
+  std::vector<double>::const_iterator params = state.custom_rate_parameters_[0].begin();
+  micm::TernaryChemicalActivationRateConstantParameters troe_params;
+  troe_params.k0_A_ = 1.0;
+  troe_params.kinf_A_ = 1.0;
+  micm::TernaryChemicalActivationRateConstant troe{ troe_params };
+  auto k = troe.calculate(state.conditions_[0], params);
+  double k0 = 1.0;
+  double kinf = 1.0;
+  EXPECT_NEAR(k, k0 / (1.0 + 42.4 * k0 / kinf) * pow(0.6, 1.0 / (1 + pow(log10(42.2 * k0 / kinf), 2))), 0.001);
+}
+
+TEST(TernaryChemicalActivationRateConstant, CalculateWithAllArugments)
+{
+  micm::State<micm::Matrix> state{ 0, 0, 1 };
+  double temperature = 301.24;
+  state.conditions_[0].temperature_ = temperature;  // [K]
+  state.conditions_[0].air_density_ = 42.2;         // [mol mol-1]
+  std::vector<double>::const_iterator params = state.custom_rate_parameters_[0].begin();
+  micm::TernaryChemicalActivationRateConstant troe{ micm::TernaryChemicalActivationRateConstantParameters{ .k0_A_ = 1.2,
+                                                                                                           .k0_B_ = 2.3,
+                                                                                                           .k0_C_ = 302.3,
+                                                                                                           .kinf_A_ = 2.6,
+                                                                                                           .kinf_B_ = -3.1,
+                                                                                                           .kinf_C_ = 402.1,
+                                                                                                           .Fc_ = 0.9,
+                                                                                                           .N_ = 1.2 } };
+  auto k = troe.calculate(state.conditions_[0], params);
+  double k0 = 1.2 * exp(302.3 / temperature) * pow(temperature / 300.0, 2.3);
+  double kinf = 2.6 * exp(402.1 / temperature) * pow(temperature / 300.0, -3.1);
+  EXPECT_NEAR(k, k0 / (1.0 + 42.2 * k0 / kinf) * pow(0.9, 1.0 / (1.0 + 1.0 / 1.2 * pow(log10(42.2 * k0 / kinf), 2))), 0.001);
+}

test/unit/process/test_troe_rate_constant.cpp

@@ -8,7 +8,7 @@ TEST(TroeRateConstant, CalculateWithMinimalArugments)
 {
   micm::State<micm::Matrix> state{ 0, 0, 1 };
   state.conditions_[0].temperature_ = 301.24;  // [K]
-  state.conditions_[0].air_density_ = 1.0;     // [mol mol-1]
+  state.conditions_[0].air_density_ = 42.2;     // [mol mol-1]
   std::vector<double>::const_iterator params = state.custom_rate_parameters_[0].begin();
   micm::TroeRateConstantParameters troe_params;
   troe_params.k0_A_ = 1.0;
@@ -17,15 +17,15 @@ TEST(TroeRateConstant, CalculateWithMinimalArugments)
   auto k = troe.calculate(state.conditions_[0], params);
   double k0 = 1.0;
   double kinf = 1.0;
-  EXPECT_NEAR(k, k0 / (1.0 + k0 / kinf) * pow(0.6, 1.0 / (1 + pow(log10(k0 / kinf), 2))), 0.001);
+  EXPECT_NEAR(k, 42.2 * k0 / (1.0 + 42.2 * k0 / kinf) * pow(0.6, 1.0 / (1 + pow(log10(42.2 * k0 / kinf), 2))), 0.001);
 }
 
 TEST(TroeRateConstant, CalculateWithAllArugments)
 {
   micm::State<micm::Matrix> state{ 0, 0, 1 };
   double temperature = 301.24;
   state.conditions_[0].temperature_ = temperature;  // [K]
-  state.conditions_[0].air_density_ = 1.0;          // [mol mol-1]
+  state.conditions_[0].air_density_ = 42.2;          // [mol mol-1]
   std::vector<double>::const_iterator params = state.custom_rate_parameters_[0].begin();
   micm::TroeRateConstant troe{ micm::TroeRateConstantParameters{ .k0_A_ = 1.2,
                                                                  .k0_B_ = 2.3,
@@ -38,5 +38,5 @@ TEST(TroeRateConstant, CalculateWithAllArugments)
   auto k = troe.calculate(state.conditions_[0], params);
   double k0 = 1.2 * exp(302.3 / temperature) * pow(temperature / 300.0, 2.3);
   double kinf = 2.6 * exp(402.1 / temperature) * pow(temperature / 300.0, -3.1);
-  EXPECT_NEAR(k, k0 / (1.0 + k0 / kinf) * pow(0.9, 1.0 / (1.0 + 1.0 / 1.2 * pow(log10(k0 / kinf), 2))), 0.001);
+  EXPECT_NEAR(k, 42.2 * k0 / (1.0 + 42.2 * k0 / kinf) * pow(0.9, 1.0 / (1.0 + 1.0 / 1.2 * pow(log10(42.2 * k0 / kinf), 2))), 0.001);
 }