WarpX
DepositCharge.H
Go to the documentation of this file.
1 /* Copyright 2019-2021 Axel Huebl, Andrew Myers
2  *
3  * This file is part of WarpX.
4  *
5  * License: BSD-3-Clause-LBNL
6  */
7 #ifndef ABLASTR_DEPOSIT_CHARGE_H_
8 #define ABLASTR_DEPOSIT_CHARGE_H_
9 
13 #include "ablastr/utils/TextMsg.H"
14 
15 #include <AMReX.H>
16 
17 #include <optional>
18 
19 
21 {
22 
49 template< typename T_PC >
50 static void
51 deposit_charge (typename T_PC::ParIterType& pti,
52  typename T_PC::RealVector const& wp,
53  amrex::Real const charge,
54  int const * const ion_lev,
55  amrex::MultiFab* rho,
56  amrex::FArrayBox& local_rho,
57  int const particle_shape,
58  std::array<amrex::Real, 3> const & dx,
59  std::array<amrex::Real, 3> const & xyzmin,
60  int const n_rz_azimuthal_modes = 0,
61  std::optional<amrex::IntVect> num_rho_deposition_guards = std::nullopt,
62  std::optional<int> depos_lev = std::nullopt,
63  std::optional<amrex::IntVect> rel_ref_ratio = std::nullopt,
64  long const offset = 0,
65  std::optional<long> np_to_deposit = std::nullopt,
66  int const icomp = 0, int const nc = 1,
67  bool const do_device_synchronize = true)
68 {
69  // deposition guards
70  amrex::IntVect ng_rho = rho->nGrowVect();
71  if (num_rho_deposition_guards.has_value()) {
72  ng_rho = num_rho_deposition_guards.value();
73  }
75  "num_rho_deposition_guards are larger than allocated!");
76  // particle shape
77  auto const[nox, noy, noz] = std::array<int, 3>{particle_shape, particle_shape, particle_shape};
78 
79  // used for MR when we want to deposit for a subset of the particles on the level in the
80  // current box; with offset, we start at a later particle index
81  if (!np_to_deposit.has_value()) {
82  np_to_deposit = pti.numParticles();
83  }
84  ABLASTR_ALWAYS_ASSERT_WITH_MESSAGE(np_to_deposit.value() + offset <= pti.numParticles(),
85  "np_to_deposit + offset are out-of-bounds for particle iterator");
86 
87  int const lev = pti.GetLevel();
88  if (!depos_lev.has_value()) {
89  depos_lev = lev;
90  }
91  ABLASTR_ALWAYS_ASSERT_WITH_MESSAGE((depos_lev.value() == (lev-1)) ||
92  (depos_lev.value() == (lev )),
93  "Deposition buffers only work for lev or lev-1");
94  if (!rel_ref_ratio.has_value()) {
95  ABLASTR_ALWAYS_ASSERT_WITH_MESSAGE(lev == depos_lev,
96  "rel_ref_ratio must be set if lev != depos_lev");
97  rel_ref_ratio = amrex::IntVect(AMREX_D_DECL(1, 1, 1));
98  }
99 
100  // If no particles, do not do anything
101  if (np_to_deposit == 0) { return; }
102 
103  // Extract deposition order and check that particles shape fits within the guard cells.
104  // NOTE: In specific situations where the staggering of rho and the charge deposition algorithm
105  // are not trivial, this check might be too strict and we might need to relax it, as currently
106  // done for the current deposition.
107 
108 #if defined(WARPX_DIM_1D_Z)
111  const amrex::IntVect shape_extent = amrex::IntVect(static_cast<int>(noz/2+1));
112 #elif defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ)
114  const amrex::IntVect shape_extent = amrex::IntVect(static_cast<int>(nox/2+1),
115  static_cast<int>(noz/2+1));
116 #elif defined(WARPX_DIM_3D)
117  const amrex::IntVect shape_extent = amrex::IntVect(static_cast<int>(nox/2+1),
118  static_cast<int>(noy/2+1),
119  static_cast<int>(noz/2+1));
120 #endif
121 
122  // On CPU: particles deposit on tile arrays, which have a small number of guard cells ng_rho
123  // On GPU: particles deposit directly on the rho array, which usually have a larger number of guard cells
124 #ifndef AMREX_USE_GPU
125  const amrex::IntVect range = ng_rho - shape_extent;
126 #else
127  const amrex::IntVect range = rho->nGrowVect() - shape_extent;
128 #endif
129  amrex::ignore_unused(range); // for release builds
131  amrex::numParticlesOutOfRange(pti, range) == 0,
132  "Particles shape does not fit within tile (CPU) or guard cells (GPU) used for charge deposition");
133 
134  ABLASTR_PROFILE_VAR_NS("ablastr::particles::deposit_charge::ChargeDeposition", blp_ppc_chd, do_device_synchronize);
135  ABLASTR_PROFILE_VAR_NS("ablastr::particles::deposit_charge::Accumulate", blp_accumulate, do_device_synchronize);
136 
137  // Get tile box where charge is deposited.
138  // The tile box is different when depositing in the buffers (depos_lev<lev)
139  // or when depositing inside the level (depos_lev=lev)
140  amrex::Box tilebox;
141  if (lev == depos_lev) {
142  tilebox = pti.tilebox();
143  } else {
144  tilebox = amrex::coarsen(pti.tilebox(), rel_ref_ratio.value());
145  }
146 
147 #ifndef AMREX_USE_GPU
148  // Staggered tile box
149  amrex::Box tb = amrex::convert( tilebox, rho->ixType().toIntVect() );
150 #endif
151 
152  tilebox.grow(ng_rho);
153 
154 #ifdef AMREX_USE_GPU
155  amrex::ignore_unused(local_rho);
156  // GPU, no tiling: rho_fab points to the full rho array
157  amrex::MultiFab rhoi(*rho, amrex::make_alias, icomp*nc, nc);
158  auto & rho_fab = rhoi.get(pti);
159 #else
160  tb.grow(ng_rho);
161 
162  // CPU, tiling: rho_fab points to local_rho
163  local_rho.resize(tb, nc);
164 
165  // local_rho is set to zero
166  local_rho.setVal(0.0);
167 
168  auto & rho_fab = local_rho;
169 #endif
170 
171  const auto GetPosition = GetParticlePosition<PIdx>(pti, offset);
172 
173  // Indices of the lower bound
174  const amrex::Dim3 lo = lbound(tilebox);
175 
176  ABLASTR_PROFILE_VAR_START(blp_ppc_chd, do_device_synchronize);
177 
178  if (nox == 1){
179  doChargeDepositionShapeN<1>(GetPosition, wp.dataPtr()+offset, ion_lev,
180  rho_fab, np_to_deposit.value(), dx, xyzmin, lo, charge,
181  n_rz_azimuthal_modes);
182  } else if (nox == 2){
183  doChargeDepositionShapeN<2>(GetPosition, wp.dataPtr()+offset, ion_lev,
184  rho_fab, np_to_deposit.value(), dx, xyzmin, lo, charge,
185  n_rz_azimuthal_modes);
186  } else if (nox == 3){
187  doChargeDepositionShapeN<3>(GetPosition, wp.dataPtr()+offset, ion_lev,
188  rho_fab, np_to_deposit.value(), dx, xyzmin, lo, charge,
189  n_rz_azimuthal_modes);
190  } else if (nox == 4){
191  doChargeDepositionShapeN<4>(GetPosition, wp.dataPtr()+offset, ion_lev,
192  rho_fab, np_to_deposit.value(), dx, xyzmin, lo, charge,
193  n_rz_azimuthal_modes);
194  }
195  ABLASTR_PROFILE_VAR_STOP(blp_ppc_chd, do_device_synchronize);
196 
197 #ifndef AMREX_USE_GPU
198  // CPU, tiling: atomicAdd local_rho into rho
199  ABLASTR_PROFILE_VAR_START(blp_accumulate, do_device_synchronize);
200  (*rho)[pti].lockAdd(local_rho, tb, tb, 0, icomp*nc, nc);
201  ABLASTR_PROFILE_VAR_STOP(blp_accumulate, do_device_synchronize);
202 #endif
203 }
204 
205 } // namespace ablastr::particles
206 
207 #endif // ABLASTR_DEPOSIT_CHARGE_H_
#define AMREX_D_DECL(a, b, c)
#define ABLASTR_PROFILE_VAR_START(vname, sync)
Definition: ProfilerWrapper.H:58
#define ABLASTR_PROFILE_VAR_NS(fname, vname, sync)
Definition: ProfilerWrapper.H:57
#define ABLASTR_PROFILE_VAR_STOP(vname, sync)
Definition: ProfilerWrapper.H:59
#define ABLASTR_ALWAYS_ASSERT_WITH_MESSAGE(EX, MSG)
Definition: TextMsg.H:75
void setVal(Real const &x, const Box &bx, int dcomp, int ncomp) noexcept
AMREX_GPU_HOST_DEVICE Box & grow(int i) noexcept
void resize(const Box &b, int N=1, Arena *ar=nullptr)
IntVect nGrowVect() const noexcept
IndexType ixType() const noexcept
const FArrayBox & get(const MFIter &mfi) const noexcept
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE IntVect toIntVect() const noexcept
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE bool allLE(const IntVect &rhs) const noexcept
Definition: DepositCharge.H:21
static void deposit_charge(typename T_PC::ParIterType &pti, typename T_PC::RealVector const &wp, amrex::Real const charge, int const *const ion_lev, amrex::MultiFab *rho, amrex::FArrayBox &local_rho, int const particle_shape, std::array< amrex::Real, 3 > const &dx, std::array< amrex::Real, 3 > const &xyzmin, int const n_rz_azimuthal_modes=0, std::optional< amrex::IntVect > num_rho_deposition_guards=std::nullopt, std::optional< int > depos_lev=std::nullopt, std::optional< amrex::IntVect > rel_ref_ratio=std::nullopt, long const offset=0, std::optional< long > np_to_deposit=std::nullopt, int const icomp=0, int const nc=1, bool const do_device_synchronize=true)
Definition: DepositCharge.H:51
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE Box convert(const Box &b, const IntVect &typ) noexcept
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE Dim3 lbound(Array4< T > const &a) noexcept
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void ignore_unused(const Ts &...)
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE Box coarsen(const Box &b, int ref_ratio) noexcept
int numParticlesOutOfRange(Iterator const &pti, int nGrow)
int nc
Definition: check_interp_points_and_weights.py:130
tuple dx
lab frame
Definition: stencil.py:429