Merge branch 'bartgol/eamxx/field-utils' into next (PR #7857)

This design pattern uses a type-agnostic wrapper for scalars,
so that we don't need to expose the templated method implementation
to the customers.

[non-BFB] For runs using the IC perturb feature

Author: bartgol

components/eamxx/src/control/atmosphere_driver.cpp

@@ -1408,13 +1408,11 @@ void AtmosphereDriver::set_initial_conditions ()
       seed = ic_pl.get<int>("perturbation_random_seed", 0);
     }
     m_atm_logger->info("      For IC perturbation, random seed: "+std::to_string(seed));
-    std::mt19937_64 engine(seed);
 
     // Get perturbation limit. Defines a range [1-perturbation_limit, 1+perturbation_limit]
     // for which the perturbation value will be randomly generated from. Create a uniform
     // distribution for this range.
     const auto perturbation_limit = ic_pl.get<Real>("perturbation_limit", 0.001);
-    std::uniform_real_distribution<Real> pdf(1-perturbation_limit, 1+perturbation_limit);
 
     // Define a level mask using reference pressure and the perturbation_minimum_pressure parameter.
     // This mask dictates which levels we apply a perturbation.
@@ -1423,13 +1421,21 @@ void AtmosphereDriver::set_initial_conditions ()
     const auto hybm_h = gll_grid->get_geometry_data("hybm").get_view<const Real*, Host>();
     constexpr auto ps0 = physics::Constants<Real>::P0;
     const auto min_pressure = ic_pl.get<Real>("perturbation_minimum_pressure", 1050.0);
-    auto pressure_mask = [&] (const int ilev) {
+
+    const auto& pmask_lt = gll_grid->get_vertical_layout(true);
+    const auto nondim = ekat::units::Units::nondimensional();
+    FieldIdentifier pmask_fid("lev_mask",pmask_lt,nondim,gll_grid->name(),DataType::IntType);
+    Field pressure_mask(pmask_fid,true);
+    auto pmask_h = pressure_mask.get_view<int*,Host>();
+    for (int ilev=0; ilev<pmask_lt.dim(0); ++ilev) {
       const auto pref = (hyam_h(ilev)*ps0 + hybm_h(ilev)*ps0)/100; // Reference pressure ps0 is in Pa, convert to millibar
-      return pref > min_pressure;
-    };
+      pmask_h(ilev) = static_cast<int>(pref > min_pressure);
+    }
+    pressure_mask.sync_to_dev();
 
     // Loop through fields and apply perturbation.
     const std::string gll_grid_name = gll_grid->name();
+    auto dofs_gids = gll_grid->get_dofs_gids();
     for (size_t f=0; f<perturbed_fields.size(); ++f) {
       const auto fname = perturbed_fields[f];
       EKAT_REQUIRE_MSG(ekat::contains(m_fields_inited[gll_grid_name], fname),
@@ -1438,7 +1444,7 @@ void AtmosphereDriver::set_initial_conditions ()
                        "  - Grid:  "+gll_grid_name+"\n");
 
       auto field = m_field_mgr->get_field(fname, gll_grid_name);
-      perturb(field, engine, pdf, seed, pressure_mask, gll_grid->get_dofs_gids());
+      perturb(field, perturbation_limit, seed, pressure_mask, dofs_gids);
     }
 
     m_atm_logger->info("    [EAMxx] Adding random perturbation to ICs ... done!");

components/eamxx/src/dynamics/homme/tests/dyn_grid_io.cpp

@@ -116,14 +116,13 @@ TEST_CASE("dyn_grid_io")
   std::vector<std::string> fnames = {"field_1", "field_2", "field_3"};
 
   // Randomize dyn fields, then remap to ctrl fields
-  std::uniform_real_distribution<Real> pdf(0.01,100.0);
-  auto engine = setup_random_test(&comm);
+  auto seed = get_random_test_seed(&comm);
   auto dyn2ctrl = gm->create_remapper(dyn_grid,phys_grid);
   for (const auto& fn : fnames) {
     auto fd = fm->get_field(fn,dyn_grid->name());
     auto fc = fm_ctrl->get_field(fn,phys_grid->name());
     dyn2ctrl->register_field(fd,fc);
-    randomize(fd,engine,pdf);
+    randomize_uniform(fd,seed++);
 
     // Init phys field to something obviously wrong
     auto fp = fm->get_field(fn,phys_grid->name());

components/eamxx/src/physics/cosp/eamxx_cosp.cpp

@@ -6,6 +6,7 @@
 #include "share/property_checks/field_within_interval_check.hpp"
 #include "share/field/field_utils.hpp"
 
+#include <ekat_team_policy_utils.hpp>
 #include <ekat_assert.hpp>
 #include <ekat_units.hpp>
 

components/eamxx/src/physics/iop_forcing/eamxx_iop_forcing_process_interface.cpp

@@ -461,17 +461,17 @@ void IOPForcing::run_impl (const double dt)
     view_1d<Pack> qv_mean, t_mean;
     view_2d<Pack> horiz_winds_mean;
     if (iop_nudge_tq){
-      horiz_contraction<Real>(m_helper_fields.at("qv_mean"), get_field_out("qv"),
-                              m_helper_fields.at("horiz_mean_weights"), &m_comm);
+      horiz_contraction(m_helper_fields.at("qv_mean"), get_field_out("qv"),
+                        m_helper_fields.at("horiz_mean_weights"), true, &m_comm);
       qv_mean = m_helper_fields.at("qv_mean").get_view<Pack*>();
 
-      horiz_contraction<Real>(m_helper_fields.at("t_mean"), get_field_out("T_mid"),
-                              m_helper_fields.at("horiz_mean_weights"), &m_comm);
+      horiz_contraction(m_helper_fields.at("t_mean"), get_field_out("T_mid"),
+                        m_helper_fields.at("horiz_mean_weights"), true, &m_comm);
       t_mean = m_helper_fields.at("t_mean").get_view<Pack*>();
     }
     if (iop_nudge_uv){
-      horiz_contraction<Real>(m_helper_fields.at("horiz_winds_mean"), get_field_out("horiz_winds"),
-                              m_helper_fields.at("horiz_mean_weights"), &m_comm);
+      horiz_contraction(m_helper_fields.at("horiz_winds_mean"), get_field_out("horiz_winds"),
+                        m_helper_fields.at("horiz_mean_weights"), true, &m_comm);
       horiz_winds_mean = m_helper_fields.at("horiz_winds_mean").get_view<Pack**>();
     }
 

components/eamxx/src/share/algorithm/tests/data_interpolation_tests.cpp

@@ -215,12 +215,8 @@ void run_tests (const std::shared_ptr<const AbstractGrid>& grid,
       // Compute the rel error
       diff[i].update(fields[i],1,1);
       diff[i].update(expected[i],-1,1);
-      diff[i].scale(1.0 / frobenius_norm<Real>(expected[i]));
-      if (frobenius_norm<Real>(diff[i])>=tol) {
-        print_field_hyperslab(expected[i]);
-        print_field_hyperslab(fields[i]);
-      }
-      REQUIRE (frobenius_norm<Real>(diff[i])<tol);
+      diff[i].scale(1.0 / frobenius_norm(expected[i]).as<Real>());
+      REQUIRE (frobenius_norm(diff[i]).as<Real>()<tol);
     }
   }
 }

components/eamxx/src/share/algorithm/tests/time_interpolation_tests.cpp

@@ -242,8 +242,8 @@ bool views_are_approx_equal(const Field& f0, const Field& f1, const Real tol)
   // simply.
   auto ft = f0.clone();
   ft.update(f1,1.0,-1.0);
-  auto d_min = field_min<Real>(ft);
-  auto d_max = field_max<Real>(ft);
+  auto d_min = field_min(ft).as<Real>();
+  auto d_max = field_max(ft).as<Real>();
   if (std::abs(d_min) > tol or std::abs(d_max) > tol) {
     printf("The two copies of (%16s) are NOT approx equal within a tolerance of %e.\n     The min and max errors are %e and %e respectively.\n",f0.name().c_str(),tol,d_min,d_max);
     return false;
@@ -315,21 +315,13 @@ std::shared_ptr<const GridsManager> get_gm (const ekat::Comm& comm, const int nc
 }
 /*-----------------------------------------------------------------------------------------------*/
 /* Create a fields manager for the test */
-std::shared_ptr<FieldManager> get_fm (const std::shared_ptr<const AbstractGrid>& grid, const util::TimeStamp& t0, const int seed)
+std::shared_ptr<FieldManager>
+get_fm (const std::shared_ptr<const AbstractGrid>& grid, const util::TimeStamp& t0, int seed)
 {
   using FL  = FieldLayout;
   using FID = FieldIdentifier;
   using namespace ShortFieldTagsNames;
 
-  // Random number generation stuff
-  // NOTES
-  //  - Use integers, so we can check answers without risk of
-  //    non bfb diffs due to different order of sums.
-  //  - Uniform_int_distribution returns an int, and the randomize
-  //    util checks that return type matches the Field data type.
-  //    So wrap the int pdf in a lambda, that does the cast.
-  std::mt19937_64 engine(seed);
-
   const int nlcols = grid->get_num_local_dofs();
   const int nlevs  = grid->get_num_vertical_levels();
 
@@ -343,13 +335,16 @@ std::shared_ptr<FieldManager> get_fm (const std::shared_ptr<const AbstractGrid>&
   auto fm = std::make_shared<FieldManager>(grid,RepoState::Closed);
 
   const auto units = ekat::units::Units::nondimensional();
+  std::vector<Real> values;
+  for (int i=1; i<=100; ++i)
+    values.push_back(static_cast<Real>(i));
   for (const auto& fl : layouts) {
     int gl_size = fl.size();
     grid->get_comm().all_reduce(&gl_size,1,MPI_SUM);
     FID fid("f_"+std::to_string(gl_size),fl,units,grid->name());
     Field f(fid);
     f.allocate_view();
-    randomize (f,engine,my_pdf);
+    randomize_discrete (f,seed++,values);
     f.get_header().get_tracking().update_time_stamp(t0);
     fm->add_field(f);
   }
@@ -364,7 +359,7 @@ std::vector<std::string> create_test_data_files(
 		const ekat::Comm& comm,
 		const std::shared_ptr<const GridsManager>& gm,
 		const util::TimeStamp& t0,
-	       	const int seed)
+	  const int seed)
 {
   // We initialize a local field manager to use for output
   auto fm = get_fm(gm->get_grid("point_grid"), t0, seed);

components/eamxx/src/share/data_managers/tests/field_mgr_tests.cpp

@@ -257,12 +257,10 @@ TEST_CASE("tracers_group", "") {
   REQUIRE ((0<=sub_idx_c1 and sub_idx_c1<=1));
 
   // Now fill tracer field with random values
-  auto engine = setup_random_test(&comm);
-  using RPDF = std::uniform_real_distribution<Real>;
-  RPDF pdf(0.0,1.0);
+  int seed = get_random_test_seed(&comm);
 
-  randomize(T1,engine,pdf);
-  randomize(T2,engine,pdf);
+  randomize_uniform(T1,seed++);
+  randomize_uniform(T2,seed++);
 
   // Check that the same values are in all individual tracers
   T1.sync_to_host();
@@ -300,7 +298,7 @@ TEST_CASE("tracers_group", "") {
   // Check that changing sub tracers change tracers group and individual tracers
   T1.deep_copy(0.0);
   auto& sub_T1 = subtracers.m_monolithic_field;
-  randomize(*sub_T1,engine,pdf);
+  randomize_uniform(*sub_T1,seed++);
   sub_T1->sync_to_host();
   auto sub_T1_h = sub_T1->get_strided_view<Real***,Host>();
 

components/eamxx/src/share/diagnostics/horiz_avg.cpp

@@ -51,7 +51,7 @@ void HorizAvgDiag::initialize_impl(const RunType /*run_type*/) {
   m_diagnostic_output.allocate_view();
 
   // scale the area field
-  auto total_area = field_sum<Real>(m_scaled_area, &m_comm);
+  auto total_area = field_sum(m_scaled_area, &m_comm).as<Real>();
   m_scaled_area.scale(sp(1.0) / total_area);
 
   if (f.get_header().has_extra_data("mask_data")) {
@@ -65,11 +65,7 @@ void HorizAvgDiag::compute_diagnostic_impl() {
   const auto &f = get_fields_in().front();
   const auto &d = m_diagnostic_output;
   // Call the horiz_contraction impl that will take care of everything
-  if (f.get_header().has_extra_data("mask_data")) {
-    horiz_contraction<Real>(d, f, m_scaled_area, &m_comm, m_dummy_field);
-  } else {
-    horiz_contraction<Real>(d, f, m_scaled_area, &m_comm);
-  }
+  horiz_contraction(d, f, m_scaled_area, true, &m_comm, m_dummy_field);
 }
 
 }  // namespace scream

components/eamxx/src/share/diagnostics/tests/aerocom_cld_test.cpp

@@ -104,10 +104,8 @@ TEST_CASE("aerocom_cld") {
   ni.allocate_view();
   ni.get_header().get_tracking().update_time_stamp(t0);
 
-  // Construct random number generator stuff
-  using RPDF = std::uniform_real_distribution<Real>;
-  RPDF pdf(0, 0.05);
-  auto engine = scream::setup_random_test();
+  // Random number generator seed
+  int seed = get_random_test_seed();
 
   // Construct the Diagnostics
   std::map<std::string, std::shared_ptr<AtmosphereDiagnostic>> diags;
@@ -125,17 +123,17 @@ TEST_CASE("aerocom_cld") {
   constexpr int ntests = 3;
   for(int itest = 0; itest < ntests; ++itest) {
     // Randomize everything to add ensure resiliency
-    randomize(tm, engine, pdf);
-    randomize(pd, engine, pdf);
-    randomize(pm, engine, pdf);
-    randomize(qv, engine, pdf);
-    randomize(qc, engine, pdf);
-    randomize(qi, engine, pdf);
-    randomize(ec, engine, pdf);
-    randomize(ei, engine, pdf);
-    randomize(cd, engine, pdf);
-    randomize(nc, engine, pdf);
-    randomize(ni, engine, pdf);
+    randomize_uniform(tm, seed++, 0, 0.05);
+    randomize_uniform(pd, seed++, 0, 0.05);
+    randomize_uniform(pm, seed++, 0, 0.05);
+    randomize_uniform(qv, seed++, 0, 0.05);
+    randomize_uniform(qc, seed++, 0, 0.05);
+    randomize_uniform(qi, seed++, 0, 0.05);
+    randomize_uniform(ec, seed++, 0, 0.05);
+    randomize_uniform(ei, seed++, 0, 0.05);
+    randomize_uniform(cd, seed++, 0, 0.05);
+    randomize_uniform(nc, seed++, 0, 0.05);
+    randomize_uniform(ni, seed++, 0, 0.05);
 
     // Create and set up the diagnostic
     params.set<std::string>("aero_com_cld_kind", "Top");

components/eamxx/src/share/diagnostics/tests/aodvis_test.cpp

@@ -66,10 +66,8 @@ TEST_CASE("aodvis") {
   sunlit.allocate_view();
   sunlit.get_header().get_tracking().update_time_stamp(t0);
 
-  // Construct random number generator stuff
-  using RPDF = std::uniform_real_distribution<Real>;
-  RPDF pdf(0, 0.005);
-  auto engine = scream::setup_random_test();
+  // Random number generator seed
+  int seed = get_random_test_seed(&comm);
 
   // Construct the Diagnostics
   std::map<std::string, std::shared_ptr<AtmosphereDiagnostic>> diags;
@@ -79,10 +77,10 @@ TEST_CASE("aodvis") {
   constexpr int ntests = 5;
   for(int itest = 0; itest < ntests; ++itest) {
     // Randomize tau
-    randomize(tau, engine, pdf);
+    randomize_uniform(tau, seed++, 0, 0.005);
 
     // Randomize sunlit
-    randomize(sunlit, engine, pdf);
+    randomize_uniform(sunlit, seed++, 0, 0.005);
 
     // Create and set up the diagnostic
     ekat::ParameterList params;