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