het_diags from mo_chm_diags

src/mam4xx/CMakeLists.txt

@@ -42,6 +42,7 @@ install(FILES
         ndrop.hpp
         mo_photo.hpp
         lin_strat_chem.hpp
+        mo_chm_diags.hpp
         DESTINATION include/mam4xx)
 
 add_library(mam4xx aero_modes.cpp)

src/mam4xx/mam4.hpp

@@ -22,6 +22,7 @@
 #include <mam4xx/hetfrz.hpp>
 #include <mam4xx/lin_strat_chem.hpp>
 #include <mam4xx/mam4_types.hpp>
+#include <mam4xx/mo_chm_diags.hpp>
 #include <mam4xx/mo_photo.hpp>
 #include <mam4xx/ndrop.hpp>
 #include <mam4xx/nucleate_ice.hpp>

src/mam4xx/mo_chm_diags.hpp

@@ -0,0 +1,69 @@
+#ifndef MAM4XX_MO_CHM_DIAGS_HPP
+#define MAM4XX_MO_CHM_DIAGS_HPP
+
+#include <haero/math.hpp>
+#include <mam4xx/gas_chem.hpp>
+#include <mam4xx/mam4_types.hpp>
+#include <mam4xx/utils.hpp>
+
+namespace mam4 {
+
+namespace mo_chm_diags {
+
+using Real = haero::Real;
+using View1D = DeviceType::view_1d<Real>;
+
+// FIXME: bad constants
+constexpr Real S_molwgt = 32.066;
+constexpr Real rearth = 6.37122e6;
+constexpr Real rgrav =
+    1.0 / 9.80616; // reciprocal of acceleration of gravity ~ m/s^2
+constexpr Real avogadro = haero::Constants::avogadro;
+constexpr const int gas_pcnst = gas_chemistry::gas_pcnst;
+constexpr const int pver = mam4::nlev; // number of vertical levels
+
+KOKKOS_INLINE_FUNCTION
+void het_diags(
+    const ThreadTeam &team,
+    const ColumnView het_rates[gas_pcnst], //[pver][gas_pcnst], //in
+    const ColumnView mmr[gas_pcnst],       //[pver][gas_pcnst],
+    const ColumnView &pdel,                //[pver],
+    const Real &wght,
+    View1D wrk_wd, //[gas_pcnst], //output
+    // Real noy_wk, //output //this isn't actually used in this function?
+    View1D sox_wk, // output
+    // Real nhx_wk, //output //this isn't actually used in this function?
+    const Real adv_mass[gas_pcnst], // constant from elsewhere
+    const Real sox_species[3]) {
+  // change to pass values for a single col
+  //===========
+  // output integrated wet deposition field
+  //===========
+
+  sox_wk[0] = 0;
+  Kokkos::parallel_for(
+      Kokkos::TeamThreadRange(team, gas_pcnst), KOKKOS_LAMBDA(int mm) {
+        //
+        // compute vertical integral
+        //
+        wrk_wd(mm) = 0;
+
+        for (int kk = 0; kk < pver; kk++) {
+          wrk_wd(mm) += het_rates[mm](kk) * mmr[mm](kk) *
+                        pdel(kk); // parallel_reduce in the future?
+        }
+
+        wrk_wd(mm) *= rgrav * wght * haero::square(rearth);
+      });
+
+  for (int mm = 0; mm < gas_pcnst; mm++) {
+    for (int i = 0; i < 3; i++) { // FIXME: bad constant (len of sox species)
+      if (sox_species[i] == mm)
+        sox_wk[0] += wrk_wd[mm] * S_molwgt / adv_mass[mm];
+    }
+  }
+} // het_diags
+
+} // namespace mo_chm_diags
+} // namespace mam4
+#endif

src/validation/CMakeLists.txt

@@ -30,3 +30,4 @@ add_subdirectory(ndrop)
 add_subdirectory(water_uptake)
 add_subdirectory(mo_photo)
 add_subdirectory(lin_strat_chem)
+add_subdirectory(mo_chm_diags)

src/validation/mo_chm_diags/CMakeLists.txt

@@ -0,0 +1,55 @@
+set(MO_CHM_DIAGS_VALIDATION_DIR ${MAM_X_VALIDATION_DIR}/mo_chm_diags)
+set(MO_CHM_DIAGS_VALIDATION_SCRIPTS_DIR ${MAM_X_VALIDATION_DIR}/scripts)
+
+
+# These subdirectories contain Skywalker drivers for MAM4 parameterizations.
+# Include directory for .mod files.
+
+include_directories(${PROJECT_BINARY_DIR}/validation)
+
+# We use a single driver for all mo_chm_diags-related parameterizations.
+add_executable(mo_chm_diags_driver
+               mo_chm_diags_driver.cpp
+               het_diags.cpp
+               )
+
+target_link_libraries(mo_chm_diags_driver skywalker;validation;${HAERO_LIBRARIES})
+
+# Copy some Python scripts from mam_x_validation to our binary directory.
+foreach(script
+        compare_mam4xx_mam4.py)
+  configure_file(
+    ${MO_CHM_DIAGS_VALIDATION_SCRIPTS_DIR}/${script}
+    ${CMAKE_CURRENT_BINARY_DIR}/${script}
+    COPYONLY
+  )
+endforeach()
+# stand_calc_sum_wght_ts_355
+# Run the driver in several configurations to produce datasets.
+set(TEST_LIST
+    het_diags_ts_355
+    )
+# # matching the tests and errors, just for convenience
+
+set(DEFAULT_TOL 1e-13)
+set(ERROR_THRESHOLDS
+   9e-2
+   )
+foreach(input tol IN ZIP_LISTS TEST_LIST ERROR_THRESHOLDS)
+  # copy the baseline file into place.
+
+  configure_file(
+    ${MO_CHM_DIAGS_VALIDATION_DIR}/mam_${input}.py
+    ${CMAKE_CURRENT_BINARY_DIR}/mam_${input}.py
+    COPYONLY
+  )
+
+  # add a test to run the skywalker driver
+  add_test(run_${input} mo_chm_diags_driver ${MO_CHM_DIAGS_VALIDATION_DIR}/${input}.yaml)
+
+  # add a test to validate mam4xx's results against the baseline.
+  # Select a threshold error slightly bigger than the largest relative error for the threshold error.
+  # compare_mam4xx_mam4.py <module1.py> <module2.py> <check_norms> <threshold error>
+  add_test(validate_${input} python3 compare_mam4xx_mam4.py mam4xx_${input}.py mam_${input}.py True ${tol})
+  set_tests_properties(validate_${input} PROPERTIES DEPENDS run_${input})
+endforeach()

src/validation/mo_chm_diags/het_diags.cpp

@@ -0,0 +1,109 @@
+// 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 <skywalker.hpp>
+#include <validation.hpp>
+
+using namespace skywalker;
+using namespace mam4;
+using namespace haero;
+using namespace mo_chm_diags;
+
+// constexpr const int gas_pcnst = gas_chemistry::gas_pcnst;
+
+void het_diags(Ensemble *ensemble) {
+  ensemble->process([=](const Input &input, Output &output) {
+    using View1DHost = typename HostType::view_1d<Real>;
+    using View1D = typename DeviceType::view_1d<Real>;
+    using ColumnView = haero::ColumnView;
+
+    const auto het_rates_in = input.get_array("het_rates");
+    const auto mmr_in = input.get_array("mmr");
+    const auto pdel_in = input.get_array("pdel");
+    const auto wght_in = input.get_array("wght");
+
+    Real wght = wght_in[0];
+
+    ColumnView het_rates[gas_pcnst];
+    ColumnView mmr[gas_pcnst];
+    View1DHost het_rates_host[gas_pcnst];
+    View1DHost mmr_host[gas_pcnst];
+
+    for (int mm = 0; mm < gas_pcnst; ++mm) {
+      het_rates[mm] = haero::testing::create_column_view(pver);
+      mmr[mm] = haero::testing::create_column_view(pver);
+
+      het_rates_host[mm] = View1DHost("het_rates_host", pver);
+      mmr_host[mm] = View1DHost("mmr_host", pver);
+    }
+
+    int count = 0;
+    // for (int kk = 0; kk < pver; ++kk) {
+    //   for (int mm = 0; mm < gas_pcnst; ++mm) {
+    for (int mm = 0; mm < gas_pcnst; ++mm) {
+      for (int kk = 0; kk < pver; ++kk) {
+        het_rates_host[mm](kk) = het_rates_in[count];
+        mmr_host[mm](kk) = mmr_in[count];
+        count++;
+      }
+    }
+
+    // transfer data to GPU.
+    for (int mm = 0; mm < gas_pcnst; ++mm) {
+      Kokkos::deep_copy(het_rates[mm], het_rates_host[mm]);
+      Kokkos::deep_copy(mmr[mm], mmr_host[mm]);
+    }
+
+    ColumnView pdel;
+    auto pdel_host =
+        View1DHost((Real *)pdel_in.data(), pver); // puts data into host
+    pdel = haero::testing::create_column_view(pver);
+    Kokkos::deep_copy(pdel, pdel_host);
+
+    std::vector<Real> vector0(gas_pcnst, 0);
+    std::vector<Real> single_vector0(1, 0);
+
+    auto wrk_wd_host = View1DHost(vector0.data(), gas_pcnst);
+    const auto wrk_wd = View1D("wrk_wd", gas_pcnst);
+    Kokkos::deep_copy(wrk_wd, wrk_wd_host);
+
+    auto sox_wk_host = View1DHost(single_vector0.data(), 1);
+    const auto sox_wk = View1D("sox_wk", 1);
+    Kokkos::deep_copy(sox_wk, sox_wk_host);
+
+    const Real sox_species[3] = {4, -1, 3};
+    const Real adv_mass[gas_pcnst] = {
+        47.998200,     34.013600,  98.078400,     64.064800, 62.132400,
+        12.011000,     115.107340, 12.011000,     12.011000, 12.011000,
+        135.064039,    58.442468,  250092.672000, 1.007400,  115.107340,
+        12.011000,     58.442468,  250092.672000, 1.007400,  135.064039,
+        58.442468,     115.107340, 12.011000,     12.011000, 12.011000,
+        250092.672000, 1.007400,   12.011000,     12.011000, 250092.672000,
+        1.007400};
+
+    auto team_policy = ThreadTeamPolicy(1u, Kokkos::AUTO);
+    Kokkos::parallel_for(
+        team_policy, KOKKOS_LAMBDA(const ThreadTeam &team) {
+          mo_chm_diags::het_diags(team, het_rates, mmr, pdel, wght, wrk_wd,
+                                  sox_wk, adv_mass, sox_species);
+        });
+
+    std::vector<Real> wrk_wd_mm(gas_pcnst);
+    Kokkos::deep_copy(wrk_wd_host, wrk_wd);
+    for (int mm = 0; mm < gas_pcnst; mm++) {
+      wrk_wd_mm[mm] = wrk_wd_host(mm);
+    }
+
+    Kokkos::deep_copy(sox_wk_host, sox_wk);
+
+    std::vector<Real> sox_wk_out(1);
+
+    output.set("wrk_wd_mm", wrk_wd_mm);
+    output.set("sox_wk", sox_wk_host[0]);
+  });
+}

src/validation/mo_chm_diags/mo_chm_diags_driver.cpp

@@ -0,0 +1,67 @@
+// mam4xx: Copyright (c) 2022,
+// Battelle Memorial Institute and
+// National Technology & Engineering Solutions of Sandia, LLC (NTESS)
+// SPDX-License-Identifier: BSD-3-Clause
+
+#include <mam4xx/mo_chm_diags.hpp>
+
+#include <iostream>
+#include <skywalker.hpp>
+#include <validation.hpp>
+
+void usage() {
+  std::cerr << "mo_chm_diags_driver: a Skywalker driver for validating the "
+               "MAM4 rename parameterizations."
+            << std::endl;
+  std::cerr << "mo_chm_diags_driver: usage:" << std::endl;
+  std::cerr << "mo_chm_diags_driver <input.yaml>" << std::endl;
+  exit(0);
+}
+
+using namespace skywalker;
+using namespace mam4;
+
+// Parameterizations used by the mo_chm_diags() process.
+void het_diags(Ensemble *ensemble);
+int main(int argc, char **argv) {
+  if (argc == 1) {
+    usage();
+  }
+  validation::initialize(argc, argv);
+  std::string input_file = argv[1];
+  std::string output_file = validation::output_name(input_file);
+  std::cout << argv[0] << ": reading " << input_file << std::endl;
+
+  // Load the ensemble. Any error encountered is fatal.
+  Ensemble *ensemble = skywalker::load_ensemble(input_file, "mam4xx");
+
+  // the settings.
+  Settings settings = ensemble->settings();
+  if (!settings.has("function")) {
+    std::cerr << "No function specified in mam4xx.settings!" << std::endl;
+    exit(1);
+  }
+
+  // Dispatch to the requested function.
+  auto func_name = settings.get("function");
+  try {
+    if (func_name == "het_diags") {
+      het_diags(ensemble);
+    } else {
+      std::cerr << "Error: Function name '" << func_name
+                << "' does not have an implemented test!" << std::endl;
+      exit(1);
+    }
+
+  } catch (std::exception &e) {
+    std::cerr << argv[0] << ": Error: " << e.what() << std::endl;
+  }
+
+  // Write out a Python module.
+  std::cout << argv[0] << ": writing " << output_file << std::endl;
+  ensemble->write(output_file);
+
+  // Clean up.
+  delete ensemble;
+  validation::finalize();
+}