From cd712b2e6113cc33192ee56eb01390b4cce5b6fd Mon Sep 17 00:00:00 2001 From: James Overfelt Date: Mon, 28 Aug 2023 11:40:58 -0600 Subject: [PATCH 1/3] Add wetdep compute_tendencies test. --- src/mam4xx/wet_dep.hpp | 51 ++--- src/validation/mam_x_validation | 2 +- src/validation/wetdep/CMakeLists.txt | 2 + src/validation/wetdep/compute_tendencies.cpp | 203 +++++++++++++++++++ src/validation/wetdep/wetdep_driver.cpp | 8 +- 5 files changed, 239 insertions(+), 27 deletions(-) create mode 100644 src/validation/wetdep/compute_tendencies.cpp diff --git a/src/mam4xx/wet_dep.hpp b/src/mam4xx/wet_dep.hpp index 46bb4affa..e963b69da 100644 --- a/src/mam4xx/wet_dep.hpp +++ b/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, @@ -1209,7 +1211,7 @@ class WetDeposition { Kokkos::View prain; Kokkos::View conicw; Kokkos::View totcond; - Kokkos::View> scratch; + Kokkos::View> scratch; Real scavimptblnum[aero_model::nimptblgrow_total][AeroConfig::num_modes()]; Real scavimptblvol[aero_model::nimptblgrow_total][AeroConfig::num_modes()]; }; @@ -1229,6 +1231,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]; @@ -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; @@ -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(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; diff --git a/src/validation/mam_x_validation b/src/validation/mam_x_validation index 914d59270..05cafe429 160000 --- a/src/validation/mam_x_validation +++ b/src/validation/mam_x_validation @@ -1 +1 @@ -Subproject commit 914d59270c1150689b7d23d8776d8e705ecc84c3 +Subproject commit 05cafe429fa67cb3911f07e7083b7b9728c0efb5 diff --git a/src/validation/wetdep/CMakeLists.txt b/src/validation/wetdep/CMakeLists.txt index 731b251f6..e17e8bb33 100644 --- a/src/validation/wetdep/CMakeLists.txt +++ b/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( diff --git a/src/validation/wetdep/compute_tendencies.cpp b/src/validation/wetdep/compute_tendencies.cpp new file mode 100644 index 000000000..a5421e392 --- /dev/null +++ b/src/validation/wetdep/compute_tendencies.cpp @@ -0,0 +1,203 @@ +// 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 +#include +#include +#include +#include +#include + +using namespace skywalker; +using namespace mam4; + +namespace { +void get_input(const Input &input, const std::string &name, const int size, + std::vector &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 &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> matrix(rows, std::vector(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 &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) { + // 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); + + int mmtoo_prevap_resusp[gas_pcnst]; + { + const int resusp[gas_pcnst] = {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, + 0, 0, 0, 0, 0, 31, 33, 34, 32, 29, + 30, 35, -2, 31, 34, 30, 35, -2, 29, 30, + 31, 32, 33, 34, 35, -2, 33, 32, 35, -2}; + for (int i = 0; i < gas_pcnst; ++i) + mmtoo_prevap_resusp[i] = resusp[i] - 1; + } + + 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 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); + }); +} diff --git a/src/validation/wetdep/wetdep_driver.cpp b/src/validation/wetdep/wetdep_driver.cpp index 29030fae3..27dff4d03 100644 --- a/src/validation/wetdep/wetdep_driver.cpp +++ b/src/validation/wetdep/wetdep_driver.cpp @@ -24,6 +24,7 @@ void test_calculate_cloudy_volume(std::unique_ptr &ensemble); void test_local_precip_production(std::unique_ptr &ensemble); void test_wetdep_resusp(std::unique_ptr &ensemble); void test_wetdepa_v2(std::unique_ptr &ensemble); +void test_compute_tendencies(std::unique_ptr &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; From fe90c41426076c97e49f5450a56ac382e0ba8324 Mon Sep 17 00:00:00 2001 From: James Overfelt Date: Mon, 28 Aug 2023 12:13:08 -0600 Subject: [PATCH 2/3] Cuda fixes. --- src/mam4xx/wet_dep.hpp | 7 ++++--- src/validation/wetdep/compute_tendencies.cpp | 10 ---------- 2 files changed, 4 insertions(+), 13 deletions(-) diff --git a/src/mam4xx/wet_dep.hpp b/src/mam4xx/wet_dep.hpp index e963b69da..8e64b6641 100644 --- a/src/mam4xx/wet_dep.hpp +++ b/src/mam4xx/wet_dep.hpp @@ -1179,8 +1179,6 @@ class WetDeposition { // 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"; } @@ -1264,6 +1262,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 @@ -1372,7 +1373,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 diff --git a/src/validation/wetdep/compute_tendencies.cpp b/src/validation/wetdep/compute_tendencies.cpp index a5421e392..856eb8643 100644 --- a/src/validation/wetdep/compute_tendencies.cpp +++ b/src/validation/wetdep/compute_tendencies.cpp @@ -70,16 +70,6 @@ void test_compute_tendencies(std::unique_ptr &ensemble) { const int gas_pcnst = 40; EKAT_ASSERT(input.get("dt") == 3600); - int mmtoo_prevap_resusp[gas_pcnst]; - { - const int resusp[gas_pcnst] = {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, - 0, 0, 0, 0, 0, 31, 33, 34, 32, 29, - 30, 35, -2, 31, 34, 30, 35, -2, 29, 30, - 31, 32, 33, 34, 35, -2, 33, 32, 35, -2}; - for (int i = 0; i < gas_pcnst; ++i) - mmtoo_prevap_resusp[i] = resusp[i] - 1; - } - Atmosphere atmosphere = validation::create_atmosphere(nlev, pblh); Surface surface = validation::create_surface(); mam4::Prognostics prognostics = validation::create_prognostics(nlev); From 716c975f238caae74c1908f65da46fab7b15bbd5 Mon Sep 17 00:00:00 2001 From: James Overfelt Date: Mon, 28 Aug 2023 12:13:50 -0600 Subject: [PATCH 3/3] Clang format. --- src/mam4xx/wet_dep.hpp | 1 - 1 file changed, 1 deletion(-) diff --git a/src/mam4xx/wet_dep.hpp b/src/mam4xx/wet_dep.hpp index 8e64b6641..af1752fad 100644 --- a/src/mam4xx/wet_dep.hpp +++ b/src/mam4xx/wet_dep.hpp @@ -1178,7 +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; - }; const char *name() const { return "MAM4 Wet Deposition"; }