Merge branch 'ambrad/hommexx/sl-fix-threading' (PR #7012)

Hommexx/SL: Fix a threading issue.

Fix a threading error in EAMxx simulations on CPU in the case that each element
has more than one thread. Add some team_barriers and rearrange a section of code
to permit more team_barriers. Clean up some compiler warnings.

[BFB]

Author: ambrad

components/homme/src/share/cxx/ComposeTransportImpl.hpp

@@ -356,9 +356,6 @@ struct ComposeTransportImpl {
     RNlev yp(pack2real(yps));
     const auto f = [&] (const int i, const int j, const int k) {
       if (k == 0) return;
-      const auto& xkm1 = x(i,j,k-1);
-      const auto& xk   = x(i,j,k  ); // also the interpolation point
-      const auto& xkp1 = x(i,j,k+1);
       yp(i,j,k) = approx_derivative(x(i,j,k-1), x(i,j,k), x(i,j,k+1),
                                     y(i,j,k-1), y(i,j,k), y(i,j,k+1));
     };
@@ -377,14 +374,19 @@ struct ComposeTransportImpl {
   {
     const auto ttr = Kokkos::TeamThreadRange(kv.team, NP*NP);
     const auto tvr = Kokkos::ThreadVectorRange(kv.team, NUM_LEV);
-    const auto f = [&] (const int idx) {
+    const auto f1 = [&] (const int idx) {
       const int i = idx / NP, j = idx % NP;
       const auto r = [&] (const int k, Real& dps, const bool final) {
         assert(k != 0 || dps == 0);
         if (final) edds(i,j,k) = dps;
         dps += divdps(i,j,k);
       };
       Dispatch<>::parallel_scan(kv.team, num_phys_lev, r);
+    };
+    Kokkos::parallel_for(ttr, f1);
+    kv.team_barrier();
+    const auto f2 = [&] (const int idx) {
+      const int i = idx / NP, j = idx % NP;
       const int kend = num_phys_lev - 1;
       const Real dps = edds(i,j,kend) + divdps(i,j,kend);
       assert(hybrid_bi(0)[0] == 0);
@@ -393,11 +395,15 @@ struct ComposeTransportImpl {
         if (kp == 0) edd(i,j,kp)[0] = 0;
       };
       Kokkos::parallel_for(tvr, s);
-      assert(edds(i,j,0) == 0);
-      const int bottom = num_phys_lev;
-      edds(i,j,bottom) = 0; // benign write race
     };
-    Kokkos::parallel_for(ttr, f);
+    Kokkos::parallel_for(ttr, f2);
+    kv.team_barrier();
+    const int bottom = num_phys_lev;
+    const auto f3 = [&] (const int idx) {
+      const int i = idx / NP, j = idx % NP;
+      Kokkos::single(Kokkos::PerThread(kv.team), [&] () { edds(i,j,bottom) = 0; });
+    };
+    Kokkos::parallel_for(ttr, f3);
   }
 };
 

components/homme/src/share/cxx/ComposeTransportImplEnhancedTrajectoryTests.cpp

@@ -25,9 +25,9 @@ get_test_team_policy (const int nelem, const int nlev, const int ncol=NP*NP) {
 }
 
 struct TestData {
+  const ComposeTransportImpl& cti;
   std::mt19937_64 engine;
   static const Real eps;
-  const ComposeTransportImpl& cti;
 
   TestData (const CTI& cti_, const int seed = 0)
     : cti(cti_), engine(seed == 0 ? std::random_device()() : seed)
@@ -103,7 +103,7 @@ void fillcols (const int n, const Real* const h, const RelnV::HostMirror& a) {
   assert(n <= a.extent_int(2));
   for (int i = 0; i < a.extent_int(0); ++i)
     for (int j = 0; j < a.extent_int(1); ++j)
-      for (size_t k = 0; k < n; ++k)
+      for (int k = 0; k < n; ++k)
         a(i,j,k) = h[k];  
 }
 

components/homme/src/share/cxx/ComposeTransportImplTrajectory.cpp

@@ -75,7 +75,6 @@ reconstruct_and_limit_dp (const KernelVariables& kv, const CSNlev& dprefp, const
     };
     Kokkos::Real2 sums;
     Dispatch<>::parallel_reduce(kv.team, tvr, g1, sums);
-    kv.team_barrier();
     const Real nmass = sums.v[0];
     if (nmass == 0) return;
     // Compensate for clipping.
@@ -157,7 +156,6 @@ KOKKOS_FUNCTION static void calc_vertically_lagrangian_levels (
   assert(hybrid_bi(0)[0] == 0);
 
   const auto ttr = TeamThreadRange(kv.team, NP*NP);
-  const auto tvr = ThreadVectorRange(kv.team, NUM_LEV);
 
   // Reconstruct an approximation to endpoint eta_dot_dpdn on Eulerian levels.
   const auto& divdp = wrk1a;
@@ -183,6 +181,7 @@ KOKKOS_FUNCTION static void calc_vertically_lagrangian_levels (
     kv.team_barrier();
     RNlevp edds(cti::pack2real(edd)), divdps(cti::pack2real(divdp));
     cti::calc_eta_dot_dpdn(kv, hybrid_bi, divdps, edd, edds);
+    kv.team_barrier();
   }
 
   // Use p0 as the reference coordinate system. p0 differs from p1 by B(eta)
@@ -203,6 +202,7 @@ KOKKOS_FUNCTION static void calc_vertically_lagrangian_levels (
   // time.
   const auto& ptp0 = dprecon;
   cti::approx_derivative(kv, pref, *eta_dot_dpdn[1], ptp0);
+  kv.team_barrier();
 
   {
     const auto& edd = *eta_dot_dpdn[0];
@@ -240,14 +240,16 @@ KOKKOS_FUNCTION static void calc_vertically_lagrangian_levels (
     if (static_cast<int>(cti::num_lev_pack) ==
         static_cast<int>(cti::max_num_lev_pack)) {
       // Re-zero eta_dot_dpdn at bottom.
+      kv.team_barrier();
       RNlevp edds(cti::pack2real(edd));
       const auto f = [&] (const int idx) {
         const int i = idx / NP, j = idx % NP;
         const int bottom = cti::num_phys_lev;
-        edds(i,j,bottom) = 0;
+        Kokkos::single(Kokkos::PerThread(kv.team), [&] () { edds(i,j,bottom) = 0; });
       };
       parallel_for(ttr, f);
     }
+    kv.team_barrier();
   }
 
   reconstruct_and_limit_dp(kv, dp3d, dt, dp_tol, *eta_dot_dpdn[0], dprecon);