Branched NO + RO2 reaction

include/micm/configure/solver_config.hpp

@@ -11,6 +11,7 @@
 #include <fstream>
 #include <iostream>
 #include <micm/process/arrhenius_rate_constant.hpp>
+#include <micm/process/branched_rate_constant.hpp>
 #include <micm/process/photolysis_rate_constant.hpp>
 #include <micm/process/process.hpp>
 #include <micm/process/ternary_chemical_activation_rate_constant.hpp>
@@ -91,6 +92,7 @@ namespace micm
     // Read from reaction configure
     std::vector<PhotolysisRateConstant> photolysis_rate_arr_;
     std::vector<ArrheniusRateConstant> arrhenius_rate_arr_;
+    std::vector<BranchedRateConstant> branched_rate_arr_;
     std::vector<TroeRateConstant> troe_rate_arr_;
     std::vector<TernaryChemicalActivationRateConstant> ternary_rate_arr_;
     std::vector<TunnelingRateConstant> tunneling_rate_arr_;
@@ -256,6 +258,10 @@ namespace micm
         {
           status = ParseArrhenius(object);
         }
+        else if (type == "BRANCHED" || type == "WENNBERG_NO_RO2")
+        {
+          status = ParseBranched(object);
+        }
         else if (type == "TERNARY_CHEMICAL_ACTIVATION")
         {
           status = ParseTernaryChemicalActivation(object);
@@ -459,6 +465,48 @@ namespace micm
       return ConfigParseStatus::Success;
     }
 
+    ConfigParseStatus ParseBranched(const json& object)
+    {
+      const std::string REACTANTS = "reactants";
+      const std::string ALKOXY_PRODUCTS = "alkoxy products";
+      const std::string NITRATE_PRODUCTS = "nitrate products";
+      const std::string X = "X";
+      const std::string Y = "Y";
+      const std::string A0 = "a0";
+      const std::string N = "n";
+
+      // Check required json objects exist
+      for (const auto& key : { REACTANTS, ALKOXY_PRODUCTS, NITRATE_PRODUCTS, X, Y, A0, N })
+      {
+        if (!ValidateJsonWithKey(object, key))
+          return ConfigParseStatus::RequiredKeyNotFound;
+      }
+
+      auto reactants = ParseReactants(object[REACTANTS]);
+      auto alkoxy_products = ParseProducts(object[ALKOXY_PRODUCTS]);
+      auto nitrate_products = ParseProducts(object[NITRATE_PRODUCTS]);
+
+      BranchedRateConstantParameters parameters;
+      parameters.X_ = object[X].get<double>();
+      parameters.Y_ = object[Y].get<double>();
+      parameters.a0_ = object[A0].get<double>();
+      parameters.n_ = object[N].get<int>();
+
+      // Alkoxy branch
+      parameters.branch_ = BranchedRateConstantParameters::Branch::Alkoxy;
+      branched_rate_arr_.push_back(BranchedRateConstant(parameters));
+      std::unique_ptr<BranchedRateConstant> rate_ptr = std::make_unique<BranchedRateConstant>(parameters);
+      processes_.push_back(Process(reactants, alkoxy_products, std::move(rate_ptr), gas_phase_));
+
+      // Nitrate branch
+      parameters.branch_ = BranchedRateConstantParameters::Branch::Nitrate;
+      branched_rate_arr_.push_back(BranchedRateConstant(parameters));
+      rate_ptr = std::make_unique<BranchedRateConstant>(parameters);
+      processes_.push_back(Process(reactants, nitrate_products, std::move(rate_ptr), gas_phase_));
+
+      return ConfigParseStatus::Success;
+    }
+
     ConfigParseStatus ParseTroe(const json& object)
     {
       const std::string REACTANTS = "reactants";

include/micm/process/branched_rate_constant.hpp

@@ -0,0 +1,111 @@
+// Copyright (C) 2023 National Center for Atmospheric Research,
+//
+// SPDX-License-Identifier: Apache-2.0
+#pragma once
+
+#include <cmath>
+#include <micm/process/rate_constant.hpp>
+#include <micm/util/constants.hpp>
+
+namespace micm
+{
+
+  struct BranchedRateConstantParameters
+  {
+    enum class Branch
+    {
+      Alkoxy,
+      Nitrate
+    };
+    /// @brief reaction branch
+    Branch branch_;
+    /// @brief pre-exponential factor
+    double X_;
+    /// @brief exponential factor
+    double Y_;
+    /// @brief branching factor
+    double a0_;
+    /// @brief number of heavy atoms in the RO2 reacting species (excluding the peroxy moiety)
+    int n_;
+  };
+
+  /// @brief A Branched rate constant
+  class BranchedRateConstant : public RateConstant
+  {
+   public:
+    const BranchedRateConstantParameters parameters_;
+    const double k0_;
+    const double z_;
+
+   public:
+    /// @brief Default constructor
+    BranchedRateConstant();
+
+    /// @brief An explicit constructor
+    /// @param parameters A set of branched rate constant parameters
+    BranchedRateConstant(const BranchedRateConstantParameters& 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;
+
+    /// @brief Calculate A(T,[M],n)
+    /// @param temperature Temperature in [K]
+    /// @param air_number_density Number density of air in [mol m-3]
+    double A(const double& temperature, const double& air_number_density) const;
+  };
+
+  inline BranchedRateConstant::BranchedRateConstant()
+      : parameters_(),
+        k0_(),
+        z_()
+  {
+  }
+
+  inline BranchedRateConstant::BranchedRateConstant(const BranchedRateConstantParameters& parameters)
+      : parameters_(parameters),
+        k0_(2.0e-22 * AVOGADRO_CONSTANT * 1.0e-6 * std::exp(parameters_.n_)),
+        z_(A(293.0, 2.45e19 / AVOGADRO_CONSTANT * 1.0e6) * (1.0 - parameters_.a0_) / parameters_.a0_)
+  {
+  }
+
+  inline std::unique_ptr<RateConstant> BranchedRateConstant::clone() const
+  {
+    return std::unique_ptr<RateConstant>{ new BranchedRateConstant{ *this } };
+  }
+
+  inline double BranchedRateConstant::calculate(
+      const Conditions& conditions,
+      const std::vector<double>::const_iterator& custom_parameters) const
+  {
+    return calculate(conditions.temperature_, conditions.air_density_);
+  }
+
+  inline double BranchedRateConstant::calculate(const double& temperature, const double& air_number_density) const
+  {
+    double pre = parameters_.X_ * std::exp(-parameters_.Y_ / temperature);
+    double Atmn = A(temperature, air_number_density);
+    if (parameters_.branch_ == BranchedRateConstantParameters::Branch::Alkoxy)
+      return pre * (z_ / (z_ + Atmn));
+    return pre * (Atmn / (Atmn + z_));
+  }
+
+  inline double BranchedRateConstant::A(const double& temperature, const double& air_number_density) const
+  {
+    double a = k0_ * air_number_density;
+    double b = 0.43 * std::pow(temperature / 298.0, -8);
+    return a / (1.0 + a / b) * std::pow(0.41, 1.0 / (1.0 + std::pow(std::log10(a / b), 2)));
+  }
+}  // namespace micm

include/micm/util/constants.hpp

@@ -1 +1,7 @@
-static constexpr double BOLTZMANN_CONSTANT = 1.3806505e-23;  // J K^{-1}
+// Copyright (C) 2023 National Center for Atmospheric Research,
+//
+// SPDX-License-Identifier: Apache-2.0
+#pragma once
+
+static constexpr double BOLTZMANN_CONSTANT = 1.3806505e-23;  // J K^{-1}
+static constexpr double AVOGADRO_CONSTANT = 6.02214076e23;   // # mol^{-1}

test/unit/configure/process/CMakeLists.txt

@@ -6,5 +6,6 @@ include(test_util)
 ################################################################################
 # Tests
 
+create_standard_test(NAME branched_config SOURCES test_branched_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)

test/unit/configure/process/test_branched_config.cpp

@@ -0,0 +1,103 @@
+#include <gtest/gtest.h>
+
+#include <micm/configure/solver_config.hpp>
+
+TEST(BranchedConfig, DetectsInvalidConfig)
+{
+  micm::SolverConfig solver_config;
+
+  // Read and parse the configure files
+  micm::ConfigParseStatus status = solver_config.ReadAndParse("./unit_configs/process/branched/missing_reactants");
+  EXPECT_EQ(micm::ConfigParseStatus::RequiredKeyNotFound, status);
+  status = solver_config.ReadAndParse("./unit_configs/process/branched/missing_alkoxy_products");
+  EXPECT_EQ(micm::ConfigParseStatus::RequiredKeyNotFound, status);
+  status = solver_config.ReadAndParse("./unit_configs/process/branched/missing_nitrate_products");
+  EXPECT_EQ(micm::ConfigParseStatus::RequiredKeyNotFound, status);
+}
+
+TEST(BranchedConfig, ParseConfig)
+{
+  micm::SolverConfig solver_config;
+
+  micm::ConfigParseStatus status = solver_config.ReadAndParse("./unit_configs/process/branched/valid");
+  EXPECT_EQ(micm::ConfigParseStatus::Success, status);
+
+  micm::SolverParameters solver_params = solver_config.GetSolverParams();
+
+  auto& process_vector = solver_params.processes_;
+  EXPECT_EQ(process_vector.size(), 4);
+
+  // 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::BranchedRateConstant* branched_rate_constant =
+        dynamic_cast<micm::BranchedRateConstant*>(process_vector[0].rate_constant_.get());
+    auto& params = branched_rate_constant->parameters_;
+    EXPECT_EQ(params.X_, 12.3);
+    EXPECT_EQ(params.Y_, 42.3);
+    EXPECT_EQ(params.a0_, 1.0e-5);
+    EXPECT_EQ(params.n_, 3);
+    EXPECT_EQ(params.branch_, micm::BranchedRateConstantParameters::Branch::Alkoxy);
+  }
+  {
+    EXPECT_EQ(process_vector[1].reactants_.size(), 3);
+    EXPECT_EQ(process_vector[1].reactants_[0].name_, "foo");
+    EXPECT_EQ(process_vector[1].reactants_[1].name_, "quz");
+    EXPECT_EQ(process_vector[1].reactants_[2].name_, "quz");
+    EXPECT_EQ(process_vector[1].products_.size(), 1);
+    EXPECT_EQ(process_vector[1].products_[0].first.name_, "quz");
+    EXPECT_EQ(process_vector[1].products_[0].second, 1.0);
+    micm::BranchedRateConstant* branched_rate_constant =
+        dynamic_cast<micm::BranchedRateConstant*>(process_vector[1].rate_constant_.get());
+    auto& params = branched_rate_constant->parameters_;
+    EXPECT_EQ(params.X_, 12.3);
+    EXPECT_EQ(params.Y_, 42.3);
+    EXPECT_EQ(params.a0_, 1.0e-5);
+    EXPECT_EQ(params.n_, 3);
+    EXPECT_EQ(params.branch_, micm::BranchedRateConstantParameters::Branch::Nitrate);
+  }
+
+  // second reaction
+  {
+    EXPECT_EQ(process_vector[2].reactants_.size(), 2);
+    EXPECT_EQ(process_vector[2].reactants_[0].name_, "bar");
+    EXPECT_EQ(process_vector[2].reactants_[1].name_, "baz");
+    EXPECT_EQ(process_vector[2].products_.size(), 1);
+    EXPECT_EQ(process_vector[2].products_[0].first.name_, "baz");
+    EXPECT_EQ(process_vector[2].products_[0].second, 1.0);
+    micm::BranchedRateConstant* branched_rate_constant =
+        dynamic_cast<micm::BranchedRateConstant*>(process_vector[2].rate_constant_.get());
+    auto& params = branched_rate_constant->parameters_;
+    EXPECT_EQ(params.X_, 0.32);
+    EXPECT_EQ(params.Y_, 2.3e8);
+    EXPECT_EQ(params.a0_, 0.423);
+    EXPECT_EQ(params.n_, 6);
+    EXPECT_EQ(params.branch_, micm::BranchedRateConstantParameters::Branch::Alkoxy);
+  }
+  {
+    EXPECT_EQ(process_vector[3].reactants_.size(), 2);
+    EXPECT_EQ(process_vector[3].reactants_[0].name_, "bar");
+    EXPECT_EQ(process_vector[3].reactants_[1].name_, "baz");
+    EXPECT_EQ(process_vector[3].products_.size(), 2);
+    EXPECT_EQ(process_vector[3].products_[0].first.name_, "bar");
+    EXPECT_EQ(process_vector[3].products_[0].second, 0.5);
+    EXPECT_EQ(process_vector[3].products_[1].first.name_, "foo");
+    EXPECT_EQ(process_vector[3].products_[1].second, 1.0);
+    micm::BranchedRateConstant* branched_rate_constant =
+        dynamic_cast<micm::BranchedRateConstant*>(process_vector[3].rate_constant_.get());
+    auto& params = branched_rate_constant->parameters_;
+    EXPECT_EQ(params.X_, 0.32);
+    EXPECT_EQ(params.Y_, 2.3e8);
+    EXPECT_EQ(params.a0_, 0.423);
+    EXPECT_EQ(params.n_, 6);
+    EXPECT_EQ(params.branch_, micm::BranchedRateConstantParameters::Branch::Nitrate);
+  }
+}

test/unit/process/CMakeLists.txt

@@ -7,6 +7,7 @@ include(test_util)
 # Tests
 
 create_standard_test(NAME arrhenius_rate_constant SOURCES test_arrhenius_rate_constant.cpp)
+create_standard_test(NAME branched_rate_constant SOURCES test_branched_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)

test/unit/process/test_branched_rate_constant.cpp

@@ -0,0 +1,48 @@
+#include <gtest/gtest.h>
+
+#include <micm/process/branched_rate_constant.hpp>
+#include <micm/solver/state.hpp>
+#include <micm/system/system.hpp>
+#include <micm/util/constants.hpp>
+
+TEST(BranchedRateConstant, CalculateAlkoxyBranchWithAllArugments)
+{
+  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::BranchedRateConstant branched{ micm::BranchedRateConstantParameters{
+      .branch_ = micm::BranchedRateConstantParameters::Branch::Alkoxy, .X_ = 1.2, .Y_ = 204.3, .a0_ = 1.0e-3, .n_ = 2 } };
+  auto k = branched.calculate(state.conditions_[0], params);
+  double air_dens_n_cm3 = 42.2 * AVOGADRO_CONSTANT * 1.0e-6;
+  double a = 2.0e-22 * std::exp(2) * 2.45e19;
+  double b = 0.43 * std::pow((293.0 / 298.0), -8.0);
+  double z = a / (1.0 + a / b) * std::pow(0.41, 1.0 / (1.0 + std::pow(std::log10(a / b), 2))) * (1.0 - 1.0e-3) / 1.0e-3;
+  a = 2.0e-22 * std::exp(2) * air_dens_n_cm3;
+  b = 0.43 * std::pow((temperature / 298.0), -8.0);
+  double A = a / (1.0 + a / b) * std::pow(0.41, 1.0 / (1.0 + std::pow(std::log10(a / b), 2)));
+
+  EXPECT_NEAR(k, 1.2 * std::exp(-204.3 / temperature) * (z / (z + A)), 1.0e-8);
+}
+
+TEST(BranchedRateConstant, CalculateNitrateBranchWithAllArugments)
+{
+  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::BranchedRateConstant branched{ micm::BranchedRateConstantParameters{
+      .branch_ = micm::BranchedRateConstantParameters::Branch::Nitrate, .X_ = 1.2, .Y_ = 204.3, .a0_ = 1.0e-3, .n_ = 2 } };
+  auto k = branched.calculate(state.conditions_[0], params);
+  double air_dens_n_cm3 = 42.2 * AVOGADRO_CONSTANT * 1.0e-6;
+  double a = 2.0e-22 * std::exp(2) * 2.45e19;
+  double b = 0.43 * std::pow((293.0 / 298.0), -8.0);
+  double z = a / (1.0 + a / b) * std::pow(0.41, 1.0 / (1.0 + std::pow(std::log10(a / b), 2))) * (1.0 - 1.0e-3) / 1.0e-3;
+  a = 2.0e-22 * std::exp(2) * air_dens_n_cm3;
+  b = 0.43 * std::pow((temperature / 298.0), -8.0);
+  double A = a / (1.0 + a / b) * std::pow(0.41, 1.0 / (1.0 + std::pow(std::log10(a / b), 2)));
+
+  EXPECT_NEAR(k, 1.2 * std::exp(-204.3 / temperature) * (A / (z + A)), 1.0e-8);
+}

test/unit/unit_configs/process/branched/missing_alkoxy_products/reactions.json

@@ -0,0 +1,20 @@
+{
+  "camp-data": [
+    {
+      "name": "Branched missing alkoxy products",
+      "type": "MECHANISM",
+      "reactions": [
+        {
+          "type": "BRANCHED",
+          "reactants": {
+            "bar": { },
+            "baz": { "qty": 2 }
+          },
+          "nitrate products": {
+            "foo": { }
+          }
+        }
+      ]
+    }
+  ]
+}