WarpX
Loading...
Searching...
No Matches
MassMatricesDeposition.H
Go to the documentation of this file.
1/* Copyright 2019 Axel Huebl, David Grote, Maxence Thevenet
2 * Remi Lehe, Weiqun Zhang, Michael Rowan
3 *
4 * This file is part of WarpX.
5 *
6 * License: BSD-3-Clause-LBNL
7 */
8#ifndef WARPX_MASS_MATRICES_DEPOSITION_H_
9#define WARPX_MASS_MATRICES_DEPOSITION_H_
10
11#include "Particles/Deposition/SharedDepositionUtils.H"
12#include "Particles/Pusher/GetAndSetPosition.H"
13#include "Particles/Pusher/UpdatePosition.H"
14#include "Particles/Gather/FieldGather.H"
15#include "Particles/ShapeFactors.H"
16#include "Utils/TextMsg.H"
17#include "Utils/WarpXAlgorithmSelection.H"
18#include "Utils/WarpXConst.H"
19#ifdef WARPX_DIM_RZ
20# include "Utils/WarpX_Complex.H"
21#endif
22
23#include <AMReX.H>
24#include <AMReX_Arena.H>
25#include <AMReX_Array4.H>
26#include <AMReX_Dim3.H>
27#include <AMReX_REAL.H>
28
41AMREX_GPU_HOST_DEVICE AMREX_INLINE
42void setMassMatricesKernels ( const amrex::ParticleReal qs,
43 const amrex::ParticleReal ms,
44 const amrex::ParticleReal dt,
45 const amrex::ParticleReal rhop,
46 const amrex::ParticleReal uxp,
47 const amrex::ParticleReal uyp,
48 const amrex::ParticleReal uzp,
49 const amrex::ParticleReal Bxp,
50 const amrex::ParticleReal Byp,
51 const amrex::ParticleReal Bzp,
52 amrex::ParticleReal& fpxx,
53 amrex::ParticleReal& fpxy,
54 amrex::ParticleReal& fpxz,
55 amrex::ParticleReal& fpyx,
56 amrex::ParticleReal& fpyy,
57 amrex::ParticleReal& fpyz,
58 amrex::ParticleReal& fpzx,
59 amrex::ParticleReal& fpzy,
60 amrex::ParticleReal& fpzz )
61{
62 using namespace amrex::literals;
63
64 constexpr amrex::ParticleReal inv_c2 = 1._prt/(PhysConst::c*PhysConst::c);
65
66 // Convert B on particle to normalized cyclotron units with dt/2.0
67 const amrex::ParticleReal gamma_bar = std::sqrt(1._prt + (uxp*uxp + uyp*uyp + uzp*uzp)*inv_c2);
68 const amrex::ParticleReal alpha = qs/ms*0.5_prt*dt/gamma_bar;
69 const amrex::ParticleReal bxp = alpha*Bxp;
70 const amrex::ParticleReal byp = alpha*Byp;
71 const amrex::ParticleReal bzp = alpha*Bzp;
72
73 // Compute Mass Matrix kernels (non-relativistic for now)
74 amrex::ParticleReal bpsq = bxp*bxp + byp*byp + bzp*bzp;
75 amrex::ParticleReal arogp = alpha*rhop/(1.0_prt + bpsq);
76
77 fpxx = arogp*(bxp*bxp + 1.0_rt);
78 fpxy = arogp*(bxp*byp + bzp);
79 fpxz = arogp*(bxp*bzp - byp);
80
81 fpyx = arogp*(byp*bxp - bzp);
82 fpyy = arogp*(byp*byp + 1.0_rt);
83 fpyz = arogp*(byp*bzp + bxp);
84
85 fpzx = arogp*(bzp*bxp + byp);
86 fpzy = arogp*(bzp*byp - bxp);
87 fpzz = arogp*(bzp*bzp + 1.0_rt);
88
89}
90
109template <int depos_order, bool full_mass_matrices>
110AMREX_GPU_HOST_DEVICE AMREX_INLINE
111void doDirectJandSigmaDepositionKernel ( [[maybe_unused]] const amrex::ParticleReal xp,
112 [[maybe_unused]] const amrex::ParticleReal yp,
113 [[maybe_unused]] const amrex::ParticleReal zp,
114 const amrex::ParticleReal wqx,
115 const amrex::ParticleReal wqy,
116 const amrex::ParticleReal wqz,
117 const amrex::ParticleReal fpxx,
118 [[maybe_unused]] const amrex::ParticleReal fpxy,
119 [[maybe_unused]] const amrex::ParticleReal fpxz,
120 [[maybe_unused]] const amrex::ParticleReal fpyx,
121 const amrex::ParticleReal fpyy,
122 [[maybe_unused]] const amrex::ParticleReal fpyz,
123 [[maybe_unused]] const amrex::ParticleReal fpzx,
124 [[maybe_unused]] const amrex::ParticleReal fpzy,
125 const amrex::ParticleReal fpzz,
126 amrex::Array4<amrex::Real> const& jx_arr,
127 amrex::Array4<amrex::Real> const& jy_arr,
128 amrex::Array4<amrex::Real> const& jz_arr,
129 amrex::Array4<amrex::Real> const& Sxx_arr,
130 [[maybe_unused]] amrex::Array4<amrex::Real> const& Sxy_arr,
131 [[maybe_unused]] amrex::Array4<amrex::Real> const& Sxz_arr,
132 [[maybe_unused]] amrex::Array4<amrex::Real> const& Syx_arr,
133 amrex::Array4<amrex::Real> const& Syy_arr,
134 [[maybe_unused]] amrex::Array4<amrex::Real> const& Syz_arr,
135 [[maybe_unused]] amrex::Array4<amrex::Real> const& Szx_arr,
136 [[maybe_unused]] amrex::Array4<amrex::Real> const& Szy_arr,
137 amrex::Array4<amrex::Real> const& Szz_arr,
138 const amrex::IntVect& jx_type,
139 const amrex::IntVect& jy_type,
140 const amrex::IntVect& jz_type,
141 const amrex::XDim3& dinv,
142 const amrex::XDim3& xyzmin,
143 const amrex::Dim3 lo )
144{
145 using namespace amrex::literals;
146
147 constexpr int NODE = amrex::IndexType::NODE;
148 constexpr int CELL = amrex::IndexType::CELL;
149
150 // MassMatrices index shift parameter
151 amrex::IntVect shift = amrex::IntVect::TheZeroVector();
152
153 // --- Compute shape factors
154 Compute_shape_factor< depos_order > const compute_shape_factor;
155#if !defined(WARPX_DIM_1D_Z)
156 // x direction
157 // Get particle position after 1/2 push back in position
158 // Keep these double to avoid bug in single precision
159 const double xmid = (xp - xyzmin.x)*dinv.x;
160
161 // j_j[xyz] leftmost grid point in x that the particle touches for the centering of each current
162 // sx_j[xyz] shape factor along x for the centering of each current
163 // There are only two possible centerings, node or cell centered, so at most only two shape factor
164 // arrays will be needed.
165 // Keep these double to avoid bug in single precision
166 double sx_node[depos_order + 1] = {0.};
167 double sx_cell[depos_order + 1] = {0.};
168 int j_node = 0;
169 int j_cell = 0;
170 if (jx_type[0] == NODE || jy_type[0] == NODE || jz_type[0] == NODE) {
171 j_node = compute_shape_factor(sx_node, xmid);
172 }
173 if (jx_type[0] == CELL || jy_type[0] == CELL || jz_type[0] == CELL) {
174 j_cell = compute_shape_factor(sx_cell, xmid - 0.5);
175 }
176
177 // Set the index shift parameter
178 if (j_node==j_cell) { shift[0] = 1; }
179
180 amrex::Real sx_jx[depos_order + 1] = {0._rt};
181 amrex::Real sx_jy[depos_order + 1] = {0._rt};
182 amrex::Real sx_jz[depos_order + 1] = {0._rt};
183 for (int ix=0; ix<=depos_order; ix++)
184 {
185 sx_jx[ix] = ((jx_type[0] == NODE) ? amrex::Real(sx_node[ix]) : amrex::Real(sx_cell[ix]));
186 sx_jy[ix] = ((jy_type[0] == NODE) ? amrex::Real(sx_node[ix]) : amrex::Real(sx_cell[ix]));
187 sx_jz[ix] = ((jz_type[0] == NODE) ? amrex::Real(sx_node[ix]) : amrex::Real(sx_cell[ix]));
188 }
189
190 int const j_jx = ((jx_type[0] == NODE) ? j_node : j_cell);
191 int const j_jy = ((jy_type[0] == NODE) ? j_node : j_cell);
192 int const j_jz = ((jz_type[0] == NODE) ? j_node : j_cell);
193#endif
194
195#if defined(WARPX_DIM_3D)
196 // y direction
197 // Keep these double to avoid bug in single precision
198 const double ymid = (yp - xyzmin.y)*dinv.y;
199 double sy_node[depos_order + 1] = {0.};
200 double sy_cell[depos_order + 1] = {0.};
201 int k_node = 0;
202 int k_cell = 0;
203 if (jx_type[1] == NODE || jy_type[1] == NODE || jz_type[1] == NODE) {
204 k_node = compute_shape_factor(sy_node, ymid);
205 }
206 if (jx_type[1] == CELL || jy_type[1] == CELL || jz_type[1] == CELL) {
207 k_cell = compute_shape_factor(sy_cell, ymid - 0.5);
208 }
209
210 // Set the index shift parameter
211 if (k_node==k_cell) { shift[1] = 1; }
212
213 amrex::Real sy_jx[depos_order + 1] = {0._rt};
214 amrex::Real sy_jy[depos_order + 1] = {0._rt};
215 amrex::Real sy_jz[depos_order + 1] = {0._rt};
216 for (int iy=0; iy<=depos_order; iy++)
217 {
218 sy_jx[iy] = ((jx_type[1] == NODE) ? amrex::Real(sy_node[iy]) : amrex::Real(sy_cell[iy]));
219 sy_jy[iy] = ((jy_type[1] == NODE) ? amrex::Real(sy_node[iy]) : amrex::Real(sy_cell[iy]));
220 sy_jz[iy] = ((jz_type[1] == NODE) ? amrex::Real(sy_node[iy]) : amrex::Real(sy_cell[iy]));
221 }
222 int const k_jx = ((jx_type[1] == NODE) ? k_node : k_cell);
223 int const k_jy = ((jy_type[1] == NODE) ? k_node : k_cell);
224 int const k_jz = ((jz_type[1] == NODE) ? k_node : k_cell);
225#endif
226
227#if !defined(WARPX_DIM_RCYLINDER) && !defined(WARPX_DIM_RSPHERE)
228 // z direction
229 // Keep these double to avoid bug in single precision
230 constexpr int zdir = WARPX_ZINDEX;
231 const double zmid = (zp - xyzmin.z)*dinv.z;
232 double sz_node[depos_order + 1] = {0.};
233 double sz_cell[depos_order + 1] = {0.};
234 int l_node = 0;
235 int l_cell = 0;
236 if (jx_type[zdir] == NODE || jy_type[zdir] == NODE || jz_type[zdir] == NODE) {
237 l_node = compute_shape_factor(sz_node, zmid);
238 }
239 if (jx_type[zdir] == CELL || jy_type[zdir] == CELL || jz_type[zdir] == CELL) {
240 l_cell = compute_shape_factor(sz_cell, zmid - 0.5);
241 }
242 amrex::Real sz_jx[depos_order + 1] = {0._rt};
243 amrex::Real sz_jy[depos_order + 1] = {0._rt};
244 amrex::Real sz_jz[depos_order + 1] = {0._rt};
245 for (int iz=0; iz<=depos_order; iz++)
246 {
247 sz_jx[iz] = ((jx_type[zdir] == NODE) ? amrex::Real(sz_node[iz]) : amrex::Real(sz_cell[iz]));
248 sz_jy[iz] = ((jy_type[zdir] == NODE) ? amrex::Real(sz_node[iz]) : amrex::Real(sz_cell[iz]));
249 sz_jz[iz] = ((jz_type[zdir] == NODE) ? amrex::Real(sz_node[iz]) : amrex::Real(sz_cell[iz]));
250 }
251 int const l_jx = ((jx_type[zdir] == NODE) ? l_node : l_cell);
252 int const l_jy = ((jy_type[zdir] == NODE) ? l_node : l_cell);
253 int const l_jz = ((jz_type[zdir] == NODE) ? l_node : l_cell);
254
255 // Set the index shift parameter
256 if (l_node==l_cell) { shift[zdir] = 1; }
257
258#endif
259
260 // Compute index offset needed when x and y comps have different location on grid
261 amrex::IntVect offset_xy, offset_xz, offset_yz;
262 for (int dir=0; dir<AMREX_SPACEDIM; dir++) {
263 offset_xy[dir] = (jx_type[dir] + jy_type[dir]) % 2;
264 offset_xz[dir] = (jx_type[dir] + jz_type[dir]) % 2;
265 offset_yz[dir] = (jy_type[dir] + jz_type[dir]) % 2;
266 }
267
268 // Deposit J and mass matrices
269#if defined(WARPX_DIM_1D_Z)
270 for (int iz=0; iz<=depos_order; iz++){
271 amrex::Gpu::Atomic::AddNoRet(
272 &jx_arr(lo.x+l_jx+iz, 0, 0, 0),
273 sz_jx[iz]*wqx);
274 amrex::Gpu::Atomic::AddNoRet(
275 &jy_arr(lo.x+l_jy+iz, 0, 0, 0),
276 sz_jy[iz]*wqy);
277 amrex::Gpu::Atomic::AddNoRet(
278 &jz_arr(lo.x+l_jz+iz, 0, 0, 0),
279 sz_jz[iz]*wqz);
280
281 // Deposit mass matrices
282 if constexpr (full_mass_matrices) {
283 for (int aa=0; aa<=depos_order; aa++){
284 // Deposit mass matrices for X-current
285 int Nc = depos_order + aa - iz;
286 amrex::Gpu::Atomic::AddNoRet(
287 &Sxx_arr(lo.x+l_jx+iz, 0, 0, Nc),
288 sz_jx[iz]*sz_jx[aa]*fpxx);
289 Nc = depos_order + shift[0]*offset_xy[0] + aa - iz;
290 amrex::Gpu::Atomic::AddNoRet(
291 &Sxy_arr(lo.x+l_jx+iz, 0, 0, Nc),
292 sz_jx[iz]*sz_jy[aa]*fpxy);
293 Nc = depos_order + shift[0]*offset_xz[0] + aa - iz;
294 amrex::Gpu::Atomic::AddNoRet(
295 &Sxz_arr(lo.x+l_jx+iz, 0, 0, Nc),
296 sz_jx[iz]*sz_jz[aa]*fpxz);
297
298 // Deposit mass matrices for Y-current
299 Nc = depos_order + shift[0]*offset_xy[0] + aa - iz;
300 amrex::Gpu::Atomic::AddNoRet(
301 &Syx_arr(lo.x+l_jy+iz, 0, 0, Nc),
302 sz_jy[iz]*sz_jx[aa]*fpyx);
303 Nc = depos_order + aa - iz;
304 amrex::Gpu::Atomic::AddNoRet(
305 &Syy_arr(lo.x+l_jy+iz, 0, 0, Nc),
306 sz_jy[iz]*sz_jy[aa]*fpyy);
307 Nc = depos_order + shift[0]*offset_yz[0] + aa - iz;
308 amrex::Gpu::Atomic::AddNoRet(
309 &Syz_arr(lo.x+l_jy+iz, 0, 0, Nc),
310 sz_jy[iz]*sz_jz[aa]*fpyz);
311
312 // Deposit mass matrices for Z-current
313 Nc = depos_order + 1 - shift[0]*offset_xz[0] + aa - iz;
314 amrex::Gpu::Atomic::AddNoRet(
315 &Szx_arr(lo.x+l_jz+iz, 0, 0, Nc),
316 sz_jz[iz]*sz_jx[aa]*fpzx);
317 Nc = depos_order + 1 - shift[0]*offset_yz[0] + aa - iz;
318 amrex::Gpu::Atomic::AddNoRet(
319 &Szy_arr(lo.x+l_jz+iz, 0, 0, Nc),
320 sz_jz[iz]*sz_jy[aa]*fpzy);
321 Nc = depos_order + aa - iz;
322 amrex::Gpu::Atomic::AddNoRet(
323 &Szz_arr(lo.x+l_jz+iz, 0, 0, Nc),
324 sz_jz[iz]*sz_jz[aa]*fpzz);
325 }
326 }
327 else { // Deposit mass matrices for diagonal PC only
328 amrex::Gpu::Atomic::AddNoRet(
329 &Sxx_arr(lo.x+l_jx+iz, 0, 0, 0),
330 sz_jx[iz]*fpxx);
331 amrex::Gpu::Atomic::AddNoRet(
332 &Syy_arr(lo.x+l_jy+iz, 0, 0, 0),
333 sz_jy[iz]*fpyy);
334 amrex::Gpu::Atomic::AddNoRet(
335 &Szz_arr(lo.x+l_jz+iz, 0, 0, 0),
336 sz_jz[iz]*fpzz);
337 }
338 }
339#elif defined(WARPX_DIM_RCYLINDER) || defined(WARPX_DIM_RSPHERE)
340 for (int ix=0; ix<=depos_order; ix++){
341 amrex::Gpu::Atomic::AddNoRet(
342 &jx_arr(lo.x+j_jx+ix, 0, 0, 0),
343 sx_jx[ix]*wqx);
344 amrex::Gpu::Atomic::AddNoRet(
345 &jy_arr(lo.x+j_jy+ix, 0, 0, 0),
346 sx_jy[ix]*wqy);
347 amrex::Gpu::Atomic::AddNoRet(
348 &jz_arr(lo.x+j_jz+ix, 0, 0, 0),
349 sx_jz[ix]*wqz);
350 //
351 amrex::Gpu::Atomic::AddNoRet(
352 &Sxx_arr(lo.x+j_jx+ix, 0, 0, 0),
353 sx_jx[ix]*sx_jx[ix]*fpxx);
354 amrex::Gpu::Atomic::AddNoRet(
355 &Syy_arr(lo.x+j_jy+ix, 0, 0, 0),
356 sx_jy[ix]*sx_jy[ix]*fpyy);
357 amrex::Gpu::Atomic::AddNoRet(
358 &Szz_arr(lo.x+j_jz+ix, 0, 0, 0),
359 sx_jz[ix]*sx_jz[ix]*fpzz);
360 }
361#elif defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ)
362 const int base_offset = 1 + 2*depos_order;
363 for (int iz=0; iz<=depos_order; iz++){
364 for (int ix=0; ix<=depos_order; ix++){
365 const amrex::Real weight_Jx = sx_jx[ix]*sz_jx[iz];
366 const amrex::Real weight_Jy = sx_jy[ix]*sz_jy[iz];
367 const amrex::Real weight_Jz = sx_jz[ix]*sz_jz[iz];
368 amrex::Gpu::Atomic::AddNoRet(
369 &jx_arr(lo.x+j_jx+ix, lo.y+l_jx+iz, 0, 0),
370 weight_Jx*wqx);
371 amrex::Gpu::Atomic::AddNoRet(
372 &jy_arr(lo.x+j_jy+ix, lo.y+l_jy+iz, 0, 0),
373 weight_Jy*wqy);
374 amrex::Gpu::Atomic::AddNoRet(
375 &jz_arr(lo.x+j_jz+ix, lo.y+l_jz+iz, 0, 0),
376 weight_Jz*wqz);
377
378 // Deposit mass matrices
379 if constexpr (full_mass_matrices) {
380 for (int bb=0; bb<=depos_order; bb++){
381 for (int aa=0; aa<=depos_order; aa++){
382 const amrex::Real weight_Ex = sx_jx[aa]*sz_jx[bb];
383 const amrex::Real weight_Ey = sx_jy[aa]*sz_jy[bb];
384 const amrex::Real weight_Ez = sx_jz[aa]*sz_jz[bb];
385
386 // Deposit mass matrices for X-current
387 int offset = base_offset;
388 int Nc = depos_order + aa - ix
389 + (depos_order + bb - iz)*offset;
390 amrex::Gpu::Atomic::AddNoRet(
391 &Sxx_arr(lo.x+j_jx+ix, lo.y+l_jx+iz, 0, Nc),
392 weight_Jx*weight_Ex*fpxx);
393 offset = base_offset + offset_xy[0];
394 Nc = depos_order + 1 - shift[0]*offset_xy[0] + aa - ix
395 + (depos_order + shift[1]*offset_xy[1] + bb - iz)*offset;
396 amrex::Gpu::Atomic::AddNoRet(
397 &Sxy_arr(lo.x+j_jx+ix, lo.y+l_jx+iz, 0, Nc),
398 weight_Jx*weight_Ey*fpxy);
399 offset = base_offset + offset_xz[0];
400 Nc = depos_order + 1 - shift[0]*offset_xz[0] + aa - ix
401 + (depos_order + shift[1]*offset_xz[1] + bb - iz)*offset;
402 amrex::Gpu::Atomic::AddNoRet(
403 &Sxz_arr(lo.x+j_jx+ix, lo.y+l_jx+iz, 0, Nc),
404 weight_Jx*weight_Ez*fpxz);
405
406 // Deposit mass matrices for Y-current
407 offset = base_offset;
408 Nc = depos_order + aa - ix
409 + (depos_order + bb - iz)*offset;
410 amrex::Gpu::Atomic::AddNoRet(
411 &Syy_arr(lo.x+j_jy+ix, lo.y+l_jy+iz, 0, Nc),
412 weight_Jy*weight_Ey*fpyy);
413 offset = base_offset + offset_xy[0];
414 Nc = depos_order + shift[0]*offset_xy[0] + aa - ix
415 + (depos_order + shift[1]*offset_xy[1] + bb - iz)*offset;
416 amrex::Gpu::Atomic::AddNoRet(
417 &Syx_arr(lo.x+j_jy+ix, lo.y+l_jy+iz, 0, Nc),
418 weight_Jy*weight_Ex*fpyx);
419 offset = base_offset + offset_yz[0];
420 Nc = depos_order + shift[0]*offset_yz[0] + aa - ix
421 + (depos_order + shift[1]*offset_yz[1] + bb - iz)*offset;
422 amrex::Gpu::Atomic::AddNoRet(
423 &Syz_arr(lo.x+j_jy+ix, lo.y+l_jy+iz, 0, Nc),
424 weight_Jy*weight_Ez*fpyz);
425
426 // Deposit mass matrices for Z-current
427 offset = base_offset;
428 Nc = depos_order + aa - ix
429 + (depos_order + bb - iz)*offset;
430 amrex::Gpu::Atomic::AddNoRet(
431 &Szz_arr(lo.x+j_jz+ix, lo.y+l_jz+iz, 0, Nc),
432 weight_Jz*weight_Ez*fpzz);
433 offset = base_offset + offset_xz[0];
434 Nc = depos_order + shift[0]*offset_xz[0] + aa - ix
435 + (depos_order + 1 - shift[1]*offset_xz[1] + bb - iz)*offset;
436 amrex::Gpu::Atomic::AddNoRet(
437 &Szx_arr(lo.x+j_jz+ix, lo.y+l_jz+iz, 0, Nc),
438 weight_Jz*weight_Ex*fpzx);
439 offset = base_offset + offset_yz[0];
440 Nc = depos_order + shift[0]*offset_yz[0] + aa - ix
441 + (depos_order + 1 - shift[1]*offset_yz[1] + bb - iz)*offset;
442 amrex::Gpu::Atomic::AddNoRet(
443 &Szy_arr(lo.x+j_jz+ix, lo.y+l_jz+iz, 0, Nc),
444 weight_Jz*weight_Ey*fpzy);
445 }
446 }
447 }
448 else { // Deposit mass matrices for diagonal PC only
449 amrex::Gpu::Atomic::AddNoRet(
450 &Syy_arr(lo.x+j_jy+ix, lo.y+l_jy+iz, 0, 0),
451 weight_Jy*fpyy);
452 amrex::Gpu::Atomic::AddNoRet(
453 &Sxx_arr(lo.x+j_jx+ix, lo.y+l_jx+iz, 0, 0),
454 weight_Jx*fpxx);
455 amrex::Gpu::Atomic::AddNoRet(
456 &Szz_arr(lo.x+j_jz+ix, lo.y+l_jz+iz, 0, 0),
457 weight_Jz*fpzz);
458 }
459 }
460 }
461#elif defined(WARPX_DIM_3D)
462 for (int iz=0; iz<=depos_order; iz++){
463 for (int iy=0; iy<=depos_order; iy++){
464 for (int ix=0; ix<=depos_order; ix++){
465 const amrex::Real weight_Jx = sx_jx[ix]*sy_jx[iy]*sz_jx[iz];
466 const amrex::Real weight_Jy = sx_jy[ix]*sy_jy[iy]*sz_jy[iz];
467 const amrex::Real weight_Jz = sx_jz[ix]*sy_jz[iy]*sz_jz[iz];
468 amrex::Gpu::Atomic::AddNoRet(
469 &jx_arr(lo.x+j_jx+ix, lo.y+k_jx+iy, lo.z+l_jx+iz),
470 weight_Jx*wqx);
471 amrex::Gpu::Atomic::AddNoRet(
472 &jy_arr(lo.x+j_jy+ix, lo.y+k_jy+iy, lo.z+l_jy+iz),
473 weight_Jy*wqy);
474 amrex::Gpu::Atomic::AddNoRet(
475 &jz_arr(lo.x+j_jz+ix, lo.y+k_jz+iy, lo.z+l_jz+iz),
476 weight_Jz*wqz);
477
478 // Deposit mass matrices
479 if constexpr (full_mass_matrices) {
480 // Should not be here. Full mass matrices not yet implemented in 3D
481 }
482 else { // Deposit mass matrices for diagonal PC only
483 amrex::Gpu::Atomic::AddNoRet(
484 &Sxx_arr(lo.x+j_jx+ix, lo.y+k_jx+iy, lo.z+l_jx+iz, 0),
485 weight_Jx*fpxx);
486 amrex::Gpu::Atomic::AddNoRet(
487 &Syy_arr(lo.x+j_jy+ix, lo.y+k_jy+iy, lo.z+l_jy+iz, 0),
488 weight_Jy*fpyy);
489 amrex::Gpu::Atomic::AddNoRet(
490 &Szz_arr(lo.x+j_jz+ix, lo.y+k_jz+iy, lo.z+l_jz+iz, 0),
491 weight_Jz*fpzz);
492 }
493 }
494 }
495 }
496#endif
497}
498
521template <int depos_order, bool full_mass_matrices>
522void doDirectJandSigmaDeposition ( const GetParticlePosition<PIdx>& GetPosition,
523 const amrex::ParticleReal* wp,
524 const amrex::ParticleReal* uxp_n,
525 const amrex::ParticleReal* uyp_n,
526 const amrex::ParticleReal* uzp_n,
527 const amrex::ParticleReal* uxp_nph,
528 const amrex::ParticleReal* uyp_nph,
529 const amrex::ParticleReal* uzp_nph,
530 amrex::FArrayBox& jx_fab,
531 amrex::FArrayBox& jy_fab,
532 amrex::FArrayBox& jz_fab,
533 amrex::Array4<amrex::Real> const& Sxx_arr,
534 amrex::Array4<amrex::Real> const& Sxy_arr,
535 amrex::Array4<amrex::Real> const& Sxz_arr,
536 amrex::Array4<amrex::Real> const& Syx_arr,
537 amrex::Array4<amrex::Real> const& Syy_arr,
538 amrex::Array4<amrex::Real> const& Syz_arr,
539 amrex::Array4<amrex::Real> const& Szx_arr,
540 amrex::Array4<amrex::Real> const& Szy_arr,
541 amrex::Array4<amrex::Real> const& Szz_arr,
542 const amrex::Array4<amrex::Real const>& Bx_arr,
543 const amrex::Array4<amrex::Real const>& By_arr,
544 const amrex::Array4<amrex::Real const>& Bz_arr,
545 const amrex::IndexType Bx_type,
546 const amrex::IndexType By_type,
547 const amrex::IndexType Bz_type,
548 const long np_to_deposit,
549 const amrex::Real dt,
550 const amrex::XDim3& dinv,
551 const amrex::XDim3& xyzmin,
552 const amrex::Dim3 lo,
553 const amrex::Real qs,
554 const amrex::Real ms )
555{
556 using namespace amrex::literals;
557
558 const amrex::Real invvol = dinv.x*dinv.y*dinv.z;
559
560 amrex::Array4<amrex::Real> const& jx_arr = jx_fab.array();
561 amrex::Array4<amrex::Real> const& jy_arr = jy_fab.array();
562 amrex::Array4<amrex::Real> const& jz_arr = jz_fab.array();
563 amrex::IntVect const jx_type = jx_fab.box().type();
564 amrex::IntVect const jy_type = jy_fab.box().type();
565 amrex::IntVect const jz_type = jz_fab.box().type();
566
567 // Loop over particles and deposit into jx_fab, jy_fab and jz_fab
568 amrex::ParallelFor(
569 np_to_deposit,
570 [=] AMREX_GPU_DEVICE (long ip) {
571 amrex::ParticleReal xp_nph, yp_nph, zp_nph;
572 GetPosition(ip, xp_nph, yp_nph, zp_nph);
573
574 // Compute magnetic field on particle
575 amrex::ParticleReal Bxp = 0.0;
576 amrex::ParticleReal Byp = 0.0;
577 amrex::ParticleReal Bzp = 0.0;
578 const int depos_order_perp = 1;
579 const int depos_order_para = 1;
580 const int n_rz_azimuthal_modes = 0;
581 doDirectGatherVectorField<depos_order_perp,depos_order_para>(
582 xp_nph, yp_nph, zp_nph,
583 Bxp, Byp, Bzp,
584 Bx_arr, By_arr, Bz_arr,
585 Bx_type, By_type, Bz_type,
586 dinv, xyzmin, lo, n_rz_azimuthal_modes );
587
588 // Compute inverse Lorentz factor, the average of gamma at time levels n and n+1
589 const amrex::ParticleReal gaminv = GetImplicitGammaInverse(uxp_n[ip], uyp_n[ip], uzp_n[ip],
590 uxp_nph[ip], uyp_nph[ip], uzp_nph[ip]);
591
592 // Compute current density kernels to deposit
593 const amrex::Real rhop = qs*wp[ip]*invvol*gaminv;
594 amrex::Real wqx = rhop*uxp_nph[ip];
595 amrex::Real wqy = rhop*uyp_nph[ip];
596 amrex::Real wqz = rhop*uzp_nph[ip];
597
598 // Set the Mass Matrices kernels
599 amrex::ParticleReal fpxx, fpxy, fpxz;
600 amrex::ParticleReal fpyx, fpyy, fpyz;
601 amrex::ParticleReal fpzx, fpzy, fpzz;
602 setMassMatricesKernels( qs, ms, dt, rhop,
603 uxp_nph[ip], uyp_nph[ip], uzp_nph[ip],
604 Bxp, Byp, Bzp,
605 fpxx, fpxy, fpxz,
606 fpyx, fpyy, fpyz,
607 fpzx, fpzy, fpzz );
608
609 doDirectJandSigmaDepositionKernel<depos_order,full_mass_matrices>(
610 xp_nph, yp_nph, zp_nph,
611 wqx, wqy, wqz,
612 fpxx, fpxy, fpxz,
613 fpyx, fpyy, fpyz,
614 fpzx, fpzy, fpzz,
615 jx_arr, jy_arr, jz_arr,
616 Sxx_arr, Sxy_arr, Sxz_arr,
617 Syx_arr, Syy_arr, Syz_arr,
618 Szx_arr, Szy_arr, Szz_arr,
619 jx_type, jy_type, jz_type,
620 dinv, xyzmin, lo );
621
622 }
623 );
624}
625
651template <int depos_order, bool full_mass_matrices>
652AMREX_GPU_HOST_DEVICE AMREX_INLINE
653void doVillasenorJandSigmaDepositionKernel ( [[maybe_unused]] const amrex::ParticleReal xp_old,
654 [[maybe_unused]] const amrex::ParticleReal yp_old,
655 [[maybe_unused]] const amrex::ParticleReal zp_old,
656 [[maybe_unused]] const amrex::ParticleReal xp_new,
657 [[maybe_unused]] const amrex::ParticleReal yp_new,
658 [[maybe_unused]] const amrex::ParticleReal zp_new,
659 const amrex::ParticleReal wq_invvol,
660 [[maybe_unused]] const amrex::ParticleReal uxp_mid,
661 [[maybe_unused]] const amrex::ParticleReal uyp_mid,
662 [[maybe_unused]] const amrex::ParticleReal uzp_mid,
663 [[maybe_unused]] const amrex::ParticleReal gaminv,
664 const amrex::ParticleReal fpxx,
665 [[maybe_unused]] const amrex::ParticleReal fpxy,
666 [[maybe_unused]] const amrex::ParticleReal fpxz,
667 [[maybe_unused]] const amrex::ParticleReal fpyx,
668 const amrex::ParticleReal fpyy,
669 [[maybe_unused]] const amrex::ParticleReal fpyz,
670 [[maybe_unused]] const amrex::ParticleReal fpzx,
671 [[maybe_unused]] const amrex::ParticleReal fpzy,
672 const amrex::ParticleReal fpzz,
673 amrex::Array4<amrex::Real> const& Jx_arr,
674 amrex::Array4<amrex::Real> const& Jy_arr,
675 amrex::Array4<amrex::Real> const& Jz_arr,
676 [[maybe_unused]] int max_crossings,
677 amrex::Array4<amrex::Real> const& Sxx_arr,
678 [[maybe_unused]] amrex::Array4<amrex::Real> const& Sxy_arr,
679 [[maybe_unused]] amrex::Array4<amrex::Real> const& Sxz_arr,
680 [[maybe_unused]] amrex::Array4<amrex::Real> const& Syx_arr,
681 amrex::Array4<amrex::Real> const& Syy_arr,
682 [[maybe_unused]] amrex::Array4<amrex::Real> const& Syz_arr,
683 [[maybe_unused]] amrex::Array4<amrex::Real> const& Szx_arr,
684 [[maybe_unused]] amrex::Array4<amrex::Real> const& Szy_arr,
685 amrex::Array4<amrex::Real> const& Szz_arr,
686 const amrex::Real dt,
687 const amrex::XDim3& dinv,
688 const amrex::XDim3& xyzmin,
689 const amrex::GpuArray<amrex::GpuArray<double,2>, AMREX_SPACEDIM> & domain_double,
690 const amrex::GpuArray<amrex::GpuArray<bool,2>, AMREX_SPACEDIM> & do_cropping,
691 const amrex::Dim3 lo )
692{
693
694 using namespace amrex::literals;
695
696#if (AMREX_SPACEDIM > 1)
697 amrex::Real constexpr one_third = 1.0_rt / 3.0_rt;
698 amrex::Real constexpr one_sixth = 1.0_rt / 6.0_rt;
699#endif
700
701 // computes current and old position in grid units
702#if defined(WARPX_DIM_RZ) || defined(WARPX_DIM_RCYLINDER)
703 amrex::Real const xp_mid = (xp_new + xp_old)*0.5_rt;
704 amrex::Real const yp_mid = (yp_new + yp_old)*0.5_rt;
705 amrex::Real const rp_new = std::sqrt(xp_new*xp_new + yp_new*yp_new);
706 amrex::Real const rp_old = std::sqrt(xp_old*xp_old + yp_old*yp_old);
707 amrex::Real const rp_mid = (rp_new + rp_old)/2._rt;
708 amrex::Real const costheta_mid = (rp_mid > 0._rt ? xp_mid/rp_mid : 1._rt);
709 amrex::Real const sintheta_mid = (rp_mid > 0._rt ? yp_mid/rp_mid : 0._rt);
710
711 // Keep these double to avoid bug in single precision
712 double x_new = (rp_new - xyzmin.x)*dinv.x;
713 double const x_old = (rp_old - xyzmin.x)*dinv.x;
714 amrex::Real const vx = (rp_new - rp_old)/dt;
715 amrex::Real const vy = (-uxp_mid*sintheta_mid + uyp_mid*costheta_mid)*gaminv;
716#if defined(WARPX_DIM_RCYLINDER)
717 amrex::Real const vz = uzp_mid*gaminv;
718#endif
719#elif defined(WARPX_DIM_RSPHERE)
720 amrex::Real const xp_mid = (xp_new + xp_old)*0.5_rt;
721 amrex::Real const yp_mid = (yp_new + yp_old)*0.5_rt;
722 amrex::Real const zp_mid = (zp_new + zp_old)*0.5_rt;
723 amrex::Real const rpxy_new = std::sqrt(xp_new*xp_new + yp_new*yp_new);
724 amrex::Real const rp_new = std::sqrt(xp_new*xp_new + yp_new*yp_new + zp_new*zp_new);
725 amrex::Real const rpxy_old = std::sqrt(xp_old*xp_old + yp_old*yp_old);
726 amrex::Real const rp_old = std::sqrt(xp_old*xp_old + yp_old*yp_old + zp_old*zp_old);
727 amrex::Real const rpxy_mid = (rpxy_new + rpxy_old)*0.5_rt;
728 amrex::Real const rp_mid = (rp_new + rp_old)*0.5_rt;
729 amrex::Real const costheta_mid = (rpxy_mid > 0._rt ? xp_mid/rpxy_mid : 1._rt);
730 amrex::Real const sintheta_mid = (rpxy_mid > 0._rt ? yp_mid/rpxy_mid : 0._rt);
731 amrex::Real const cosphi_mid = (rp_mid > 0._rt ? rpxy_mid/rp_mid : 1._rt);
732 amrex::Real const sinphi_mid = (rp_mid > 0._rt ? zp_mid/rp_mid : 0._rt);
733
734 // Keep these double to avoid bug in single precision
735 double x_new = (rp_new - xyzmin.x)*dinv.x;
736 double const x_old = (rp_old - xyzmin.x)*dinv.x;
737 amrex::Real const vx = (rp_new - rp_old)/dt;
738 amrex::Real const vy = (-uxp_mid*sintheta_mid + uyp_mid*costheta_mid)*gaminv;
739 amrex::Real const vz = (-uxp_mid*costheta_mid*cosphi_mid - uyp_mid*sintheta_mid*cosphi_mid + uzp_mid*sinphi_mid)*gaminv;
740#elif defined(WARPX_DIM_XZ)
741 // Keep these double to avoid bug in single precision
742 double x_new = (xp_new - xyzmin.x)*dinv.x;
743 double const x_old = (xp_old - xyzmin.x)*dinv.x;
744 amrex::Real const vx = (xp_new - xp_old)/dt;
745 amrex::Real const vy = uyp_mid*gaminv;
746#elif defined(WARPX_DIM_1D_Z)
747 amrex::Real const vx = uxp_mid*gaminv;
748 amrex::Real const vy = uyp_mid*gaminv;
749#elif defined(WARPX_DIM_3D)
750 // Keep these double to avoid bug in single precision
751 double x_new = (xp_new - xyzmin.x)*dinv.x;
752 double const x_old = (xp_old - xyzmin.x)*dinv.x;
753 double y_new = (yp_new - xyzmin.y)*dinv.y;
754 double const y_old = (yp_old - xyzmin.y)*dinv.y;
755 amrex::Real const vx = (xp_new - xp_old)/dt;
756 amrex::Real const vy = (yp_new - yp_old)/dt;
757#endif
758
759#if !defined(WARPX_DIM_RCYLINDER) && !defined(WARPX_DIM_RSPHERE)
760 // Keep these double to avoid bug in single precision
761 double z_new = (zp_new - xyzmin.z)*dinv.z;
762 double const z_old = (zp_old - xyzmin.z)*dinv.z;
763 amrex::Real const vz = (zp_new - zp_old)/dt;
764#endif
765
766 // Define velocity kernels to deposit
767 amrex::Real const wqx = wq_invvol*vx;
768 amrex::Real const wqy = wq_invvol*vy;
769 amrex::Real const wqz = wq_invvol*vz;
770
771 // Compute total change in particle position (always do before cropping)
772#if !defined(WARPX_DIM_1D_Z)
773 const double dxp = x_new - x_old;
774#endif
775#if defined(WARPX_DIM_3D)
776 const double dyp = y_new - y_old;
777#endif
778#if !defined(WARPX_DIM_RCYLINDER) && !defined(WARPX_DIM_RSPHERE)
779 const double dzp = z_new - z_old;
780#endif
781
782 // Crop particle orbits at absorbing domain boundaries
783#if defined(WARPX_DIM_3D)
784 ParticleUtils::crop_at_boundary(x_old, y_old, z_old, x_new, y_new, z_new,
785 domain_double[0][0], domain_double[0][1], do_cropping[0]);
786 ParticleUtils::crop_at_boundary(y_old, z_old, x_old, y_new, z_new, x_new,
787 domain_double[1][0], domain_double[1][1], do_cropping[1]);
788 ParticleUtils::crop_at_boundary(z_old, x_old, y_old, z_new, x_new, y_new,
789 domain_double[2][0], domain_double[2][1], do_cropping[2]);
790#elif defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ)
791 ParticleUtils::crop_at_boundary(x_old, z_old, x_new, z_new,
792 domain_double[0][0], domain_double[0][1], do_cropping[0]);
793 ParticleUtils::crop_at_boundary(z_old, x_old, z_new, x_new,
794 domain_double[1][0], domain_double[1][1], do_cropping[1]);
795#elif defined(WARPX_DIM_1D_Z)
796 ParticleUtils::crop_at_boundary(z_new, domain_double[0][0], domain_double[0][1], do_cropping[0]);
797#elif defined(WARPX_DIM_RCYLINDER) || defined(WARPX_DIM_RSPHERE)
798 ParticleUtils::crop_at_boundary(x_new, domain_double[0][0], domain_double[0][1], do_cropping[0]);
799#endif
800
801 // 1) Determine the number of segments.
802 // 2) Loop over segments and deposit current.
803
804 // cell crossings are defined at cell edges if depos_order is odd
805 // cell crossings are defined at cell centers if depos_order is even
806
807 int num_segments = 1;
808 double shift = 0.0;
809 if ( (depos_order % 2) == 0 ) { shift = 0.5; }
810
811#if defined(WARPX_DIM_3D)
812
813 // compute cell crossings in X-direction
814 const auto i_old = static_cast<int>(x_old-shift);
815 const auto i_new = static_cast<int>(x_new-shift);
816 const int cell_crossings_x = std::abs(i_new-i_old);
817 num_segments += cell_crossings_x;
818
819 // compute cell crossings in Y-direction
820 const auto j_old = static_cast<int>(y_old-shift);
821 const auto j_new = static_cast<int>(y_new-shift);
822 const int cell_crossings_y = std::abs(j_new-j_old);
823 num_segments += cell_crossings_y;
824
825 // compute cell crossings in Z-direction
826 const auto k_old = static_cast<int>(z_old-shift);
827 const auto k_new = static_cast<int>(z_new-shift);
828 const int cell_crossings_z = std::abs(k_new-k_old);
829 num_segments += cell_crossings_z;
830
831 // Compute initial particle cell locations in each direction
832 // used to find the position at cell crossings.
833 // Keep these double to avoid bug in single precision
834 const auto dirX_sign = static_cast<double>(dxp < 0. ? -1. : 1.);
835 const auto dirY_sign = static_cast<double>(dyp < 0. ? -1. : 1.);
836 const auto dirZ_sign = static_cast<double>(dzp < 0. ? -1. : 1.);
837 double Xcell = 0., Ycell = 0., Zcell = 0.;
838 if (num_segments > 1) {
839 Xcell = static_cast<double>(i_old) + shift + 0.5*(1.-dirX_sign);
840 Ycell = static_cast<double>(j_old) + shift + 0.5*(1.-dirY_sign);
841 Zcell = static_cast<double>(k_old) + shift + 0.5*(1.-dirZ_sign);
842 }
843
844 // loop over the number of segments and deposit
845 const Compute_shape_factor< depos_order-1 > compute_shape_factor_cell;
846 const Compute_shape_factor_pair< depos_order > compute_shape_factors_node;
847 double dxp_seg, dyp_seg, dzp_seg;
848 double x0_new, y0_new, z0_new;
849 double x0_old = x_old;
850 double y0_old = y_old;
851 double z0_old = z_old;
852
853 for (int ns=0; ns<num_segments; ns++) {
854
855 if (ns == num_segments-1) { // final segment
856
857 x0_new = x_new;
858 y0_new = y_new;
859 z0_new = z_new;
860 dxp_seg = x0_new - x0_old;
861 dyp_seg = y0_new - y0_old;
862 dzp_seg = z0_new - z0_old;
863
864 }
865 else {
866
867 x0_new = Xcell + dirX_sign;
868 y0_new = Ycell + dirY_sign;
869 z0_new = Zcell + dirZ_sign;
870 dxp_seg = x0_new - x0_old;
871 dyp_seg = y0_new - y0_old;
872 dzp_seg = z0_new - z0_old;
873
874 if ( (dyp == 0. || std::abs(dxp_seg) < std::abs(dxp/dyp*dyp_seg))
875 && (dzp == 0. || std::abs(dxp_seg) < std::abs(dxp/dzp*dzp_seg)) ) {
876 Xcell = x0_new;
877 dyp_seg = dyp/dxp*dxp_seg;
878 dzp_seg = dzp/dxp*dxp_seg;
879 y0_new = y0_old + dyp_seg;
880 z0_new = z0_old + dzp_seg;
881 }
882 else if (dzp == 0. || std::abs(dyp_seg) < std::abs(dyp/dzp*dzp_seg)) {
883 Ycell = y0_new;
884 dxp_seg = dxp/dyp*dyp_seg;
885 dzp_seg = dzp/dyp*dyp_seg;
886 x0_new = x0_old + dxp_seg;
887 z0_new = z0_old + dzp_seg;
888 }
889 else {
890 Zcell = z0_new;
891 dxp_seg = dxp/dzp*dzp_seg;
892 dyp_seg = dyp/dzp*dzp_seg;
893 x0_new = x0_old + dxp_seg;
894 y0_new = y0_old + dyp_seg;
895 }
896
897 }
898
899 // Compute the segment factors (each equal to dt_seg/dt for nonzero dxp, dyp, or dzp)
900 const auto seg_factor_x = static_cast<amrex::Real>(dxp == 0. ? 1._rt : dxp_seg/dxp);
901 const auto seg_factor_y = static_cast<amrex::Real>(dyp == 0. ? 1._rt : dyp_seg/dyp);
902 const auto seg_factor_z = static_cast<amrex::Real>(dzp == 0. ? 1._rt : dzp_seg/dzp);
903
904 // Compute cell-based weights using the average segment position
905 // Keep these double to avoid bug in single precision
906 double sx_cell[depos_order] = {0.};
907 double sy_cell[depos_order] = {0.};
908 double sz_cell[depos_order] = {0.};
909 double const x0_bar = (x0_new + x0_old)/2.0;
910 double const y0_bar = (y0_new + y0_old)/2.0;
911 double const z0_bar = (z0_new + z0_old)/2.0;
912 const int i0_cell = compute_shape_factor_cell( sx_cell, x0_bar-0.5 );
913 const int j0_cell = compute_shape_factor_cell( sy_cell, y0_bar-0.5 );
914 const int k0_cell = compute_shape_factor_cell( sz_cell, z0_bar-0.5 );
915
916 if constexpr (depos_order >= 3) { // higher-order correction to the cell-based weights
917 const Compute_shape_factor_pair<depos_order-1> compute_shape_factors_cell;
918 double sx_old_cell[depos_order] = {0.};
919 double sx_new_cell[depos_order] = {0.};
920 double sy_old_cell[depos_order] = {0.};
921 double sy_new_cell[depos_order] = {0.};
922 double sz_old_cell[depos_order] = {0.};
923 double sz_new_cell[depos_order] = {0.};
924 const int i0_cell_2 = compute_shape_factors_cell( sx_old_cell, sx_new_cell, x0_old-0.5, x0_new-0.5 );
925 const int j0_cell_2 = compute_shape_factors_cell( sy_old_cell, sy_new_cell, y0_old-0.5, y0_new-0.5 );
926 const int k0_cell_2 = compute_shape_factors_cell( sz_old_cell, sz_new_cell, z0_old-0.5, z0_new-0.5 );
927 amrex::ignore_unused(i0_cell_2, j0_cell_2, k0_cell_2);
928 for (int m=0; m<depos_order; m++) {
929 sx_cell[m] = (4.0*sx_cell[m] + sx_old_cell[m] + sx_new_cell[m])/6.0;
930 sy_cell[m] = (4.0*sy_cell[m] + sy_old_cell[m] + sy_new_cell[m])/6.0;
931 sz_cell[m] = (4.0*sz_cell[m] + sz_old_cell[m] + sz_new_cell[m])/6.0;
932 }
933 }
934
935 // Compute node-based weights using the old and new segment positions
936 // Keep these double to avoid bug in single precision
937 double sx_old_node[depos_order+1] = {0.};
938 double sx_new_node[depos_order+1] = {0.};
939 double sy_old_node[depos_order+1] = {0.};
940 double sy_new_node[depos_order+1] = {0.};
941 double sz_old_node[depos_order+1] = {0.};
942 double sz_new_node[depos_order+1] = {0.};
943 const int i0_node = compute_shape_factors_node( sx_old_node, sx_new_node, x0_old, x0_new );
944 const int j0_node = compute_shape_factors_node( sy_old_node, sy_new_node, y0_old, y0_new );
945 const int k0_node = compute_shape_factors_node( sz_old_node, sz_new_node, z0_old, z0_new );
946
947 // deposit Jx and Sxx for this segment
948 amrex::Real weight;
949 for (int i=0; i<=depos_order-1; i++) {
950 for (int j=0; j<=depos_order; j++) {
951 for (int k=0; k<=depos_order; k++) {
952 weight = sx_cell[i]*( sy_old_node[j]*sz_old_node[k]*one_third
953 + sy_old_node[j]*sz_new_node[k]*one_sixth
954 + sy_new_node[j]*sz_old_node[k]*one_sixth
955 + sy_new_node[j]*sz_new_node[k]*one_third )*seg_factor_x;
956 amrex::Gpu::Atomic::AddNoRet( &Jx_arr(lo.x+i0_cell+i, lo.y+j0_node+j, lo.z+k0_node+k), wqx*weight);
957 amrex::Gpu::Atomic::AddNoRet( &Sxx_arr(lo.x+i0_cell+i, lo.y+j0_node+j, lo.z+k0_node+k, 0), fpxx*weight*weight);
958 }
959 }
960 }
961
962 // deposit Jy and Syy or this segment
963 for (int i=0; i<=depos_order; i++) {
964 for (int j=0; j<=depos_order-1; j++) {
965 for (int k=0; k<=depos_order; k++) {
966 weight = sy_cell[j]*( sx_old_node[i]*sz_old_node[k]*one_third
967 + sx_old_node[i]*sz_new_node[k]*one_sixth
968 + sx_new_node[i]*sz_old_node[k]*one_sixth
969 + sx_new_node[i]*sz_new_node[k]*one_third )*seg_factor_y;
970 amrex::Gpu::Atomic::AddNoRet( &Jy_arr(lo.x+i0_node+i, lo.y+j0_cell+j, lo.z+k0_node+k), wqy*weight);
971 amrex::Gpu::Atomic::AddNoRet( &Syy_arr(lo.x+i0_node+i, lo.y+j0_cell+j, lo.z+k0_node+k, 0), fpyy*weight*weight);
972 }
973 }
974 }
975
976 // deposit Jz and Sz for this segment
977 for (int i=0; i<=depos_order; i++) {
978 for (int j=0; j<=depos_order; j++) {
979 for (int k=0; k<=depos_order-1; k++) {
980 weight = sz_cell[k]*( sx_old_node[i]*sy_old_node[j]*one_third
981 + sx_old_node[i]*sy_new_node[j]*one_sixth
982 + sx_new_node[i]*sy_old_node[j]*one_sixth
983 + sx_new_node[i]*sy_new_node[j]*one_third )*seg_factor_z;
984 amrex::Gpu::Atomic::AddNoRet( &Jz_arr(lo.x+i0_node+i, lo.y+j0_node+j, lo.z+k0_cell+k), wqz*weight);
985 amrex::Gpu::Atomic::AddNoRet( &Szz_arr(lo.x+i0_node+i, lo.y+j0_node+j, lo.z+k0_cell+k, 0), fpzz*weight*weight);
986 }
987 }
988 }
989
990 // update old segment values
991 if (ns < num_segments-1) {
992 x0_old = x0_new;
993 y0_old = y0_new;
994 z0_old = z0_new;
995 }
996
997 } // end loop over segments
998
999#elif defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ)
1000
1001 // compute cell crossings in X-direction
1002 const auto i_old = static_cast<int>(x_old-shift);
1003 const auto i_new = static_cast<int>(x_new-shift);
1004 const int cell_crossings_x = std::abs(i_new-i_old);
1005 num_segments += cell_crossings_x;
1006
1007 // compute cell crossings in Z-direction
1008 const auto k_old = static_cast<int>(z_old-shift);
1009 const auto k_new = static_cast<int>(z_new-shift);
1010 const int cell_crossings_z = std::abs(k_new-k_old);
1011 num_segments += cell_crossings_z;
1012
1013 // Compute initial particle cell locations in each direction
1014 // used to find the position at cell crossings.
1015 // Keep these double to avoid bug in single precision
1016 const auto dirX_sign = static_cast<double>(dxp < 0. ? -1. : 1.);
1017 const auto dirZ_sign = static_cast<double>(dzp < 0. ? -1. : 1.);
1018 double Xcell = 0., Zcell = 0.;
1019 if (num_segments > 1) {
1020 Xcell = static_cast<double>(i_old) + shift + 0.5*(1.-dirX_sign);
1021 Zcell = static_cast<double>(k_old) + shift + 0.5*(1.-dirZ_sign);
1022 }
1023
1024 // loop over the number of segments and deposit
1025 const Compute_shape_factor< depos_order-1 > compute_shape_factor_cell;
1026 const Compute_shape_factor_pair< depos_order > compute_shape_factors_node;
1027 double dxp_seg, dzp_seg;
1028 double x0_new, z0_new;
1029 double x0_old = x_old;
1030 double z0_old = z_old;
1031
1032 constexpr int num_segments_max = 1 + 4*AMREX_SPACEDIM;
1033 AMREX_ALWAYS_ASSERT_WITH_MESSAGE( num_segments <= num_segments_max,
1034 "Error: num_segments must be less than or equal to 1 + 4*AMREX_SPACEDIM.");
1035
1036 // Save the start index and interpolation weights for each segment
1037 int i0_cell[num_segments_max];
1038 int i0_node[num_segments_max];
1039 int k0_cell[num_segments_max];
1040 int k0_node[num_segments_max];
1041 amrex::Real weight_cellX_nodeZ[num_segments_max][depos_order][depos_order+1];
1042 amrex::Real weight_nodeX_cellZ[num_segments_max][depos_order+1][depos_order];
1043 amrex::Real weight_nodeX_nodeZ[num_segments_max][depos_order+1][depos_order+1];
1044
1045 const auto i_mid = static_cast<int>(0.5*(x_new+x_old)-shift);
1046 const auto k_mid = static_cast<int>(0.5*(z_new+z_old)-shift);
1047 int SegNumX[num_segments_max];
1048 int SegNumZ[num_segments_max];
1049
1050 for (int ns=0; ns<num_segments; ns++) {
1051
1052 if (ns == num_segments-1) { // final segment
1053
1054 x0_new = x_new;
1055 z0_new = z_new;
1056 dxp_seg = x0_new - x0_old;
1057 dzp_seg = z0_new - z0_old;
1058
1059 }
1060 else {
1061
1062 x0_new = Xcell + dirX_sign;
1063 z0_new = Zcell + dirZ_sign;
1064 dxp_seg = x0_new - x0_old;
1065 dzp_seg = z0_new - z0_old;
1066
1067 if (dzp == 0. || std::abs(dxp_seg) < std::abs(dxp/dzp*dzp_seg)) {
1068 Xcell = x0_new;
1069 dzp_seg = dzp/dxp*dxp_seg;
1070 z0_new = z0_old + dzp_seg;
1071 }
1072 else {
1073 Zcell = z0_new;
1074 dxp_seg = dxp/dzp*dzp_seg;
1075 x0_new = x0_old + dxp_seg;
1076 }
1077
1078 }
1079
1080 // Compute the segment factors (each equal to dt_seg/dt for nonzero dxp, or dzp)
1081 const auto seg_factor_x = static_cast<amrex::Real>(dxp == 0. ? 1._rt : dxp_seg/dxp);
1082 const auto seg_factor_z = static_cast<amrex::Real>(dzp == 0. ? 1._rt : dzp_seg/dzp);
1083
1084 // Compute cell-based weights using the average segment position
1085 // Keep these double to avoid bug in single precision
1086 double sx_cell[depos_order] = {0.};
1087 double sz_cell[depos_order] = {0.};
1088 double const x0_bar = (x0_new + x0_old)/2.0;
1089 double const z0_bar = (z0_new + z0_old)/2.0;
1090 i0_cell[ns] = compute_shape_factor_cell( sx_cell, x0_bar-0.5 );
1091 k0_cell[ns] = compute_shape_factor_cell( sz_cell, z0_bar-0.5 );
1092
1093 // Set the segment number for the mass matrix component calc
1094 if constexpr (full_mass_matrices) {
1095 const auto i0_mid = static_cast<int>(x0_bar-shift);
1096 const auto k0_mid = static_cast<int>(z0_bar-shift);
1097 SegNumX[ns] = 1 + i0_mid - i_mid;
1098 SegNumZ[ns] = 1 + k0_mid - k_mid;
1099 }
1100
1101 if constexpr (depos_order >= 3) { // higher-order correction to the cell-based weights
1102 const Compute_shape_factor_pair<depos_order-1> compute_shape_factors_cell;
1103 double sx_old_cell[depos_order] = {0.};
1104 double sx_new_cell[depos_order] = {0.};
1105 double sz_old_cell[depos_order] = {0.};
1106 double sz_new_cell[depos_order] = {0.};
1107 const int i0_cell_2 = compute_shape_factors_cell( sx_old_cell, sx_new_cell, x0_old-0.5, x0_new-0.5 );
1108 const int k0_cell_2 = compute_shape_factors_cell( sz_old_cell, sz_new_cell, z0_old-0.5, z0_new-0.5 );
1109 amrex::ignore_unused(i0_cell_2, k0_cell_2);
1110 for (int m=0; m<depos_order; m++) {
1111 sx_cell[m] = (4.0*sx_cell[m] + sx_old_cell[m] + sx_new_cell[m])/6.0;
1112 sz_cell[m] = (4.0*sz_cell[m] + sz_old_cell[m] + sz_new_cell[m])/6.0;
1113 }
1114 }
1115
1116 // Compute node-based weights using the old and new segment positions
1117 // Keep these double to avoid bug in single precision
1118 double sx_old_node[depos_order+1] = {0.};
1119 double sx_new_node[depos_order+1] = {0.};
1120 double sz_old_node[depos_order+1] = {0.};
1121 double sz_new_node[depos_order+1] = {0.};
1122 i0_node[ns] = compute_shape_factors_node( sx_old_node, sx_new_node, x0_old, x0_new );
1123 k0_node[ns] = compute_shape_factors_node( sz_old_node, sz_new_node, z0_old, z0_new );
1124
1125 // deposit Jx and Sx for this segment
1126 amrex::Real weight;
1127 for (int i=0; i<=depos_order-1; i++) {
1128 for (int k=0; k<=depos_order; k++) {
1129 const int i_J = lo.x + i0_cell[ns] + i;
1130 const int k_J = lo.y + k0_node[ns] + k;
1131 weight = sx_cell[i]*(sz_old_node[k] + sz_new_node[k])/2.0_rt*seg_factor_x;
1132 amrex::Gpu::Atomic::AddNoRet(&Jx_arr(i_J, k_J, 0, 0), wqx*weight);
1133 if constexpr (full_mass_matrices) { weight_cellX_nodeZ[ns][i][k] = weight; }
1134 else {
1135 amrex::Gpu::Atomic::AddNoRet(&Sxx_arr(i_J, k_J, 0, 0), fpxx*weight*weight);
1136 }
1137 }
1138 }
1139
1140 // deposit out-of-plane Jy and Sy for this segment
1141 const auto seg_factor_y = std::min(seg_factor_x,seg_factor_z);
1142 for (int i=0; i<=depos_order; i++) {
1143 for (int k=0; k<=depos_order; k++) {
1144 const int i_J = lo.x + i0_node[ns] + i;
1145 const int k_J = lo.y + k0_node[ns] + k;
1146 weight = ( sx_old_node[i]*sz_old_node[k]*one_third
1147 + sx_old_node[i]*sz_new_node[k]*one_sixth
1148 + sx_new_node[i]*sz_old_node[k]*one_sixth
1149 + sx_new_node[i]*sz_new_node[k]*one_third )*seg_factor_y;
1150 amrex::Gpu::Atomic::AddNoRet(&Jy_arr(i_J, k_J, 0, 0), wqy*weight);
1151 if constexpr (full_mass_matrices) { weight_nodeX_nodeZ[ns][i][k] = weight; }
1152 else {
1153 amrex::Gpu::Atomic::AddNoRet(&Syy_arr(i_J, k_J, 0, 0), fpyy*weight*weight);
1154 }
1155 }
1156 }
1157
1158 // deposit Jz and Szz for this segment
1159 for (int i=0; i<=depos_order; i++) {
1160 for (int k=0; k<=depos_order-1; k++) {
1161 const int i_J = lo.x + i0_node[ns] + i;
1162 const int k_J = lo.y + k0_cell[ns] + k;
1163 weight = sz_cell[k]*(sx_old_node[i] + sx_new_node[i])/2.0_rt*seg_factor_z;
1164 amrex::Gpu::Atomic::AddNoRet(&Jz_arr(i_J, k_J, 0, 0), wqz*weight);
1165 if constexpr (full_mass_matrices) { weight_nodeX_cellZ[ns][i][k] = weight; }
1166 else {
1167 amrex::Gpu::Atomic::AddNoRet(&Szz_arr(i_J, k_J, 0, 0), fpzz*weight*weight);
1168 }
1169 }
1170 }
1171
1172 // update old segment values
1173 if (ns < num_segments-1) {
1174 x0_old = x0_new;
1175 z0_old = z0_new;
1176 }
1177
1178 } // end loop over segments
1179
1180 if constexpr (full_mass_matrices) {
1181
1182 // Loop over segments and deposit full mass matrices
1183 for (int ns=0; ns<num_segments; ns++) {
1184
1185 // Deposit Sxx, Sxz, and Sxy for this segment
1186 for (int i=0; i<=depos_order-1; i++) {
1187 for (int k=0; k<=depos_order; k++) {
1188 const int i_J = lo.x + i0_cell[ns] + i;
1189 const int k_J = lo.y + k0_node[ns] + k;
1190 const amrex::Real weight_J = weight_cellX_nodeZ[ns][i][k];
1191 for (int ms=0; ms<num_segments; ms++) {
1192 const int SegShiftX = max_crossings + SegNumX[ms] - SegNumX[ns];
1193 const int SegShiftZ = max_crossings + SegNumZ[ms] - SegNumZ[ns];
1194 // Deposit Sxx
1195 const int Ncomp_xx0 = 1 + 2*(depos_order-1) + 2*max_crossings;
1196 for (int iE=0; iE<=depos_order-1; iE++) {
1197 for (int kE=0; kE<=depos_order; kE++) {
1198 const amrex::Real weight_E = weight_cellX_nodeZ[ms][iE][kE];
1199 const int comp_xx = depos_order-1 - i + iE + SegShiftX
1200 + Ncomp_xx0*(depos_order - k + kE + SegShiftZ);
1201 amrex::Gpu::Atomic::AddNoRet(&Sxx_arr(i_J, k_J, 0, comp_xx), fpxx*weight_J*weight_E);
1202 }
1203 }
1204 // Deposit Sxz
1205 const int Ncomp_xz0 = 2*depos_order + 2*max_crossings;
1206 for (int iE=0; iE<=depos_order; iE++) {
1207 for (int kE=0; kE<=depos_order-1; kE++) {
1208 const amrex::Real weight_E = weight_nodeX_cellZ[ms][iE][kE];
1209 const int comp_xz = depos_order-1 - i + iE + SegShiftX
1210 + Ncomp_xz0*(depos_order - k + kE + SegShiftZ);
1211 amrex::Gpu::Atomic::AddNoRet(&Sxz_arr(i_J, k_J, 0, comp_xz), fpxz*weight_J*weight_E);
1212 }
1213 }
1214 // Deposit Sxy
1215 const int Ncomp_xy0 = 2*depos_order + 2*max_crossings;
1216 for (int iE=0; iE<=depos_order; iE++) {
1217 for (int kE=0; kE<=depos_order; kE++) {
1218 const amrex::Real weight_E = weight_nodeX_nodeZ[ms][iE][kE];
1219 const int comp_xy = depos_order-1 - i + iE + SegShiftX
1220 + Ncomp_xy0*(depos_order - k + kE + SegShiftZ);
1221 amrex::Gpu::Atomic::AddNoRet(&Sxy_arr(i_J, k_J, 0, comp_xy), fpxy*weight_J*weight_E);
1222 }
1223 }
1224
1225 }
1226 }
1227 }
1228
1229 // Deposit Szx, Szz, and Szy for this segment
1230 for (int i=0; i<=depos_order; i++) {
1231 for (int k=0; k<=depos_order-1; k++) {
1232 const int i_J = lo.x + i0_node[ns] + i;
1233 const int k_J = lo.y + k0_cell[ns] + k;
1234 const amrex::Real weight_J = weight_nodeX_cellZ[ns][i][k];
1235 for (int ms=0; ms<num_segments; ms++) {
1236 const int SegShiftX = max_crossings + SegNumX[ms] - SegNumX[ns];
1237 const int SegShiftZ = max_crossings + SegNumZ[ms] - SegNumZ[ns];
1238 // Deposit Szx
1239 const int Ncomp_zx0 = 2*depos_order + 2*max_crossings;
1240 for (int iE=0; iE<=depos_order-1; iE++) {
1241 for (int kE=0; kE<=depos_order; kE++) {
1242 const amrex::Real weight_E = weight_cellX_nodeZ[ms][iE][kE];
1243 const int comp_zx = depos_order - i + iE + SegShiftX
1244 + Ncomp_zx0*(depos_order-1 - k + kE + SegShiftZ);
1245 amrex::Gpu::Atomic::AddNoRet( &Szx_arr(i_J, k_J, 0, comp_zx), fpzx*weight_J*weight_E);
1246 }
1247 }
1248 // Deposit Szz
1249 const int Ncomp_zz0 = 1 + 2*depos_order + 2*max_crossings;
1250 for (int iE=0; iE<=depos_order; iE++) {
1251 for (int kE=0; kE<=depos_order-1; kE++) {
1252 const amrex::Real weight_E = weight_nodeX_cellZ[ms][iE][kE];
1253 const int comp_zz = depos_order - i + iE + SegShiftX
1254 + Ncomp_zz0*(depos_order-1 - k + kE + SegShiftZ);
1255 amrex::Gpu::Atomic::AddNoRet( &Szz_arr(i_J, k_J, 0, comp_zz), fpzz*weight_J*weight_E);
1256 }
1257 }
1258 // Deposit Szy
1259 const int Ncomp_zy0 = 1 + 2*depos_order + 2*max_crossings;
1260 for (int iE=0; iE<=depos_order; iE++) {
1261 for (int kE=0; kE<=depos_order; kE++) {
1262 const amrex::Real weight_E = weight_nodeX_nodeZ[ms][iE][kE];
1263 const int comp_zy = depos_order - i + iE + SegShiftX
1264 + Ncomp_zy0*(depos_order-1 - k + kE + SegShiftZ);
1265 amrex::Gpu::Atomic::AddNoRet( &Szy_arr(i_J, k_J, 0, comp_zy), fpzy*weight_J*weight_E);
1266 }
1267 }
1268
1269 }
1270 }
1271 }
1272
1273 // Deposit Syx, Syz, and Syy for this segment
1274 for (int i=0; i<=depos_order; i++) {
1275 for (int k=0; k<=depos_order; k++) {
1276 const int i_J = lo.x + i0_node[ns] + i;
1277 const int k_J = lo.y + k0_node[ns] + k;
1278 const amrex::Real weight_J = weight_nodeX_nodeZ[ns][i][k];
1279 for (int ms=0; ms<num_segments; ms++) {
1280 const int SegShiftX = max_crossings + SegNumX[ms] - SegNumX[ns];
1281 const int SegShiftZ = max_crossings + SegNumZ[ms] - SegNumZ[ns];
1282 // Deposit Syx
1283 const int Ncomp_yx0 = 2*depos_order + 2*max_crossings;
1284 for (int iE=0; iE<=depos_order-1; iE++) {
1285 for (int kE=0; kE<=depos_order; kE++) {
1286 const amrex::Real weight_E = weight_cellX_nodeZ[ms][iE][kE];
1287 const int comp_yx = depos_order - i + iE + SegShiftX
1288 + Ncomp_yx0*(depos_order - k + kE + SegShiftZ);
1289 amrex::Gpu::Atomic::AddNoRet( &Syx_arr(i_J, k_J, 0, comp_yx), fpyx*weight_J*weight_E);
1290 }
1291 }
1292 // Deposit Syz
1293 const int Ncomp_yz0 = 1 + 2*depos_order + 2*max_crossings;
1294 for (int iE=0; iE<=depos_order; iE++) {
1295 for (int kE=0; kE<=depos_order-1; kE++) {
1296 const amrex::Real weight_E = weight_nodeX_cellZ[ms][iE][kE];
1297 const int comp_yz = depos_order - i + iE + SegShiftX
1298 + Ncomp_yz0*(depos_order - k + kE + SegShiftZ);
1299 amrex::Gpu::Atomic::AddNoRet( &Syz_arr(i_J, k_J, 0, comp_yz), fpyz*weight_J*weight_E);
1300 }
1301 }
1302 // Deposit Syy
1303 const int Ncomp_yy0 = 1 + 2*depos_order + 2*max_crossings;
1304 for (int iE=0; iE<=depos_order; iE++) {
1305 for (int kE=0; kE<=depos_order; kE++) {
1306 const amrex::Real weight_E = weight_nodeX_nodeZ[ms][iE][kE];
1307 const int comp_yy = depos_order - i + iE + SegShiftX
1308 + Ncomp_yy0*(depos_order - k + kE + SegShiftZ);
1309 amrex::Gpu::Atomic::AddNoRet( &Syy_arr(i_J, k_J, 0, comp_yy), fpyy*weight_J*weight_E);
1310 }
1311 }
1312
1313 }
1314 }
1315 }
1316
1317 }
1318
1319 }
1320
1321#elif defined(WARPX_DIM_RCYLINDER) || defined(WARPX_DIM_RSPHERE)
1322
1323 // compute cell crossings in X-direction
1324 const auto i_old = static_cast<int>(x_old-shift);
1325 const auto i_new = static_cast<int>(x_new-shift);
1326 const int cell_crossings_x = std::abs(i_new-i_old);
1327 num_segments += cell_crossings_x;
1328
1329 // Compute the initial cell location used to find the cell crossings.
1330 // Keep these double to avoid bug in single precision
1331 const auto dirX_sign = static_cast<double>(dxp < 0. ? -1. : 1.);
1332 double Xcell = static_cast<double>(i_old) + shift + 0.5*(1.-dirX_sign);
1333
1334 // loop over the number of segments and deposit
1335 const Compute_shape_factor< depos_order-1 > compute_shape_factor_cell;
1336 const Compute_shape_factor_pair< depos_order > compute_shape_factors_node;
1337 double dxp_seg;
1338 double x0_new;
1339 double x0_old = x_old;
1340
1341 for (int ns=0; ns<num_segments; ns++) {
1342
1343 if (ns == num_segments-1) { // final segment
1344 x0_new = x_new;
1345 dxp_seg = x0_new - x0_old;
1346 }
1347 else {
1348 Xcell = Xcell + dirX_sign;
1349 x0_new = Xcell;
1350 dxp_seg = x0_new - x0_old;
1351 }
1352
1353 // Compute the segment factor (equal to dt_seg/dt for nonzero dxp)
1354 const auto seg_factor = static_cast<amrex::Real>(dxp == 0. ? 1._rt : dxp_seg/dxp);
1355
1356 // Compute cell-based weights using the average segment position
1357 // Keep these double to avoid bug in single precision
1358 double sx_cell[depos_order] = {0.};
1359 double const x0_bar = (x0_new + x0_old)/2.0;
1360 const int i0_cell = compute_shape_factor_cell( sx_cell, x0_bar-0.5 );
1361
1362 if constexpr (depos_order >= 3) { // higher-order correction to the cell-based weights
1363 const Compute_shape_factor_pair<depos_order-1> compute_shape_factors_cell;
1364 double sx_old_cell[depos_order] = {0.};
1365 double sx_new_cell[depos_order] = {0.};
1366 const int i0_cell_2 = compute_shape_factors_cell( sx_old_cell, sx_new_cell, x0_old-0.5, x0_new-0.5 );
1367 amrex::ignore_unused(i0_cell_2);
1368 for (int m=0; m<depos_order; m++) {
1369 sx_cell[m] = (4.0*sx_cell[m] + sx_old_cell[m] + sx_new_cell[m])/6.0;
1370 }
1371 }
1372
1373 // Compute node-based weights using the old and new segment positions
1374 // Keep these double to avoid bug in single precision
1375 double sx_old_node[depos_order+1] = {0.};
1376 double sx_new_node[depos_order+1] = {0.};
1377 const int i0_node = compute_shape_factors_node( sx_old_node, sx_new_node, x0_old, x0_new );
1378
1379 // deposit out-of-plane Jy, Jz, Syy, and Szz for this segment
1380 for (int i=0; i<=depos_order; i++) {
1381 const amrex::Real weight = 0.5_rt*(sx_old_node[i] + sx_new_node[i])*seg_factor;
1382 amrex::Gpu::Atomic::AddNoRet( &Jy_arr(lo.x+i0_node+i, 0, 0), wqy*weight);
1383 amrex::Gpu::Atomic::AddNoRet( &Jz_arr(lo.x+i0_node+i, 0, 0), wqz*weight);
1384 //
1385 amrex::Gpu::Atomic::AddNoRet( &Syy_arr(lo.x+i0_node+i, 0, 0), fpyy*weight*weight);
1386 amrex::Gpu::Atomic::AddNoRet( &Szz_arr(lo.x+i0_node+i, 0, 0), fpzz*weight*weight);
1387 }
1388
1389 // deposit Jx and Sxx for this segment
1390 for (int i=0; i<=depos_order-1; i++) {
1391 const amrex::Real weight = sx_cell[i]*seg_factor;
1392 amrex::Gpu::Atomic::AddNoRet( &Jx_arr(lo.x+i0_cell+i, 0, 0), wqx*weight);
1393 //
1394 amrex::Gpu::Atomic::AddNoRet( &Sxx_arr(lo.x+i0_cell+i, 0, 0), fpxx*weight*weight);
1395 }
1396
1397 // update old segment values
1398 if (ns < num_segments-1) {
1399 x0_old = x0_new;
1400 }
1401
1402 }
1403
1404#elif defined(WARPX_DIM_1D_Z)
1405
1406 // compute cell crossings in Z-direction
1407 const auto k_old = static_cast<int>(z_old-shift);
1408 const auto k_new = static_cast<int>(z_new-shift);
1409 const int cell_crossings_z = std::abs(k_new-k_old);
1410 num_segments += cell_crossings_z;
1411
1412 // Compute initial particle cell location used to find cell crossings.
1413 // Keep these double to avoid bug in single precision
1414 const auto dirZ_sign = static_cast<double>(dzp < 0. ? -1. : 1.);
1415 double Zcell = static_cast<double>(k_old) + shift + 0.5*(1.-dirZ_sign);
1416
1417 // loop over the number of segments and deposit
1418 const Compute_shape_factor< depos_order-1 > compute_shape_factor_cell;
1419 const Compute_shape_factor_pair< depos_order > compute_shape_factors_node;
1420 double dzp_seg;
1421 double z0_new;
1422 double z0_old = z_old;
1423
1424 constexpr int num_segments_max = 1 + 4*AMREX_SPACEDIM;
1425 AMREX_ALWAYS_ASSERT_WITH_MESSAGE( num_segments <= num_segments_max,
1426 "Error: num_segments must be less than or equal to 1 + 4*AMREX_SPACEDIM.");
1427
1428 // Save the start index and interpolation weights for each segment
1429 int k0_cell[num_segments_max];
1430 int k0_node[num_segments_max];
1431 amrex::Real weight_cell[num_segments_max][depos_order];
1432 amrex::Real weight_node[num_segments_max][depos_order+1];
1433
1434 const auto k_mid = static_cast<int>(0.5*(z_new+z_old)-shift);
1435 int SegNum[num_segments_max];
1436
1437 for (int ns=0; ns<num_segments; ns++) {
1438
1439 if (ns == num_segments-1) { // final segment
1440 z0_new = z_new;
1441 dzp_seg = z0_new - z0_old;
1442 }
1443 else {
1444 Zcell = Zcell + dirZ_sign;
1445 z0_new = Zcell;
1446 dzp_seg = z0_new - z0_old;
1447 }
1448
1449 // Compute the segment factor (equal to dt_seg/dt for nonzero dzp)
1450 const auto seg_factor = static_cast<amrex::Real>(dzp == 0. ? 1._rt : dzp_seg/dzp);
1451
1452 // Compute cell-based weights using the average segment position
1453 // Keep these double to avoid bug in single precision
1454 double sz_cell[depos_order] = {0.};
1455 double const z0_bar = (z0_new + z0_old)/2.0;
1456 k0_cell[ns] = compute_shape_factor_cell( sz_cell, z0_bar-0.5 );
1457
1458 // Set the segment number for the mass matrix component calc
1459 if constexpr (full_mass_matrices) {
1460 const auto k0_mid = static_cast<int>(z0_bar-shift);
1461 SegNum[ns] = 1 + k0_mid - k_mid;
1462 }
1463
1464 if constexpr (depos_order >= 3) { // higher-order correction to the cell-based weights
1465 const Compute_shape_factor_pair<depos_order-1> compute_shape_factors_cell;
1466 double sz_old_cell[depos_order] = {0.};
1467 double sz_new_cell[depos_order] = {0.};
1468 const int k0_cell_2 = compute_shape_factors_cell( sz_old_cell, sz_new_cell, z0_old-0.5, z0_new-0.5 );
1469 amrex::ignore_unused(k0_cell_2);
1470 for (int m=0; m<depos_order; m++) {
1471 sz_cell[m] = (4.0*sz_cell[m] + sz_old_cell[m] + sz_new_cell[m])/6.0;
1472 }
1473 }
1474
1475 // Compute node-based weights using the old and new segment positions
1476 // Keep these double to avoid bug in single precision
1477 double sz_old_node[depos_order+1] = {0.};
1478 double sz_new_node[depos_order+1] = {0.};
1479 k0_node[ns] = compute_shape_factors_node( sz_old_node, sz_new_node, z0_old, z0_new );
1480
1481 // deposit out-of-plane Jx, Jy, Sx, and Sy for this segment
1482 for (int k=0; k<=depos_order; k++) {
1483 const amrex::Real weight = 0.5_rt*(sz_old_node[k] + sz_new_node[k])*seg_factor;
1484 const int k_J = lo.x + k0_node[ns] + k;
1485 amrex::Gpu::Atomic::AddNoRet(&Jx_arr(k_J, 0, 0), wqx*weight);
1486 amrex::Gpu::Atomic::AddNoRet(&Jy_arr(k_J, 0, 0), wqy*weight);
1487 if constexpr (full_mass_matrices) { weight_node[ns][k] = weight; }
1488 else {
1489 amrex::Gpu::Atomic::AddNoRet(&Sxx_arr(k_J, 0, 0), fpxx*weight*weight);
1490 amrex::Gpu::Atomic::AddNoRet(&Syy_arr(k_J, 0, 0), fpyy*weight*weight);
1491 }
1492 }
1493
1494 // deposit Jz and Szz for this segment
1495 for (int k=0; k<=depos_order-1; k++) {
1496 const amrex::Real weight = sz_cell[k]*seg_factor;
1497 const int k_J = lo.x + k0_cell[ns] + k;
1498 amrex::Gpu::Atomic::AddNoRet(&Jz_arr(k_J, 0, 0), wqz*weight);
1499 if constexpr (full_mass_matrices) { weight_cell[ns][k] = weight; }
1500 else {
1501 amrex::Gpu::Atomic::AddNoRet(&Szz_arr(k_J, 0, 0), fpzz*weight*weight);
1502 }
1503 }
1504
1505 // update old segment values
1506 if (ns < num_segments-1) {
1507 z0_old = z0_new;
1508 }
1509
1510 }
1511
1512 if constexpr (full_mass_matrices) {
1513
1514 // Loop over segments and deposit full mass matrices
1515 for (int ns=0; ns<num_segments; ns++) {
1516
1517 // Deposit Sxx, Sxy, Sxz, Syx, Syy, and Syz for this segment
1518 for (int k=0; k<=depos_order; k++) {
1519
1520 const int k_J = lo.x + k0_node[ns] + k;
1521 const amrex::Real weight_J = weight_node[ns][k];
1522 for (int ms=0; ms<num_segments; ms++) {
1523 const int SegShift = max_crossings + SegNum[ms] - SegNum[ns];
1524 for (int kE=0; kE<=depos_order; kE++) {
1525 const amrex::Real weight_E = weight_node[ms][kE];
1526 const int comp_yy = depos_order - k + kE + SegShift;
1527 amrex::Gpu::Atomic::AddNoRet(&Sxx_arr(k_J, 0, comp_yy), fpxx*weight_J*weight_E);
1528 amrex::Gpu::Atomic::AddNoRet(&Syy_arr(k_J, 0, comp_yy), fpyy*weight_J*weight_E);
1529 amrex::Gpu::Atomic::AddNoRet(&Sxy_arr(k_J, 0, comp_yy), fpxy*weight_J*weight_E);
1530 amrex::Gpu::Atomic::AddNoRet(&Syx_arr(k_J, 0, comp_yy), fpyx*weight_J*weight_E);
1531 }
1532 for (int kE=0; kE<=depos_order-1; kE++) {
1533 const amrex::Real weight_E = weight_cell[ms][kE];
1534 const int comp_yz = depos_order - k + kE + SegShift;
1535 amrex::Gpu::Atomic::AddNoRet(&Sxz_arr(k_J, 0, comp_yz), fpxz*weight_J*weight_E);
1536 amrex::Gpu::Atomic::AddNoRet(&Syz_arr(k_J, 0, comp_yz), fpyz*weight_J*weight_E);
1537 }
1538 }
1539
1540 }
1541
1542 // Deposit Szx, Szy, and Szz for this segment
1543 for (int k=0; k<=depos_order-1; k++) {
1544
1545 const int k_J = lo.x + k0_cell[ns] + k;
1546 const amrex::Real weight_J = weight_cell[ns][k];
1547 for (int ms=0; ms<num_segments; ms++) {
1548 const int SegShift = max_crossings + SegNum[ms] - SegNum[ns];
1549 for (int kE=0; kE<=depos_order-1; kE++) {
1550 const amrex::Real weight_E = weight_cell[ms][kE];
1551 const int comp_zz = depos_order-1 - k + kE + SegShift;
1552 amrex::Gpu::Atomic::AddNoRet(&Szz_arr(k_J, 0, comp_zz), fpzz*weight_J*weight_E);
1553 }
1554 for (int kE=0; kE<=depos_order; kE++) {
1555 const amrex::Real weight_E = weight_node[ms][kE];
1556 const int comp_zy = depos_order-1 - k + kE + SegShift;
1557 amrex::Gpu::Atomic::AddNoRet(&Szx_arr(k_J, 0, comp_zy), fpzx*weight_J*weight_E);
1558 amrex::Gpu::Atomic::AddNoRet(&Szy_arr(k_J, 0, comp_zy), fpzy*weight_J*weight_E);
1559 }
1560 }
1561
1562 }
1563
1564 }
1565
1566 }
1567
1568#endif
1569}
1570
1598template <int depos_order, bool full_mass_matrices>
1599void doVillasenorJandSigmaDeposition ( [[maybe_unused]] const amrex::ParticleReal* xp_n_data,
1600 [[maybe_unused]] const amrex::ParticleReal* yp_n_data,
1601 [[maybe_unused]] const amrex::ParticleReal* zp_n_data,
1602 const GetParticlePosition<PIdx>& GetPosition,
1603 const amrex::ParticleReal* wp,
1604 const amrex::ParticleReal* uxp_n,
1605 const amrex::ParticleReal* uyp_n,
1606 const amrex::ParticleReal* uzp_n,
1607 const amrex::ParticleReal* uxp_nph,
1608 const amrex::ParticleReal* uyp_nph,
1609 const amrex::ParticleReal* uzp_nph,
1610 amrex::Array4<amrex::Real> const& Jx_arr,
1611 amrex::Array4<amrex::Real> const& Jy_arr,
1612 amrex::Array4<amrex::Real> const& Jz_arr,
1613 const int max_crossings,
1614 amrex::Array4<amrex::Real> const& Sxx_arr,
1615 amrex::Array4<amrex::Real> const& Sxy_arr,
1616 amrex::Array4<amrex::Real> const& Sxz_arr,
1617 amrex::Array4<amrex::Real> const& Syx_arr,
1618 amrex::Array4<amrex::Real> const& Syy_arr,
1619 amrex::Array4<amrex::Real> const& Syz_arr,
1620 amrex::Array4<amrex::Real> const& Szx_arr,
1621 amrex::Array4<amrex::Real> const& Szy_arr,
1622 amrex::Array4<amrex::Real> const& Szz_arr,
1623 const amrex::Array4<amrex::Real const>& Bx_arr,
1624 const amrex::Array4<amrex::Real const>& By_arr,
1625 const amrex::Array4<amrex::Real const>& Bz_arr,
1626 const amrex::IndexType Bx_type,
1627 const amrex::IndexType By_type,
1628 const amrex::IndexType Bz_type,
1629 const long np_to_deposit,
1630 const amrex::Real dt,
1631 const amrex::XDim3& dinv,
1632 const amrex::XDim3& xyzmin,
1633 const amrex::GpuArray<amrex::GpuArray<double,2>, AMREX_SPACEDIM> & domain_double,
1634 const amrex::GpuArray<amrex::GpuArray<bool,2>, AMREX_SPACEDIM> & do_cropping,
1635 const amrex::Dim3 lo,
1636 const amrex::Real qs,
1637 const amrex::Real ms )
1638{
1639 using namespace amrex::literals;
1640
1641 const amrex::Real invvol = dinv.x*dinv.y*dinv.z;
1642
1643 // Loop over particles and deposit into Jx_arr, Jy_arr and Jz_arr
1644 amrex::ParallelFor(
1645 np_to_deposit,
1646 [=] AMREX_GPU_DEVICE (long const ip) {
1647
1648 // Skip particles with zero weight.
1649 // This should only be the case for particles that will be suborbited.
1650 if (wp[ip] == 0.) { return; }
1651
1652 amrex::ParticleReal xp_nph, yp_nph, zp_nph;
1653 GetPosition(ip, xp_nph, yp_nph, zp_nph);
1654
1655 // Compute magnetic field on particle
1656 amrex::ParticleReal Bxp = 0.0;
1657 amrex::ParticleReal Byp = 0.0;
1658 amrex::ParticleReal Bzp = 0.0;
1659 const int depos_order_perp = 1;
1660 const int depos_order_para = 1;
1661 const int n_rz_azimuthal_modes = 0;
1662 doDirectGatherVectorField<depos_order_perp,depos_order_para>(
1663 xp_nph, yp_nph, zp_nph,
1664 Bxp, Byp, Bzp,
1665 Bx_arr, By_arr, Bz_arr,
1666 Bx_type, By_type, Bz_type,
1667 dinv, xyzmin, lo, n_rz_azimuthal_modes );
1668
1669 // Compute inverse Lorentz factor, the average of gamma at time levels n and n+1
1670 const amrex::ParticleReal gaminv = GetImplicitGammaInverse(uxp_n[ip], uyp_n[ip], uzp_n[ip],
1671 uxp_nph[ip], uyp_nph[ip], uzp_nph[ip]);
1672
1673 // Compute current density kernels to deposit
1674 const amrex::Real wq_invvol = qs*wp[ip]*invvol;
1675 const amrex::Real rhop = wq_invvol*gaminv;
1676
1677 // Set the Mass Matrices kernels
1678 amrex::ParticleReal fpxx, fpxy, fpxz;
1679 amrex::ParticleReal fpyx, fpyy, fpyz;
1680 amrex::ParticleReal fpzx, fpzy, fpzz;
1681 setMassMatricesKernels( qs, ms, dt, rhop,
1682 uxp_nph[ip], uyp_nph[ip], uzp_nph[ip],
1683 Bxp, Byp, Bzp,
1684 fpxx, fpxy, fpxz,
1685 fpyx, fpyy, fpyz,
1686 fpzx, fpzy, fpzz );
1687
1688 amrex::ParticleReal const xp_n = (xp_n_data ? xp_n_data[ip] : 0._prt);
1689 amrex::ParticleReal const yp_n = (yp_n_data ? yp_n_data[ip] : 0._prt);
1690 amrex::ParticleReal const zp_n = (zp_n_data ? zp_n_data[ip] : 0._prt);
1691
1692 // Compute position at time n + 1
1693 amrex::ParticleReal const xp_np1 = 2._prt*xp_nph - xp_n;
1694 amrex::ParticleReal const yp_np1 = 2._prt*yp_nph - yp_n;
1695 amrex::ParticleReal const zp_np1 = 2._prt*zp_nph - zp_n;
1696
1697 doVillasenorJandSigmaDepositionKernel<depos_order,full_mass_matrices>(
1698 xp_n, yp_n, zp_n,
1699 xp_np1, yp_np1, zp_np1,
1700 wq_invvol,
1701 uxp_nph[ip], uyp_nph[ip], uzp_nph[ip],
1702 gaminv,
1703 fpxx, fpxy, fpxz,
1704 fpyx, fpyy, fpyz,
1705 fpzx, fpzy, fpzz,
1706 Jx_arr, Jy_arr, Jz_arr,
1707 max_crossings,
1708 Sxx_arr, Sxy_arr, Sxz_arr,
1709 Syx_arr, Syy_arr, Syz_arr,
1710 Szx_arr, Szy_arr, Szz_arr,
1711 dt, dinv, xyzmin, domain_double, do_cropping, lo );
1712
1713 });
1714}
1715
1716#endif // WARPX_MASSMATRICESDEPOSITION_H_
AMReX_Arena.H
AMReX_Array4.H
AMREX_ALWAYS_ASSERT_WITH_MESSAGE
#define AMREX_ALWAYS_ASSERT_WITH_MESSAGE(EX, MSG)
AMReX_Dim3.H
AMREX_INLINE
#define AMREX_INLINE
AMREX_GPU_DEVICE
#define AMREX_GPU_DEVICE
AMREX_GPU_HOST_DEVICE
#define AMREX_GPU_HOST_DEVICE
offset
Array4< int const > offset
AMReX_REAL.H
doDirectGatherVectorField
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void doDirectGatherVectorField(const amrex::ParticleReal xp, const amrex::ParticleReal yp, const amrex::ParticleReal zp, amrex::ParticleReal &Fxp, amrex::ParticleReal &Fyp, amrex::ParticleReal &Fzp, amrex::Array4< amrex::Real const > const &Fx_arr, amrex::Array4< amrex::Real const > const &Fy_arr, amrex::Array4< amrex::Real const > const &Fz_arr, const amrex::IndexType Fx_type, const amrex::IndexType Fy_type, const amrex::IndexType Fz_type, const amrex::XDim3 &dinv, const amrex::XDim3 &xyzmin, const amrex::Dim3 &lo, const int n_rz_azimuthal_modes)
Gather vector field F for a single particle.
Definition FieldGather.H:36
setMassMatricesKernels
AMREX_GPU_HOST_DEVICE AMREX_INLINE void setMassMatricesKernels(const amrex::ParticleReal qs, const amrex::ParticleReal ms, const amrex::ParticleReal dt, const amrex::ParticleReal rhop, const amrex::ParticleReal uxp, const amrex::ParticleReal uyp, const amrex::ParticleReal uzp, const amrex::ParticleReal Bxp, const amrex::ParticleReal Byp, const amrex::ParticleReal Bzp, amrex::ParticleReal &fpxx, amrex::ParticleReal &fpxy, amrex::ParticleReal &fpxz, amrex::ParticleReal &fpyx, amrex::ParticleReal &fpyy, amrex::ParticleReal &fpyz, amrex::ParticleReal &fpzx, amrex::ParticleReal &fpzy, amrex::ParticleReal &fpzz)
Set the mass matrices kernels for thread thread_num.
Definition MassMatricesDeposition.H:42
doDirectJandSigmaDepositionKernel
AMREX_GPU_HOST_DEVICE AMREX_INLINE void doDirectJandSigmaDepositionKernel(const amrex::ParticleReal xp, const amrex::ParticleReal yp, const amrex::ParticleReal zp, const amrex::ParticleReal wqx, const amrex::ParticleReal wqy, const amrex::ParticleReal wqz, const amrex::ParticleReal fpxx, const amrex::ParticleReal fpxy, const amrex::ParticleReal fpxz, const amrex::ParticleReal fpyx, const amrex::ParticleReal fpyy, const amrex::ParticleReal fpyz, const amrex::ParticleReal fpzx, const amrex::ParticleReal fpzy, const amrex::ParticleReal fpzz, amrex::Array4< amrex::Real > const &jx_arr, amrex::Array4< amrex::Real > const &jy_arr, amrex::Array4< amrex::Real > const &jz_arr, amrex::Array4< amrex::Real > const &Sxx_arr, amrex::Array4< amrex::Real > const &Sxy_arr, amrex::Array4< amrex::Real > const &Sxz_arr, amrex::Array4< amrex::Real > const &Syx_arr, amrex::Array4< amrex::Real > const &Syy_arr, amrex::Array4< amrex::Real > const &Syz_arr, amrex::Array4< amrex::Real > const &Szx_arr, amrex::Array4< amrex::Real > const &Szy_arr, amrex::Array4< amrex::Real > const &Szz_arr, const amrex::IntVect &jx_type, const amrex::IntVect &jy_type, const amrex::IntVect &jz_type, const amrex::XDim3 &dinv, const amrex::XDim3 &xyzmin, const amrex::Dim3 lo)
Kernel for the direct deposition of J and S (mass matrices) for thread thread_num.
Definition MassMatricesDeposition.H:111
doVillasenorJandSigmaDeposition
void doVillasenorJandSigmaDeposition(const amrex::ParticleReal *xp_n_data, const amrex::ParticleReal *yp_n_data, const amrex::ParticleReal *zp_n_data, const GetParticlePosition< PIdx > &GetPosition, const amrex::ParticleReal *wp, const amrex::ParticleReal *uxp_n, const amrex::ParticleReal *uyp_n, const amrex::ParticleReal *uzp_n, const amrex::ParticleReal *uxp_nph, const amrex::ParticleReal *uyp_nph, const amrex::ParticleReal *uzp_nph, amrex::Array4< amrex::Real > const &Jx_arr, amrex::Array4< amrex::Real > const &Jy_arr, amrex::Array4< amrex::Real > const &Jz_arr, const int max_crossings, amrex::Array4< amrex::Real > const &Sxx_arr, amrex::Array4< amrex::Real > const &Sxy_arr, amrex::Array4< amrex::Real > const &Sxz_arr, amrex::Array4< amrex::Real > const &Syx_arr, amrex::Array4< amrex::Real > const &Syy_arr, amrex::Array4< amrex::Real > const &Syz_arr, amrex::Array4< amrex::Real > const &Szx_arr, amrex::Array4< amrex::Real > const &Szy_arr, amrex::Array4< amrex::Real > const &Szz_arr, const amrex::Array4< amrex::Real const > &Bx_arr, const amrex::Array4< amrex::Real const > &By_arr, const amrex::Array4< amrex::Real const > &Bz_arr, const amrex::IndexType Bx_type, const amrex::IndexType By_type, const amrex::IndexType Bz_type, const long np_to_deposit, const amrex::Real dt, const amrex::XDim3 &dinv, const amrex::XDim3 &xyzmin, const amrex::GpuArray< amrex::GpuArray< double, 2 >, 3 > &domain_double, const amrex::GpuArray< amrex::GpuArray< bool, 2 >, 3 > &do_cropping, const amrex::Dim3 lo, const amrex::Real qs, const amrex::Real ms)
Villasenor and Buneman deposition of J and mass matrices for thread thread_num.
Definition MassMatricesDeposition.H:1599
doDirectJandSigmaDeposition
void doDirectJandSigmaDeposition(const GetParticlePosition< PIdx > &GetPosition, const amrex::ParticleReal *wp, const amrex::ParticleReal *uxp_n, const amrex::ParticleReal *uyp_n, const amrex::ParticleReal *uzp_n, const amrex::ParticleReal *uxp_nph, const amrex::ParticleReal *uyp_nph, const amrex::ParticleReal *uzp_nph, amrex::FArrayBox &jx_fab, amrex::FArrayBox &jy_fab, amrex::FArrayBox &jz_fab, amrex::Array4< amrex::Real > const &Sxx_arr, amrex::Array4< amrex::Real > const &Sxy_arr, amrex::Array4< amrex::Real > const &Sxz_arr, amrex::Array4< amrex::Real > const &Syx_arr, amrex::Array4< amrex::Real > const &Syy_arr, amrex::Array4< amrex::Real > const &Syz_arr, amrex::Array4< amrex::Real > const &Szx_arr, amrex::Array4< amrex::Real > const &Szy_arr, amrex::Array4< amrex::Real > const &Szz_arr, const amrex::Array4< amrex::Real const > &Bx_arr, const amrex::Array4< amrex::Real const > &By_arr, const amrex::Array4< amrex::Real const > &Bz_arr, const amrex::IndexType Bx_type, const amrex::IndexType By_type, const amrex::IndexType Bz_type, const long np_to_deposit, const amrex::Real dt, const amrex::XDim3 &dinv, const amrex::XDim3 &xyzmin, const amrex::Dim3 lo, const amrex::Real qs, const amrex::Real ms)
direct deposition of J and mass matrices for thread thread_num
Definition MassMatricesDeposition.H:522
doVillasenorJandSigmaDepositionKernel
AMREX_GPU_HOST_DEVICE AMREX_INLINE void doVillasenorJandSigmaDepositionKernel(const amrex::ParticleReal xp_old, const amrex::ParticleReal yp_old, const amrex::ParticleReal zp_old, const amrex::ParticleReal xp_new, const amrex::ParticleReal yp_new, const amrex::ParticleReal zp_new, const amrex::ParticleReal wq_invvol, const amrex::ParticleReal uxp_mid, const amrex::ParticleReal uyp_mid, const amrex::ParticleReal uzp_mid, const amrex::ParticleReal gaminv, const amrex::ParticleReal fpxx, const amrex::ParticleReal fpxy, const amrex::ParticleReal fpxz, const amrex::ParticleReal fpyx, const amrex::ParticleReal fpyy, const amrex::ParticleReal fpyz, const amrex::ParticleReal fpzx, const amrex::ParticleReal fpzy, const amrex::ParticleReal fpzz, amrex::Array4< amrex::Real > const &Jx_arr, amrex::Array4< amrex::Real > const &Jy_arr, amrex::Array4< amrex::Real > const &Jz_arr, int max_crossings, amrex::Array4< amrex::Real > const &Sxx_arr, amrex::Array4< amrex::Real > const &Sxy_arr, amrex::Array4< amrex::Real > const &Sxz_arr, amrex::Array4< amrex::Real > const &Syx_arr, amrex::Array4< amrex::Real > const &Syy_arr, amrex::Array4< amrex::Real > const &Syz_arr, amrex::Array4< amrex::Real > const &Szx_arr, amrex::Array4< amrex::Real > const &Szy_arr, amrex::Array4< amrex::Real > const &Szz_arr, const amrex::Real dt, const amrex::XDim3 &dinv, const amrex::XDim3 &xyzmin, const amrex::GpuArray< amrex::GpuArray< double, 2 >, 3 > &domain_double, const amrex::GpuArray< amrex::GpuArray< bool, 2 >, 3 > &do_cropping, const amrex::Dim3 lo)
Kernel for the Villasenor deposition of J and S (mass matrices) for thread thread_num.
Definition MassMatricesDeposition.H:653
RigidAdvanceMode::vz
@ vz
Definition RigidInjectedParticleContainer.H:27
PhysicalSpecies::alpha
@ alpha
Definition SpeciesPhysicalProperties.H:18
GetImplicitGammaInverse
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::ParticleReal GetImplicitGammaInverse(const amrex::ParticleReal uxp_n, const amrex::ParticleReal uyp_n, const amrex::ParticleReal uzp_n, const amrex::ParticleReal uxp_nph, const amrex::ParticleReal uyp_nph, const amrex::ParticleReal uzp_nph) noexcept
Compute the inverse Lorentz factor for the position update in the implicit methods,...
Definition UpdatePosition.H:77
TextMsg.H
WarpX_Complex.H
amrex::BaseFab::box
const Box & box() const noexcept
amrex::BaseFab::array
Array4< T const > array() const noexcept
amrex::BoxND::type
__host__ __device__ IntVectND< dim > type() const noexcept
amrex::IntVectND< 3 >::TheZeroVector
__host__ static __device__ constexpr IntVectND< dim > TheZeroVector() noexcept
amrex::Real
amrex_real Real
amrex::ParticleReal
amrex_particle_real ParticleReal
ParticleUtils::crop_at_boundary
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void crop_at_boundary(double &x1, double xmin, double xmax, amrex::GpuArray< bool, 2 > const &do_cropping)
Definition ParticleUtils.H:256
ablastr::constant::SI::c
constexpr auto c
vacuum speed of light [m/s]
Definition constant.H:153
amrex::Gpu::Atomic::AddNoRet
__host__ __device__ AMREX_FORCE_INLINE void AddNoRet(T *sum, T value) noexcept
amrex::ignore_unused
__host__ __device__ void ignore_unused(const Ts &...)
amrex::ParallelFor
std::enable_if_t< std::is_integral_v< T > > ParallelFor(TypeList< CTOs... > ctos, std::array< int, sizeof...(CTOs)> const &runtime_options, T N, F &&f)
amrex::shift
__host__ __device__ BoxND< dim > shift(const BoxND< dim > &b, int dir, int nzones) noexcept
amrex::IndexType
IndexTypeND< 3 > IndexType
amrex::Array4
ArrayND< T, 4, true > Array4
amrex::IntVect
IntVectND< 3 > IntVect
Compute_shape_factor_pair
Definition ShapeFactors.H:168
Compute_shape_factor
Definition ShapeFactors.H:29
GetParticlePosition
Functor that can be used to extract the positions of the macroparticles inside a ParallelFor kernel.
Definition GetAndSetPosition.H:75
Generated by doxygen 1.13.2