Merge branch 'tcclevenger/eamxx/add_tracer_with_turbulence_advection' (PR #6789)

Processes can choose for a tracer to be advected only by dynamics (not by SHOC) by passing
turbulence_advected=false to add_tracer(). Guards set up to ensure all tracers requests are
consistent among processes.

[non-BFB for MAMxx]

Author: tcclevenger

components/eamxx/src/control/atmosphere_driver.cpp

@@ -58,9 +58,7 @@ namespace control {
  *     Note: at this stage, atm procs that act on non-ref grid(s) should be able to create their
  *           remappers. The AD will *not* take care of remapping inputs/outputs of the process.
  *  4) Register all fields and all groups from all atm procs inside the field managers, and proceed
- *     to allocate fields. Each field manager (there is one FM per grid) will take care of
- *     accommodating all requests for packing as well as (if possible) bundling of groups.
- *     For more details, see the documentation in the share/field/field_request.hpp header.
+ *     to allocate fields. For more details, see the documentation in the share/field/field_request.hpp header.
  *  5) Set all the fields into the atm procs. Before this point, all the atm procs had were the
  *     FieldIdentifiers for their input/output fields and FieldGroupInfo for their input/output
  *     field groups. Now, we pass actual Field and FieldGroup objects to them, where both the
@@ -405,7 +403,7 @@ void AtmosphereDriver::reset_accumulated_fields ()
     }
 
     auto accum_group = m_field_mgr->get_field_group("ACCUMULATED", grid_name);
-    for (auto f_it : accum_group.m_fields) {
+    for (auto f_it : accum_group.m_individual_fields) {
       auto& track = f_it.second->get_header().get_tracking();
       f_it.second->deep_copy(zero);
       track.set_accum_start_time(m_current_ts);
@@ -535,8 +533,54 @@ void AtmosphereDriver::create_fields()
   m_field_mgr = std::make_shared<field_mgr_type>(m_grids_manager);
   m_field_mgr->registration_begins();
 
-  // By now, the processes should have fully built the ids of their
-  // required/computed fields and groups. Let them register them in the FM
+  // Before registering fields, check that Field Requests for tracers are compatible
+  {
+    // Create map from tracer name to a vector which contains the field requests for that tracer.
+    std::map<std::string, std::set<FieldRequest>> tracer_requests;
+    auto gather_tracer_requests = [&] (FieldRequest req) {
+      if (not ekat::contains(req.groups, "tracers")) return;
+
+      std::string fname = req.fid.name();
+      if (tracer_requests.find(fname) == tracer_requests.end()) {
+        tracer_requests[fname] = {req};
+      } else {
+        tracer_requests[fname].emplace(req);
+      }
+    };
+    for (const auto& req : m_atm_process_group->get_required_field_requests()){
+      gather_tracer_requests(req);
+    }
+    for (const auto& req : m_atm_process_group->get_computed_field_requests()) {
+      gather_tracer_requests(req);
+    }
+
+    // Go through the map entry for each tracer and check that every one
+    // has the same request for turbulence advection.
+    for (auto fr : tracer_requests) {
+      const auto reqs = fr.second;
+
+      std::set<bool> turb_advect_types;
+      for (auto req : reqs) {
+        turb_advect_types.emplace(ekat::contains(req.groups, "turbulence_advected_tracers"));
+      }
+
+      if (turb_advect_types.size()!=1) {
+        std::ostringstream ss;
+        ss << "Error! Incompatible tracer request. Turbulence advection requests not consistent among processes.\n"
+              "  - Tracer name: " + fr.first + "\n"
+              "  - Requests (process name, grid name, is tracers turbulence advected):\n";
+        for (auto req : reqs) {
+          const auto grid_name = req.fid.get_grid_name();
+          const bool turb_advect = ekat::contains(req.groups, "turbulence_advected_tracers");
+          ss << "    - (" + req.calling_process + ", " + grid_name + ", " + (turb_advect ? "true" : "false") + ")\n";
+        }
+        EKAT_ERROR_MSG(ss.str());
+      }
+    }
+  }
+
+  // Register required/computed fields. By now, the processes should have
+  // fully built the ids of their required/computed fields and groups
   for (const auto& req : m_atm_process_group->get_required_field_requests()) {
     m_field_mgr->register_field(req);
   }
@@ -599,8 +643,8 @@ void AtmosphereDriver::create_fields()
     m_field_mgr->add_to_group(fid, "RESTART");
   }
   for (const auto& g : m_atm_process_group->get_groups_in()) {
-    if (g.m_bundle) {
-      m_field_mgr->add_to_group(g.m_bundle->get_header().get_identifier(), "RESTART");
+    if (g.m_monolithic_field) {
+      m_field_mgr->add_to_group(g.m_monolithic_field->get_header().get_identifier(), "RESTART");
     } else {
       for (const auto& fn : g.m_info->m_fields_names) {
         m_field_mgr->add_to_group(fn, g.grid_name(), "RESTART");
@@ -1128,19 +1172,19 @@ void AtmosphereDriver::set_initial_conditions ()
   // ...then the input groups
   m_atm_logger->debug("    [EAMxx] Processing input groups ...");
   for (const auto& g : m_atm_process_group->get_groups_in()) {
-    if (g.m_bundle) {
-      process_ic_field(*g.m_bundle);
+    if (g.m_monolithic_field) {
+      process_ic_field(*g.m_monolithic_field);
     }
-    for (auto it : g.m_fields) {
+    for (auto it : g.m_individual_fields) {
       process_ic_field(*it.second);
     }
   }
   m_atm_logger->debug("    [EAMxx] Processing input groups ... done!");
 
-  // Some fields might be the subfield of a group's bundled field. In that case,
-  // we only need to init one: either the bundled field, or all the individual subfields.
+  // Some fields might be the subfield of a group's monolithic field. In that case,
+  // we only need to init one: either the monolithic field, or all the individual subfields.
   // So loop over the fields that appear to require loading from file, and remove
-  // them from the list if they are the subfield of a bundled field already inited
+  // them from the list if they are the subfield of a groups monolithic field already inited
   // (perhaps via initialize_constant_field, or copied from another field).
   for (auto& it1 : ic_fields_names) {
     const auto& grid_name =  it1.first;
@@ -1244,17 +1288,17 @@ void AtmosphereDriver::set_initial_conditions ()
   }
   m_atm_logger->debug("    [EAMxx] Processing fields to copy ... done!");
 
-  // It is possible to have a bundled group G1=(f1,f2,f3),
+  // It is possible to have a monolithically allocated group G1=(f1,f2,f3),
   // where the IC are read from file for f1, f2, and f3. In that case,
-  // the time stamp for the bundled G1 has not be inited, but the data
+  // the time stamp for the monolithic field of G1 has not be inited, but the data
   // is valid (all entries have been inited). Let's fix that.
   m_atm_logger->debug("    [EAMxx] Processing subfields ...");
   for (const auto& g : m_atm_process_group->get_groups_in()) {
-    if (g.m_bundle) {
-      auto& track = g.m_bundle->get_header().get_tracking();
+    if (g.m_monolithic_field) {
+      auto& track = g.m_monolithic_field->get_header().get_tracking();
       if (not track.get_time_stamp().is_valid()) {
-        // The bundled field has not been inited. Check if all the subfields
-        // have been inited. If so, init the timestamp of the bundled field too.
+        // The groups monolithic field has not been inited. Check if all the subfields
+        // have been inited. If so, init the timestamp of the monlithic field too.
         const auto& children = track.get_children();
         bool all_inited = children.size()>0; // If no children, then something is off, so mark as not good
         for (auto wp : children) {
@@ -1632,7 +1676,7 @@ void AtmosphereDriver::run (const int dt) {
     }
 
     auto rescale_group = m_field_mgr->get_field_group("DIVIDE_BY_DT", gname);
-    for (auto f_it : rescale_group.m_fields) {
+    for (auto f_it : rescale_group.m_individual_fields) {
       f_it.second->scale(Real(1) / dt);
     }
   }

components/eamxx/src/control/tests/dummy_atm_proc.hpp

@@ -93,8 +93,8 @@ class DummyProcess : public scream::AtmosphereProcess {
       });
     } else if (m_name=="Group to Group") {
       const auto& g = get_group_out("The Group");
-      const auto view_B = g.m_fields.at("B")->get_view<Real**>();
-      const auto view_C = g.m_fields.at("C")->get_view<Real**>();
+      const auto view_B = g.m_individual_fields.at("B")->get_view<Real**>();
+      const auto view_C = g.m_individual_fields.at("C")->get_view<Real**>();
 
       Kokkos::parallel_for(policy,KOKKOS_LAMBDA(const int idx) {
         const int icol = idx / nlevs;
@@ -105,8 +105,8 @@ class DummyProcess : public scream::AtmosphereProcess {
       });
     } else {
       const auto& g = get_group_in("The Group");
-      const auto view_B = g.m_fields.at("B")->get_view<const Real**>();
-      const auto view_C = g.m_fields.at("C")->get_view<const Real**>();
+      const auto view_B = g.m_individual_fields.at("B")->get_view<const Real**>();
+      const auto view_C = g.m_individual_fields.at("C")->get_view<const Real**>();
       const auto view_A = get_field_out("A").get_view<Real**>();
 
       Kokkos::parallel_for(policy,KOKKOS_LAMBDA(const int idx) {
@@ -128,7 +128,7 @@ class DummyProcess : public scream::AtmosphereProcess {
 
   std::string m_name;
 
-  DummyType     m_dummy_type; 
+  DummyType     m_dummy_type;
 };
 
 } // namespace scream

components/eamxx/src/dynamics/homme/eamxx_homme_fv_phys.cpp

@@ -111,7 +111,7 @@ void HommeDynamics::fv_phys_dyn_to_fv_phys (const bool restart) {
     t.T_mid = Homme::ExecView<Real***>("T_mid_tmp", nelem, npg, npacks*N);
     t.horiz_winds = Homme::ExecView<Real****>("horiz_winds_tmp", nelem, npg, 2, npacks*N);
     // Really need just the first tracer.
-    const auto qsize = get_group_out("tracers", pgn).m_bundle->get_view<Real***>().extent_int(1);
+    const auto qsize = get_group_out("tracers", pgn).m_monolithic_field->get_view<Real***>().extent_int(1);
     t.tracers = Homme::ExecView<Real****>("tracers_tmp", nelem, npg, qsize, npacks*N);
     remap_dyn_to_fv_phys(&t);
     assert(ncols == nelem*npg);
@@ -138,7 +138,7 @@ void HommeDynamics::fv_phys_dyn_to_fv_phys (const bool restart) {
       auto f = get_field_out(n,pgn);
       f.get_header().get_tracking().update_time_stamp(timestamp());
     }
-    auto Q = get_group_out("tracers",pgn).m_bundle;
+    auto Q = get_group_out("tracers",pgn).m_monolithic_field;
     Q->get_header().get_tracking().update_time_stamp(timestamp());
   }
   update_pressure(m_phys_grid);
@@ -169,7 +169,7 @@ void HommeDynamics::remap_dyn_to_fv_phys (GllFvRemapTmp* t) const {
   const auto npg = m_phys_grid_pgN*m_phys_grid_pgN;
   const auto& gn = m_phys_grid->name();
   const auto nlev = get_field_out("T_mid", gn).get_view<Real**>().extent_int(1);
-  const auto nq = get_group_out("tracers").m_bundle->get_view<Real***>().extent_int(1);
+  const auto nq = get_group_out("tracers").m_monolithic_field->get_view<Real***>().extent_int(1);
   assert(get_field_out("T_mid", gn).get_view<Real**>().extent_int(0) == nelem*npg);
   assert(get_field_out("horiz_winds", gn).get_view<Real***>().extent_int(1) == 2);
 
@@ -189,7 +189,7 @@ void HommeDynamics::remap_dyn_to_fv_phys (GllFvRemapTmp* t) const {
     t ? t->horiz_winds.data() : get_field_out("horiz_winds", gn).get_view<Real***>().data(),
     nelem, npg, 2, nlev);
   const auto q = Homme::GllFvRemap::Phys3T(
-    t ? t->tracers.data() : get_group_out("tracers", gn).m_bundle->get_view<Real***>().data(),
+    t ? t->tracers.data() : get_group_out("tracers", gn).m_monolithic_field->get_view<Real***>().data(),
     nelem, npg, nq, nlev);
   const auto dp = Homme::GllFvRemap::Phys2T(
     get_field_out("pseudo_density", gn).get_view<Real**>().data(),
@@ -209,7 +209,7 @@ void HommeDynamics::remap_fv_phys_to_dyn () const {
   const auto npg = m_phys_grid_pgN*m_phys_grid_pgN;
   const auto& gn = m_phys_grid->name();
   const auto nlev = m_helper_fields.at("FT_phys").get_view<const Real**>().extent_int(1);
-  const auto nq = get_group_in("tracers", gn).m_bundle->get_view<const Real***>().extent_int(1);
+  const auto nq = get_group_in("tracers", gn).m_monolithic_field->get_view<const Real***>().extent_int(1);
   assert(m_helper_fields.at("FT_phys").get_view<const Real**>().extent_int(0) == nelem*npg);
 
   const auto uv_ndim = m_helper_fields.at("FM_phys").get_view<const Real***>().extent_int(1);
@@ -222,7 +222,7 @@ void HommeDynamics::remap_fv_phys_to_dyn () const {
     m_helper_fields.at("FM_phys").get_view<const Real***>().data(),
     nelem, npg, uv_ndim, nlev);
   const auto q = Homme::GllFvRemap::CPhys3T(
-    get_group_in("tracers", gn).m_bundle->get_view<const Real***>().data(),
+    get_group_in("tracers", gn).m_monolithic_field->get_view<const Real***>().data(),
     nelem, npg, nq, nlev);
 
   gfr.run_fv_phys_to_dyn(time_idx, T, uv, q);

components/eamxx/src/dynamics/homme/eamxx_homme_process_interface.cpp

@@ -180,7 +180,7 @@ void HommeDynamics::set_grids (const std::shared_ptr<const GridsManager> grids_m
   add_field<Computed>("omega",              pg_scalar3d_mid, Pa/s,  pgn,N);
 
   add_tracer<Updated >("qv", m_phys_grid, kg/kg, N);
-  add_group<Updated>("tracers",pgn,N, true);
+  add_group<Updated>("tracers",pgn,N, MonolithicAlloc::Required);
 
   if (fv_phys_active()) {
     // [CGLL ICs in pg2] Read CGLL IC data even though our in/out format is
@@ -196,7 +196,7 @@ void HommeDynamics::set_grids (const std::shared_ptr<const GridsManager> grids_m
     add_field<Required>("T_mid",         rg_scalar3d_mid,K,     rgn,N);
     add_field<Required>("ps",            rg_scalar2d    ,Pa,    rgn);
     add_field<Required>("phis",          rg_scalar2d    ,m2/s2, rgn);
-    add_group<Required>("tracers",rgn,N, true);
+    add_group<Required>("tracers",rgn,N, MonolithicAlloc::Required);
     fv_phys_rrtmgp_active_gases_init(grids_manager);
     // This is needed for the dp_ref init in initialize_homme_state.
     add_field<Computed>("pseudo_density",rg_scalar3d_mid,Pa,    rgn,N);
@@ -397,7 +397,7 @@ void HommeDynamics::initialize_impl (const RunType run_type)
     // ftype!=FORCING_0:
     //  1) remap Q_pgn->FQ_dyn
     // Remap Q directly into FQ, tendency computed in pre_process step
-    m_p2d_remapper->register_field(*get_group_out("Q",pgn).m_bundle,m_helper_fields.at("FQ_dyn"));
+    m_p2d_remapper->register_field(*get_group_out("Q",pgn).m_monolithic_field,m_helper_fields.at("FQ_dyn"));
     m_p2d_remapper->register_field(m_helper_fields.at("FT_phys"),m_helper_fields.at("FT_dyn"));
 
     // FM has 3 components on dyn grid, but only 2 on phys grid
@@ -412,7 +412,7 @@ void HommeDynamics::initialize_impl (const RunType run_type)
     m_d2p_remapper->register_field(get_internal_field("v_dyn"),get_field_out("horiz_winds"));
     m_d2p_remapper->register_field(get_internal_field("dp3d_dyn"), get_field_out("pseudo_density"));
     m_d2p_remapper->register_field(get_internal_field("ps_dyn"), get_field_out("ps"));
-    m_d2p_remapper->register_field(m_helper_fields.at("Q_dyn"),*get_group_out("Q",pgn).m_bundle);
+    m_d2p_remapper->register_field(m_helper_fields.at("Q_dyn"),*get_group_out("Q",pgn).m_monolithic_field);
     m_d2p_remapper->register_field(m_helper_fields.at("omega_dyn"), get_field_out("omega"));
 
     m_p2d_remapper->registration_ends();
@@ -463,7 +463,7 @@ void HommeDynamics::initialize_impl (const RunType run_type)
   using Interval = FieldWithinIntervalCheck;
   using LowerBound = FieldLowerBoundCheck;
 
-  add_postcondition_check<LowerBound>(*get_group_out("Q",pgn).m_bundle,m_phys_grid,0,true);
+  add_postcondition_check<LowerBound>(*get_group_out("Q",pgn).m_monolithic_field,m_phys_grid,0,true);
   add_postcondition_check<Interval>(get_field_out("T_mid",pgn),m_phys_grid,100.0, 500.0,false);
   add_postcondition_check<Interval>(get_field_out("horiz_winds",pgn),m_phys_grid,-400.0, 400.0,false);
   add_postcondition_check<Interval>(get_field_out("ps"),m_phys_grid,30000.0, 120000.0,false);
@@ -686,7 +686,7 @@ void HommeDynamics::homme_post_process (const double dt) {
   const auto dp_dry_view    = get_field_out("pseudo_density_dry").get_view<Pack**>();
   const auto p_dry_int_view = get_field_out("p_dry_int").get_view<Pack**>();
   const auto p_dry_mid_view = get_field_out("p_dry_mid").get_view<Pack**>();
-  const auto Q_view   = get_group_out("Q",pgn).m_bundle->get_view<Pack***>();
+  const auto Q_view   = get_group_out("Q",pgn).m_monolithic_field->get_view<Pack***>();
 
   const auto T_view  = get_field_out("T_mid").get_view<Pack**>();
   const auto T_prev_view = m_helper_fields.at("FT_phys").get_view<Pack**>();
@@ -1003,7 +1003,7 @@ void HommeDynamics::restart_homme_state () {
   auto qv_prev_ref = std::make_shared<Field>();
   auto Q_dyn = m_helper_fields.at("Q_dyn");
   if (params.ftype==Homme::ForcingAlg::FORCING_2) {
-    auto Q_old = *get_group_in("Q",pgn).m_bundle;
+    auto Q_old = *get_group_in("Q",pgn).m_monolithic_field;
     m_ic_remapper->register_field(Q_old,Q_dyn);
 
     // Grab qv_ref_old from Q_old
@@ -1106,7 +1106,7 @@ void HommeDynamics::initialize_homme_state () {
   m_ic_remapper->register_field(get_field_in("ps",rgn),get_internal_field("ps_dyn"));
   m_ic_remapper->register_field(get_field_in("phis",rgn),m_helper_fields.at("phis_dyn"));
   m_ic_remapper->register_field(get_field_in("T_mid",rgn),get_internal_field("vtheta_dp_dyn"));
-  m_ic_remapper->register_field(*get_group_in("tracers",rgn).m_bundle,m_helper_fields.at("Q_dyn"));
+  m_ic_remapper->register_field(*get_group_in("tracers",rgn).m_monolithic_field,m_helper_fields.at("Q_dyn"));
   m_ic_remapper->registration_ends();
   m_ic_remapper->remap_fwd();
 

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

@@ -29,7 +29,7 @@ void IOPForcing::set_grids(const std::shared_ptr<const GridsManager> grids_manag
   add_field<Updated>("T_mid", scalar3d_mid, K, grid_name, pack_size);
 
   add_tracer<Updated>("qv", m_grid, kg/kg, pack_size);
-  add_group<Updated>("tracers", grid_name, pack_size, true);
+  add_group<Updated>("tracers", grid_name, pack_size, MonolithicAlloc::Required);
 
   // Sanity check that iop data manager is setup by driver
   EKAT_REQUIRE_MSG(m_iop_data_manager,
@@ -64,8 +64,8 @@ set_computed_group_impl (const FieldGroup& group)
   EKAT_REQUIRE_MSG(name=="tracers",
     "Error! IOPForcing was not expecting a field group called '" << name << "\n");
 
-  EKAT_REQUIRE_MSG(group.m_info->m_bundled,
-      "Error! IOPForcing expects bundled fields for tracers.\n");
+  EKAT_REQUIRE_MSG(group.m_info->m_monolithic_allocation,
+      "Error! IOPForcing expects a monolithic allocation for tracers.\n");
 
   m_num_tracers = group.m_info->size();
 }
@@ -343,7 +343,7 @@ void IOPForcing::run_impl (const double dt)
   const auto horiz_winds = get_field_out("horiz_winds").get_view<Pack***>();
   const auto T_mid = get_field_out("T_mid").get_view<Pack**>();
   const auto qv = get_field_out("qv").get_view<Pack**>();
-  const auto Q = get_group_out("tracers").m_bundle->get_view<Pack***>();
+  const auto Q = get_group_out("tracers").m_monolithic_field->get_view<Pack***>();
 
   // Load data from IOP files, if necessary
   m_iop_data_manager->read_iop_file_data(timestamp());

components/eamxx/src/physics/mam/eamxx_mam_generic_process_interface.cpp

@@ -140,22 +140,24 @@ void MAMGenericInterface::add_tracers_interstitial_aerosol() {
   // ---------------------------------------------------------------------
   // These variables are "Updated" or inputs/outputs for the process
   // ---------------------------------------------------------------------
-  // NOTE: Cloud borne aerosols are not updated in this process but are included
-  // to create data structures.
+  // NOTE:
+  //   - Cloud borne aerosols are not updated in this process but are included
+  //     to create data structures.
+  //   - For interstitial aerosols, we have dynamics advect, but not turbulence.
 
   // interstitial and cloudborne aerosol tracers of interest: mass (q) and
   // number (n) mixing ratios
   for(int mode = 0; mode < mam_coupling::num_aero_modes(); ++mode) {
     // interstitial aerosol tracers of interest: number (n) mixing ratios
     const std::string int_nmr_field_name =
         mam_coupling::int_aero_nmr_field_name(mode);
-    add_tracer<Updated>(int_nmr_field_name, grid_, n_unit);
+    add_tracer<Updated>(int_nmr_field_name, grid_, n_unit, 1, TracerAdvection::DynamicsOnly);
     for(int a = 0; a < mam_coupling::num_aero_species(); ++a) {
       // (interstitial) aerosol tracers of interest: mass (q) mixing ratios
       const std::string int_mmr_field_name =
           mam_coupling::int_aero_mmr_field_name(mode, a);
       if(not int_mmr_field_name.empty()) {
-        add_tracer<Updated>(int_mmr_field_name, grid_, q_unit);
+        add_tracer<Updated>(int_mmr_field_name, grid_, q_unit, 1, TracerAdvection::DynamicsOnly);
       }
     }  // end for loop num species
   }    // end for loop for num modes