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}
39
80template <class PC, class F, std::enable_if_t<IsParticleContainer<PC>::value, int> foo = 0>
81auto
82ReduceSum (PC const& pc, F&& f)
83 -> decltype(particle_detail::call_f(f, typename PC::ConstPTDType(), int()))
84{
85 return ReduceSum(pc, 0, pc.finestLevel(), std::forward<F>(f));
86}
87
128template <class PC, class F, std::enable_if_t<IsParticleContainer<PC>::value, int> foo = 0>
129auto
130ReduceSum (PC const& pc, int lev, F&& f)
131 -> decltype(particle_detail::call_f(f, typename PC::ConstPTDType(), int()))
132{
133 return ReduceSum(pc, lev, lev, std::forward<F>(f));
134}
135
177template <class PC, class F, std::enable_if_t<IsParticleContainer<PC>::value, int> foo = 0>
178auto
179ReduceSum (PC const& pc, int lev_min, int lev_max, F const& f)
180 -> decltype(particle_detail::call_f(f, typename PC::ConstPTDType(), int()))
181{
182 using value_type = decltype(particle_detail::call_f(f, typename PC::ConstPTDType(), int()));
183 value_type sm = 0;
184
185#ifdef AMREX_USE_GPU
186 if (Gpu::inLaunchRegion())
187 {
188 ReduceOps<ReduceOpSum> reduce_op;
189 ReduceData<value_type> reduce_data(reduce_op);
190 using ReduceTuple = typename decltype(reduce_data)::Type;
191
192 for (int lev = lev_min; lev <= lev_max; ++lev)
193 {
194 const auto& plev = pc.GetParticles(lev);
195 for (const auto& kv : plev)
196 {
197 const auto& tile = plev.at(kv.first);
198 const auto np = tile.numParticles();
199 const auto& ptd = tile.getConstParticleTileData();
200 reduce_op.eval(np, reduce_data,
201 [=] AMREX_GPU_DEVICE (const int i) -> ReduceTuple {
202 return particle_detail::call_f(f, ptd, i);
203 });
204 }
205 }
206
207 ReduceTuple hv = reduce_data.value(reduce_op);
208 sm = amrex::get<0>(hv);
209 }
210 else
211#endif
212 {
213 for (int lev = lev_min; lev <= lev_max; ++lev)
214 {
215 const auto& plev = pc.GetParticles(lev);
216 Vector<std::pair<int, int> > grid_tile_ids;
217 Vector<const typename PC::ParticleTileType*> ptile_ptrs;
218 for (auto& kv : plev)
219 {
220 grid_tile_ids.push_back(kv.first);
221 ptile_ptrs.push_back(&(kv.second));
222 }
223
224 // We avoid ParIter here to allow reductions after grids change but prior to calling Redistribute()
225#ifdef AMREX_USE_OMP
226#pragma omp parallel for if (!system::regtest_reduction) reduction(+:sm)
227#endif
228 for (int pmap_it = 0; pmap_it < static_cast<int>(ptile_ptrs.size()); ++pmap_it)
229 {
230 const auto& tile = plev.at(grid_tile_ids[pmap_it]);
231 const auto np = tile.numParticles();
232 const auto& ptd = tile.getConstParticleTileData();
233 for (int i = 0; i < np; ++i) {
234 sm += particle_detail::call_f(f, ptd, i);
235 }
236 }
237 }
238 }
239
240 return sm;
241}
242
282template <class PC, class F, std::enable_if_t<IsParticleContainer<PC>::value, int> foo = 0>
283auto
284ReduceMax (PC const& pc, F&& f)
285 -> decltype(particle_detail::call_f(f, typename PC::ConstPTDType(), int()))
286{
287 return ReduceMax(pc, 0, pc.finestLevel(), std::forward<F>(f));
288}
289
331template <class PC, class F, std::enable_if_t<IsParticleContainer<PC>::value, int> foo = 0>
332auto
333ReduceMax (PC const& pc, int lev, F&& f)
334 -> decltype(particle_detail::call_f(f, typename PC::ConstPTDType(), int()))
335{
336 return ReduceMax(pc, lev, lev, std::forward<F>(f));
337}
338
380template <class PC, class F, std::enable_if_t<IsParticleContainer<PC>::value, int> foo = 0>
381auto
382ReduceMax (PC const& pc, int lev_min, int lev_max, F const& f)
383 -> decltype(particle_detail::call_f(f, typename PC::ConstPTDType(), int()))
384{
385 using value_type = decltype(particle_detail::call_f(f, typename PC::ConstPTDType(), int()));
386 constexpr value_type value_lowest = std::numeric_limits<value_type>::lowest();
387 value_type r = value_lowest;
388
389#ifdef AMREX_USE_GPU
390 if (Gpu::inLaunchRegion())
391 {
392 ReduceOps<ReduceOpMax> reduce_op;
393 ReduceData<value_type> reduce_data(reduce_op);
394 using ReduceTuple = typename decltype(reduce_data)::Type;
395
396 for (int lev = lev_min; lev <= lev_max; ++lev)
397 {
398 const auto& plev = pc.GetParticles(lev);
399 for (const auto& kv : plev)
400 {
401 const auto& tile = plev.at(kv.first);
402 const auto np = tile.numParticles();
403 const auto& ptd = tile.getConstParticleTileData();
404 reduce_op.eval(np, reduce_data,
405 [=] AMREX_GPU_DEVICE (const int i) -> ReduceTuple {
406 return particle_detail::call_f(f, ptd, i);
407 });
408 }
409 }
410
411 ReduceTuple hv = reduce_data.value(reduce_op);
412 r = amrex::get<0>(hv);
413 }
414 else
415#endif
416 {
417 for (int lev = lev_min; lev <= lev_max; ++lev)
418 {
419 const auto& plev = pc.GetParticles(lev);
420 Vector<std::pair<int, int> > grid_tile_ids;
421 Vector<const typename PC::ParticleTileType*> ptile_ptrs;
422 for (auto& kv : plev)
423 {
424 grid_tile_ids.push_back(kv.first);
425 ptile_ptrs.push_back(&(kv.second));
426 }
427
428 // We avoid ParIter here to allow reductions after grids change but prior to calling Redistribute()
429#ifdef AMREX_USE_OMP
430#pragma omp parallel for if (!system::regtest_reduction) reduction(max:r)
431#endif
432 for (int pmap_it = 0; pmap_it < static_cast<int>(ptile_ptrs.size()); ++pmap_it)
433 {
434 const auto& tile = plev.at(grid_tile_ids[pmap_it]);
435 const auto np = tile.numParticles();
436 const auto& ptd = tile.getConstParticleTileData();
437 for (int i = 0; i < np; ++i) {
438 r = std::max(r, particle_detail::call_f(f, ptd, i));
439 }
440 }
441 }
442 }
443
444 return r;
445}
446
486template <class PC, class F, std::enable_if_t<IsParticleContainer<PC>::value, int> foo = 0>
487auto
488ReduceMin (PC const& pc, F&& f)
489 -> decltype(particle_detail::call_f(f, typename PC::ConstPTDType(), int()))
490{
491 return ReduceMin(pc, 0, pc.finestLevel(), std::forward<F>(f));
492}
493
534template <class PC, class F, std::enable_if_t<IsParticleContainer<PC>::value, int> foo = 0>
535auto
536ReduceMin (PC const& pc, int lev, F&& f)
537 -> decltype(particle_detail::call_f(f, typename PC::ConstPTDType(), int()))
538{
539 return ReduceMin(pc, lev, lev, std::forward<F>(f));
540}
541
583template <class PC, class F, std::enable_if_t<IsParticleContainer<PC>::value, int> foo = 0>
584auto
585ReduceMin (PC const& pc, int lev_min, int lev_max, F const& f)
586 -> decltype(particle_detail::call_f(f, typename PC::ConstPTDType(), int()))
587{
588 using value_type = decltype(particle_detail::call_f(f, typename PC::ConstPTDType(), int()));
589 constexpr value_type value_max = std::numeric_limits<value_type>::max();
590 value_type r = value_max;
591
592#ifdef AMREX_USE_GPU
593 if (Gpu::inLaunchRegion())
594 {
595 ReduceOps<ReduceOpMin> reduce_op;
596 ReduceData<value_type> reduce_data(reduce_op);
597 using ReduceTuple = typename decltype(reduce_data)::Type;
598
599 for (int lev = lev_min; lev <= lev_max; ++lev)
600 {
601 const auto& plev = pc.GetParticles(lev);
602 for (const auto& kv : plev)
603 {
604 const auto& tile = plev.at(kv.first);
605 const auto np = tile.numParticles();
606 const auto& ptd = tile.getConstParticleTileData();
607 reduce_op.eval(np, reduce_data,
608 [=] AMREX_GPU_DEVICE (const int i) -> ReduceTuple {
609 return particle_detail::call_f(f, ptd, i);
610 });
611 }
612 }
613
614 ReduceTuple hv = reduce_data.value(reduce_op);
615 r = amrex::get<0>(hv);
616 }
617 else
618#endif
619 {
620 for (int lev = lev_min; lev <= lev_max; ++lev)
621 {
622 const auto& plev = pc.GetParticles(lev);
623 Vector<std::pair<int, int> > grid_tile_ids;
624 Vector<const typename PC::ParticleTileType*> ptile_ptrs;
625 for (auto& kv : plev)
626 {
627 grid_tile_ids.push_back(kv.first);
628 ptile_ptrs.push_back(&(kv.second));
629 }
630
631 // We avoid ParIter here to allow reductions after grids change but prior to calling Redistribute()
632#ifdef AMREX_USE_OMP
633#pragma omp parallel for if (!system::regtest_reduction) reduction(min:r)
634#endif
635 for (int pmap_it = 0; pmap_it < static_cast<int>(ptile_ptrs.size()); ++pmap_it)
636 {
637 const auto& tile = plev.at(grid_tile_ids[pmap_it]);
638 const auto np = tile.numParticles();
639 const auto& ptd = tile.getConstParticleTileData();
640 for (int i = 0; i < np; ++i) {
641 r = std::min(r, particle_detail::call_f(f, ptd, i));
642 }
643 }
644 }
645 }
646
647 return r;
648}
649
689template <class PC, class F, std::enable_if_t<IsParticleContainer<PC>::value, int> foo = 0>
690bool
691ReduceLogicalAnd (PC const& pc, F&& f)
692{
693 return ReduceLogicalAnd(pc, 0, pc.finestLevel(), std::forward<F>(f));
694}
695
736template <class PC, class F, std::enable_if_t<IsParticleContainer<PC>::value, int> foo = 0>
737bool
738ReduceLogicalAnd (PC const& pc, int lev, F&& f)
739{
740 return ReduceLogicalAnd(pc, lev, lev, std::forward<F>(f));
741}
742
784template <class PC, class F, std::enable_if_t<IsParticleContainer<PC>::value, int> foo = 0>
785bool
786ReduceLogicalAnd (PC const& pc, int lev_min, int lev_max, F const& f)
787{
788 int r = true;
789
790#ifdef AMREX_USE_GPU
791 if (Gpu::inLaunchRegion())
792 {
793 ReduceOps<ReduceOpLogicalAnd> reduce_op;
794 ReduceData<int> reduce_data(reduce_op);
795 using ReduceTuple = typename decltype(reduce_data)::Type;
796
797 for (int lev = lev_min; lev <= lev_max; ++lev)
798 {
799 const auto& plev = pc.GetParticles(lev);
800 for (const auto& kv : plev)
801 {
802 const auto& tile = plev.at(kv.first);
803 const auto np = tile.numParticles();
804 const auto& ptd = tile.getConstParticleTileData();
805 reduce_op.eval(np, reduce_data,
806 [=] AMREX_GPU_DEVICE (const int i) -> ReduceTuple {
807 return particle_detail::call_f(f, ptd, i);
808 });
809 }
810 }
811
812 ReduceTuple hv = reduce_data.value(reduce_op);
813 r = amrex::get<0>(hv);
814 }
815 else
816#endif
817 {
818 for (int lev = lev_min; lev <= lev_max; ++lev)
819 {
820 const auto& plev = pc.GetParticles(lev);
821 Vector<std::pair<int, int> > grid_tile_ids;
822 Vector<const typename PC::ParticleTileType*> ptile_ptrs;
823 for (auto& kv : plev)
824 {
825 grid_tile_ids.push_back(kv.first);
826 ptile_ptrs.push_back(&(kv.second));
827 }
828
829 // We avoid ParIter here to allow reductions after grids change but prior to calling Redistribute()
830#ifdef AMREX_USE_OMP
831#pragma omp parallel for if (!system::regtest_reduction) reduction(&&:r)
832#endif
833 for (int pmap_it = 0; pmap_it < static_cast<int>(ptile_ptrs.size()); ++pmap_it)
834 {
835 const auto& tile = plev.at(grid_tile_ids[pmap_it]);
836 const auto np = tile.numParticles();
837 const auto& ptd = tile.getConstParticleTileData();
838 for (int i = 0; i < np; ++i) {
839 r = r && particle_detail::call_f(f, ptd, i);
840 }
841 }
842 }
843 }
844
845 return r;
846}
847
887template <class PC, class F, std::enable_if_t<IsParticleContainer<PC>::value, int> foo = 0>
888bool
889ReduceLogicalOr (PC const& pc, F&& f)
890{
891 return ReduceLogicalOr(pc, 0, pc.finestLevel(), std::forward<F>(f));
892}
893
934template <class PC, class F, std::enable_if_t<IsParticleContainer<PC>::value, int> foo = 0>
935bool
936ReduceLogicalOr (PC const& pc, int lev, F&& f)
937{
938 return ReduceLogicalOr(pc, lev, lev, std::forward<F>(f));
939}
940
982template <class PC, class F, std::enable_if_t<IsParticleContainer<PC>::value, int> foo = 0>
983bool
984ReduceLogicalOr (PC const& pc, int lev_min, int lev_max, F const& f)
985{
986 int r = false;
987
988#ifdef AMREX_USE_GPU
989 if (Gpu::inLaunchRegion())
990 {
991 ReduceOps<ReduceOpLogicalOr> reduce_op;
992 ReduceData<int> reduce_data(reduce_op);
993 using ReduceTuple = typename decltype(reduce_data)::Type;
994
995 for (int lev = lev_min; lev <= lev_max; ++lev)
996 {
997 const auto& plev = pc.GetParticles(lev);
998 for (const auto& kv : plev)
999 {
1000 const auto& tile = plev.at(kv.first);
1001 const auto np = tile.numParticles();
1002 const auto& ptd = tile.getConstParticleTileData();
1003 reduce_op.eval(np, reduce_data,
1004 [=] AMREX_GPU_DEVICE (const int i) -> ReduceTuple
1005 {
1006 return particle_detail::call_f(f, ptd, i);
1007 });
1008 }
1009 }
1010
1011 ReduceTuple hv = reduce_data.value(reduce_op);
1012 r = amrex::get<0>(hv);
1013 }
1014 else
1015#endif
1016 {
1017 for (int lev = lev_min; lev <= lev_max; ++lev)
1018 {
1019 const auto& plev = pc.GetParticles(lev);
1020 Vector<std::pair<int, int> > grid_tile_ids;
1021 Vector<const typename PC::ParticleTileType*> ptile_ptrs;
1022 for (auto& kv : plev)
1023 {
1024 grid_tile_ids.push_back(kv.first);
1025 ptile_ptrs.push_back(&(kv.second));
1026 }
1027
1028 // We avoid ParIter here to allow reductions after grids change but prior to calling Redistribute()
1029#ifdef AMREX_USE_OMP
1030#pragma omp parallel for if (!system::regtest_reduction) reduction(||:r)
1031#endif
1032 for (int pmap_it = 0; pmap_it < static_cast<int>(ptile_ptrs.size()); ++pmap_it)
1033 {
1034 const auto& tile = plev.at(grid_tile_ids[pmap_it]);
1035 const auto np = tile.numParticles();
1036 const auto& ptd = tile.getConstParticleTileData();
1037 for (int i = 0; i < np; ++i) {
1038 r = r || particle_detail::call_f(f, ptd, i);
1039 }
1040 }
1041 }
1042 }
1043
1044 return r;
1045}
1046
1109template <class RD, class PC, class F, class ReduceOps,
1110 std::enable_if_t<IsParticleContainer<PC>::value, int> foo = 0>
1111typename RD::Type
1112ParticleReduce (PC const& pc, F&& f, ReduceOps& reduce_ops)
1113{
1114 return ParticleReduce<RD>(pc, 0, pc.finestLevel(), std::forward<F>(f), reduce_ops);
1115}
1116
1180template <class RD, class PC, class F, class ReduceOps,
1181 std::enable_if_t<IsParticleContainer<PC>::value, int> foo = 0>
1182typename RD::Type
1183ParticleReduce (PC const& pc, int lev, F&& f, ReduceOps& reduce_ops)
1184{
1185 return ParticleReduce<RD>(pc, lev, lev, std::forward<F>(f), reduce_ops);
1186}
1187
1252template <class RD, class PC, class F, class ReduceOps,
1253 std::enable_if_t<IsParticleContainer<PC>::value, int> foo = 0>
1254typename RD::Type
1255ParticleReduce (PC const& pc, int lev_min, int lev_max, F const& f, ReduceOps& reduce_ops)
1256{
1257 RD reduce_data(reduce_ops);
1258 for (int lev = lev_min; lev <= lev_max; ++lev) {
1259 const auto& plev = pc.GetParticles(lev);
1260 Vector<std::pair<int, int> > grid_tile_ids;
1261 Vector<const typename PC::ParticleTileType*> ptile_ptrs;
1262 for (auto& kv : plev)
1263 {
1264 grid_tile_ids.push_back(kv.first);
1265 ptile_ptrs.push_back(&(kv.second));
1266 }
1267
1268 // We avoid ParIter here to allow reductions after grids change but prior to calling Redistribute()
1269#if !defined(AMREX_USE_GPU) && defined(AMREX_USE_OMP)
1270#pragma omp parallel for
1271#endif
1272 for (int pmap_it = 0; pmap_it < static_cast<int>(ptile_ptrs.size()); ++pmap_it)
1273 {
1274 const auto& tile = plev.at(grid_tile_ids[pmap_it]);
1275 const auto np = tile.numParticles();
1276 const auto& ptd = tile.getConstParticleTileData();
1277 reduce_ops.eval(np, reduce_data,
1278 [=] AMREX_GPU_DEVICE (const int i) noexcept
1279 {
1280 return particle_detail::call_f(f, ptd, i);
1281 });
1282 }
1283 }
1284 return reduce_data.value(reduce_ops);
1285}
1286}
1287#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:257
amrex::ReduceData::value
Type value()
Definition AMReX_Reduce.H:289
amrex::ReduceOps
Definition AMReX_Reduce.H:377
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:446
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:1112
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