Add wetdep compute_tendencies test.

src/mam4xx/wet_dep.hpp

@@ -976,8 +976,8 @@ void wetdepa_v2(const Real deltat, const Real pdel, const Real cmfdqr,
       haero::max(small_value_12, cldc * conicw + (cmfdqr + dlf) * deltat);
   fracp = utils::min_max_bound(0.0, 1.0, fracp) * cldc;
 
-  Real srcc; // tendency for convective rain scavenging [kg/kg/s]
-  Real finc; // fraction of rem. rate by conv. rain [fraction]
+  Real srcc = 0; // tendency for convective rain scavenging [kg/kg/s]
+  Real finc = 0; // fraction of rem. rate by conv. rain [fraction]
   // 2 is for convective:
   wetdep_scavenging(2, is_strat_cloudborne, deltat, fracp, precabc, cldvcu,
                     scavcoef, sol_factb, sol_factic, tracer_incu, tracer_mean,
@@ -1053,16 +1053,18 @@ void wetdepa_v2(const Real deltat, const Real pdel, const Real cmfdqr,
   }
 
   // ****************** update scavengingfor output ***************
-  Real scavt_ik;   // scavenging tend at current  [kg/kg/s]
-  Real iscavt_ik;  // incloud scavenging tends at current  [kg/kg/s]
-  Real icscavt_ik; // incloud, convective scavenging tends at current  [kg/kg/s]
-  Real isscavt_ik; // incloud, stratiform scavenging tends at current  [kg/kg/s]
-  Real bcscavt_ik; // below cloud, convective scavenging tends at current
-                   // [kg/kg/s]
-  Real bsscavt_ik; // below cloud, stratiform scavenging tends at current
-                   // [kg/kg/s]
-  Real rcscavt_ik; // resuspension, convective tends at current  [kg/kg/s]
-  Real rsscavt_ik; // resuspension, stratiform tends at current  [kg/kg/s]
+  Real scavt_ik = 0;  // scavenging tend at current  [kg/kg/s]
+  Real iscavt_ik = 0; // incloud scavenging tends at current  [kg/kg/s]
+  Real icscavt_ik =
+      0; // incloud, convective scavenging tends at current  [kg/kg/s]
+  Real isscavt_ik =
+      0; // incloud, stratiform scavenging tends at current  [kg/kg/s]
+  Real bcscavt_ik = 0; // below cloud, convective scavenging tends at current
+                       // [kg/kg/s]
+  Real bsscavt_ik = 0; // below cloud, stratiform scavenging tends at current
+                       // [kg/kg/s]
+  Real rcscavt_ik = 0; // resuspension, convective tends at current  [kg/kg/s]
+  Real rsscavt_ik = 0; // resuspension, stratiform tends at current  [kg/kg/s]
   update_scavenging(mam_prevap_resusp_optcc, pdel, omsm, srcc, srcs, srct, fins,
                     finc, fracev_st, fracev_cu, resusp_c, resusp_s, precs,
                     evaps, cmfdqr, evapc, scavt_ik, iscavt_ik, icscavt_ik,
@@ -1176,9 +1178,6 @@ class WetDeposition {
     //   raindrop number the 130 thru 230 all use the new prevap_resusp code
     //   block in subr wetdepa_v2
     int mam_prevap_resusp_optcc = 0;
-
-    // change mode order as mmode_loop_aa loops in a different order
-    static constexpr int mode_order_change[4] = {0, 1, 3, 2};
   };
 
   const char *name() const { return "MAM4 Wet Deposition"; }
@@ -1209,7 +1208,7 @@ class WetDeposition {
   Kokkos::View<Real *> prain;
   Kokkos::View<Real *> conicw;
   Kokkos::View<Real *> totcond;
-  Kokkos::View<Real[1], Kokkos::MemoryTraits<Kokkos::Atomic>> scratch;
+  Kokkos::View<Real *, Kokkos::MemoryTraits<Kokkos::Atomic>> scratch;
   Real scavimptblnum[aero_model::nimptblgrow_total][AeroConfig::num_modes()];
   Real scavimptblvol[aero_model::nimptblgrow_total][AeroConfig::num_modes()];
 };
@@ -1229,6 +1228,7 @@ inline void WetDeposition::init(const AeroConfig &aero_config,
   Kokkos::resize(prain, nlev);
   Kokkos::resize(conicw, nlev);
   Kokkos::resize(totcond, nlev);
+  Kokkos::resize(scratch, 1);
 
   const int num_modes = AeroConfig::num_modes();
   Real dgnum_amode[num_modes];
@@ -1261,6 +1261,9 @@ void WetDeposition::compute_tendencies(
   const int nlev = config_.nlev;
   static constexpr int gas_pcnst = 40;
 
+  // change mode order as mmode_loop_aa loops in a different order
+  static constexpr int mode_order_change[4] = {0, 1, 3, 2};
+
   //   0 = no resuspension
   // 130 = non-linear resuspension of aerosol mass   based on scavenged aerosol
   // mass 230 = non-linear resuspension of aerosol number based on raindrop
@@ -1346,7 +1349,7 @@ void WetDeposition::compute_tendencies(
   // calc_sfc_flux does.
   // TODO: Determine a better way to implement calc_sfc_flux without
   // team_fence() to avoid locking if the team size is not nlev;
-  EKAT_KERNEL_REQUIRE(team.team_size() == nlev);
+  EKAT_KERNEL_REQUIRE(team.team_size() == nlev || team.team_size() == 1);
   Kokkos::parallel_for(
       Kokkos::TeamThreadRange(team, nlev), KOKKOS_CLASS_LAMBDA(int k) {
         aerdepwetis[k] = 0;
@@ -1369,7 +1372,7 @@ void WetDeposition::compute_tendencies(
           // then aitken = 1,
           // then pcarbon - 3,
           // then coarse = 2
-          const int imode = Config::mode_order_change[mtmp];
+          const int imode = mode_order_change[mtmp];
 
           // loop over interstitial (1) and cloud-borne (2) forms
           // BSINGH (09/12/2014):Do cloudborne first for unified convection
@@ -1439,7 +1442,6 @@ void WetDeposition::compute_tendencies(
                       f_act_conv = f_act_conv_coarse_nacl;
                   }
                 }
-
                 if (lphase == 1) {
                   // traces reflects changes from modal_aero_calcsize and is the
                   // "most current" q
@@ -1458,15 +1460,16 @@ void WetDeposition::compute_tendencies(
                   const int k_p1 =
                       static_cast<int>(haero::min(k + 1, nlev - 1));
 
-                  Real fracis;  // fraction of species not scavenged [fraction]
-                  Real scavt;   // scavenging tend [kg/kg/s]
-                  Real iscavt;  // incloud scavenging tends [kg/kg/s]
-                  Real icscavt; // incloud, convective [kg/kg/s]
-                  Real isscavt; // incloud, stratiform [kg/kg/s]
-                  Real bcscavt; // below cloud, convective [kg/kg/s]
-                  Real bsscavt; // below cloud, stratiform [kg/kg/s]
-                  Real rcscavt; // resuspension, convective [kg/kg/s]
-                  Real rsscavt; // resuspension, stratiform [kg/kg/s]
+                  Real fracis =
+                      0; // fraction of species not scavenged [fraction]
+                  Real scavt = 0;   // scavenging tend [kg/kg/s]
+                  Real iscavt = 0;  // incloud scavenging tends [kg/kg/s]
+                  Real icscavt = 0; // incloud, convective [kg/kg/s]
+                  Real isscavt = 0; // incloud, stratiform [kg/kg/s]
+                  Real bcscavt = 0; // below cloud, convective [kg/kg/s]
+                  Real bsscavt = 0; // below cloud, stratiform [kg/kg/s]
+                  Real rcscavt = 0; // resuspension, convective [kg/kg/s]
+                  Real rsscavt = 0; // resuspension, stratiform [kg/kg/s]
                   // is_strat_cloudborne = true if tracer is
                   // stratiform-cloudborne aerosol; else false
                   const bool is_strat_cloudborne = false;

src/validation/wetdep/CMakeLists.txt

@@ -18,6 +18,7 @@ add_executable(wetdep_driver
                local_precip_production.cpp
                wetdep_resusp.cpp
                wetdepa_v2.cpp
+               compute_tendencies.cpp
               )
 target_link_libraries(wetdep_driver skywalker;validation;${HAERO_LIBRARIES})
 
@@ -43,6 +44,7 @@ foreach (TEST_NAME
 	 wetdepa_v2_0
 	 wetdepa_v2_130
 	 wetdepa_v2_230
+	 compute_tendencies
 	)
   set(SKYWALKER_RESULT ${TEST_NAME}_output_ts_355)
   configure_file(

src/validation/wetdep/compute_tendencies.cpp

@@ -0,0 +1,193 @@
+// mam4xx: Copyright (c) 2022,
+// Battelle Memorial Institute and
+// National Technology & Engineering Solutions of Sandia, LLC (NTESS)
+// SPDX-License-Identifier: BSD-3-Clause
+
+#include "Kokkos_Core.hpp"
+#include <catch2/catch.hpp>
+#include <iomanip>
+#include <iostream>
+#include <mam4xx/wet_dep.hpp>
+#include <skywalker.hpp>
+#include <validation.hpp>
+
+using namespace skywalker;
+using namespace mam4;
+
+namespace {
+void get_input(const Input &input, const std::string &name, const int size,
+               std::vector<Real> &host, ColumnView &dev) {
+  host = input.get_array(name);
+  dev = mam4::validation::create_column_view(size);
+
+  EKAT_ASSERT(host.size() == size);
+  auto host_view = Kokkos::create_mirror_view(dev);
+  for (int n = 0; n < size; ++n)
+    host_view[n] = host[n];
+  Kokkos::deep_copy(dev, host_view);
+}
+void get_input(const Input &input, const std::string &name, const int rows,
+               const int cols, std::vector<Real> &host,
+               Diagnostics::ColumnTracerView &dev) {
+  host = input.get_array(name);
+  ColumnView col_view = mam4::validation::create_column_view(rows * cols);
+  dev = Diagnostics::ColumnTracerView(col_view.data(), rows, cols);
+  EKAT_ASSERT(host.size() == rows * cols);
+  {
+    std::vector<std::vector<Real>> matrix(rows, std::vector<Real>(cols));
+    for (int j = 0, n = 0; j < cols; ++j)
+      for (int i = 0; i < rows; ++i, ++n)
+        matrix[i][j] = host[n];
+    auto host_view = Kokkos::create_mirror_view(dev);
+    for (int i = 0; i < rows; ++i)
+      for (int j = 0; j < cols; ++j)
+        host_view(i, j) = matrix[i][j];
+    Kokkos::deep_copy(dev, host_view);
+  }
+}
+void set_output(Output &output, const std::string &name, const int rows,
+                const int cols, std::vector<Real> &host,
+                const Diagnostics::ColumnTracerView &dev) {
+  host.resize(rows * cols);
+  auto host_view = Kokkos::create_mirror_view(dev);
+  Kokkos::deep_copy(host_view, dev);
+  for (int j = 0, n = 0; j < cols; ++j)
+    for (int i = 0; i < rows; ++i, ++n) {
+      host[n] = host_view(i, j);
+    }
+  output.set(name, host);
+}
+} // namespace
+void test_compute_tendencies(std::unique_ptr<Ensemble> &ensemble) {
+  // We don't need any settings for this particular test.
+  // Settings settings = ensemble->settings();
+  // Run the ensemble.
+  ensemble->process([=](const Input &input, Output &output) {
+    const int nlev = 72;
+    const Real t = 0;
+    const Real dt = 36000;
+    const Real pblh = 1000;
+    const int gas_pcnst = 40;
+    EKAT_ASSERT(input.get("dt") == 3600);
+
+    Atmosphere atmosphere = validation::create_atmosphere(nlev, pblh);
+    Surface surface = validation::create_surface();
+    mam4::Prognostics prognostics = validation::create_prognostics(nlev);
+    mam4::Diagnostics diagnostics = validation::create_diagnostics(nlev);
+    mam4::Tendencies tendencies = validation::create_tendencies(nlev);
+    mam4::AeroConfig aero_config;
+    mam4::WetDeposition::Config wetdep_config;
+
+    mam4::WetDeposition wetdep;
+    wetdep.init(aero_config, wetdep_config);
+
+    diagnostics.deep_convective_cloud_fraction =
+        mam4::validation::create_column_view(nlev);
+    diagnostics.shallow_convective_cloud_fraction =
+        mam4::validation::create_column_view(nlev);
+    diagnostics.deep_convective_cloud_condensate =
+        mam4::validation::create_column_view(nlev);
+    diagnostics.shallow_convective_cloud_condensate =
+        mam4::validation::create_column_view(nlev);
+    diagnostics.deep_convective_precipitation_production =
+        mam4::validation::create_column_view(nlev);
+    diagnostics.shallow_convective_precipitation_production =
+        mam4::validation::create_column_view(nlev);
+    diagnostics.deep_convective_precipitation_evaporation =
+        mam4::validation::create_column_view(nlev);
+    diagnostics.shallow_convective_precipitation_evaporation =
+        mam4::validation::create_column_view(nlev);
+    diagnostics.total_convective_detrainment =
+        mam4::validation::create_column_view(nlev);
+    diagnostics.shallow_convective_detrainment =
+        mam4::validation::create_column_view(nlev);
+    diagnostics.shallow_convective_ratio =
+        mam4::validation::create_column_view(nlev);
+    diagnostics.evaporation_of_falling_precipitation =
+        mam4::validation::create_column_view(nlev);
+    diagnostics.aerosol_wet_deposition_interstitial =
+        mam4::validation::create_column_view(nlev);
+    diagnostics.aerosol_wet_deposition_cloud_water =
+        mam4::validation::create_column_view(nlev);
+    for (int i = 0; i < AeroConfig::num_modes(); ++i)
+      diagnostics.wet_geometric_mean_diameter_i[i] =
+          mam4::validation::create_column_view(nlev);
+    auto mixing_ratio = mam4::validation::create_column_view(nlev * gas_pcnst);
+    diagnostics.tracer_mixing_ratio =
+        Diagnostics::ColumnTracerView(mixing_ratio.data(), nlev, gas_pcnst);
+    auto mixing_ratio_dt =
+        mam4::validation::create_column_view(nlev * gas_pcnst);
+    diagnostics.d_tracer_mixing_ratio_dt =
+        Diagnostics::ColumnTracerView(mixing_ratio_dt.data(), nlev, gas_pcnst);
+    Kokkos::parallel_for(
+        "init_column_views", nlev, KOKKOS_LAMBDA(int i) {
+          diagnostics.deep_convective_cloud_fraction[i] = 0;
+          diagnostics.shallow_convective_cloud_fraction[i] = 0;
+          diagnostics.deep_convective_cloud_condensate[i] = 0;
+          diagnostics.shallow_convective_cloud_condensate[i] = 0;
+          diagnostics.deep_convective_precipitation_production[i] = 0;
+          diagnostics.shallow_convective_precipitation_production[i] = 0;
+          diagnostics.deep_convective_precipitation_evaporation[i] = 0;
+          diagnostics.shallow_convective_precipitation_evaporation[i] = 0;
+          diagnostics.total_convective_detrainment[i] = 0;
+          diagnostics.shallow_convective_detrainment[i] = 0;
+          diagnostics.shallow_convective_ratio[i] = 0;
+          diagnostics.evaporation_of_falling_precipitation[i] = 0;
+          diagnostics.aerosol_wet_deposition_interstitial[i] = 0;
+          diagnostics.aerosol_wet_deposition_cloud_water[i] = 0;
+          for (int j = 0; j < gas_pcnst; ++j) {
+            diagnostics.tracer_mixing_ratio(i, j) = 0;
+            diagnostics.d_tracer_mixing_ratio_dt(i, j) = 0;
+          }
+        });
+    Kokkos::fence();
+    std::vector<Real> temperature_host, pdel_host, pmid_host, cldfrac_host,
+        icwmr_host, rprd_host, evapc_host, dqdt_host, dp_host, dgn_awet_host,
+        evapr_host;
+    ColumnView temperature_dev, pmid_dev, pdel_dev;
+
+    get_input(input, "state_t", nlev, temperature_host, temperature_dev);
+    atmosphere.temperature = temperature_dev;
+    get_input(input, "state_pmid", nlev, pmid_host, pmid_dev);
+    atmosphere.pressure = pmid_dev;
+    get_input(input, "state_pdel", nlev, pdel_host, pdel_dev);
+    atmosphere.hydrostatic_dp = pdel_dev;
+    get_input(input, "dp_frac", nlev, cldfrac_host,
+              diagnostics.deep_convective_cloud_fraction);
+    get_input(input, "cldfrac", nlev, cldfrac_host,
+              diagnostics.total_convective_detrainment);
+    get_input(input, "icwmrdp", nlev, icwmr_host,
+              diagnostics.deep_convective_cloud_condensate);
+    get_input(input, "rprddp", nlev, rprd_host,
+              diagnostics.deep_convective_precipitation_production);
+    get_input(input, "evapcdp", nlev, evapc_host,
+              diagnostics.deep_convective_precipitation_evaporation);
+    get_input(input, "qnew", nlev, gas_pcnst, dp_host,
+              diagnostics.tracer_mixing_ratio);
+    get_input(input, "evapr", nlev, evapr_host,
+              diagnostics.evaporation_of_falling_precipitation);
+    ColumnView awet_tmp;
+    get_input(input, "dgn_awet", AeroConfig::num_modes(), dgn_awet_host,
+              awet_tmp);
+    for (int i = 0; i < AeroConfig::num_modes(); ++i) {
+      auto host_view = Kokkos::create_mirror_view(
+          diagnostics.wet_geometric_mean_diameter_i[i]);
+      for (int j = 0; j < nlev; ++j)
+        host_view[j] = dgn_awet_host[i];
+      Kokkos::deep_copy(diagnostics.wet_geometric_mean_diameter_i[i],
+                        host_view);
+    }
+    // NOTE: we haven't parallelized wetdep over vertical levels because of
+    // NOTE: data dependencies, so we run this serially
+    auto team_policy = haero::ThreadTeamPolicy(1u, 1u);
+    Kokkos::parallel_for(
+        team_policy, KOKKOS_LAMBDA(const ThreadTeam &team) {
+          wetdep.compute_tendencies(aero_config, team, t, dt, atmosphere,
+                                    surface, prognostics, diagnostics,
+                                    tendencies);
+        });
+    Kokkos::fence();
+    set_output(output, "dqdt", nlev, gas_pcnst, dqdt_host,
+               diagnostics.d_tracer_mixing_ratio_dt);
+  });
+}

src/validation/wetdep/wetdep_driver.cpp

@@ -24,6 +24,7 @@ void test_calculate_cloudy_volume(std::unique_ptr<Ensemble> &ensemble);
 void test_local_precip_production(std::unique_ptr<Ensemble> &ensemble);
 void test_wetdep_resusp(std::unique_ptr<Ensemble> &ensemble);
 void test_wetdepa_v2(std::unique_ptr<Ensemble> &ensemble);
+void test_compute_tendencies(std::unique_ptr<Ensemble> &ensemble);
 
 void usage(const std::string &prog_name) {
   std::cerr << prog_name << ": usage:" << std::endl;
@@ -55,14 +56,15 @@ int main(int argc, char **argv) {
       name != "wetdep_resusp_noprecip" && name != "wetdep_scavenging" &&
       name != "compute_evap_frac" && name != "rain_mix_ratio" &&
       name != "calculate_cloudy_volume" && name != "local_precip_production" &&
-      name != "wetdep_resusp" && name != "wetdepa_v2") {
+      name != "wetdep_resusp" && name != "wetdepa_v2" &&
+      name != "compute_tendencies") {
     std::cerr << "Invalid name: " << name << std::endl;
     std::cerr << "Currently the only valid name are: "
               << "wetdep_clddiag, update_scavenging, wetdep_prevap, "
                  "wetdep_resusp_nonlinear, wetdep_resusp_noprecip, "
                  "wetdep_scavenging, compute_evap_frac, rain_mix_ratio, "
                  "calculate_cloudy_volume, local_precip_production, "
-                 "wetdep_resusp, wetdepa_v2."
+                 "wetdep_resusp, wetdepa_v2 compute_tendencies."
               << std::endl;
     exit(1);
   }
@@ -92,6 +94,8 @@ int main(int argc, char **argv) {
       test_wetdep_resusp(ensemble);
     } else if (name == "wetdepa_v2") {
       test_wetdepa_v2(ensemble);
+    } else if (name == "compute_tendencies") {
+      test_compute_tendencies(ensemble);
     }
     // Write out a Python module.
     std::cout << argv[0] << ": writing " << output_file << std::endl;