Tunneling reaction rate constant (#156)

* add tunneling rate constant and tests * add tunneling config parsing and tests * Auto-format code using Clang-Format --------- Co-authored-by: GitHub Actions <actions@github.com>
NCAR · Jul 27, 2023 · 5f2891b · 5f2891b
1 parent 01f3cd6
commit 5f2891b
Show file tree

Hide file tree

Showing 18 changed files with 378 additions and 17 deletions.
diff --git a/include/micm/configure/solver_config.hpp b/include/micm/configure/solver_config.hpp
@@ -15,6 +15,7 @@
 #include <micm/process/process.hpp>
 #include <micm/process/ternary_chemical_activation_rate_constant.hpp>
 #include <micm/process/troe_rate_constant.hpp>
+#include <micm/process/tunneling_rate_constant.hpp>
 #include <micm/system/phase.hpp>
 #include <micm/system/property.hpp>
 #include <micm/system/species.hpp>
@@ -92,6 +93,7 @@ namespace micm
     std::vector<ArrheniusRateConstant> arrhenius_rate_arr_;
     std::vector<TroeRateConstant> troe_rate_arr_;
     std::vector<TernaryChemicalActivationRateConstant> ternary_rate_arr_;
+    std::vector<TunnelingRateConstant> tunneling_rate_arr_;
     std::vector<Species> emission_arr_;
     std::vector<Species> first_order_loss_arr_;
 
@@ -262,6 +264,10 @@ namespace micm
         {
           status = ParseTroe(object);
         }
+        else if (type == "TUNNELING" || type == "WENNBERG_TUNNELING")
+        {
+          status = ParseTunneling(object);
+        }
         else if (type == "EMISSION")
         {
           status = ParseEmission(object);
@@ -570,6 +576,44 @@ namespace micm
       return ConfigParseStatus::Success;
     }
 
+    ConfigParseStatus ParseTunneling(const json& object)
+    {
+      const std::string REACTANTS = "reactants";
+      const std::string PRODUCTS = "products";
+
+      // Check required json objects exist
+      for (const auto& key : { REACTANTS, PRODUCTS })
+      {
+        if (!ValidateJsonWithKey(object, key))
+          return ConfigParseStatus::RequiredKeyNotFound;
+      }
+
+      auto reactants = ParseReactants(object[REACTANTS]);
+      auto products = ParseProducts(object[PRODUCTS]);
+
+      TunnelingRateConstantParameters parameters;
+      if (object.contains("A"))
+      {
+        parameters.A_ = object["A"].get<double>();
+      }
+      if (object.contains("B"))
+      {
+        parameters.B_ = object["B"].get<double>();
+      }
+      if (object.contains("C"))
+      {
+        parameters.C_ = object["C"].get<double>();
+      }
+
+      tunneling_rate_arr_.push_back(TunnelingRateConstant(parameters));
+
+      std::unique_ptr<TunnelingRateConstant> rate_ptr = std::make_unique<TunnelingRateConstant>(parameters);
+
+      processes_.push_back(Process(reactants, products, std::move(rate_ptr), gas_phase_));
+
+      return ConfigParseStatus::Success;
+    }
+
     ConfigParseStatus ParseEmission(const json& object)
     {
       std::vector<std::string> required_keys = { "species" };

diff --git a/include/micm/process/tunneling_rate_constant.hpp b/include/micm/process/tunneling_rate_constant.hpp
@@ -0,0 +1,79 @@
+// 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 TunnelingRateConstantParameters
+  {
+    /// @brief Pre-exponential factor [(mol m−3)^(−(𝑛−1)) s−1]
+    double A_ = 1.0;
+    /// @brief Linear temperature-dependent parameter [K]
+    double B_ = 0.0;
+    /// @brief Cubed temperature-dependent parameter [K^3]
+    double C_ = 0.0;
+  };
+
+  /// @brief Rate constant for tunneling reactions
+  class TunnelingRateConstant : public RateConstant
+  {
+   public:
+    const TunnelingRateConstantParameters parameters_;
+
+   public:
+    /// @brief Default constructor
+    TunnelingRateConstant();
+
+    /// @brief An explicit constructor
+    /// @param parameters A set of troe rate constants
+    TunnelingRateConstant(const TunnelingRateConstantParameters& 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]
+    /// @return the calculated rate constant
+    double calculate(const double& temperature) const;
+  };
+
+  inline TunnelingRateConstant::TunnelingRateConstant()
+      : parameters_()
+  {
+  }
+
+  inline TunnelingRateConstant::TunnelingRateConstant(const TunnelingRateConstantParameters& parameters)
+      : parameters_(parameters)
+  {
+  }
+
+  inline std::unique_ptr<RateConstant> TunnelingRateConstant::clone() const
+  {
+    return std::unique_ptr<RateConstant>{ new TunnelingRateConstant{ *this } };
+  }
+
+  inline double TunnelingRateConstant::calculate(
+      const Conditions& conditions,
+      const std::vector<double>::const_iterator& custom_parameters) const
+  {
+    return calculate(conditions.temperature_);
+  }
+
+  inline double TunnelingRateConstant::calculate(const double& temperature) const
+  {
+    return parameters_.A_ * std::exp(-parameters_.B_ / temperature + parameters_.C_ / std::pow(temperature, 3));
+  }
+
+}  // namespace micm
diff --git a/test/unit/configure/process/CMakeLists.txt b/test/unit/configure/process/CMakeLists.txt
@@ -6,4 +6,5 @@ include(test_util)
 ################################################################################
 # Tests
 
-create_standard_test(NAME ternary_chemical_activation_config SOURCES test_ternary_chemical_activation_config.cpp)
+create_standard_test(NAME ternary_chemical_activation_config SOURCES test_ternary_chemical_activation_config.cpp)
+create_standard_test(NAME tunneling_config SOURCES test_tunneling_config.cpp)
diff --git a/test/unit/configure/process/test_tunneling_config.cpp b/test/unit/configure/process/test_tunneling_config.cpp
@@ -0,0 +1,63 @@
+#include <gtest/gtest.h>
+
+#include <micm/configure/solver_config.hpp>
+
+TEST(TunnelingConfig, DetectsInvalidConfig)
+{
+  micm::SolverConfig solver_config;
+
+  // Read and parse the configure files
+  micm::ConfigParseStatus status = solver_config.ReadAndParse("./unit_configs/process/tunneling/missing_reactants");
+  EXPECT_EQ(micm::ConfigParseStatus::RequiredKeyNotFound, status);
+  status = solver_config.ReadAndParse("./unit_configs/process/tunneling/missing_products");
+  EXPECT_EQ(micm::ConfigParseStatus::RequiredKeyNotFound, status);
+}
+
+TEST(TunnelingConfig, ParseConfig)
+{
+  micm::SolverConfig solver_config;
+
+  micm::ConfigParseStatus status = solver_config.ReadAndParse("./unit_configs/process/tunneling/valid");
+  EXPECT_EQ(micm::ConfigParseStatus::Success, status);
+
+  micm::SolverParameters solver_params = solver_config.GetSolverParams();
+
+  auto& process_vector = solver_params.processes_;
+
+  // first reaction
+  {
+    EXPECT_EQ(process_vector[0].reactants_.size(), 3);
+    EXPECT_EQ(process_vector[0].reactants_[0].name_, "foo");
+    EXPECT_EQ(process_vector[0].reactants_[1].name_, "quz");
+    EXPECT_EQ(process_vector[0].reactants_[2].name_, "quz");
+    EXPECT_EQ(process_vector[0].products_.size(), 2);
+    EXPECT_EQ(process_vector[0].products_[0].first.name_, "bar");
+    EXPECT_EQ(process_vector[0].products_[0].second, 1.0);
+    EXPECT_EQ(process_vector[0].products_[1].first.name_, "baz");
+    EXPECT_EQ(process_vector[0].products_[1].second, 3.2);
+    micm::TunnelingRateConstant* tunneling_rate_constant =
+        dynamic_cast<micm::TunnelingRateConstant*>(process_vector[0].rate_constant_.get());
+    auto& params = tunneling_rate_constant->parameters_;
+    EXPECT_EQ(params.A_, 1.0);
+    EXPECT_EQ(params.B_, 0.0);
+    EXPECT_EQ(params.C_, 0.0);
+  }
+
+  // second reaction
+  {
+    EXPECT_EQ(process_vector[1].reactants_.size(), 2);
+    EXPECT_EQ(process_vector[1].reactants_[0].name_, "bar");
+    EXPECT_EQ(process_vector[1].reactants_[1].name_, "baz");
+    EXPECT_EQ(process_vector[1].products_.size(), 2);
+    EXPECT_EQ(process_vector[1].products_[0].first.name_, "bar");
+    EXPECT_EQ(process_vector[1].products_[0].second, 0.5);
+    EXPECT_EQ(process_vector[1].products_[1].first.name_, "foo");
+    EXPECT_EQ(process_vector[1].products_[1].second, 1.0);
+    micm::TunnelingRateConstant* tunneling_rate_constant =
+        dynamic_cast<micm::TunnelingRateConstant*>(process_vector[1].rate_constant_.get());
+    auto& params = tunneling_rate_constant->parameters_;
+    EXPECT_EQ(params.A_, 32.1);
+    EXPECT_EQ(params.B_, -2.3);
+    EXPECT_EQ(params.C_, 102.3);
+  }
+}
diff --git a/test/unit/process/CMakeLists.txt b/test/unit/process/CMakeLists.txt
@@ -10,6 +10,7 @@ create_standard_test(NAME arrhenius_rate_constant SOURCES test_arrhenius_rate_co
 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 tunneling_rate_constant SOURCES test_tunneling_rate_constant.cpp)
 create_standard_test(NAME process_set SOURCES test_process_set.cpp)
 
 # GPU tests

diff --git a/test/unit/process/test_ternary_chemical_activation_rate_constant.cpp b/test/unit/process/test_ternary_chemical_activation_rate_constant.cpp
@@ -10,11 +10,11 @@ TEST(TernaryChemicalActivationRateConstant, CalculateWithMinimalArugments)
   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);
+  micm::TernaryChemicalActivationRateConstantParameters ternary_params;
+  ternary_params.k0_A_ = 1.0;
+  ternary_params.kinf_A_ = 1.0;
+  micm::TernaryChemicalActivationRateConstant ternary{ ternary_params };
+  auto k = ternary.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);
@@ -27,15 +27,16 @@ TEST(TernaryChemicalActivationRateConstant, CalculateWithAllArugments)
   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);
+  micm::TernaryChemicalActivationRateConstant ternary{ 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 = ternary.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);

diff --git a/test/unit/process/test_tunneling_rate_constant.cpp b/test/unit/process/test_tunneling_rate_constant.cpp
@@ -0,0 +1,27 @@
+#include <gtest/gtest.h>
+
+#include <micm/process/tunneling_rate_constant.hpp>
+#include <micm/solver/state.hpp>
+#include <micm/system/system.hpp>
+
+TEST(TunnelingRateConstant, CalculateWithMinimalArugments)
+{
+  micm::State<micm::Matrix> state{ 0, 0, 1 };
+  state.conditions_[0].temperature_ = 301.24;  // [K]
+  std::vector<double>::const_iterator params = state.custom_rate_parameters_[0].begin();
+  micm::TunnelingRateConstantParameters tunneling_params;
+  micm::TunnelingRateConstant tunneling{ tunneling_params };
+  auto k = tunneling.calculate(state.conditions_[0], params);
+  EXPECT_NEAR(k, 1.0, 1.0e-8);
+}
+
+TEST(TunnelingRateConstant, CalculateWithAllArugments)
+{
+  micm::State<micm::Matrix> state{ 0, 0, 1 };
+  double temperature = 301.24;
+  state.conditions_[0].temperature_ = temperature;  // [K]
+  std::vector<double>::const_iterator params = state.custom_rate_parameters_[0].begin();
+  micm::TunnelingRateConstant tunneling{ micm::TunnelingRateConstantParameters{ .A_ = 1.2, .B_ = 2.3, .C_ = 302.3 } };
+  auto k = tunneling.calculate(state.conditions_[0], params);
+  EXPECT_NEAR(k, 1.2 * std::exp(-2.3 / temperature) * std::exp(302.3 / std::pow(temperature, 3)), 1.0e-8);
+}
diff --git a/test/unit/unit_configs/process/ternary_chemical_activation/missing_products/reactions.json b/test/unit/unit_configs/process/ternary_chemical_activation/missing_products/reactions.json
@@ -1,7 +1,7 @@
 {
   "camp-data": [
     {
-      "name": "Ternary missing reactants",
+      "name": "Ternary missing products",
       "type": "MECHANISM",
       "reactions": [
         {

diff --git a/test/unit/unit_configs/process/ternary_chemical_activation/valid/reactions.json b/test/unit/unit_configs/process/ternary_chemical_activation/valid/reactions.json
@@ -1,7 +1,7 @@
 {
   "camp-data": [
     {
-      "name": "Ternary missing reactants",
+      "name": "Valid ternary chemical activation reactions",
       "type": "MECHANISM",
       "reactions": [
         {

diff --git a/test/unit/unit_configs/process/tunneling/missing_products/reactions.json b/test/unit/unit_configs/process/tunneling/missing_products/reactions.json
@@ -0,0 +1,17 @@
+{
+  "camp-data": [
+    {
+      "name": "Tunneling missing products",
+      "type": "MECHANISM",
+      "reactions": [
+        {
+          "type": "TUNNELING",
+          "reactants": {
+            "bar": { },
+            "baz": { "qty": 2 }
+          }
+        }
+      ]
+    }
+  ]
+}
diff --git a/test/unit/unit_configs/process/tunneling/missing_products/species.json b/test/unit/unit_configs/process/tunneling/missing_products/species.json
@@ -0,0 +1,16 @@
+{
+  "camp-data": [
+    {
+      "name": "foo",
+      "type": "CHEM_SPEC"
+    },
+    {
+      "name": "bar",
+      "type": "CHEM_SPEC"
+    },
+    {
+      "name": "baz",
+      "type": "CHEM_SPEC"
+    }
+  ]
+}
diff --git a/test/unit/unit_configs/process/tunneling/missing_products/tolerance.json b/test/unit/unit_configs/process/tunneling/missing_products/tolerance.json
@@ -0,0 +1,8 @@
+{
+  "camp-data": [
+    {
+      "type": "RELATIVE_TOLERANCE",
+      "value": 1.0e-4
+    }
+  ]
+}
diff --git a/test/unit/unit_configs/process/tunneling/missing_reactants/reactions.json b/test/unit/unit_configs/process/tunneling/missing_reactants/reactions.json
@@ -0,0 +1,17 @@
+{
+  "camp-data": [
+    {
+      "name": "Tunneling missing reactants",
+      "type": "MECHANISM",
+      "reactions": [
+        {
+          "type": "TUNNELING",
+          "products": {
+            "bar": { "yield": 1.0 },
+            "baz": { "yield": 3.2 }
+          }
+        }
+      ]
+    }
+  ]
+}