Block-Structured AMR Software Framework
 
Loading...
Searching...
No Matches
AMReX_ParticleReduce.H
Go to the documentation of this file.
1#ifndef AMREX_PARTICLEREDUCE_H_
2#define AMREX_PARTICLEREDUCE_H_
3#include <AMReX_Config.H>
4
5#include <AMReX_IntVect.H>
6#include <AMReX_Box.H>
7#include <AMReX_Gpu.H>
8#include <AMReX_Print.H>
9#include <AMReX_GpuUtility.H>
10#include <AMReX_TypeTraits.H>
11#include <AMReX_ParticleUtil.H>
12#include <AMReX_Vector.H>
13
14#include <limits>
15#include <type_traits>
16
17namespace amrex {
18
19namespace particle_detail {
20
21template <typename F, typename T_ParticleType, int NAR, int NAI>
22AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
23auto call_f (F const& f,
24 const ConstParticleTileData<T_ParticleType, NAR, NAI>& p,
25 const int i) noexcept
26{
27 if constexpr ( ! T_ParticleType::is_soa_particle &&
28 IsCallable<F, T_ParticleType const&>::value) {
29 return f(p.m_aos[i]);
30 } else if constexpr (IsCallable<F, decltype(p.getSuperParticle(i))>::value) {
31 return f(p.getSuperParticle(i));
32 } else {
33 return f(p, i);
34 }
35}
36}
37
78template <class PC, class F, std::enable_if_t<IsParticleContainer<PC>::value, int> foo = 0>
79auto
80ReduceSum (PC const& pc, F&& f)
81 -> decltype(particle_detail::call_f(f, typename PC::ConstPTDType(), int()))
82{
83 return ReduceSum(pc, 0, pc.finestLevel(), std::forward<F>(f));
84}
85
126template <class PC, class F, std::enable_if_t<IsParticleContainer<PC>::value, int> foo = 0>
127auto
128ReduceSum (PC const& pc, int lev, F&& f)
129 -> decltype(particle_detail::call_f(f, typename PC::ConstPTDType(), int()))
130{
131 return ReduceSum(pc, lev, lev, std::forward<F>(f));
132}
133
175template <class PC, class F, std::enable_if_t<IsParticleContainer<PC>::value, int> foo = 0>
176auto
177ReduceSum (PC const& pc, int lev_min, int lev_max, F const& f)
178 -> decltype(particle_detail::call_f(f, typename PC::ConstPTDType(), int()))
179{
180 using value_type = decltype(particle_detail::call_f(f, typename PC::ConstPTDType(), int()));
181 value_type sm = 0;
182
183#ifdef AMREX_USE_GPU
184 if (Gpu::inLaunchRegion())
185 {
186 ReduceOps<ReduceOpSum> reduce_op;
187 ReduceData<value_type> reduce_data(reduce_op);
188 using ReduceTuple = typename decltype(reduce_data)::Type;
189
190 for (int lev = lev_min; lev <= lev_max; ++lev)
191 {
192 const auto& plev = pc.GetParticles(lev);
193 for (const auto& kv : plev)
194 {
195 const auto& tile = plev.at(kv.first);
196 const auto np = tile.numParticles();
197 const auto& ptd = tile.getConstParticleTileData();
198 reduce_op.eval(np, reduce_data,
199 [=] AMREX_GPU_DEVICE (const int i) -> ReduceTuple {
200 return particle_detail::call_f(f, ptd, i);
201 });
202 }
203 }
204
205 ReduceTuple hv = reduce_data.value(reduce_op);
206 sm = amrex::get<0>(hv);
207 }
208 else
209#endif
210 {
211 for (int lev = lev_min; lev <= lev_max; ++lev)
212 {
213 const auto& plev = pc.GetParticles(lev);
214 Vector<std::pair<int, int> > grid_tile_ids;
215 Vector<const typename PC::ParticleTileType*> ptile_ptrs;
216 for (auto& kv : plev)
217 {
218 grid_tile_ids.push_back(kv.first);
219 ptile_ptrs.push_back(&(kv.second));
220 }
221
222 // We avoid ParIter here to allow reductions after grids change but prior to calling Redistribute()
223#ifdef AMREX_USE_OMP
224#pragma omp parallel for if (!system::regtest_reduction) reduction(+:sm)
225#endif
226 for (int pmap_it = 0; pmap_it < static_cast<int>(ptile_ptrs.size()); ++pmap_it)
227 {
228 const auto& tile = plev.at(grid_tile_ids[pmap_it]);
229 const auto np = tile.numParticles();
230 const auto& ptd = tile.getConstParticleTileData();
231 for (int i = 0; i < np; ++i) {
232 sm += particle_detail::call_f(f, ptd, i);
233 }
234 }
235 }
236 }
237
238 return sm;
239}
240
280template <class PC, class F, std::enable_if_t<IsParticleContainer<PC>::value, int> foo = 0>
281auto
282ReduceMax (PC const& pc, F&& f)
283 -> decltype(particle_detail::call_f(f, typename PC::ConstPTDType(), int()))
284{
285 return ReduceMax(pc, 0, pc.finestLevel(), std::forward<F>(f));
286}
287
329template <class PC, class F, std::enable_if_t<IsParticleContainer<PC>::value, int> foo = 0>
330auto
331ReduceMax (PC const& pc, int lev, F&& f)
332 -> decltype(particle_detail::call_f(f, typename PC::ConstPTDType(), int()))
333{
334 return ReduceMax(pc, lev, lev, std::forward<F>(f));
335}
336
378template <class PC, class F, std::enable_if_t<IsParticleContainer<PC>::value, int> foo = 0>
379auto
380ReduceMax (PC const& pc, int lev_min, int lev_max, F const& f)
381 -> decltype(particle_detail::call_f(f, typename PC::ConstPTDType(), int()))
382{
383 using value_type = decltype(particle_detail::call_f(f, typename PC::ConstPTDType(), int()));
384 constexpr value_type value_lowest = std::numeric_limits<value_type>::lowest();
385 value_type r = value_lowest;
386
387#ifdef AMREX_USE_GPU
388 if (Gpu::inLaunchRegion())
389 {
390 ReduceOps<ReduceOpMax> reduce_op;
391 ReduceData<value_type> reduce_data(reduce_op);
392 using ReduceTuple = typename decltype(reduce_data)::Type;
393
394 for (int lev = lev_min; lev <= lev_max; ++lev)
395 {
396 const auto& plev = pc.GetParticles(lev);
397 for (const auto& kv : plev)
398 {
399 const auto& tile = plev.at(kv.first);
400 const auto np = tile.numParticles();
401 const auto& ptd = tile.getConstParticleTileData();
402 reduce_op.eval(np, reduce_data,
403 [=] AMREX_GPU_DEVICE (const int i) -> ReduceTuple {
404 return particle_detail::call_f(f, ptd, i);
405 });
406 }
407 }
408
409 ReduceTuple hv = reduce_data.value(reduce_op);
410 r = amrex::get<0>(hv);
411 }
412 else
413#endif
414 {
415 for (int lev = lev_min; lev <= lev_max; ++lev)
416 {
417 const auto& plev = pc.GetParticles(lev);
418 Vector<std::pair<int, int> > grid_tile_ids;
419 Vector<const typename PC::ParticleTileType*> ptile_ptrs;
420 for (auto& kv : plev)
421 {
422 grid_tile_ids.push_back(kv.first);
423 ptile_ptrs.push_back(&(kv.second));
424 }
425
426 // We avoid ParIter here to allow reductions after grids change but prior to calling Redistribute()
427#ifdef AMREX_USE_OMP
428#pragma omp parallel for if (!system::regtest_reduction) reduction(max:r)
429#endif
430 for (int pmap_it = 0; pmap_it < static_cast<int>(ptile_ptrs.size()); ++pmap_it)
431 {
432 const auto& tile = plev.at(grid_tile_ids[pmap_it]);
433 const auto np = tile.numParticles();
434 const auto& ptd = tile.getConstParticleTileData();
435 for (int i = 0; i < np; ++i) {
436 r = std::max(r, particle_detail::call_f(f, ptd, i));
437 }
438 }
439 }
440 }
441
442 return r;
443}
444
484template <class PC, class F, std::enable_if_t<IsParticleContainer<PC>::value, int> foo = 0>
485auto
486ReduceMin (PC const& pc, F&& f)
487 -> decltype(particle_detail::call_f(f, typename PC::ConstPTDType(), int()))
488{
489 return ReduceMin(pc, 0, pc.finestLevel(), std::forward<F>(f));
490}
491
532template <class PC, class F, std::enable_if_t<IsParticleContainer<PC>::value, int> foo = 0>
533auto
534ReduceMin (PC const& pc, int lev, F&& f)
535 -> decltype(particle_detail::call_f(f, typename PC::ConstPTDType(), int()))
536{
537 return ReduceMin(pc, lev, lev, std::forward<F>(f));
538}
539
581template <class PC, class F, std::enable_if_t<IsParticleContainer<PC>::value, int> foo = 0>
582auto
583ReduceMin (PC const& pc, int lev_min, int lev_max, F const& f)
584 -> decltype(particle_detail::call_f(f, typename PC::ConstPTDType(), int()))
585{
586 using value_type = decltype(particle_detail::call_f(f, typename PC::ConstPTDType(), int()));
587 constexpr value_type value_max = std::numeric_limits<value_type>::max();
588 value_type r = value_max;
589
590#ifdef AMREX_USE_GPU
591 if (Gpu::inLaunchRegion())
592 {
593 ReduceOps<ReduceOpMin> reduce_op;
594 ReduceData<value_type> reduce_data(reduce_op);
595 using ReduceTuple = typename decltype(reduce_data)::Type;
596
597 for (int lev = lev_min; lev <= lev_max; ++lev)
598 {
599 const auto& plev = pc.GetParticles(lev);
600 for (const auto& kv : plev)
601 {
602 const auto& tile = plev.at(kv.first);
603 const auto np = tile.numParticles();
604 const auto& ptd = tile.getConstParticleTileData();
605 reduce_op.eval(np, reduce_data,
606 [=] AMREX_GPU_DEVICE (const int i) -> ReduceTuple {
607 return particle_detail::call_f(f, ptd, i);
608 });
609 }
610 }
611
612 ReduceTuple hv = reduce_data.value(reduce_op);
613 r = amrex::get<0>(hv);
614 }
615 else
616#endif
617 {
618 for (int lev = lev_min; lev <= lev_max; ++lev)
619 {
620 const auto& plev = pc.GetParticles(lev);
621 Vector<std::pair<int, int> > grid_tile_ids;
622 Vector<const typename PC::ParticleTileType*> ptile_ptrs;
623 for (auto& kv : plev)
624 {
625 grid_tile_ids.push_back(kv.first);
626 ptile_ptrs.push_back(&(kv.second));
627 }
628
629 // We avoid ParIter here to allow reductions after grids change but prior to calling Redistribute()
630#ifdef AMREX_USE_OMP
631#pragma omp parallel for if (!system::regtest_reduction) reduction(min:r)
632#endif
633 for (int pmap_it = 0; pmap_it < static_cast<int>(ptile_ptrs.size()); ++pmap_it)
634 {
635 const auto& tile = plev.at(grid_tile_ids[pmap_it]);
636 const auto np = tile.numParticles();
637 const auto& ptd = tile.getConstParticleTileData();
638 for (int i = 0; i < np; ++i) {
639 r = std::min(r, particle_detail::call_f(f, ptd, i));
640 }
641 }
642 }
643 }
644
645 return r;
646}
647
687template <class PC, class F, std::enable_if_t<IsParticleContainer<PC>::value, int> foo = 0>
688bool
689ReduceLogicalAnd (PC const& pc, F&& f)
690{
691 return ReduceLogicalAnd(pc, 0, pc.finestLevel(), std::forward<F>(f));
692}
693
734template <class PC, class F, std::enable_if_t<IsParticleContainer<PC>::value, int> foo = 0>
735bool
736ReduceLogicalAnd (PC const& pc, int lev, F&& f)
737{
738 return ReduceLogicalAnd(pc, lev, lev, std::forward<F>(f));
739}
740
782template <class PC, class F, std::enable_if_t<IsParticleContainer<PC>::value, int> foo = 0>
783bool
784ReduceLogicalAnd (PC const& pc, int lev_min, int lev_max, F const& f)
785{
786 int r = true;
787
788#ifdef AMREX_USE_GPU
789 if (Gpu::inLaunchRegion())
790 {
791 ReduceOps<ReduceOpLogicalAnd> reduce_op;
792 ReduceData<int> reduce_data(reduce_op);
793 using ReduceTuple = typename decltype(reduce_data)::Type;
794
795 for (int lev = lev_min; lev <= lev_max; ++lev)
796 {
797 const auto& plev = pc.GetParticles(lev);
798 for (const auto& kv : plev)
799 {
800 const auto& tile = plev.at(kv.first);
801 const auto np = tile.numParticles();
802 const auto& ptd = tile.getConstParticleTileData();
803 reduce_op.eval(np, reduce_data,
804 [=] AMREX_GPU_DEVICE (const int i) -> ReduceTuple {
805 return particle_detail::call_f(f, ptd, i);
806 });
807 }
808 }
809
810 ReduceTuple hv = reduce_data.value(reduce_op);
811 r = amrex::get<0>(hv);
812 }
813 else
814#endif
815 {
816 for (int lev = lev_min; lev <= lev_max; ++lev)
817 {
818 const auto& plev = pc.GetParticles(lev);
819 Vector<std::pair<int, int> > grid_tile_ids;
820 Vector<const typename PC::ParticleTileType*> ptile_ptrs;
821 for (auto& kv : plev)
822 {
823 grid_tile_ids.push_back(kv.first);
824 ptile_ptrs.push_back(&(kv.second));
825 }
826
827 // We avoid ParIter here to allow reductions after grids change but prior to calling Redistribute()
828#ifdef AMREX_USE_OMP
829#pragma omp parallel for if (!system::regtest_reduction) reduction(&&:r)
830#endif
831 for (int pmap_it = 0; pmap_it < static_cast<int>(ptile_ptrs.size()); ++pmap_it)
832 {
833 const auto& tile = plev.at(grid_tile_ids[pmap_it]);
834 const auto np = tile.numParticles();
835 const auto& ptd = tile.getConstParticleTileData();
836 for (int i = 0; i < np; ++i) {
837 r = r && particle_detail::call_f(f, ptd, i);
838 }
839 }
840 }
841 }
842
843 return r;
844}
845
885template <class PC, class F, std::enable_if_t<IsParticleContainer<PC>::value, int> foo = 0>
886bool
887ReduceLogicalOr (PC const& pc, F&& f)
888{
889 return ReduceLogicalOr(pc, 0, pc.finestLevel(), std::forward<F>(f));
890}
891
932template <class PC, class F, std::enable_if_t<IsParticleContainer<PC>::value, int> foo = 0>
933bool
934ReduceLogicalOr (PC const& pc, int lev, F&& f)
935{
936 return ReduceLogicalOr(pc, lev, lev, std::forward<F>(f));
937}
938
980template <class PC, class F, std::enable_if_t<IsParticleContainer<PC>::value, int> foo = 0>
981bool
982ReduceLogicalOr (PC const& pc, int lev_min, int lev_max, F const& f)
983{
984 int r = false;
985
986#ifdef AMREX_USE_GPU
987 if (Gpu::inLaunchRegion())
988 {
989 ReduceOps<ReduceOpLogicalOr> reduce_op;
990 ReduceData<int> reduce_data(reduce_op);
991 using ReduceTuple = typename decltype(reduce_data)::Type;
992
993 for (int lev = lev_min; lev <= lev_max; ++lev)
994 {
995 const auto& plev = pc.GetParticles(lev);
996 for (const auto& kv : plev)
997 {
998 const auto& tile = plev.at(kv.first);
999 const auto np = tile.numParticles();
1000 const auto& ptd = tile.getConstParticleTileData();
1001 reduce_op.eval(np, reduce_data,
1002 [=] AMREX_GPU_DEVICE (const int i) -> ReduceTuple
1003 {
1004 return particle_detail::call_f(f, ptd, i);
1005 });
1006 }
1007 }
1008
1009 ReduceTuple hv = reduce_data.value(reduce_op);
1010 r = amrex::get<0>(hv);
1011 }
1012 else
1013#endif
1014 {
1015 for (int lev = lev_min; lev <= lev_max; ++lev)
1016 {
1017 const auto& plev = pc.GetParticles(lev);
1018 Vector<std::pair<int, int> > grid_tile_ids;
1019 Vector<const typename PC::ParticleTileType*> ptile_ptrs;
1020 for (auto& kv : plev)
1021 {
1022 grid_tile_ids.push_back(kv.first);
1023 ptile_ptrs.push_back(&(kv.second));
1024 }
1025
1026 // We avoid ParIter here to allow reductions after grids change but prior to calling Redistribute()
1027#ifdef AMREX_USE_OMP
1028#pragma omp parallel for if (!system::regtest_reduction) reduction(||:r)
1029#endif
1030 for (int pmap_it = 0; pmap_it < static_cast<int>(ptile_ptrs.size()); ++pmap_it)
1031 {
1032 const auto& tile = plev.at(grid_tile_ids[pmap_it]);
1033 const auto np = tile.numParticles();
1034 const auto& ptd = tile.getConstParticleTileData();
1035 for (int i = 0; i < np; ++i) {
1036 r = r || particle_detail::call_f(f, ptd, i);
1037 }
1038 }
1039 }
1040 }
1041
1042 return r;
1043}
1044
1107template <class RD, class PC, class F, class ReduceOps,
1108 std::enable_if_t<IsParticleContainer<PC>::value, int> foo = 0>
1109typename RD::Type
1110ParticleReduce (PC const& pc, F&& f, ReduceOps& reduce_ops)
1111{
1112 return ParticleReduce<RD>(pc, 0, pc.finestLevel(), std::forward<F>(f), reduce_ops);
1113}
1114
1178template <class RD, class PC, class F, class ReduceOps,
1179 std::enable_if_t<IsParticleContainer<PC>::value, int> foo = 0>
1180typename RD::Type
1181ParticleReduce (PC const& pc, int lev, F&& f, ReduceOps& reduce_ops)
1182{
1183 return ParticleReduce<RD>(pc, lev, lev, std::forward<F>(f), reduce_ops);
1184}
1185
1250template <class RD, class PC, class F, class ReduceOps,
1251 std::enable_if_t<IsParticleContainer<PC>::value, int> foo = 0>
1252typename RD::Type
1253ParticleReduce (PC const& pc, int lev_min, int lev_max, F const& f, ReduceOps& reduce_ops)
1254{
1255 RD reduce_data(reduce_ops);
1256 for (int lev = lev_min; lev <= lev_max; ++lev) {
1257 const auto& plev = pc.GetParticles(lev);
1258 Vector<std::pair<int, int> > grid_tile_ids;
1259 Vector<const typename PC::ParticleTileType*> ptile_ptrs;
1260 for (auto& kv : plev)
1261 {
1262 grid_tile_ids.push_back(kv.first);
1263 ptile_ptrs.push_back(&(kv.second));
1264 }
1265
1266 // We avoid ParIter here to allow reductions after grids change but prior to calling Redistribute()
1267#if !defined(AMREX_USE_GPU) && defined(AMREX_USE_OMP)
1268#pragma omp parallel for
1269#endif
1270 for (int pmap_it = 0; pmap_it < static_cast<int>(ptile_ptrs.size()); ++pmap_it)
1271 {
1272 const auto& tile = plev.at(grid_tile_ids[pmap_it]);
1273 const auto np = tile.numParticles();
1274 const auto& ptd = tile.getConstParticleTileData();
1275 reduce_ops.eval(np, reduce_data,
1276 [=] AMREX_GPU_DEVICE (const int i) noexcept
1277 {
1278 return particle_detail::call_f(f, ptd, i);
1279 });
1280 }
1281 }
1282 return reduce_data.value(reduce_ops);
1283}
1284}
1285#endif
AMReX_Box.H
AMREX_FORCE_INLINE
#define AMREX_FORCE_INLINE
Definition AMReX_Extension.H:119
AMREX_GPU_DEVICE
#define AMREX_GPU_DEVICE
Definition AMReX_GpuQualifiers.H:18
AMREX_GPU_HOST_DEVICE
#define AMREX_GPU_HOST_DEVICE
Definition AMReX_GpuQualifiers.H:20
AMReX_GpuUtility.H
AMReX_Gpu.H
AMReX_IntVect.H
AMReX_ParticleUtil.H
AMReX_Print.H
AMReX_TypeTraits.H
AMReX_Vector.H
amrex::ReduceData
Definition AMReX_Reduce.H:249
amrex::ReduceData::value
Type value()
Definition AMReX_Reduce.H:281
amrex::ReduceOps
Definition AMReX_Reduce.H:364
amrex::ReduceOps::eval
std::enable_if_t< IsFabArray< MF >::value > eval(MF const &mf, IntVect const &nghost, D &reduce_data, F &&f)
Definition AMReX_Reduce.H:433
amrex::Vector
This class is a thin wrapper around std::vector. Unlike vector, Vector::operator[] provides bound che...
Definition AMReX_Vector.H:28
amrex::Vector::size
Long size() const noexcept
Definition AMReX_Vector.H:53
amrex::Gpu::inLaunchRegion
bool inLaunchRegion() noexcept
Definition AMReX_GpuControl.H:92
amrex::particle_detail::call_f
__host__ __device__ auto call_f(F const &f, const PTDType< T_ParticleType, NAR, NAI > &p, const int i, Array4< T > const &fabarr, GpuArray< Real, 3 > const &plo, GpuArray< Real, 3 > const &dxi) noexcept
Definition AMReX_ParticleMesh.H:16
amrex
Definition AMReX_Amr.cpp:49
amrex::ReduceMax
FAB::value_type ReduceMax(FabArray< FAB > const &fa, int nghost, F &&f)
Definition AMReX_FabArrayUtility.H:531
amrex::ReduceLogicalAnd
bool ReduceLogicalAnd(FabArray< FAB > const &fa, int nghost, F &&f)
Definition AMReX_FabArrayUtility.H:758
amrex::ReduceLogicalOr
bool ReduceLogicalOr(FabArray< FAB > const &fa, int nghost, F &&f)
Definition AMReX_FabArrayUtility.H:905
amrex::ParticleReduce
RD::Type ParticleReduce(PC const &pc, F &&f, ReduceOps &reduce_ops)
A general reduction method for the particles in a ParticleContainer that can run on either CPUs or GP...
Definition AMReX_ParticleReduce.H:1110
amrex::ReduceMin
FAB::value_type ReduceMin(FabArray< FAB > const &fa, int nghost, F &&f)
Definition AMReX_FabArrayUtility.H:305
amrex::ReduceSum
FAB::value_type ReduceSum(FabArray< FAB > const &fa, int nghost, F &&f)
Definition AMReX_FabArrayUtility.H:16
amrex::ConstParticleTileData
Definition AMReX_ParticleTile.H:513
amrex::IsCallable
Test if a given type T is callable with arguments of type Args...
Definition AMReX_TypeTraits.H:209