Merge pull request #227 from eagles-project/oscar/mo_usrrxt

Oscar/mo usrrxt

src/mam4xx/gas_chem.hpp

@@ -24,6 +24,51 @@ const Real rel_err = 1.0e-3;
 const Real high_rel_err = 1.0e-4;
 const int max_time_steps = 1000;
 
+KOKKOS_INLINE_FUNCTION
+void usrrxt(Real rxt[rxntot], // inout
+            const Real temperature, const Real invariants[nfs], const Real mtot,
+            const int usr_HO2_HO2_ndx, const int usr_DMS_OH_ndx,
+            const int usr_SO2_OH_ndx, const int inv_h2o_ndx) {
+
+  /*-----------------------------------------------------------------
+   ... ho2 + ho2 --> h2o2
+   note: this rate involves the water vapor number density
+  -----------------------------------------------------------------*/
+  const Real one = 1.0;
+  if (usr_HO2_HO2_ndx > 0) {
+    // BAD CONSTANT
+    const Real ko = 3.5e-13 * haero::exp(430.0 / temperature);
+    const Real kinf = 1.7e-33 * mtot * haero::exp(1000. / temperature);
+    const Real fc = one + 1.4e-21 * invariants[inv_h2o_ndx] *
+                              haero::exp(2200. / temperature);
+    rxt[usr_HO2_HO2_ndx] = (ko + kinf) * fc;
+  }
+
+  /*-----------------------------------------------------------------
+       ... DMS + OH  --> .5 * SO2
+   -----------------------------------------------------------------*/
+  if (usr_DMS_OH_ndx > 0) {
+    // BAD CONSTANT
+    const Real ko =
+        one + 5.5e-31 * haero::exp(7460. / temperature) * mtot * 0.21;
+    rxt[usr_DMS_OH_ndx] =
+        1.7e-42 * haero::exp(7810. / temperature) * mtot * 0.21 / ko;
+  }
+
+  /*-----------------------------------------------------------------
+         ... SO2 + OH  --> SO4  (REFERENCE?? - not Liao)
+  -----------------------------------------------------------------*/
+  if (usr_SO2_OH_ndx > 0) {
+    // BAD CONSTANT
+    const Real fc = 3.0e-31 * haero::pow(300. / temperature, 3.3);
+    const Real ko = fc * mtot / (one + fc * mtot / 1.5e-12);
+    rxt[usr_SO2_OH_ndx] =
+        ko * haero::pow(0.6, one / (one + haero::square(haero::log10(
+                                              fc * mtot / 1.5e-12))));
+  }
+
+} // usrrxt
+
 // initialize the solver (error tolerance)
 KOKKOS_INLINE_FUNCTION
 void imp_slv_inti(Real epsilon[clscnt4]) {

src/validation/gas_chem/CMakeLists.txt

@@ -18,6 +18,7 @@ add_executable(gas_chem_driver
                imp_sol.cpp
                adjrxt.cpp
                setrxt.cpp
+               usrrxt.cpp
                )
 target_link_libraries(gas_chem_driver skywalker;validation;${HAERO_LIBRARIES})
 
@@ -42,6 +43,8 @@ set(TEST_LIST
     imp_sol_ts_355
     adjrxt_ts_1400
     setrxt_ts_1400
+    usrrxt_merged
+    usrrxt_ts_1416
     )
 
 set(DEFAULT_TOL 1e-12)
@@ -56,6 +59,8 @@ set(ERROR_THRESHOLDS
    ${DEFAULT_TOL}
    ${DEFAULT_TOL}
    ${DEFAULT_TOL}
+   ${DEFAULT_TOL}
+   ${DEFAULT_TOL}
    )
 
 foreach(input tol IN ZIP_LISTS TEST_LIST ERROR_THRESHOLDS)

src/validation/gas_chem/gas_chem_driver.cpp

@@ -33,6 +33,7 @@ void newton_raphson_iter(Ensemble *ensemble);
 void imp_sol(Ensemble *ensemble);
 void adjrxt(Ensemble *ensemble);
 void setrxt(Ensemble *ensemble);
+void usrrxt(Ensemble *ensemble);
 
 int main(int argc, char **argv) {
   if (argc == 1) {
@@ -72,6 +73,12 @@ int main(int argc, char **argv) {
       adjrxt(ensemble);
     } else if (func_name == "setrxt") {
       setrxt(ensemble);
+    } else if (func_name == "usrrxt") {
+      usrrxt(ensemble);
+    } else {
+      std::cerr << "Error: Function name '" << func_name
+                << "' does not have an implemented test!" << std::endl;
+      exit(1);
     }
 
   } catch (std::exception &e) {

src/validation/gas_chem/usrrxt.cpp

@@ -0,0 +1,36 @@
+// mam4xx: Copyright (c) 2022,
+// Battelle Memorial Institute and
+// National Technology & Engineering Solutions of Sandia, LLC (NTESS)
+// SPDX-License-Identifier: BSD-3-Clause
+
+#include <mam4xx/mam4.hpp>
+
+#include <mam4xx/aero_config.hpp>
+#include <mam4xx/gas_chem.hpp>
+#include <skywalker.hpp>
+#include <validation.hpp>
+
+using namespace skywalker;
+using namespace mam4;
+using namespace gas_chemistry;
+
+void usrrxt(Ensemble *ensemble) {
+
+  ensemble->process([=](const Input &input, Output &output) {
+    const Real temperature = input.get_array("temp")[0];
+    const Real mtot = input.get_array("mtot")[0];
+
+    auto rxt = input.get_array("rxt");
+    const auto invariants = input.get_array("invariants");
+    const int usr_HO2_HO2_ndx = int(input.get_array("usr_HO2_HO2_ndx")[0]) - 1;
+    const int usr_DMS_OH_ndx = int(input.get_array("usr_DMS_OH_ndx")[0]) - 1;
+    const int usr_SO2_OH_ndx = int(input.get_array("usr_SO2_OH_ndx")[0]) - 1;
+    const int inv_h2o_ndx = int(input.get_array("inv_h2o_ndx")[0]) - 1;
+
+    usrrxt(rxt.data(), // inout
+           temperature, invariants.data(), mtot, usr_HO2_HO2_ndx,
+           usr_DMS_OH_ndx, usr_SO2_OH_ndx, inv_h2o_ndx);
+
+    output.set("rxt", rxt);
+  });
+}