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
20namespace particle_detail {
21
22template <typename F, typename T_ParticleType, int NAR, int NAI>
23AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
24auto call_f (F const& f,
25 const ConstParticleTileData<T_ParticleType, NAR, NAI>& p,
26 const int i) noexcept
27{
28 if constexpr ( ! T_ParticleType::is_soa_particle &&
29 IsCallable<F, T_ParticleType const&>::value) {
30 return f(p.m_aos[i]);
31 } else if constexpr (IsCallable<F, decltype(p.getSuperParticle(i))>::value) {
32 return f(p.getSuperParticle(i));
33 } else {
34 return f(p, i);
35 }
36}
37
38template <typename F, template<class, class> class PTDType, class RType, class IType>
39AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
40auto call_f (F const& f,
41 const PTDType<RType, IType>& p,
42 const int i) noexcept
43{
44 return f(p, i);
45}
46}
48
89template <class PC, class F, std::enable_if_t<IsParticleContainer<PC>::value, int> foo = 0>
90auto
91ReduceSum (PC const& pc, F&& f)
92 -> decltype(particle_detail::call_f(f, typename PC::ConstPTDType(), int()))
93{
94 return ReduceSum(pc, 0, pc.finestLevel(), std::forward<F>(f));
95}
96
137template <class PC, class F, std::enable_if_t<IsParticleContainer<PC>::value, int> foo = 0>
138auto
139ReduceSum (PC const& pc, int lev, F&& f)
140 -> decltype(particle_detail::call_f(f, typename PC::ConstPTDType(), int()))
141{
142 return ReduceSum(pc, lev, lev, std::forward<F>(f));
143}
144
186template <class PC, class F, std::enable_if_t<IsParticleContainer<PC>::value, int> foo = 0>
187auto
188ReduceSum (PC const& pc, int lev_min, int lev_max, F const& f)
189 -> decltype(particle_detail::call_f(f, typename PC::ConstPTDType(), int()))
190{
191 using value_type = decltype(particle_detail::call_f(f, typename PC::ConstPTDType(), int()));
192 value_type sm = 0;
193
194#ifdef AMREX_USE_GPU
195 if (Gpu::inLaunchRegion())
196 {
197 ReduceOps<ReduceOpSum> reduce_op;
198 ReduceData<value_type> reduce_data(reduce_op);
199 using ReduceTuple = typename decltype(reduce_data)::Type;
200
201 for (int lev = lev_min; lev <= lev_max; ++lev)
202 {
203 const auto& plev = pc.GetParticles(lev);
204 for (const auto& kv : plev)
205 {
206 const auto& tile = plev.at(kv.first);
207 const auto np = tile.numParticles();
208 const auto& ptd = tile.getConstParticleTileData();
209 reduce_op.eval(np, reduce_data,
210 [=] AMREX_GPU_DEVICE (const int i) -> ReduceTuple {
211 return particle_detail::call_f(f, ptd, i);
212 });
213 }
214 }
215
216 ReduceTuple hv = reduce_data.value(reduce_op);
217 sm = amrex::get<0>(hv);
218 }
219 else
220#endif
221 {
222 for (int lev = lev_min; lev <= lev_max; ++lev)
223 {
224 const auto& plev = pc.GetParticles(lev);
225 Vector<std::pair<int, int> > grid_tile_ids;
226 Vector<const typename PC::ParticleTileType*> ptile_ptrs;
227 for (auto& kv : plev)
228 {
229 grid_tile_ids.push_back(kv.first);
230 ptile_ptrs.push_back(&(kv.second));
231 }
232
233 // We avoid ParIter here to allow reductions after grids change but prior to calling Redistribute()
234#ifdef AMREX_USE_OMP
235#pragma omp parallel for if (!system::regtest_reduction) reduction(+:sm)
236#endif
237 for (int pmap_it = 0; pmap_it < static_cast<int>(ptile_ptrs.size()); ++pmap_it)
238 {
239 const auto& tile = plev.at(grid_tile_ids[pmap_it]);
240 const auto np = tile.numParticles();
241 const auto& ptd = tile.getConstParticleTileData();
242 for (int i = 0; i < np; ++i) {
243 sm += particle_detail::call_f(f, ptd, i);
244 }
245 }
246 }
247 }
248
249 return sm;
250}
251
291template <class PC, class F, std::enable_if_t<IsParticleContainer<PC>::value, int> foo = 0>
292auto
293ReduceMax (PC const& pc, F&& f)
294 -> decltype(particle_detail::call_f(f, typename PC::ConstPTDType(), int()))
295{
296 return ReduceMax(pc, 0, pc.finestLevel(), std::forward<F>(f));
297}
298
340template <class PC, class F, std::enable_if_t<IsParticleContainer<PC>::value, int> foo = 0>
341auto
342ReduceMax (PC const& pc, int lev, F&& f)
343 -> decltype(particle_detail::call_f(f, typename PC::ConstPTDType(), int()))
344{
345 return ReduceMax(pc, lev, lev, std::forward<F>(f));
346}
347
389template <class PC, class F, std::enable_if_t<IsParticleContainer<PC>::value, int> foo = 0>
390auto
391ReduceMax (PC const& pc, int lev_min, int lev_max, F const& f)
392 -> decltype(particle_detail::call_f(f, typename PC::ConstPTDType(), int()))
393{
394 using value_type = decltype(particle_detail::call_f(f, typename PC::ConstPTDType(), int()));
395 constexpr value_type value_lowest = std::numeric_limits<value_type>::lowest();
396 value_type r = value_lowest;
397
398#ifdef AMREX_USE_GPU
399 if (Gpu::inLaunchRegion())
400 {
401 ReduceOps<ReduceOpMax> reduce_op;
402 ReduceData<value_type> reduce_data(reduce_op);
403 using ReduceTuple = typename decltype(reduce_data)::Type;
404
405 for (int lev = lev_min; lev <= lev_max; ++lev)
406 {
407 const auto& plev = pc.GetParticles(lev);
408 for (const auto& kv : plev)
409 {
410 const auto& tile = plev.at(kv.first);
411 const auto np = tile.numParticles();
412 const auto& ptd = tile.getConstParticleTileData();
413 reduce_op.eval(np, reduce_data,
414 [=] AMREX_GPU_DEVICE (const int i) -> ReduceTuple {
415 return particle_detail::call_f(f, ptd, i);
416 });
417 }
418 }
419
420 ReduceTuple hv = reduce_data.value(reduce_op);
421 r = amrex::get<0>(hv);
422 }
423 else
424#endif
425 {
426 for (int lev = lev_min; lev <= lev_max; ++lev)
427 {
428 const auto& plev = pc.GetParticles(lev);
429 Vector<std::pair<int, int> > grid_tile_ids;
430 Vector<const typename PC::ParticleTileType*> ptile_ptrs;
431 for (auto& kv : plev)
432 {
433 grid_tile_ids.push_back(kv.first);
434 ptile_ptrs.push_back(&(kv.second));
435 }
436
437 // We avoid ParIter here to allow reductions after grids change but prior to calling Redistribute()
438#ifdef AMREX_USE_OMP
439#pragma omp parallel for if (!system::regtest_reduction) reduction(max:r)
440#endif
441 for (int pmap_it = 0; pmap_it < static_cast<int>(ptile_ptrs.size()); ++pmap_it)
442 {
443 const auto& tile = plev.at(grid_tile_ids[pmap_it]);
444 const auto np = tile.numParticles();
445 const auto& ptd = tile.getConstParticleTileData();
446 for (int i = 0; i < np; ++i) {
447 r = std::max(r, particle_detail::call_f(f, ptd, i));
448 }
449 }
450 }
451 }
452
453 return r;
454}
455
495template <class PC, class F, std::enable_if_t<IsParticleContainer<PC>::value, int> foo = 0>
496auto
497ReduceMin (PC const& pc, F&& f)
498 -> decltype(particle_detail::call_f(f, typename PC::ConstPTDType(), int()))
499{
500 return ReduceMin(pc, 0, pc.finestLevel(), std::forward<F>(f));
501}
502
543template <class PC, class F, std::enable_if_t<IsParticleContainer<PC>::value, int> foo = 0>
544auto
545ReduceMin (PC const& pc, int lev, F&& f)
546 -> decltype(particle_detail::call_f(f, typename PC::ConstPTDType(), int()))
547{
548 return ReduceMin(pc, lev, lev, std::forward<F>(f));
549}
550
592template <class PC, class F, std::enable_if_t<IsParticleContainer<PC>::value, int> foo = 0>
593auto
594ReduceMin (PC const& pc, int lev_min, int lev_max, F const& f)
595 -> decltype(particle_detail::call_f(f, typename PC::ConstPTDType(), int()))
596{
597 using value_type = decltype(particle_detail::call_f(f, typename PC::ConstPTDType(), int()));
598 constexpr value_type value_max = std::numeric_limits<value_type>::max();
599 value_type r = value_max;
600
601#ifdef AMREX_USE_GPU
602 if (Gpu::inLaunchRegion())
603 {
604 ReduceOps<ReduceOpMin> reduce_op;
605 ReduceData<value_type> reduce_data(reduce_op);
606 using ReduceTuple = typename decltype(reduce_data)::Type;
607
608 for (int lev = lev_min; lev <= lev_max; ++lev)
609 {
610 const auto& plev = pc.GetParticles(lev);
611 for (const auto& kv : plev)
612 {
613 const auto& tile = plev.at(kv.first);
614 const auto np = tile.numParticles();
615 const auto& ptd = tile.getConstParticleTileData();
616 reduce_op.eval(np, reduce_data,
617 [=] AMREX_GPU_DEVICE (const int i) -> ReduceTuple {
618 return particle_detail::call_f(f, ptd, i);
619 });
620 }
621 }
622
623 ReduceTuple hv = reduce_data.value(reduce_op);
624 r = amrex::get<0>(hv);
625 }
626 else
627#endif
628 {
629 for (int lev = lev_min; lev <= lev_max; ++lev)
630 {
631 const auto& plev = pc.GetParticles(lev);
632 Vector<std::pair<int, int> > grid_tile_ids;
633 Vector<const typename PC::ParticleTileType*> ptile_ptrs;
634 for (auto& kv : plev)
635 {
636 grid_tile_ids.push_back(kv.first);
637 ptile_ptrs.push_back(&(kv.second));
638 }
639
640 // We avoid ParIter here to allow reductions after grids change but prior to calling Redistribute()
641#ifdef AMREX_USE_OMP
642#pragma omp parallel for if (!system::regtest_reduction) reduction(min:r)
643#endif
644 for (int pmap_it = 0; pmap_it < static_cast<int>(ptile_ptrs.size()); ++pmap_it)
645 {
646 const auto& tile = plev.at(grid_tile_ids[pmap_it]);
647 const auto np = tile.numParticles();
648 const auto& ptd = tile.getConstParticleTileData();
649 for (int i = 0; i < np; ++i) {
650 r = std::min(r, particle_detail::call_f(f, ptd, i));
651 }
652 }
653 }
654 }
655
656 return r;
657}
658
698template <class PC, class F, std::enable_if_t<IsParticleContainer<PC>::value, int> foo = 0>
699bool
700ReduceLogicalAnd (PC const& pc, F&& f)
701{
702 return ReduceLogicalAnd(pc, 0, pc.finestLevel(), std::forward<F>(f));
703}
704
745template <class PC, class F, std::enable_if_t<IsParticleContainer<PC>::value, int> foo = 0>
746bool
747ReduceLogicalAnd (PC const& pc, int lev, F&& f)
748{
749 return ReduceLogicalAnd(pc, lev, lev, std::forward<F>(f));
750}
751
793template <class PC, class F, std::enable_if_t<IsParticleContainer<PC>::value, int> foo = 0>
794bool
795ReduceLogicalAnd (PC const& pc, int lev_min, int lev_max, F const& f)
796{
797 int r = true;
798
799#ifdef AMREX_USE_GPU
800 if (Gpu::inLaunchRegion())
801 {
802 ReduceOps<ReduceOpLogicalAnd> reduce_op;
803 ReduceData<int> reduce_data(reduce_op);
804 using ReduceTuple = typename decltype(reduce_data)::Type;
805
806 for (int lev = lev_min; lev <= lev_max; ++lev)
807 {
808 const auto& plev = pc.GetParticles(lev);
809 for (const auto& kv : plev)
810 {
811 const auto& tile = plev.at(kv.first);
812 const auto np = tile.numParticles();
813 const auto& ptd = tile.getConstParticleTileData();
814 reduce_op.eval(np, reduce_data,
815 [=] AMREX_GPU_DEVICE (const int i) -> ReduceTuple {
816 return particle_detail::call_f(f, ptd, i);
817 });
818 }
819 }
820
821 ReduceTuple hv = reduce_data.value(reduce_op);
822 r = amrex::get<0>(hv);
823 }
824 else
825#endif
826 {
827 for (int lev = lev_min; lev <= lev_max; ++lev)
828 {
829 const auto& plev = pc.GetParticles(lev);
830 Vector<std::pair<int, int> > grid_tile_ids;
831 Vector<const typename PC::ParticleTileType*> ptile_ptrs;
832 for (auto& kv : plev)
833 {
834 grid_tile_ids.push_back(kv.first);
835 ptile_ptrs.push_back(&(kv.second));
836 }
837
838 // We avoid ParIter here to allow reductions after grids change but prior to calling Redistribute()
839#ifdef AMREX_USE_OMP
840#pragma omp parallel for if (!system::regtest_reduction) reduction(&&:r)
841#endif
842 for (int pmap_it = 0; pmap_it < static_cast<int>(ptile_ptrs.size()); ++pmap_it)
843 {
844 const auto& tile = plev.at(grid_tile_ids[pmap_it]);
845 const auto np = tile.numParticles();
846 const auto& ptd = tile.getConstParticleTileData();
847 for (int i = 0; i < np; ++i) {
848 r = r && particle_detail::call_f(f, ptd, i);
849 }
850 }
851 }
852 }
853
854 return r;
855}
856
896template <class PC, class F, std::enable_if_t<IsParticleContainer<PC>::value, int> foo = 0>
897bool
898ReduceLogicalOr (PC const& pc, F&& f)
899{
900 return ReduceLogicalOr(pc, 0, pc.finestLevel(), std::forward<F>(f));
901}
902
943template <class PC, class F, std::enable_if_t<IsParticleContainer<PC>::value, int> foo = 0>
944bool
945ReduceLogicalOr (PC const& pc, int lev, F&& f)
946{
947 return ReduceLogicalOr(pc, lev, lev, std::forward<F>(f));
948}
949
991template <class PC, class F, std::enable_if_t<IsParticleContainer<PC>::value, int> foo = 0>
992bool
993ReduceLogicalOr (PC const& pc, int lev_min, int lev_max, F const& f)
994{
995 int r = false;
996
997#ifdef AMREX_USE_GPU
998 if (Gpu::inLaunchRegion())
999 {
1000 ReduceOps<ReduceOpLogicalOr> reduce_op;
1001 ReduceData<int> reduce_data(reduce_op);
1002 using ReduceTuple = typename decltype(reduce_data)::Type;
1003
1004 for (int lev = lev_min; lev <= lev_max; ++lev)
1005 {
1006 const auto& plev = pc.GetParticles(lev);
1007 for (const auto& kv : plev)
1008 {
1009 const auto& tile = plev.at(kv.first);
1010 const auto np = tile.numParticles();
1011 const auto& ptd = tile.getConstParticleTileData();
1012 reduce_op.eval(np, reduce_data,
1013 [=] AMREX_GPU_DEVICE (const int i) -> ReduceTuple
1014 {
1015 return particle_detail::call_f(f, ptd, i);
1016 });
1017 }
1018 }
1019
1020 ReduceTuple hv = reduce_data.value(reduce_op);
1021 r = amrex::get<0>(hv);
1022 }
1023 else
1024#endif
1025 {
1026 for (int lev = lev_min; lev <= lev_max; ++lev)
1027 {
1028 const auto& plev = pc.GetParticles(lev);
1029 Vector<std::pair<int, int> > grid_tile_ids;
1030 Vector<const typename PC::ParticleTileType*> ptile_ptrs;
1031 for (auto& kv : plev)
1032 {
1033 grid_tile_ids.push_back(kv.first);
1034 ptile_ptrs.push_back(&(kv.second));
1035 }
1036
1037 // We avoid ParIter here to allow reductions after grids change but prior to calling Redistribute()
1038#ifdef AMREX_USE_OMP
1039#pragma omp parallel for if (!system::regtest_reduction) reduction(||:r)
1040#endif
1041 for (int pmap_it = 0; pmap_it < static_cast<int>(ptile_ptrs.size()); ++pmap_it)
1042 {
1043 const auto& tile = plev.at(grid_tile_ids[pmap_it]);
1044 const auto np = tile.numParticles();
1045 const auto& ptd = tile.getConstParticleTileData();
1046 for (int i = 0; i < np; ++i) {
1047 r = r || particle_detail::call_f(f, ptd, i);
1048 }
1049 }
1050 }
1051 }
1052
1053 return r;
1054}
1055
1118template <class RD, class PC, class F, class ReduceOps,
1119 std::enable_if_t<IsParticleContainer<PC>::value, int> foo = 0>
1120typename RD::Type
1121ParticleReduce (PC const& pc, F&& f, ReduceOps& reduce_ops)
1122{
1123 return ParticleReduce<RD>(pc, 0, pc.finestLevel(), std::forward<F>(f), reduce_ops);
1124}
1125
1189template <class RD, class PC, class F, class ReduceOps,
1190 std::enable_if_t<IsParticleContainer<PC>::value, int> foo = 0>
1191typename RD::Type
1192ParticleReduce (PC const& pc, int lev, F&& f, ReduceOps& reduce_ops)
1193{
1194 return ParticleReduce<RD>(pc, lev, lev, std::forward<F>(f), reduce_ops);
1195}
1196
1261template <class RD, class PC, class F, class ReduceOps,
1262 std::enable_if_t<IsParticleContainer<PC>::value, int> foo = 0>
1263typename RD::Type
1264ParticleReduce (PC const& pc, int lev_min, int lev_max, F const& f, ReduceOps& reduce_ops)
1265{
1266 RD reduce_data(reduce_ops);
1267 for (int lev = lev_min; lev <= lev_max; ++lev) {
1268 const auto& plev = pc.GetParticles(lev);
1269 Vector<std::pair<int, int> > grid_tile_ids;
1270 Vector<const typename PC::ParticleTileType*> ptile_ptrs;
1271 for (auto& kv : plev)
1272 {
1273 grid_tile_ids.push_back(kv.first);
1274 ptile_ptrs.push_back(&(kv.second));
1275 }
1276
1277 // We avoid ParIter here to allow reductions after grids change but prior to calling Redistribute()
1278#if !defined(AMREX_USE_GPU) && defined(AMREX_USE_OMP)
1279#pragma omp parallel for
1280#endif
1281 for (int pmap_it = 0; pmap_it < static_cast<int>(ptile_ptrs.size()); ++pmap_it)
1282 {
1283 const auto& tile = plev.at(grid_tile_ids[pmap_it]);
1284 const auto np = tile.numParticles();
1285 const auto& ptd = tile.getConstParticleTileData();
1286 reduce_ops.eval(np, reduce_data,
1287 [=] AMREX_GPU_DEVICE (const int i) noexcept
1288 {
1289 return particle_detail::call_f(f, ptd, i);
1290 });
1291 }
1292 }
1293 return reduce_data.value(reduce_ops);
1294}
1295}
1296#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:453
amrex::ReduceData::value
Type value()
Definition AMReX_Reduce.H:488
amrex::ReduceOps
Definition AMReX_Reduce.H:612
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:746
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
Definition AMReX_Amr.cpp:49
amrex::ReduceMax
FAB::value_type ReduceMax(FabArray< FAB > const &fa, int nghost, F &&f)
Definition AMReX_FabArrayUtility.H:555
amrex::ReduceLogicalAnd
bool ReduceLogicalAnd(FabArray< FAB > const &fa, int nghost, F &&f)
Definition AMReX_FabArrayUtility.H:794
amrex::ReduceLogicalOr
bool ReduceLogicalOr(FabArray< FAB > const &fa, int nghost, F &&f)
Definition AMReX_FabArrayUtility.H:949
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:1121
amrex::ReduceMin
FAB::value_type ReduceMin(FabArray< FAB > const &fa, int nghost, F &&f)
Definition AMReX_FabArrayUtility.H:317
amrex::ReduceSum
FAB::value_type ReduceSum(FabArray< FAB > const &fa, int nghost, F &&f)
Definition AMReX_FabArrayUtility.H:16