WarpX
FilterCreateTransformFromFAB.H
Go to the documentation of this file.
1 /* Copyright 2020 Luca Fedeli, Neil Zaim
2  *
3  * This file is part of WarpX.
4  *
5  * License: BSD-3-Clause-LBNL
6  */
7 
8 #ifndef WARPX_FILTER_CREATE_TRANSFORM_FROM_FAB_H_
9 #define WARPX_FILTER_CREATE_TRANSFORM_FROM_FAB_H_
10 
12 
13 #include <AMReX_REAL.H>
14 #include <AMReX_TypeTraits.H>
15 
46 template <int N, typename DstPC, typename DstTile, typename FAB, typename Index,
47  typename CreateFunc1, typename CreateFunc2, typename TransFunc,
49 Index filterCreateTransformFromFAB (DstPC& pc1, DstPC& pc2,
50  DstTile& dst1, DstTile& dst2, const amrex::Box box,
51  const FAB *src_FAB, const Index* mask,
52  const Index dst1_index, const Index dst2_index,
53  CreateFunc1&& create1, CreateFunc2&& create2,
54  TransFunc&& transform, const amrex::Geometry& geom_lev_zero) noexcept
55 {
56  using namespace amrex;
57 
58  const auto ncells = box.volume();
59  if (ncells == 0) { return 0; }
60 
61  constexpr int spacedim = AMREX_SPACEDIM;
62 
63 #if defined(WARPX_DIM_1D_Z)
64  const Real zlo_global = geom_lev_zero.ProbLo(0);
65  const Real dz = geom_lev_zero.CellSize(0);
66 #elif defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ)
67  const Real xlo_global = geom_lev_zero.ProbLo(0);
68  const Real dx = geom_lev_zero.CellSize(0);
69  const Real zlo_global = geom_lev_zero.ProbLo(1);
70  const Real dz = geom_lev_zero.CellSize(1);
71 #elif defined(WARPX_DIM_3D)
72  const Real xlo_global = geom_lev_zero.ProbLo(0);
73  const Real dx = geom_lev_zero.CellSize(0);
74  const Real ylo_global = geom_lev_zero.ProbLo(1);
75  const Real dy = geom_lev_zero.CellSize(1);
76  const Real zlo_global = geom_lev_zero.ProbLo(2);
77  const Real dz = geom_lev_zero.CellSize(2);
78 #endif
79 
80  const auto arrNumPartCreation = src_FAB->array();
81  Gpu::DeviceVector<Index> offsets(ncells);
82  auto total = amrex::Scan::ExclusiveSum(ncells, mask, offsets.data());
83  const Index num_added = N*total;
84  auto old_np1 = dst1.size();
85  auto new_np1 = std::max(dst1_index + num_added, dst1.numParticles());
86  dst1.resize(new_np1);
87 
88  auto old_np2 = dst2.size();
89  auto new_np2 = std::max(dst2_index + num_added, dst2.numParticles());
90  dst2.resize(new_np2);
91 
92  auto *p_offsets = offsets.dataPtr();
93 
94  const auto dst1_data = dst1.getParticleTileData();
95  const auto dst2_data = dst2.getParticleTileData();
96 
97  // For loop over all cells in the box. If mask is true in the given cell,
98  // we create the particles in the cell and apply a transform function to the
99  // created particles.
100  amrex::ParallelForRNG(ncells,
101  [=] AMREX_GPU_DEVICE (int i, amrex::RandomEngine const& engine) noexcept
102  {
103  if (mask[i])
104  {
105  const IntVect iv = box.atOffset(i);
106  const int j = iv[0];
107  const int k = (spacedim >= 2) ? iv[1] : 0;
108  const int l = (spacedim == 3) ? iv[2] : 0;
109 
110  // Currently all particles are created on nodes. This makes it useless
111  // to use N>1 (for now).
112 #if defined(WARPX_DIM_1D_Z)
113  Real const x = 0.0;
114  Real const y = 0.0;
115  Real const z = zlo_global + j*dz;
116 #elif defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ)
117  Real const x = xlo_global + j*dx;
118  Real const y = 0.0;
119  Real const z = zlo_global + k*dz;
120 #elif defined(WARPX_DIM_3D)
121  Real const x = xlo_global + j*dx;
122  Real const y = ylo_global + k*dy;
123  Real const z = zlo_global + l*dz;
124 #endif
125 
126  for (int n = 0; n < N; ++n)
127  {
128  create1(dst1_data, N*p_offsets[i] + dst1_index + n, engine, x, y, z);
129  create2(dst2_data, N*p_offsets[i] + dst2_index + n, engine, x, y, z);
130  }
131  transform(dst1_data, dst2_data, N*p_offsets[i] + dst1_index,
132  N*p_offsets[i] + dst2_index, N, arrNumPartCreation(j,k,l));
133  }
134  });
135 
137  0, 0,
138  pc1.getUserRealAttribs(), pc1.getUserIntAttribs(),
139  pc1.getParticleComps(), pc1.getParticleiComps(),
140  pc1.getUserRealAttribParser(),
141  pc1.getUserIntAttribParser(),
142 #ifdef WARPX_QED
143  false, // do not initialize QED quantities, since they were initialized
144  // when calling the CreateFunc functor
145  pc1.get_breit_wheeler_engine_ptr(),
146  pc1.get_quantum_sync_engine_ptr(),
147 #endif
148  pc1.getIonizationInitialLevel(),
149  old_np1, new_np1);
151  0, 0,
152  pc2.getUserRealAttribs(), pc2.getUserIntAttribs(),
153  pc2.getParticleComps(), pc2.getParticleiComps(),
154  pc2.getUserRealAttribParser(),
155  pc2.getUserIntAttribParser(),
156 #ifdef WARPX_QED
157  false, // do not initialize QED quantities, since they were initialized
158  // when calling the CreateFunc functor
159  pc2.get_breit_wheeler_engine_ptr(),
160  pc2.get_quantum_sync_engine_ptr(),
161 #endif
162  pc2.getIonizationInitialLevel(),
163  old_np2, new_np2);
164 
166  return num_added;
167 }
168 
202 template <int N, typename DstPC, typename DstTile, typename FABs, typename Index,
203  typename FilterFunc, typename CreateFunc1, typename CreateFunc2,
204  typename TransFunc>
205 Index filterCreateTransformFromFAB (DstPC& pc1, DstPC& pc2, DstTile& dst1, DstTile& dst2, const amrex::Box box,
206  const FABs& src_FABs, const Index dst1_index,
207  const Index dst2_index, FilterFunc&& filter,
208  CreateFunc1&& create1, CreateFunc2&& create2,
209  TransFunc && transform, const amrex::Geometry& geom_lev_zero) noexcept
210 {
211  using namespace amrex;
212 
213  FArrayBox NumPartCreation(box, 1);
214  // This may be unnecessary because of the Gpu::streamSynchronize() in
215  // filterCreateTransformFromFAB called below, but let us keep it for safety.
216  const Elixir tmp_eli = NumPartCreation.elixir();
217  auto arrNumPartCreation = NumPartCreation.array();
218 
219  const auto ncells = box.volume();
220  if (ncells == 0) { return 0; }
221 
223 
224  auto *p_mask = mask.dataPtr();
225 
226  // for loop over all cells in the box. We apply the filter function to each cell
227  // and store the result in arrNumPartCreation. If the result is strictly greater than
228  // 0, the mask is set to true at the given cell position.
230  [=] AMREX_GPU_DEVICE (int i, int j, int k, amrex::RandomEngine const& engine){
231  arrNumPartCreation(i,j,k) = filter(src_FABs,i,j,k,engine);
232  const IntVect iv(AMREX_D_DECL(i,j,k));
233  const auto mask_position = box.index(iv);
234  p_mask[mask_position] = (arrNumPartCreation(i,j,k) > 0);
235  });
236 
237  return filterCreateTransformFromFAB<N>(pc1, pc2, dst1, dst2, box, &NumPartCreation,
238  mask.dataPtr(), dst1_index, dst2_index,
239  std::forward<CreateFunc1>(create1),
240  std::forward<CreateFunc2>(create2),
241  std::forward<TransFunc>(transform),
242  geom_lev_zero);
243 }
244 
245 #endif // WARPX_FILTER_CREATE_TRANSFORM_FROM_FAB_H_
#define AMREX_GPU_DEVICE
Array4< int const > mask
#define AMREX_D_DECL(a, b, c)
Index filterCreateTransformFromFAB(DstPC &pc1, DstPC &pc2, DstTile &dst1, DstTile &dst2, const amrex::Box box, const FAB *src_FAB, const Index *mask, const Index dst1_index, const Index dst2_index, CreateFunc1 &&create1, CreateFunc2 &&create2, TransFunc &&transform, const amrex::Geometry &geom_lev_zero) noexcept
Apply a filter on a list of FABs, then create and apply a transform operation to the particles depend...
Definition: FilterCreateTransformFromFAB.H:49
AMREX_FORCE_INLINE Array4< Real const > array() const noexcept
Elixir elixir() noexcept
T * data() noexcept
T * dataPtr() noexcept
void DefaultInitializeRuntimeAttributes(PTile &ptile, const int n_external_attr_real, const int n_external_attr_int, const std::vector< std::string > &user_real_attribs, const std::vector< std::string > &user_int_attribs, const std::map< std::string, int > &particle_comps, const std::map< std::string, int > &particle_icomps, const std::vector< amrex::Parser * > &user_real_attrib_parser, const std::vector< amrex::Parser * > &user_int_attrib_parser, const bool do_qed_comps, BreitWheelerEngine *p_bw_engine, QuantumSynchrotronEngine *p_qs_engine, const int ionization_initial_level, int start, int stop)
Default initialize runtime attributes in a tile. This routine does not initialize the first n_externa...
Definition: DefaultInitialization.H:118
FAB
void synchronize() noexcept
T ExclusiveSum(N n, T const *in, T *out, RetSum a_ret_sum=retSum)
AMREX_ATTRIBUTE_FLATTEN_FOR void ParallelForRNG(T n, L const &f) noexcept
std::enable_if_t< B, T > EnableIf_t
i
Definition: check_interp_points_and_weights.py:174
int dz
Definition: compute_domain.py:36
tuple dx
lab frame
Definition: stencil.py:429
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE int const * dataPtr() const noexcept