add implementation and testing for MOSAIC function drh_mutual (#109)

* add coefficients for MDRH(T) calculation

* add implementation of drh_mutual

* add fix array size and typos

* add verification datasets for drh_mutual

* T -> T_K for consistency

* testing for drh_mutual

Author: zdaq12

src/core/impl/TChem_Impl_MOSAIC.hpp

@@ -0,0 +1,7 @@
+
+TChem-atm:
+  function: drh_mutual
+input:
+  fixed:
+    j_index: [         1]
+    T_K: [  0.29142360000000002E+003]

src/verification/data_sets/mosaic/drh_mutual_input_ts_0.yaml

@@ -0,0 +1,17 @@
+
+# This file was generated by PartMC-mosaic.
+
+from math import nan as nan, inf as inf
+
+# Object is just a dynamic container that stores input/output data.
+class Object(object):
+    pass
+# Settings are stored here.
+settings = Object()
+# Input is stored here.
+input = Object()
+input.j_index=[[         1],]
+input.T_K=[[  0.29142360000000002E+003],]
+# Output data is stored here.
+output = Object()
+output.drh_mutual=[[  0.80202077162854152E+002],]

src/verification/data_sets/mosaic/drh_mutual_output_ts_0.py

@@ -5,5 +5,5 @@ input:
   fixed:
     Po_298: [  0.57000000000000003E-004]
     DH: [  0.15600000000000000E+003]
-    T: [  0.29142360000000002E+003]
+    T_K: [  0.29142360000000002E+003]
     Po: [  0.29005922138061241E-316]

src/verification/data_sets/mosaic/fn_Po_input_ts_0.yaml

@@ -12,7 +12,7 @@ class Object(object):
 input = Object()
 input.Po_298=[[  0.57000000000000003E-004],]
 input.DH=[[  0.15600000000000000E+003],]
-input.T=[[  0.29142360000000002E+003],]
+input.T_K=[[  0.29142360000000002E+003],]
 input.Po=[[  0.29005922138061241E-316],]
 # Output data is stored here.
 output = Object()

src/verification/data_sets/mosaic/fn_Po_output_ts_0.py

@@ -10,6 +10,7 @@ add_executable(${DRIVER_NAME}.x
                calculate_XT.cpp
                fnlog_gamZ.cpp
                fn_Keq.cpp
+               drh_mutual.cpp
                fn_Po.cpp
                molality_0.cpp
                bin_molality.cpp
@@ -31,6 +32,7 @@ SET(TEST_LIST
      calculate_XT_input_ts_0
      fnlog_gamZ_input_ts_0
      fn_Keq_input_ts_0
+     drh_mutual_input_ts_0
      fn_Po_input_ts_0
      molality_0_input_ts_0
      bin_molality_input_ts_0
@@ -47,6 +49,7 @@ SET(TEST_BASELINE
     calculate_XT_output_ts_0
     fnlog_gamZ_output_ts_0
     fn_Keq_output_ts_0
+    drh_mutual_output_ts_0
     fn_Po_output_ts_0
     molality_0_output_ts_0
     bin_molality_output_ts_0
@@ -64,6 +67,7 @@ set(ERROR_THRESHOLDS
     ${DEFAULT_TOL} # calculate_XT_output_ts_0
     1e-6           # fnlog_gamZ_output_ts_0
     1e-6           # fn_Keq_output_ts_0
+    1e-6           # drh_mutual_output_ts_0
     1e-8           # fn_Po_output_ts_0
     1e-5           # molality_0_output_ts_0
     1e-5           # bin_molality_output_ts_0

src/verification/mosaic/CMakeLists.txt

@@ -0,0 +1,59 @@
+#include "TChem.hpp"
+#include "TChem_Impl_MOSAIC.hpp"
+#include <verification.hpp>
+#include "skywalker.hpp"
+
+using device_type = typename Tines::UseThisDevice<TChem::exec_space>::type;
+using real_type_1d_view = TChem::real_type_1d_view;
+using oridinal_type_1d_view = TChem::ordinal_type_1d_view;
+using ordinal_type = TChem::ordinal_type;
+using namespace skywalker;
+using namespace TChem;
+
+void drh_mutual(Ensemble *ensemble) {
+  ensemble->process([=](const Input &input, Output &output) {
+
+    const auto j_index_arr = input.get_array("j_index");
+    const auto T_K_arr = input.get_array("T_K");
+
+    real_type_1d_view j_index("j_index", 1);
+    verification::convert_1d_vector_to_1d_view_device(j_index_arr, j_index);
+
+    real_type_1d_view T_K("T_K", 1);
+    verification::convert_1d_vector_to_1d_view_device(T_K_arr, T_K);
+
+    // Reals or int that are defined outside of the parallel_for region are passed as const.
+    real_type_1d_view outputs_drh_mutual("outputs_drh_mutual", 1);
+
+    const auto mmd = TChem::Impl::MosaicModelData<device_type>();
+
+    std::string profile_name ="Verification_test_drh_mutual";
+    using policy_type =
+          typename TChem::UseThisTeamPolicy<TChem::exec_space>::type;
+    const auto exec_space_instance = TChem::exec_space();
+    const auto host_exec_space = TChem::host_exec_space();
+    policy_type policy(exec_space_instance, 1, Kokkos::AUTO());
+
+     Kokkos::parallel_for(
+    profile_name,
+    policy,
+    KOKKOS_LAMBDA(const typename policy_type::member_type& member) {
+
+      Real& drh_mutual = outputs_drh_mutual(0);
+
+    // Perform the adjustment calculation
+    TChem::Impl::MOSAIC<real_type, device_type>::drh_mutual(
+      mmd,
+      j_index(0)-1,
+      T_K(0),
+      drh_mutual);
+    });
+
+    const auto outputs_drh_mutual_h = Kokkos::create_mirror_view_and_copy(host_exec_space, outputs_drh_mutual);
+
+    Real drh_mutual = outputs_drh_mutual_h(0);
+
+    output.set("drh_mutual", drh_mutual);
+
+  });
+}

src/verification/mosaic/drh_mutual.cpp

@@ -15,16 +15,16 @@ void fn_Po(Ensemble *ensemble) {
 
     const auto Po_298_arr = input.get_array("Po_298");
     const auto DH_arr = input.get_array("DH");
-    const auto T_arr = input.get_array("T");
+    const auto T_K_arr = input.get_array("T_K");
 
     real_type_1d_view Po_298("Po_298", 1);
     verification::convert_1d_vector_to_1d_view_device(Po_298_arr, Po_298);
 
     real_type_1d_view DH("DH", 1);
     verification::convert_1d_vector_to_1d_view_device(DH_arr, DH);
 
-    real_type_1d_view T("T", 1);
-    verification::convert_1d_vector_to_1d_view_device(T_arr, T);
+    real_type_1d_view T_K("T_K", 1);
+    verification::convert_1d_vector_to_1d_view_device(T_K_arr, T_K);
 
     // Reals or int that are defined outside of the parallel_for region are passed as const.
     real_type_1d_view outputs_fn_Po("outputs_fn_Po", 1);
@@ -47,7 +47,7 @@ void fn_Po(Ensemble *ensemble) {
     TChem::Impl::MOSAIC<real_type, device_type>::fn_Po(
       Po_298(0),
       DH(0),
-      T(0),
+      T_K(0),
       Po);
     });
 

src/verification/mosaic/fn_Po.cpp

@@ -29,6 +29,8 @@ void fnlog_gamZ(Ensemble *ensemble);
 
 void fn_Keq(Ensemble *ensemble);
 
+void drh_mutual(Ensemble *ensemble);
+
 void fn_Po(Ensemble *ensemble);
 
 void molality_0(Ensemble *ensemble);
@@ -78,6 +80,8 @@ int main(int argc, char **argv) {
       fnlog_gamZ(ensemble);
     } else if (func_name == "fn_Keq") {
       fn_Keq(ensemble);
+    } else if (func_name == "drh_mutual") {
+      drh_mutual(ensemble);
     } else if (func_name == "fn_Po") {
       fn_Po(ensemble);
     } else if (func_name == "molality_0") {