Merge branch 'bartgol/eamxx/fix-horiz-remap-bug' (PR #7350)

Fixes an issue in coarsening remap in case not all
the fine grid cols participate to the remap.
This bug affected ARM site output.

Fixes #7346

[BFB] except for EAMxx ARM output

Author: bartgol

components/eamxx/src/diagnostics/precip_surf_mass_flux.cpp

@@ -63,27 +63,29 @@ void PrecipSurfMassFlux::compute_diagnostic_impl()
   const auto use_liq = m_type & s_liq;
   const auto use_ice = m_type & s_ice;
 
-  std::int64_t dt;
+  double dt = 0;
   if (use_ice) {
     auto mass_ice = get_field_in("precip_ice_surf_mass");
     mass_ice_d = mass_ice.get_view<const Real*>();
 
     const auto& t_start = mass_ice.get_header().get_tracking().get_accum_start_time ();
     const auto& t_now   = mass_ice.get_header().get_tracking().get_time_stamp ();
-    dt = t_now-t_start;
-  }
-  if (use_liq) {
+    dt = t_now.seconds_from(t_start);
+    if (use_liq) {
+      // Ensure liq/ice have same accumulation times
+      auto mass_liq = get_field_in("precip_liq_surf_mass");
+      mass_liq_d = mass_liq.get_view<const Real*>();
+      EKAT_REQUIRE_MSG (t_now==mass_liq.get_header().get_tracking().get_time_stamp() and
+                        t_start==mass_liq.get_header().get_tracking().get_accum_start_time(),
+          "Error! Liquid and ice precip mass fields have different accumulation time stamps!\n");
+    }
+  } else if (use_liq) {
     auto mass_liq = get_field_in("precip_liq_surf_mass");
     mass_liq_d = mass_liq.get_view<const Real*>();
 
     const auto& t_start = mass_liq.get_header().get_tracking().get_accum_start_time ();
     const auto& t_now   = mass_liq.get_header().get_tracking().get_time_stamp ();
-    if (use_ice) {
-      EKAT_REQUIRE_MSG (dt==(t_now-t_start),
-          "Error! Liquid and ice precip mass fields have different accumulation time stamps!\n");
-    } else {
-      dt = t_now-t_start;
-    }
+    dt = t_now.seconds_from(t_start);
   }
 
   if (dt==0) {

components/eamxx/src/diagnostics/tests/precip_surf_mass_flux_tests.cpp

@@ -48,7 +48,7 @@ void run(std::mt19937_64& engine)
   auto gm = create_gm(comm,ncols);
 
   // Create timestep
-  const int dt=1800;
+  const double dt=1800;
 
   // Construct random input data
   using RPDF = std::uniform_real_distribution<Real>;
@@ -120,28 +120,29 @@ void run(std::mt19937_64& engine)
   diag_liq->compute_diagnostic();
   diag_ice->compute_diagnostic();
 
-  Field preicp_total_f = diag_total->get_diagnostic().clone();
-  Field preicp_liq_f   = diag_liq->get_diagnostic().clone();
-  Field preicp_ice_f   = diag_ice->get_diagnostic().clone();
-  preicp_total_f.deep_copy(0);
-  preicp_liq_f.deep_copy(0);
-  preicp_ice_f.deep_copy(0);
-  auto precip_total_v = preicp_total_f.get_view<Real*>();
-  auto precip_liq_v   = preicp_liq_f.get_view<Real*>();
-  auto precip_ice_v   = preicp_ice_f.get_view<Real*>();
+  Field precip_total_f = diag_total->get_diagnostic().clone();
+  Field precip_liq_f   = diag_liq->get_diagnostic().clone();
+  Field precip_ice_f   = diag_ice->get_diagnostic().clone();
+  precip_total_f.deep_copy(0);
+  precip_liq_f.deep_copy(0);
+  precip_ice_f.deep_copy(0);
+  auto precip_total_v = precip_total_f.get_view<Real*>();
+  auto precip_liq_v   = precip_liq_f.get_view<Real*>();
+  auto precip_ice_v   = precip_ice_f.get_view<Real*>();
   const auto rhodt = PC::RHO_H2O*dt;
   Kokkos::parallel_for("precip_total_surf_mass_flux_test",
                        typename KT::RangePolicy(0,ncols),
                        KOKKOS_LAMBDA(const int& icol) {
     precip_liq_v(icol)   = precip_liq_surf_mass_v(icol)/rhodt;
     precip_ice_v(icol)   = precip_ice_surf_mass_v(icol)/rhodt;
-    precip_total_v(icol) = precip_liq_v(icol) + precip_ice_v(icol);
+    precip_total_v(icol) = precip_ice_v(icol);
+    precip_total_v(icol) += precip_liq_surf_mass_v(icol)/rhodt;
   });
   Kokkos::fence();
 
-  REQUIRE(views_are_equal(diag_total->get_diagnostic(),preicp_total_f));
-  REQUIRE(views_are_equal(diag_liq->get_diagnostic(),preicp_liq_f));
-  REQUIRE(views_are_equal(diag_ice->get_diagnostic(),preicp_ice_f));
+  REQUIRE(views_are_equal(diag_total->get_diagnostic(),precip_total_f));
+  REQUIRE(views_are_equal(diag_liq->get_diagnostic(),precip_liq_f));
+  REQUIRE(views_are_equal(diag_ice->get_diagnostic(),precip_ice_f));
 
   // Finalize the diagnostic
   diag_total->finalize();

components/eamxx/src/diagnostics/tests/vertical_layer_tests.cpp

@@ -134,10 +134,10 @@ void run (const std::string& diag_name, const std::string& location)
         // If interface, check value, otherwise perform int->mid averaging and check value
         auto int_val = prev_int_val + delta;
         if (location=="interfaces") {
-          REQUIRE_THAT(d_h(icol,ilev), Catch::Matchers::WithinRel(int_val,1e-5));
+          REQUIRE_THAT(d_h(icol,ilev), Catch::Matchers::WithinRel(int_val,Real(1e-5)));
         } else {
           auto mid_val = (int_val + prev_int_val) / 2;
-          REQUIRE_THAT(d_h(icol,ilev), Catch::Matchers::WithinRel(mid_val,1e-5));
+          REQUIRE_THAT(d_h(icol,ilev), Catch::Matchers::WithinRel(mid_val,Real(1e-5)));
         }
         prev_int_val = int_val;
       }

components/eamxx/src/dynamics/homme/interface/dyn_grid_mod.F90

@@ -4,13 +4,11 @@ module dyn_grid_mod
   use shr_kind_mod,       only: r8 => shr_kind_r8
   use dimensions_mod,     only: nelem, nelemd, nelemdmax, np
   use edgetype_mod,       only: EdgeBuffer_t
+  use mpi
 
   implicit none
   private
 
-! We need MPI in here, so include it
-#include <mpif.h>
-
   public :: dyn_grid_init, get_my_dyn_data, cleanup_grid_init_data
 
   type (EdgeBuffer_t) :: edge

components/eamxx/src/dynamics/homme/interface/phys_grid_mod.F90

@@ -3,6 +3,7 @@ module phys_grid_mod
   use iso_c_binding, only: c_int, c_double
   use parallel_mod,  only: abortmp, MPIinteger_t, MPIreal_t
   use kinds,         only: iulog
+  use mpi
 
   implicit none
   private
@@ -55,8 +56,6 @@ module phys_grid_mod
 
   type(pg_specs_t), target :: pg_specs (pgN_min:pgN_max)
 
-! To get MPI_IN_PLACE and MPI_DATATYPE_NULL
-#include <mpif.h>
   ! Note: in this module, we often use MPI_IN_PLACE,0,MPI_DATATYPE_NULL
   !       for the src array specs in Allgatherv calls. These special values 
   !       inform MPI that src array is aliasing the dst one, so MPI will

components/eamxx/src/share/grid/remap/horiz_interp_remapper_data.cpp

@@ -88,41 +88,10 @@ get_my_triplets (const std::string& map_file) const
 
   scorpio::release_file(map_file);
 
-  // 1.2 Dofs in grid are likely 0-based, while row/col ids in map file
-  // are likely 1-based. To match dofs, we need to offset the row/cols
-  // ids we just read in.
-  int map_file_min_row = std::numeric_limits<int>::max();
-  int map_file_min_col = std::numeric_limits<int>::max();
-  for (int id=0; id<nlweights; id++) {
-    map_file_min_row = std::min(rows[id],map_file_min_row);
-    map_file_min_col = std::min(cols[id],map_file_min_col);
-  }
-  int global_map_file_min_row, global_map_file_min_col;
-  comm.all_reduce(&map_file_min_row,&global_map_file_min_row,1,MPI_MIN);
-  comm.all_reduce(&map_file_min_col,&global_map_file_min_col,1,MPI_MIN);
-
-  gid_type row_offset = global_map_file_min_row;
-  gid_type col_offset = global_map_file_min_col;
-  if (type==InterpType::Refine) {
-    row_offset -= fine_grid->get_global_min_dof_gid();
-  } else {
-    col_offset -= fine_grid->get_global_min_dof_gid();
-  }
-  for (auto& id : rows) {
-    id -= row_offset;
-  }
-  for (auto& id : cols) {
-    id -= col_offset;
-  }
-
   // Create a grid based on the row gids I read in (may be duplicated across ranks)
-  std::vector<gid_type> unique_gids;
   const auto& gids = type==InterpType::Refine ? rows : cols;
-  for (auto gid : gids) {
-    if (not ekat::contains(unique_gids,gid)) {
-      unique_gids.push_back(gid);
-    }
-  }
+  std::set<gid_type> temp (gids.begin(),gids.end());
+  std::vector<gid_type> unique_gids (temp.begin(),temp.end());
   auto io_grid = std::make_shared<PointGrid> ("helper",unique_gids.size(),0,comm);
   auto io_grid_gids_h = io_grid->get_dofs_gids().get_view<gid_type*,Host>();
   int k = 0;
@@ -150,10 +119,18 @@ get_my_triplets (const std::string& map_file) const
   MPI_Type_create_struct (3,lengths,displacements,types,&mpi_triplet_t);
   MPI_Type_commit(&mpi_triplet_t);
 
-  // Create import-export
-  GridImportExport imp_exp (fine_grid,io_grid);
+  // Create import-export and gather my triplets
   std::map<int,std::vector<Triplet>> my_triplets_map;
-  imp_exp.gather(mpi_triplet_t,io_triplets,my_triplets_map);
+  try {
+    GridImportExport imp_exp (fine_grid,io_grid);
+    imp_exp.gather(mpi_triplet_t,io_triplets,my_triplets_map);
+  } catch (...) {
+    // Print the map file name, to help debugging
+    std::cout << "[HorizRemapperData] Something went wrong while performing a grid import-export operation.\n"
+              << " - map file name : " << map_file << "\n"
+              << " - fine grid name: " << fine_grid->name() << "\n";
+    throw;
+  }
   MPI_Type_free(&mpi_triplet_t);
 
   std::vector<Triplet> my_triplets;

components/eamxx/src/share/grid/remap/horiz_interp_remapper_data.hpp

@@ -17,13 +17,16 @@ enum class InterpType {
   Coarsen
 };
 
-// A small struct to hold horiz remap data, which can
-// be shared across multiple horiz remappers
+// A small struct to hold horiz remap data, which can be shared across multiple horiz remappers
+// NOTE: the client will call the build method, which will read the map file, and create the
+//       CRS matrix data for online interpolation.
 struct HorizRemapperData {
   using KT = KokkosTypes<DefaultDevice>;
   template<typename T>
   using view_1d = typename KT::template view_1d<T>;
 
+  // The last argument specifies the base index for gids in the map file
+  // For ncremap-type files, all indices are 1-based
   void build (const std::string& map_file,
               const std::shared_ptr<const AbstractGrid>& fine_grid,
               const ekat::Comm& comm,

components/eamxx/src/share/io/tests/CMakeLists.txt

@@ -88,6 +88,12 @@ CreateUnitTest(io_remap_test "io_remap_test.cpp"
   MPI_RANKS 1 ${SCREAM_TEST_MAX_RANKS}
 )
 
+## Test remap output when map file is sub-sampling (ARM-style)
+CreateUnitTest(io_horiz_sampling "io_horiz_sampling.cpp"
+  LIBS scream_io LABELS io remap
+  MPI_RANKS 1 ${SCREAM_TEST_MAX_RANKS}
+)
+
 ## Test single-column reader
 CreateUnitTest(io_scm_reader "io_scm_reader.cpp"
   LIBS scream_io LABELS io