Use reduction to compute min and max particle distances in NeighborPa…

…rticles test. (#3212)

This replaces a much slower atomics-based check.

## Summary

## Additional background

## Checklist

The proposed changes:
- [ ] fix a bug or incorrect behavior in AMReX
- [ ] add new capabilities to AMReX
- [ ] changes answers in the test suite to more than roundoff level
- [ ] are likely to significantly affect the results of downstream AMReX
users
- [ ] include documentation in the code and/or rst files, if appropriate

---------

Co-authored-by: Weiqun Zhang <[email protected]>

Tests/Particles/NeighborParticles/Constants.H

@@ -8,8 +8,8 @@ namespace Params
     // This is designed to represent MFiX-like conditions where the grid spacing is
     //     roughly 2.5 times the particle diameter.  In main.cpp we set grid spacing to 1
     //     so here we set cutoff to diameter = 1/2.5 --> cutoff = 0.2
-    static constexpr amrex::Real cutoff = 0.2  ;
-    static constexpr amrex::Real min_r  = 1.e-4;
+    static constexpr amrex::ParticleReal cutoff = 0.2  ;
+    static constexpr amrex::ParticleReal min_r  = 1.e-4;
 }
 
 #endif

Tests/Particles/NeighborParticles/MDParticleContainer.cpp

@@ -1,8 +1,9 @@
+#include <AMReX_Reduce.H>
+#include <AMReX_SPACE.H>
+
 #include "MDParticleContainer.H"
 #include "Constants.H"
-
 #include "CheckPair.H"
-#include <AMReX_SPACE.H>
 
 using namespace amrex;
 
@@ -148,8 +149,9 @@ std::pair<Real, Real> MDParticleContainer::minAndMaxDistance()
     const int lev = 0;
     auto& plev  = GetParticles(lev);
 
-    Real min_d = std::numeric_limits<Real>::max();
-    Real max_d = std::numeric_limits<Real>::min();
+    ReduceOps<ReduceOpMin, ReduceOpMax> reduce_op;
+    ReduceData<ParticleReal, ParticleReal> reduce_data(reduce_op);
+    using ReduceTuple = typename decltype(reduce_data)::Type;
 
     for(MFIter mfi = MakeMFIter(lev); mfi.isValid(); ++mfi)
     {
@@ -159,41 +161,40 @@ std::pair<Real, Real> MDParticleContainer::minAndMaxDistance()
 
         auto& ptile = plev[index];
         auto& aos   = ptile.GetArrayOfStructs();
-        const size_t np = aos.numParticles();
 
         auto nbor_data = m_neighbor_list[lev][index].data();
         ParticleType* pstruct = aos().dataPtr();
 
-        Gpu::DeviceScalar<Real> min_d_gpu(min_d);
-        Gpu::DeviceScalar<Real> max_d_gpu(max_d);
-
-        Real* pmin_d = min_d_gpu.dataPtr();
-        Real* pmax_d = max_d_gpu.dataPtr();
-
-        AMREX_FOR_1D ( np, i,
-        {
-            ParticleType& p1 = pstruct[i];
+        ParticleReal min_start = std::numeric_limits<ParticleReal>::max();
+        ParticleReal max_start = std::numeric_limits<ParticleReal>::lowest();
 
-            for (const auto& p2 : nbor_data.getNeighbors(i))
-            {
-                AMREX_D_TERM(Real dx = p1.pos(0) - p2.pos(0);,
-                             Real dy = p1.pos(1) - p2.pos(1);,
-                             Real dz = p1.pos(2) - p2.pos(2);)
+        reduce_op.eval(aos.numParticles(), reduce_data,
+                       [=] AMREX_GPU_DEVICE (int i) -> ReduceTuple
+                       {
+                           ParticleType& p1 = pstruct[i];
 
-                Real r2 = AMREX_D_TERM(dx*dx, + dy*dy, + dz*dz);
-                r2 = amrex::max(r2, Params::min_r*Params::min_r);
-                Real r = sqrt(r2);
+                           ParticleReal min_d = min_start;
+                           ParticleReal max_d = max_start;
 
-                Gpu::Atomic::Min(pmin_d, r);
-                Gpu::Atomic::Max(pmax_d, r);
-            }
-        });
+                           for (const auto& p2 : nbor_data.getNeighbors(i))
+                           {
+                               AMREX_D_TERM(ParticleReal dx = p1.pos(0) - p2.pos(0);,
+                                            ParticleReal dy = p1.pos(1) - p2.pos(1);,
+                                            ParticleReal dz = p1.pos(2) - p2.pos(2);)
 
-        Gpu::Device::streamSynchronize();
+                               ParticleReal r2 = AMREX_D_TERM(dx*dx, + dy*dy, + dz*dz);
+                               r2 = amrex::max(r2, Params::min_r*Params::min_r);
+                               auto r = ParticleReal(std::sqrt(r2));
 
-        min_d = std::min(min_d, min_d_gpu.dataValue());
-        max_d = std::max(max_d, max_d_gpu.dataValue());
+                               min_d = std::min(min_d, r);
+                               max_d = std::max(max_d, r);
+                           }
+                           return {min_d, max_d};
+                       });
     }
+
+    ParticleReal min_d = amrex::get<0>(reduce_data.value(reduce_op));
+    ParticleReal max_d = amrex::get<1>(reduce_data.value(reduce_op));
     ParallelDescriptor::ReduceRealMin(min_d, ParallelDescriptor::IOProcessorNumber());
     ParallelDescriptor::ReduceRealMax(max_d, ParallelDescriptor::IOProcessorNumber());