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