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 upx,
47 const amrex::ParticleReal upy,
48 const amrex::ParticleReal upz,
49 const amrex::ParticleReal Bpx,
50 const amrex::ParticleReal Bpy,
51 const amrex::ParticleReal Bpz,
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 auto inv_c2 = PhysConst::inv_c2_v<amrex::ParticleReal>;
65
66 // Convert B on particle to normalized cyclotron units with dt/2.0
67 const amrex::ParticleReal gamma_bar = std::sqrt(1._prt + (upx*upx + upy*upy + upz*upz)*inv_c2);
68 const amrex::ParticleReal alpha = qs/ms*0.5_prt*dt/gamma_bar;
69 const amrex::ParticleReal bpx = alpha*Bpx;
70 const amrex::ParticleReal bpy = alpha*Bpy;
71 const amrex::ParticleReal bpz = alpha*Bpz;
72
73 // Compute Mass Matrix kernels (non-relativistic for now)
74 const amrex::ParticleReal bpsq = bpx*bpx + bpy*bpy + bpz*bpz;
75 const amrex::ParticleReal arogp = alpha*rhop/(1.0_prt + bpsq);
76
77 fpxx = arogp*(bpx*bpx + 1.0_rt);
78 fpxy = arogp*(bpx*bpy + bpz);
79 fpxz = arogp*(bpx*bpz - bpy);
80
81 fpyx = arogp*(bpy*bpx - bpz);
82 fpyy = arogp*(bpy*bpy + 1.0_rt);
83 fpyz = arogp*(bpy*bpz + bpx);
84
85 fpzx = arogp*(bpz*bpx + bpy);
86 fpzy = arogp*(bpz*bpy - bpx);
87 fpzz = arogp*(bpz*bpz + 1.0_rt);
88
89}
90
110template <int depos_order, bool full_mass_matrices, bool deposit_J>
111AMREX_GPU_HOST_DEVICE AMREX_INLINE
112void doDirectJandSigmaDepositionKernel ( [[maybe_unused]] const amrex::ParticleReal xp,
113 [[maybe_unused]] const amrex::ParticleReal yp,
114 [[maybe_unused]] const amrex::ParticleReal zp,
115 const amrex::ParticleReal wqx,
116 const amrex::ParticleReal wqy,
117 const amrex::ParticleReal wqz,
118 const amrex::ParticleReal fpxx,
119 [[maybe_unused]] const amrex::ParticleReal fpxy,
120 [[maybe_unused]] const amrex::ParticleReal fpxz,
121 [[maybe_unused]] const amrex::ParticleReal fpyx,
122 const amrex::ParticleReal fpyy,
123 [[maybe_unused]] const amrex::ParticleReal fpyz,
124 [[maybe_unused]] const amrex::ParticleReal fpzx,
125 [[maybe_unused]] const amrex::ParticleReal fpzy,
126 const amrex::ParticleReal fpzz,
127 amrex::Array4<amrex::Real> const& jx_arr,
128 amrex::Array4<amrex::Real> const& jy_arr,
129 amrex::Array4<amrex::Real> const& jz_arr,
130 amrex::Array4<amrex::Real> const& Sxx_arr,
131 [[maybe_unused]] amrex::Array4<amrex::Real> const& Sxy_arr,
132 [[maybe_unused]] amrex::Array4<amrex::Real> const& Sxz_arr,
133 [[maybe_unused]] amrex::Array4<amrex::Real> const& Syx_arr,
134 amrex::Array4<amrex::Real> const& Syy_arr,
135 [[maybe_unused]] amrex::Array4<amrex::Real> const& Syz_arr,
136 [[maybe_unused]] amrex::Array4<amrex::Real> const& Szx_arr,
137 [[maybe_unused]] amrex::Array4<amrex::Real> const& Szy_arr,
138 amrex::Array4<amrex::Real> const& Szz_arr,
139 const amrex::IntVect& jx_type,
140 const amrex::IntVect& jy_type,
141 const amrex::IntVect& jz_type,
142 const amrex::XDim3& dinv,
143 const amrex::XDim3& xyzmin,
144 const amrex::Dim3 lo )
145{
146 using namespace amrex::literals;
147
148 constexpr int NODE = amrex::IndexType::NODE;
149 constexpr int CELL = amrex::IndexType::CELL;
150
151 // MassMatrices index shift parameter
152 amrex::IntVect shift = amrex::IntVect::TheZeroVector();
153
154 // --- Compute shape factors
155 Compute_shape_factor< depos_order > const compute_shape_factor;
156#if !defined(WARPX_DIM_1D_Z)
157 // x direction
158 // Get particle position after 1/2 push back in position
159 // Keep these double to avoid bug in single precision
160 const double xmid = (xp - xyzmin.x)*dinv.x;
161
162 // j_j[xyz] leftmost grid point in x that the particle touches for the centering of each current
163 // sx_j[xyz] shape factor along x for the centering of each current
164 // There are only two possible centerings, node or cell centered, so at most only two shape factor
165 // arrays will be needed.
166 // Keep these double to avoid bug in single precision
167 double sx_node[depos_order + 1] = {0.};
168 double sx_cell[depos_order + 1] = {0.};
169 int j_node = 0;
170 int j_cell = 0;
171 if (jx_type[0] == NODE || jy_type[0] == NODE || jz_type[0] == NODE) {
172 j_node = compute_shape_factor(sx_node, xmid);
173 }
174 if (jx_type[0] == CELL || jy_type[0] == CELL || jz_type[0] == CELL) {
175 j_cell = compute_shape_factor(sx_cell, xmid - 0.5);
176 }
177
178 // Set the index shift parameter
179 if (j_node==j_cell) { shift[0] = 1; }
180
181 amrex::Real sx_jx[depos_order + 1] = {0._rt};
182 amrex::Real sx_jy[depos_order + 1] = {0._rt};
183 amrex::Real sx_jz[depos_order + 1] = {0._rt};
184 for (int ix=0; ix<=depos_order; ix++)
185 {
186 sx_jx[ix] = ((jx_type[0] == NODE) ? amrex::Real(sx_node[ix]) : amrex::Real(sx_cell[ix]));
187 sx_jy[ix] = ((jy_type[0] == NODE) ? amrex::Real(sx_node[ix]) : amrex::Real(sx_cell[ix]));
188 sx_jz[ix] = ((jz_type[0] == NODE) ? amrex::Real(sx_node[ix]) : amrex::Real(sx_cell[ix]));
189 }
190
191 int const j_jx = ((jx_type[0] == NODE) ? j_node : j_cell);
192 int const j_jy = ((jy_type[0] == NODE) ? j_node : j_cell);
193 int const j_jz = ((jz_type[0] == NODE) ? j_node : j_cell);
194#endif
195
196#if defined(WARPX_DIM_3D)
197 // y direction
198 // Keep these double to avoid bug in single precision
199 const double ymid = (yp - xyzmin.y)*dinv.y;
200 double sy_node[depos_order + 1] = {0.};
201 double sy_cell[depos_order + 1] = {0.};
202 int k_node = 0;
203 int k_cell = 0;
204 if (jx_type[1] == NODE || jy_type[1] == NODE || jz_type[1] == NODE) {
205 k_node = compute_shape_factor(sy_node, ymid);
206 }
207 if (jx_type[1] == CELL || jy_type[1] == CELL || jz_type[1] == CELL) {
208 k_cell = compute_shape_factor(sy_cell, ymid - 0.5);
209 }
210
211 // Set the index shift parameter
212 if (k_node==k_cell) { shift[1] = 1; }
213
214 amrex::Real sy_jx[depos_order + 1] = {0._rt};
215 amrex::Real sy_jy[depos_order + 1] = {0._rt};
216 amrex::Real sy_jz[depos_order + 1] = {0._rt};
217 for (int iy=0; iy<=depos_order; iy++)
218 {
219 sy_jx[iy] = ((jx_type[1] == NODE) ? amrex::Real(sy_node[iy]) : amrex::Real(sy_cell[iy]));
220 sy_jy[iy] = ((jy_type[1] == NODE) ? amrex::Real(sy_node[iy]) : amrex::Real(sy_cell[iy]));
221 sy_jz[iy] = ((jz_type[1] == NODE) ? amrex::Real(sy_node[iy]) : amrex::Real(sy_cell[iy]));
222 }
223 int const k_jx = ((jx_type[1] == NODE) ? k_node : k_cell);
224 int const k_jy = ((jy_type[1] == NODE) ? k_node : k_cell);
225 int const k_jz = ((jz_type[1] == NODE) ? k_node : k_cell);
226#endif
227
228#if !defined(WARPX_DIM_RCYLINDER) && !defined(WARPX_DIM_RSPHERE)
229 // z direction
230 // Keep these double to avoid bug in single precision
231 constexpr int zdir = WARPX_ZINDEX;
232 const double zmid = (zp - xyzmin.z)*dinv.z;
233 double sz_node[depos_order + 1] = {0.};
234 double sz_cell[depos_order + 1] = {0.};
235 int l_node = 0;
236 int l_cell = 0;
237 if (jx_type[zdir] == NODE || jy_type[zdir] == NODE || jz_type[zdir] == NODE) {
238 l_node = compute_shape_factor(sz_node, zmid);
239 }
240 if (jx_type[zdir] == CELL || jy_type[zdir] == CELL || jz_type[zdir] == CELL) {
241 l_cell = compute_shape_factor(sz_cell, zmid - 0.5);
242 }
243 amrex::Real sz_jx[depos_order + 1] = {0._rt};
244 amrex::Real sz_jy[depos_order + 1] = {0._rt};
245 amrex::Real sz_jz[depos_order + 1] = {0._rt};
246 for (int iz=0; iz<=depos_order; iz++)
247 {
248 sz_jx[iz] = ((jx_type[zdir] == NODE) ? amrex::Real(sz_node[iz]) : amrex::Real(sz_cell[iz]));
249 sz_jy[iz] = ((jy_type[zdir] == NODE) ? amrex::Real(sz_node[iz]) : amrex::Real(sz_cell[iz]));
250 sz_jz[iz] = ((jz_type[zdir] == NODE) ? amrex::Real(sz_node[iz]) : amrex::Real(sz_cell[iz]));
251 }
252 int const l_jx = ((jx_type[zdir] == NODE) ? l_node : l_cell);
253 int const l_jy = ((jy_type[zdir] == NODE) ? l_node : l_cell);
254 int const l_jz = ((jz_type[zdir] == NODE) ? l_node : l_cell);
255
256 // Set the index shift parameter
257 if (l_node==l_cell) { shift[zdir] = 1; }
258
259#endif
260
261 // Compute index offset needed when x and y comps have different location on grid
262 amrex::IntVect offset_xy, offset_xz, offset_yz;
263 for (int dir=0; dir<AMREX_SPACEDIM; dir++) {
264 offset_xy[dir] = (jx_type[dir] + jy_type[dir]) % 2;
265 offset_xz[dir] = (jx_type[dir] + jz_type[dir]) % 2;
266 offset_yz[dir] = (jy_type[dir] + jz_type[dir]) % 2;
267 }
268
269 // Deposit J and mass matrices
270#if defined(WARPX_DIM_1D_Z)
271 for (int iz=0; iz<=depos_order; iz++){
272 if constexpr (deposit_J) {
273 amrex::Gpu::Atomic::AddNoRet(&jx_arr(lo.x+l_jx+iz, 0, 0, 0), sz_jx[iz]*wqx);
274 amrex::Gpu::Atomic::AddNoRet(&jy_arr(lo.x+l_jy+iz, 0, 0, 0), sz_jy[iz]*wqy);
275 amrex::Gpu::Atomic::AddNoRet(&jz_arr(lo.x+l_jz+iz, 0, 0, 0), sz_jz[iz]*wqz);
276 }
277
278 // Deposit mass matrices
279 if constexpr (full_mass_matrices) {
280 for (int aa=0; aa<=depos_order; aa++){
281
282 const int col_base = depos_order + aa - iz;
283 int Nc = 0;
284
285 // Reduced deposit for diagonal mass matrices
286 if (aa <= iz) {
287 Nc = col_base;
288 amrex::Gpu::Atomic::AddNoRet(
289 &Sxx_arr(lo.x+l_jx+iz, 0, 0, Nc),
290 sz_jx[iz]*sz_jx[aa]*fpxx);
291 amrex::Gpu::Atomic::AddNoRet(
292 &Syy_arr(lo.x+l_jy+iz, 0, 0, Nc),
293 sz_jy[iz]*sz_jy[aa]*fpyy);
294 amrex::Gpu::Atomic::AddNoRet(
295 &Szz_arr(lo.x+l_jz+iz, 0, 0, Nc),
296 sz_jz[iz]*sz_jz[aa]*fpzz);
297 }
298
299 // Deposit off-diagonal mass matrices for X-current
300 Nc = col_base + shift[0]*offset_xy[0];
301 amrex::Gpu::Atomic::AddNoRet(
302 &Sxy_arr(lo.x+l_jx+iz, 0, 0, Nc),
303 sz_jx[iz]*sz_jy[aa]*fpxy);
304 Nc = col_base + shift[0]*offset_xz[0];
305 amrex::Gpu::Atomic::AddNoRet(
306 &Sxz_arr(lo.x+l_jx+iz, 0, 0, Nc),
307 sz_jx[iz]*sz_jz[aa]*fpxz);
308
309 // Deposit off-diagonal mass matrices for Y-current
310 Nc = col_base + shift[0]*offset_xy[0];
311 amrex::Gpu::Atomic::AddNoRet(
312 &Syx_arr(lo.x+l_jy+iz, 0, 0, Nc),
313 sz_jy[iz]*sz_jx[aa]*fpyx);
314 Nc = col_base + shift[0]*offset_yz[0];
315 amrex::Gpu::Atomic::AddNoRet(
316 &Syz_arr(lo.x+l_jy+iz, 0, 0, Nc),
317 sz_jy[iz]*sz_jz[aa]*fpyz);
318
319 // Deposit off-diagonal mass matrices for Z-current
320 Nc = col_base + 1 - shift[0]*offset_xz[0];
321 amrex::Gpu::Atomic::AddNoRet(
322 &Szx_arr(lo.x+l_jz+iz, 0, 0, Nc),
323 sz_jz[iz]*sz_jx[aa]*fpzx);
324 Nc = col_base + 1 - shift[0]*offset_yz[0];
325 amrex::Gpu::Atomic::AddNoRet(
326 &Szy_arr(lo.x+l_jz+iz, 0, 0, Nc),
327 sz_jz[iz]*sz_jy[aa]*fpzy);
328 }
329 }
330 else { // Deposit mass matrices for diagonal PC only
331 amrex::Gpu::Atomic::AddNoRet(
332 &Sxx_arr(lo.x+l_jx+iz, 0, 0, 0),
333 sz_jx[iz]*fpxx);
334 amrex::Gpu::Atomic::AddNoRet(
335 &Syy_arr(lo.x+l_jy+iz, 0, 0, 0),
336 sz_jy[iz]*fpyy);
337 amrex::Gpu::Atomic::AddNoRet(
338 &Szz_arr(lo.x+l_jz+iz, 0, 0, 0),
339 sz_jz[iz]*fpzz);
340 }
341 }
342#elif defined(WARPX_DIM_RCYLINDER) || defined(WARPX_DIM_RSPHERE)
343 for (int ix=0; ix<=depos_order; ix++){
344 if constexpr (deposit_J) {
345 amrex::Gpu::Atomic::AddNoRet(&jx_arr(lo.x+j_jx+ix, 0, 0, 0), sx_jx[ix]*wqx);
346 amrex::Gpu::Atomic::AddNoRet(&jy_arr(lo.x+j_jy+ix, 0, 0, 0), sx_jy[ix]*wqy);
347 amrex::Gpu::Atomic::AddNoRet(&jz_arr(lo.x+j_jz+ix, 0, 0, 0), sx_jz[ix]*wqz);
348 }
349 amrex::Gpu::Atomic::AddNoRet(&Sxx_arr(lo.x+j_jx+ix, 0, 0, 0), sx_jx[ix]*sx_jx[ix]*fpxx);
350 amrex::Gpu::Atomic::AddNoRet(&Syy_arr(lo.x+j_jy+ix, 0, 0, 0), sx_jy[ix]*sx_jy[ix]*fpyy);
351 amrex::Gpu::Atomic::AddNoRet(&Szz_arr(lo.x+j_jz+ix, 0, 0, 0), sx_jz[ix]*sx_jz[ix]*fpzz);
352 }
353#elif defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ)
354 const int base_offset = 1 + 2*depos_order;
355
356 for (int iz=0; iz<=depos_order; iz++){
357
358 for (int ix=0; ix<=depos_order; ix++){
359
360 const amrex::Real weight_Jx = sx_jx[ix]*sz_jx[iz];
361 const amrex::Real weight_Jy = sx_jy[ix]*sz_jy[iz];
362 const amrex::Real weight_Jz = sx_jz[ix]*sz_jz[iz];
363
364 if constexpr (deposit_J) {
365 amrex::Gpu::Atomic::AddNoRet(&jx_arr(lo.x+j_jx+ix, lo.y+l_jx+iz, 0, 0), weight_Jx*wqx);
366 amrex::Gpu::Atomic::AddNoRet(&jy_arr(lo.x+j_jy+ix, lo.y+l_jy+iz, 0, 0), weight_Jy*wqy);
367 amrex::Gpu::Atomic::AddNoRet(&jz_arr(lo.x+j_jz+ix, lo.y+l_jz+iz, 0, 0), weight_Jz*wqz);
368 }
369
370 // Deposit mass matrices
371 if constexpr (full_mass_matrices) {
372
373 const int Ncomp0 = 1 + 2*depos_order;
374
375 for (int bb=0; bb<=depos_order; bb++){
376
377 const int row_base = depos_order + bb - iz;
378
379 for (int aa=0; aa<=depos_order; aa++){
380 const int col_base = depos_order + aa - ix;
381 const amrex::Real weight_Ex = sx_jx[aa]*sz_jx[bb];
382 const amrex::Real weight_Ey = sx_jy[aa]*sz_jy[bb];
383 const amrex::Real weight_Ez = sx_jz[aa]*sz_jz[bb];
384
385 int offset;
386 int Nc;
387
388 // Reduced deposit for diagonal mass matrices
389 if (col_base <= Ncomp0 - row_base) {
390 offset = base_offset;
391 Nc = col_base + row_base*offset;
392 amrex::Gpu::Atomic::AddNoRet(
393 &Sxx_arr(lo.x+j_jx+ix, lo.y+l_jx+iz, 0, Nc),
394 weight_Jx*weight_Ex*fpxx);
395 amrex::Gpu::Atomic::AddNoRet(
396 &Syy_arr(lo.x+j_jy+ix, lo.y+l_jy+iz, 0, Nc),
397 weight_Jy*weight_Ey*fpyy);
398 amrex::Gpu::Atomic::AddNoRet(
399 &Szz_arr(lo.x+j_jz+ix, lo.y+l_jz+iz, 0, Nc),
400 weight_Jz*weight_Ez*fpzz);
401 }
402
403 // Deposit off-diagonal mass matrices for X-current
404 offset = base_offset + offset_xy[0];
405 Nc = col_base + 1 - shift[0]*offset_xy[0]
406 + (row_base + shift[1]*offset_xy[1])*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 = col_base + 1 - shift[0]*offset_xz[0]
412 + (row_base + shift[1]*offset_xz[1])*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 off-diagonal mass matrices for Y-current
418 offset = base_offset + offset_xy[0];
419 Nc = col_base + shift[0]*offset_xy[0]
420 + (row_base + shift[1]*offset_xy[1])*offset;
421 amrex::Gpu::Atomic::AddNoRet(
422 &Syx_arr(lo.x+j_jy+ix, lo.y+l_jy+iz, 0, Nc),
423 weight_Jy*weight_Ex*fpyx);
424 offset = base_offset + offset_yz[0];
425 Nc = col_base + shift[0]*offset_yz[0]
426 + (row_base + shift[1]*offset_yz[1])*offset;
427 amrex::Gpu::Atomic::AddNoRet(
428 &Syz_arr(lo.x+j_jy+ix, lo.y+l_jy+iz, 0, Nc),
429 weight_Jy*weight_Ez*fpyz);
430
431 // Deposit off-diagonal mass matrices for Z-current
432 offset = base_offset + offset_xz[0];
433 Nc = col_base + shift[0]*offset_xz[0]
434 + (row_base + 1 - shift[1]*offset_xz[1])*offset;
435 amrex::Gpu::Atomic::AddNoRet(
436 &Szx_arr(lo.x+j_jz+ix, lo.y+l_jz+iz, 0, Nc),
437 weight_Jz*weight_Ex*fpzx);
438 offset = base_offset + offset_yz[0];
439 Nc = col_base + shift[0]*offset_yz[0]
440 + (row_base + 1 - shift[1]*offset_yz[1])*offset;
441 amrex::Gpu::Atomic::AddNoRet(
442 &Szy_arr(lo.x+j_jz+ix, lo.y+l_jz+iz, 0, Nc),
443 weight_Jz*weight_Ey*fpzy);
444 }
445
446 }
447
448 }
449 else { // Deposit mass matrices for diagonal PC only
450 amrex::Gpu::Atomic::AddNoRet(
451 &Syy_arr(lo.x+j_jy+ix, lo.y+l_jy+iz, 0, 0),
452 weight_Jy*fpyy);
453 amrex::Gpu::Atomic::AddNoRet(
454 &Sxx_arr(lo.x+j_jx+ix, lo.y+l_jx+iz, 0, 0),
455 weight_Jx*fpxx);
456 amrex::Gpu::Atomic::AddNoRet(
457 &Szz_arr(lo.x+j_jz+ix, lo.y+l_jz+iz, 0, 0),
458 weight_Jz*fpzz);
459 }
460 }
461 }
462#elif defined(WARPX_DIM_3D)
463 for (int iz=0; iz<=depos_order; iz++){
464 for (int iy=0; iy<=depos_order; iy++){
465 for (int ix=0; ix<=depos_order; ix++){
466 const amrex::Real weight_Jx = sx_jx[ix]*sy_jx[iy]*sz_jx[iz];
467 const amrex::Real weight_Jy = sx_jy[ix]*sy_jy[iy]*sz_jy[iz];
468 const amrex::Real weight_Jz = sx_jz[ix]*sy_jz[iy]*sz_jz[iz];
469
470 if constexpr (deposit_J) {
471 amrex::Gpu::Atomic::AddNoRet(&jx_arr(lo.x+j_jx+ix, lo.y+k_jx+iy, lo.z+l_jx+iz), weight_Jx*wqx);
472 amrex::Gpu::Atomic::AddNoRet(&jy_arr(lo.x+j_jy+ix, lo.y+k_jy+iy, lo.z+l_jy+iz), weight_Jy*wqy);
473 amrex::Gpu::Atomic::AddNoRet(&jz_arr(lo.x+j_jz+ix, lo.y+k_jz+iy, lo.z+l_jz+iz), weight_Jz*wqz);
474 }
475
476 // Deposit mass matrices
477 if constexpr (full_mass_matrices) {
478 // Should not be here. Full mass matrices not yet implemented in 3D
479 }
480 else { // Deposit mass matrices for diagonal PC only
481 amrex::Gpu::Atomic::AddNoRet(
482 &Sxx_arr(lo.x+j_jx+ix, lo.y+k_jx+iy, lo.z+l_jx+iz, 0),
483 weight_Jx*fpxx);
484 amrex::Gpu::Atomic::AddNoRet(
485 &Syy_arr(lo.x+j_jy+ix, lo.y+k_jy+iy, lo.z+l_jy+iz, 0),
486 weight_Jy*fpyy);
487 amrex::Gpu::Atomic::AddNoRet(
488 &Szz_arr(lo.x+j_jz+ix, lo.y+k_jz+iy, lo.z+l_jz+iz, 0),
489 weight_Jz*fpzz);
490 }
491 }
492 }
493 }
494#endif
495}
496
521template <int depos_order, bool full_mass_matrices>
522void doDirectSigmaDeposition (const GetParticlePosition<PIdx>& GetPosition,
523 [[maybe_unused]] const int* nsuborbits,
524 const amrex::ParticleReal* wp,
525 const amrex::ParticleReal* uxp_n,
526 const amrex::ParticleReal* uyp_n,
527 const amrex::ParticleReal* uzp_n,
528 const amrex::ParticleReal* uxp_nph,
529 const amrex::ParticleReal* uyp_nph,
530 const amrex::ParticleReal* uzp_nph,
531 amrex::Array4<amrex::Real> const& Sxx_arr,
532 amrex::Array4<amrex::Real> const& Sxy_arr,
533 amrex::Array4<amrex::Real> const& Sxz_arr,
534 amrex::Array4<amrex::Real> const& Syx_arr,
535 amrex::Array4<amrex::Real> const& Syy_arr,
536 amrex::Array4<amrex::Real> const& Syz_arr,
537 amrex::Array4<amrex::Real> const& Szx_arr,
538 amrex::Array4<amrex::Real> const& Szy_arr,
539 amrex::Array4<amrex::Real> const& Szz_arr,
540 const amrex::IntVect& jx_type,
541 const amrex::IntVect& jy_type,
542 const amrex::IntVect& jz_type,
543 const amrex::Array4<amrex::Real const>& Bx_arr,
544 const amrex::Array4<amrex::Real const>& By_arr,
545 const amrex::Array4<amrex::Real const>& Bz_arr,
546 const amrex::IndexType Bx_type,
547 const amrex::IndexType By_type,
548 const amrex::IndexType Bz_type,
549 const long np_to_deposit,
550 const amrex::Real dt,
551 const amrex::XDim3& dinv,
552 const amrex::XDim3& xyzmin,
553 const amrex::Dim3 lo,
554 const amrex::Real qs,
555 const amrex::Real ms)
556{
557 using namespace amrex::literals;
558
559 const amrex::Real invvol = dinv.x*dinv.y*dinv.z;
560
561 // Loop over particles and deposit mass matrices
562 amrex::ParallelFor(
563 np_to_deposit,
564 [=] AMREX_GPU_DEVICE (long ip) {
565
566 // Skip mass matrix deposition for particles with suborbits.
567 if (nsuborbits && nsuborbits[ip] > 1) { return; }
568
569 amrex::ParticleReal xp_nph, yp_nph, zp_nph;
570 GetPosition(ip, xp_nph, yp_nph, zp_nph);
571
572 // Compute magnetic field on particle
573 amrex::ParticleReal Bxp = 0.0;
574 amrex::ParticleReal Byp = 0.0;
575 amrex::ParticleReal Bzp = 0.0;
576 doDirectGatherVectorField</*depos_order_perp=*/depos_order,/*depos_order_para=*/depos_order>(
577 xp_nph, yp_nph, zp_nph,
578 Bxp, Byp, Bzp,
579 Bx_arr, By_arr, Bz_arr,
580 Bx_type, By_type, Bz_type,
581 dinv, xyzmin, lo, /*n_rz_azimuthal_modes=*/0 );
582
583 // Compute inverse Lorentz factor, the average of gamma at time levels n and n+1
584 const amrex::ParticleReal gaminv = GetImplicitGammaInverse(uxp_n[ip], uyp_n[ip], uzp_n[ip],
585 uxp_nph[ip], uyp_nph[ip], uzp_nph[ip]);
586
587 // Compute current density kernels to deposit
588 const amrex::Real rhop = qs*wp[ip]*invvol*gaminv;
589 const amrex::Real wqx = rhop*uxp_nph[ip];
590 const amrex::Real wqy = rhop*uyp_nph[ip];
591 const amrex::Real wqz = rhop*uzp_nph[ip];
592
593 // Set the Mass Matrices kernels
594 amrex::ParticleReal fpxx, fpxy, fpxz;
595 amrex::ParticleReal fpyx, fpyy, fpyz;
596 amrex::ParticleReal fpzx, fpzy, fpzz;
597 setMassMatricesKernels( qs, ms, dt, rhop,
598 uxp_nph[ip], uyp_nph[ip], uzp_nph[ip],
599 Bxp, Byp, Bzp,
600 fpxx, fpxy, fpxz,
601 fpyx, fpyy, fpyz,
602 fpzx, fpzy, fpzz );
603
604 // Pass dummy arrays for Jx, Jy, Jz (which will not be used)
605 amrex::Array4<amrex::Real> const dummy_Jx{};
606 amrex::Array4<amrex::Real> const dummy_Jy{};
607 amrex::Array4<amrex::Real> const dummy_Jz{};
608
609 doDirectJandSigmaDepositionKernel<depos_order,full_mass_matrices,/*deposit_J=*/false>(
610 xp_nph, yp_nph, zp_nph,
611 wqx, wqy, wqz,
612 fpxx, fpxy, fpxz,
613 fpyx, fpyy, fpyz,
614 fpzx, fpzy, fpzz,
615 dummy_Jx, dummy_Jy, dummy_Jz,
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, bool deposit_J>
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 if constexpr (deposit_J) {
957 amrex::Gpu::Atomic::AddNoRet( &Jx_arr(lo.x+i0_cell+i, lo.y+j0_node+j, lo.z+k0_node+k), wqx*weight);
958 }
959 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);
960 }
961 }
962 }
963
964 // deposit Jy and Syy or this segment
965 for (int i=0; i<=depos_order; i++) {
966 for (int j=0; j<=depos_order-1; j++) {
967 for (int k=0; k<=depos_order; k++) {
968 weight = sy_cell[j]*( sx_old_node[i]*sz_old_node[k]*one_third
969 + sx_old_node[i]*sz_new_node[k]*one_sixth
970 + sx_new_node[i]*sz_old_node[k]*one_sixth
971 + sx_new_node[i]*sz_new_node[k]*one_third )*seg_factor_y;
972 if constexpr (deposit_J) {
973 amrex::Gpu::Atomic::AddNoRet( &Jy_arr(lo.x+i0_node+i, lo.y+j0_cell+j, lo.z+k0_node+k), wqy*weight);
974 }
975 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);
976 }
977 }
978 }
979
980 // deposit Jz and Sz for this segment
981 for (int i=0; i<=depos_order; i++) {
982 for (int j=0; j<=depos_order; j++) {
983 for (int k=0; k<=depos_order-1; k++) {
984 weight = sz_cell[k]*( sx_old_node[i]*sy_old_node[j]*one_third
985 + sx_old_node[i]*sy_new_node[j]*one_sixth
986 + sx_new_node[i]*sy_old_node[j]*one_sixth
987 + sx_new_node[i]*sy_new_node[j]*one_third )*seg_factor_z;
988 if constexpr (deposit_J) {
989 amrex::Gpu::Atomic::AddNoRet( &Jz_arr(lo.x+i0_node+i, lo.y+j0_node+j, lo.z+k0_cell+k), wqz*weight);
990 }
991 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);
992 }
993 }
994 }
995
996 // update old segment values
997 if (ns < num_segments-1) {
998 x0_old = x0_new;
999 y0_old = y0_new;
1000 z0_old = z0_new;
1001 }
1002
1003 } // end loop over segments
1004
1005#elif defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ)
1006
1007 // compute cell crossings in X-direction
1008 const auto i_old = static_cast<int>(x_old-shift);
1009 const auto i_new = static_cast<int>(x_new-shift);
1010 const int cell_crossings_x = std::abs(i_new-i_old);
1011 num_segments += cell_crossings_x;
1012
1013 // compute cell crossings in Z-direction
1014 const auto k_old = static_cast<int>(z_old-shift);
1015 const auto k_new = static_cast<int>(z_new-shift);
1016 const int cell_crossings_z = std::abs(k_new-k_old);
1017 num_segments += cell_crossings_z;
1018
1019 // Compute initial particle cell locations in each direction
1020 // used to find the position at cell crossings.
1021 // Keep these double to avoid bug in single precision
1022 const auto dirX_sign = static_cast<double>(dxp < 0. ? -1. : 1.);
1023 const auto dirZ_sign = static_cast<double>(dzp < 0. ? -1. : 1.);
1024 double Xcell = 0., Zcell = 0.;
1025 if (num_segments > 1) {
1026 Xcell = static_cast<double>(i_old) + shift + 0.5*(1.-dirX_sign);
1027 Zcell = static_cast<double>(k_old) + shift + 0.5*(1.-dirZ_sign);
1028 }
1029
1030 // loop over the number of segments and deposit
1031 const Compute_shape_factor< depos_order-1 > compute_shape_factor_cell;
1032 const Compute_shape_factor_pair< depos_order > compute_shape_factors_node;
1033 double dxp_seg, dzp_seg;
1034 double x0_new, z0_new;
1035 double x0_old = x_old;
1036 double z0_old = z_old;
1037
1038 constexpr int num_segments_max = 1 + 4*AMREX_SPACEDIM;
1039 AMREX_ALWAYS_ASSERT_WITH_MESSAGE( num_segments <= num_segments_max,
1040 "Error: num_segments must be less than or equal to 1 + 4*AMREX_SPACEDIM.");
1041
1042 // Save the start index and interpolation weights for each segment
1043 int i0_cell[num_segments_max];
1044 int i0_node[num_segments_max];
1045 int k0_cell[num_segments_max];
1046 int k0_node[num_segments_max];
1047 amrex::Real weight_cellX_nodeZ[num_segments_max][depos_order][depos_order+1];
1048 amrex::Real weight_nodeX_cellZ[num_segments_max][depos_order+1][depos_order];
1049 amrex::Real weight_nodeX_nodeZ[num_segments_max][depos_order+1][depos_order+1];
1050
1051 const auto i_mid = static_cast<int>(0.5*(x_new+x_old)-shift);
1052 const auto k_mid = static_cast<int>(0.5*(z_new+z_old)-shift);
1053 int SegNumX[num_segments_max];
1054 int SegNumZ[num_segments_max];
1055
1056 for (int ns=0; ns<num_segments; ns++) {
1057
1058 if (ns == num_segments-1) { // final segment
1059
1060 x0_new = x_new;
1061 z0_new = z_new;
1062 dxp_seg = x0_new - x0_old;
1063 dzp_seg = z0_new - z0_old;
1064
1065 }
1066 else {
1067
1068 x0_new = Xcell + dirX_sign;
1069 z0_new = Zcell + dirZ_sign;
1070 dxp_seg = x0_new - x0_old;
1071 dzp_seg = z0_new - z0_old;
1072
1073 if (dzp == 0. || std::abs(dxp_seg) < std::abs(dxp/dzp*dzp_seg)) {
1074 Xcell = x0_new;
1075 dzp_seg = dzp/dxp*dxp_seg;
1076 z0_new = z0_old + dzp_seg;
1077 }
1078 else {
1079 Zcell = z0_new;
1080 dxp_seg = dxp/dzp*dzp_seg;
1081 x0_new = x0_old + dxp_seg;
1082 }
1083
1084 }
1085
1086 // Compute the segment factors (each equal to dt_seg/dt for nonzero dxp, or dzp)
1087 const auto seg_factor_x = static_cast<amrex::Real>(dxp == 0. ? 1._rt : dxp_seg/dxp);
1088 const auto seg_factor_z = static_cast<amrex::Real>(dzp == 0. ? 1._rt : dzp_seg/dzp);
1089
1090 // Compute cell-based weights using the average segment position
1091 // Keep these double to avoid bug in single precision
1092 double sx_cell[depos_order] = {0.};
1093 double sz_cell[depos_order] = {0.};
1094 double const x0_bar = (x0_new + x0_old)/2.0;
1095 double const z0_bar = (z0_new + z0_old)/2.0;
1096 i0_cell[ns] = compute_shape_factor_cell( sx_cell, x0_bar-0.5 );
1097 k0_cell[ns] = compute_shape_factor_cell( sz_cell, z0_bar-0.5 );
1098
1099 // Set the segment number for the mass matrix component calc
1100 if constexpr (full_mass_matrices) {
1101 const auto i0_mid = static_cast<int>(x0_bar-shift);
1102 const auto k0_mid = static_cast<int>(z0_bar-shift);
1103 SegNumX[ns] = 1 + i0_mid - i_mid;
1104 SegNumZ[ns] = 1 + k0_mid - k_mid;
1105 }
1106
1107 if constexpr (depos_order >= 3) { // higher-order correction to the cell-based weights
1108 const Compute_shape_factor_pair<depos_order-1> compute_shape_factors_cell;
1109 double sx_old_cell[depos_order] = {0.};
1110 double sx_new_cell[depos_order] = {0.};
1111 double sz_old_cell[depos_order] = {0.};
1112 double sz_new_cell[depos_order] = {0.};
1113 const int i0_cell_2 = compute_shape_factors_cell( sx_old_cell, sx_new_cell, x0_old-0.5, x0_new-0.5 );
1114 const int k0_cell_2 = compute_shape_factors_cell( sz_old_cell, sz_new_cell, z0_old-0.5, z0_new-0.5 );
1115 amrex::ignore_unused(i0_cell_2, k0_cell_2);
1116 for (int m=0; m<depos_order; m++) {
1117 sx_cell[m] = (4.0*sx_cell[m] + sx_old_cell[m] + sx_new_cell[m])/6.0;
1118 sz_cell[m] = (4.0*sz_cell[m] + sz_old_cell[m] + sz_new_cell[m])/6.0;
1119 }
1120 }
1121
1122 // Compute node-based weights using the old and new segment positions
1123 // Keep these double to avoid bug in single precision
1124 double sx_old_node[depos_order+1] = {0.};
1125 double sx_new_node[depos_order+1] = {0.};
1126 double sz_old_node[depos_order+1] = {0.};
1127 double sz_new_node[depos_order+1] = {0.};
1128 i0_node[ns] = compute_shape_factors_node( sx_old_node, sx_new_node, x0_old, x0_new );
1129 k0_node[ns] = compute_shape_factors_node( sz_old_node, sz_new_node, z0_old, z0_new );
1130
1131 // deposit Jx and Sx for this segment
1132 amrex::Real weight;
1133 for (int i=0; i<=depos_order-1; i++) {
1134 for (int k=0; k<=depos_order; k++) {
1135 const int i_J = lo.x + i0_cell[ns] + i;
1136 const int k_J = lo.y + k0_node[ns] + k;
1137 weight = sx_cell[i]*(sz_old_node[k] + sz_new_node[k])/2.0_rt*seg_factor_x;
1138 if constexpr (deposit_J) {
1139 amrex::Gpu::Atomic::AddNoRet(&Jx_arr(i_J, k_J, 0, 0), wqx*weight);
1140 }
1141 if constexpr (full_mass_matrices) { weight_cellX_nodeZ[ns][i][k] = weight; }
1142 else {
1143 amrex::Gpu::Atomic::AddNoRet(&Sxx_arr(i_J, k_J, 0, 0), fpxx*weight*weight);
1144 }
1145 }
1146 }
1147
1148 // deposit out-of-plane Jy and Sy for this segment
1149 const auto seg_factor_y = std::min(seg_factor_x,seg_factor_z);
1150 for (int i=0; i<=depos_order; i++) {
1151 for (int k=0; k<=depos_order; k++) {
1152 const int i_J = lo.x + i0_node[ns] + i;
1153 const int k_J = lo.y + k0_node[ns] + k;
1154 weight = ( sx_old_node[i]*sz_old_node[k]*one_third
1155 + sx_old_node[i]*sz_new_node[k]*one_sixth
1156 + sx_new_node[i]*sz_old_node[k]*one_sixth
1157 + sx_new_node[i]*sz_new_node[k]*one_third )*seg_factor_y;
1158 if constexpr (deposit_J) {
1159 amrex::Gpu::Atomic::AddNoRet(&Jy_arr(i_J, k_J, 0, 0), wqy*weight);
1160 }
1161 if constexpr (full_mass_matrices) { weight_nodeX_nodeZ[ns][i][k] = weight; }
1162 else {
1163 amrex::Gpu::Atomic::AddNoRet(&Syy_arr(i_J, k_J, 0, 0), fpyy*weight*weight);
1164 }
1165 }
1166 }
1167
1168 // deposit Jz and Szz for this segment
1169 for (int i=0; i<=depos_order; i++) {
1170 for (int k=0; k<=depos_order-1; k++) {
1171 const int i_J = lo.x + i0_node[ns] + i;
1172 const int k_J = lo.y + k0_cell[ns] + k;
1173 weight = sz_cell[k]*(sx_old_node[i] + sx_new_node[i])/2.0_rt*seg_factor_z;
1174 if constexpr (deposit_J) {
1175 amrex::Gpu::Atomic::AddNoRet(&Jz_arr(i_J, k_J, 0, 0), wqz*weight);
1176 }
1177 if constexpr (full_mass_matrices) { weight_nodeX_cellZ[ns][i][k] = weight; }
1178 else {
1179 amrex::Gpu::Atomic::AddNoRet(&Szz_arr(i_J, k_J, 0, 0), fpzz*weight*weight);
1180 }
1181 }
1182 }
1183
1184 // update old segment values
1185 if (ns < num_segments-1) {
1186 x0_old = x0_new;
1187 z0_old = z0_new;
1188 }
1189
1190 } // end loop over segments
1191
1192 if constexpr (full_mass_matrices) {
1193
1194 const int Ncomp_base = 2*depos_order + 2*max_crossings;
1195 const int Ncomp_xx0 = Ncomp_base - 1;
1196 const int Ncomp_xz0 = Ncomp_base;
1197 const int Ncomp_zx0 = Ncomp_base;
1198 const int Ncomp_zz0 = 1 + Ncomp_base;
1199 const int Ncomp_xy0 = Ncomp_base;
1200 const int Ncomp_yx0 = Ncomp_base;
1201 const int Ncomp_yz0 = 1 + Ncomp_base;
1202 const int Ncomp_yy0 = 1 + Ncomp_base;
1203 const int Ncomp_zy0 = 1 + Ncomp_base;
1204
1205 const int width_xx1 = depos_order + max_crossings; // (Ncomp_xx[1] - 1)/2
1206 const int width_zz1 = width_xx1 - 1; // (Ncomp_zz[1] - 1)/2
1207 const int width_yy1 = width_xx1; // (Ncomp_yy[1] - 1)/2
1208
1209 // Loop over segments and deposit full mass matrices
1210 for (int ns=0; ns<num_segments; ns++) {
1211
1212 // Deposit Sxx, Sxz, and Sxy for this segment
1213 for (int i = 0; i <= depos_order - 1; i++) {
1214 for (int k = 0; k <= depos_order; k++) {
1215 const int i_J = lo.x + i0_cell[ns] + i;
1216 const int k_J = lo.y + k0_node[ns] + k;
1217 const amrex::Real weight_J = weight_cellX_nodeZ[ns][i][k];
1218 for (int ms=0; ms<num_segments; ms++) {
1219 const int SegShiftX = max_crossings + SegNumX[ms] - SegNumX[ns];
1220 const int SegShiftZ = max_crossings + SegNumZ[ms] - SegNumZ[ns];
1221 // Deposit Sxx
1222 for (int kE = 0; kE <= depos_order; kE++) {
1223 const int row_xx = depos_order - k + kE + SegShiftZ;
1224 const int above_diag = (row_xx > width_xx1) ? 1 : 0;
1225 for (int iE = 0; iE <= depos_order - 1; iE++) {
1226 const int col_xx = depos_order - 1 - i + iE + SegShiftX;
1227 if (col_xx > Ncomp_xx0 - row_xx - above_diag) { break; } // Reduced deposit for diagonal mass matrices
1228 const int comp_xx = col_xx + Ncomp_xx0*row_xx;
1229 const amrex::Real weight_E = weight_cellX_nodeZ[ms][iE][kE];
1230 amrex::Gpu::Atomic::AddNoRet(&Sxx_arr(i_J, k_J, 0, comp_xx), fpxx*weight_J*weight_E);
1231 }
1232 }
1233 // Deposit Sxz
1234 for (int iE = 0; iE <= depos_order; iE++) {
1235 for (int kE = 0; kE <= depos_order - 1; kE++) {
1236 const amrex::Real weight_E = weight_nodeX_cellZ[ms][iE][kE];
1237 const int comp_xz = depos_order - 1 - i + iE + SegShiftX
1238 + Ncomp_xz0*(depos_order - k + kE + SegShiftZ);
1239 amrex::Gpu::Atomic::AddNoRet(&Sxz_arr(i_J, k_J, 0, comp_xz), fpxz*weight_J*weight_E);
1240 }
1241 }
1242 // Deposit Sxy
1243 for (int iE = 0; iE <= depos_order; iE++) {
1244 for (int kE = 0; kE <= depos_order; kE++) {
1245 const amrex::Real weight_E = weight_nodeX_nodeZ[ms][iE][kE];
1246 const int comp_xy = depos_order - 1 - i + iE + SegShiftX
1247 + Ncomp_xy0*(depos_order - k + kE + SegShiftZ);
1248 amrex::Gpu::Atomic::AddNoRet(&Sxy_arr(i_J, k_J, 0, comp_xy), fpxy*weight_J*weight_E);
1249 }
1250 }
1251 }
1252 }
1253 }
1254
1255 // Deposit Szx, Szz, and Szy for this segment
1256 for (int i=0; i<=depos_order; i++) {
1257 for (int k=0; k<=depos_order-1; k++) {
1258 const int i_J = lo.x + i0_node[ns] + i;
1259 const int k_J = lo.y + k0_cell[ns] + k;
1260 const amrex::Real weight_J = weight_nodeX_cellZ[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 Szx
1265 for (int iE=0; iE<=depos_order-1; iE++) {
1266 for (int kE=0; kE<=depos_order; kE++) {
1267 const amrex::Real weight_E = weight_cellX_nodeZ[ms][iE][kE];
1268 const int comp_zx = depos_order - i + iE + SegShiftX
1269 + Ncomp_zx0*(depos_order-1 - k + kE + SegShiftZ);
1270 amrex::Gpu::Atomic::AddNoRet( &Szx_arr(i_J, k_J, 0, comp_zx), fpzx*weight_J*weight_E);
1271 }
1272 }
1273 // Deposit Szz
1274 for (int kE=0; kE<=depos_order-1; kE++) {
1275 const int row_zz = depos_order - 1 - k + kE + SegShiftZ;
1276 const int above_diag = (row_zz > width_zz1) ? 1 : 0;
1277 for (int iE=0; iE<=depos_order; iE++) {
1278 const int col_zz = depos_order - i + iE + SegShiftX;
1279 if (col_zz > Ncomp_zz0 - 2 - row_zz - above_diag) { break; } // Reduced deposit for diagonal mass matrices
1280 const int comp_zz = col_zz + Ncomp_zz0*row_zz;
1281 const amrex::Real weight_E = weight_nodeX_cellZ[ms][iE][kE];
1282 amrex::Gpu::Atomic::AddNoRet( &Szz_arr(i_J, k_J, 0, comp_zz), fpzz*weight_J*weight_E);
1283 }
1284 }
1285 // Deposit Szy
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_zy = depos_order - i + iE + SegShiftX
1290 + Ncomp_zy0*(depos_order-1 - k + kE + SegShiftZ);
1291 amrex::Gpu::Atomic::AddNoRet( &Szy_arr(i_J, k_J, 0, comp_zy), fpzy*weight_J*weight_E);
1292 }
1293 }
1294 }
1295 }
1296 }
1297
1298 // Deposit Syx, Syz, and Syy for this segment
1299 for (int i=0; i<=depos_order; i++) {
1300 for (int k=0; k<=depos_order; k++) {
1301 const int i_J = lo.x + i0_node[ns] + i;
1302 const int k_J = lo.y + k0_node[ns] + k;
1303 const amrex::Real weight_J = weight_nodeX_nodeZ[ns][i][k];
1304 for (int ms=0; ms<num_segments; ms++) {
1305 const int SegShiftX = max_crossings + SegNumX[ms] - SegNumX[ns];
1306 const int SegShiftZ = max_crossings + SegNumZ[ms] - SegNumZ[ns];
1307 // Deposit Syx
1308 for (int iE=0; iE<=depos_order-1; iE++) {
1309 for (int kE=0; kE<=depos_order; kE++) {
1310 const amrex::Real weight_E = weight_cellX_nodeZ[ms][iE][kE];
1311 const int comp_yx = depos_order - i + iE + SegShiftX
1312 + Ncomp_yx0*(depos_order - k + kE + SegShiftZ);
1313 amrex::Gpu::Atomic::AddNoRet( &Syx_arr(i_J, k_J, 0, comp_yx), fpyx*weight_J*weight_E);
1314 }
1315 }
1316 // Deposit Syz
1317 for (int iE=0; iE<=depos_order; iE++) {
1318 for (int kE=0; kE<=depos_order-1; kE++) {
1319 const amrex::Real weight_E = weight_nodeX_cellZ[ms][iE][kE];
1320 const int comp_yz = depos_order - i + iE + SegShiftX
1321 + Ncomp_yz0*(depos_order - k + kE + SegShiftZ);
1322 amrex::Gpu::Atomic::AddNoRet( &Syz_arr(i_J, k_J, 0, comp_yz), fpyz*weight_J*weight_E);
1323 }
1324 }
1325 // Deposit Syy
1326 for (int kE=0; kE<=depos_order; kE++) {
1327 const int row_yy = depos_order - k + kE + SegShiftZ;
1328 const int above_diag = (row_yy > width_yy1) ? 1 : 0;
1329 for (int iE=0; iE<=depos_order; iE++) {
1330 const int col_yy = depos_order - i + iE + SegShiftX;
1331 if (col_yy > Ncomp_yy0 - 1 - row_yy - above_diag) { break; } // Reduced deposit for diagonal mass matrices
1332 const int comp_yy = col_yy + Ncomp_yy0*row_yy;
1333 const amrex::Real weight_E = weight_nodeX_nodeZ[ms][iE][kE];
1334 amrex::Gpu::Atomic::AddNoRet( &Syy_arr(i_J, k_J, 0, comp_yy), fpyy*weight_J*weight_E);
1335 }
1336 }
1337
1338 }
1339 }
1340 }
1341
1342 }
1343
1344 }
1345
1346#elif defined(WARPX_DIM_RCYLINDER) || defined(WARPX_DIM_RSPHERE)
1347
1348 // compute cell crossings in X-direction
1349 const auto i_old = static_cast<int>(x_old-shift);
1350 const auto i_new = static_cast<int>(x_new-shift);
1351 const int cell_crossings_x = std::abs(i_new-i_old);
1352 num_segments += cell_crossings_x;
1353
1354 // Compute the initial cell location used to find the cell crossings.
1355 // Keep these double to avoid bug in single precision
1356 const auto dirX_sign = static_cast<double>(dxp < 0. ? -1. : 1.);
1357 double Xcell = static_cast<double>(i_old) + shift + 0.5*(1.-dirX_sign);
1358
1359 // loop over the number of segments and deposit
1360 const Compute_shape_factor< depos_order-1 > compute_shape_factor_cell;
1361 const Compute_shape_factor_pair< depos_order > compute_shape_factors_node;
1362 double dxp_seg;
1363 double x0_new;
1364 double x0_old = x_old;
1365
1366 for (int ns=0; ns<num_segments; ns++) {
1367
1368 if (ns == num_segments-1) { // final segment
1369 x0_new = x_new;
1370 dxp_seg = x0_new - x0_old;
1371 }
1372 else {
1373 Xcell = Xcell + dirX_sign;
1374 x0_new = Xcell;
1375 dxp_seg = x0_new - x0_old;
1376 }
1377
1378 // Compute the segment factor (equal to dt_seg/dt for nonzero dxp)
1379 const auto seg_factor = static_cast<amrex::Real>(dxp == 0. ? 1._rt : dxp_seg/dxp);
1380
1381 // Compute cell-based weights using the average segment position
1382 // Keep these double to avoid bug in single precision
1383 double sx_cell[depos_order] = {0.};
1384 double const x0_bar = (x0_new + x0_old)/2.0;
1385 const int i0_cell = compute_shape_factor_cell( sx_cell, x0_bar-0.5 );
1386
1387 if constexpr (depos_order >= 3) { // higher-order correction to the cell-based weights
1388 const Compute_shape_factor_pair<depos_order-1> compute_shape_factors_cell;
1389 double sx_old_cell[depos_order] = {0.};
1390 double sx_new_cell[depos_order] = {0.};
1391 const int i0_cell_2 = compute_shape_factors_cell( sx_old_cell, sx_new_cell, x0_old-0.5, x0_new-0.5 );
1392 amrex::ignore_unused(i0_cell_2);
1393 for (int m=0; m<depos_order; m++) {
1394 sx_cell[m] = (4.0*sx_cell[m] + sx_old_cell[m] + sx_new_cell[m])/6.0;
1395 }
1396 }
1397
1398 // Compute node-based weights using the old and new segment positions
1399 // Keep these double to avoid bug in single precision
1400 double sx_old_node[depos_order+1] = {0.};
1401 double sx_new_node[depos_order+1] = {0.};
1402 const int i0_node = compute_shape_factors_node( sx_old_node, sx_new_node, x0_old, x0_new );
1403
1404 // deposit out-of-plane Jy, Jz, Syy, and Szz for this segment
1405 for (int i=0; i<=depos_order; i++) {
1406 const amrex::Real weight = 0.5_rt*(sx_old_node[i] + sx_new_node[i])*seg_factor;
1407 if constexpr (deposit_J) {
1408 amrex::Gpu::Atomic::AddNoRet( &Jy_arr(lo.x+i0_node+i, 0, 0), wqy*weight);
1409 amrex::Gpu::Atomic::AddNoRet( &Jz_arr(lo.x+i0_node+i, 0, 0), wqz*weight);
1410 }
1411 //
1412 amrex::Gpu::Atomic::AddNoRet( &Syy_arr(lo.x+i0_node+i, 0, 0, 0), fpyy*weight*weight);
1413 amrex::Gpu::Atomic::AddNoRet( &Szz_arr(lo.x+i0_node+i, 0, 0, 0), fpzz*weight*weight);
1414 }
1415
1416 // deposit Jx and Sxx for this segment
1417 for (int i=0; i<=depos_order-1; i++) {
1418 const amrex::Real weight = sx_cell[i]*seg_factor;
1419 if constexpr (deposit_J) {
1420 amrex::Gpu::Atomic::AddNoRet( &Jx_arr(lo.x+i0_cell+i, 0, 0), wqx*weight);
1421 }
1422 //
1423 amrex::Gpu::Atomic::AddNoRet( &Sxx_arr(lo.x+i0_cell+i, 0, 0, 0), fpxx*weight*weight);
1424 }
1425
1426 // update old segment values
1427 if (ns < num_segments-1) {
1428 x0_old = x0_new;
1429 }
1430
1431 }
1432
1433#elif defined(WARPX_DIM_1D_Z)
1434
1435 // compute cell crossings in Z-direction
1436 const auto k_old = static_cast<int>(z_old-shift);
1437 const auto k_new = static_cast<int>(z_new-shift);
1438 const int cell_crossings_z = std::abs(k_new-k_old);
1439 num_segments += cell_crossings_z;
1440
1441 // Compute initial particle cell location used to find cell crossings.
1442 // Keep these double to avoid bug in single precision
1443 const auto dirZ_sign = static_cast<double>(dzp < 0. ? -1. : 1.);
1444 double Zcell = static_cast<double>(k_old) + shift + 0.5*(1.-dirZ_sign);
1445
1446 // loop over the number of segments and deposit
1447 const Compute_shape_factor< depos_order-1 > compute_shape_factor_cell;
1448 const Compute_shape_factor_pair< depos_order > compute_shape_factors_node;
1449 double dzp_seg;
1450 double z0_new;
1451 double z0_old = z_old;
1452
1453 constexpr int num_segments_max = 1 + 4*AMREX_SPACEDIM;
1454 AMREX_ALWAYS_ASSERT_WITH_MESSAGE( num_segments <= num_segments_max,
1455 "Error: num_segments must be less than or equal to 1 + 4*AMREX_SPACEDIM.");
1456
1457 // Save the start index and interpolation weights for each segment
1458 int k0_cell[num_segments_max];
1459 int k0_node[num_segments_max];
1460 amrex::Real weight_cell[num_segments_max][depos_order];
1461 amrex::Real weight_node[num_segments_max][depos_order+1];
1462
1463 const auto k_mid = static_cast<int>(0.5*(z_new+z_old)-shift);
1464 int SegNum[num_segments_max];
1465
1466 for (int ns=0; ns<num_segments; ns++) {
1467
1468 if (ns == num_segments-1) { // final segment
1469 z0_new = z_new;
1470 dzp_seg = z0_new - z0_old;
1471 }
1472 else {
1473 Zcell = Zcell + dirZ_sign;
1474 z0_new = Zcell;
1475 dzp_seg = z0_new - z0_old;
1476 }
1477
1478 // Compute the segment factor (equal to dt_seg/dt for nonzero dzp)
1479 const auto seg_factor = static_cast<amrex::Real>(dzp == 0. ? 1._rt : dzp_seg/dzp);
1480
1481 // Compute cell-based weights using the average segment position
1482 // Keep these double to avoid bug in single precision
1483 double sz_cell[depos_order] = {0.};
1484 double const z0_bar = (z0_new + z0_old)/2.0;
1485 k0_cell[ns] = compute_shape_factor_cell( sz_cell, z0_bar-0.5 );
1486
1487 // Set the segment number for the mass matrix component calc
1488 if constexpr (full_mass_matrices) {
1489 const auto k0_mid = static_cast<int>(z0_bar-shift);
1490 SegNum[ns] = 1 + k0_mid - k_mid;
1491 }
1492
1493 if constexpr (depos_order >= 3) { // higher-order correction to the cell-based weights
1494 const Compute_shape_factor_pair<depos_order-1> compute_shape_factors_cell;
1495 double sz_old_cell[depos_order] = {0.};
1496 double sz_new_cell[depos_order] = {0.};
1497 const int k0_cell_2 = compute_shape_factors_cell( sz_old_cell, sz_new_cell, z0_old-0.5, z0_new-0.5 );
1498 amrex::ignore_unused(k0_cell_2);
1499 for (int m=0; m<depos_order; m++) {
1500 sz_cell[m] = (4.0*sz_cell[m] + sz_old_cell[m] + sz_new_cell[m])/6.0;
1501 }
1502 }
1503
1504 // Compute node-based weights using the old and new segment positions
1505 // Keep these double to avoid bug in single precision
1506 double sz_old_node[depos_order+1] = {0.};
1507 double sz_new_node[depos_order+1] = {0.};
1508 k0_node[ns] = compute_shape_factors_node( sz_old_node, sz_new_node, z0_old, z0_new );
1509
1510 // deposit out-of-plane Jx, Jy, Sx, and Sy for this segment
1511 for (int k=0; k<=depos_order; k++) {
1512 const amrex::Real weight = 0.5_rt*(sz_old_node[k] + sz_new_node[k])*seg_factor;
1513 const int k_J = lo.x + k0_node[ns] + k;
1514 if constexpr (deposit_J) {
1515 amrex::Gpu::Atomic::AddNoRet(&Jx_arr(k_J, 0, 0), wqx*weight);
1516 amrex::Gpu::Atomic::AddNoRet(&Jy_arr(k_J, 0, 0), wqy*weight);
1517 }
1518 if constexpr (full_mass_matrices) { weight_node[ns][k] = weight; }
1519 else {
1520 amrex::Gpu::Atomic::AddNoRet(&Sxx_arr(k_J, 0, 0, 0), fpxx*weight*weight);
1521 amrex::Gpu::Atomic::AddNoRet(&Syy_arr(k_J, 0, 0, 0), fpyy*weight*weight);
1522 }
1523 }
1524
1525 // deposit Jz and Szz for this segment
1526 for (int k=0; k<=depos_order-1; k++) {
1527 const amrex::Real weight = sz_cell[k]*seg_factor;
1528 const int k_J = lo.x + k0_cell[ns] + k;
1529 if constexpr (deposit_J) {
1530 amrex::Gpu::Atomic::AddNoRet(&Jz_arr(k_J, 0, 0), wqz*weight);
1531 }
1532 if constexpr (full_mass_matrices) { weight_cell[ns][k] = weight; }
1533 else {
1534 amrex::Gpu::Atomic::AddNoRet(&Szz_arr(k_J, 0, 0, 0), fpzz*weight*weight);
1535 }
1536 }
1537
1538 // update old segment values
1539 if (ns < num_segments-1) {
1540 z0_old = z0_new;
1541 }
1542
1543 }
1544
1545 if constexpr (full_mass_matrices) {
1546
1547 const int width_zz = depos_order - 1 + max_crossings;
1548 const int width_yy = depos_order + max_crossings;
1549
1550 // Loop over segments and deposit full mass matrices
1551 for (int ns=0; ns<num_segments; ns++) {
1552
1553 // Deposit Sxx, Sxy, Sxz, Syx, Syy, and Syz for this segment
1554 for (int k=0; k<=depos_order; k++) {
1555
1556 const int k_J = lo.x + k0_node[ns] + k;
1557 const amrex::Real weight_J = weight_node[ns][k];
1558 for (int ms=0; ms<num_segments; ms++) {
1559 const int SegShift = max_crossings + SegNum[ms] - SegNum[ns];
1560 for (int kE=0; kE<=depos_order; kE++) {
1561 const amrex::Real weight_E = weight_node[ms][kE];
1562 const int comp_yy = depos_order - k + kE + SegShift;
1563 if (comp_yy <= width_yy) { // Reduced deposit for diagonal mass matrices
1564 amrex::Gpu::Atomic::AddNoRet(&Sxx_arr(k_J, 0, 0, comp_yy), fpxx*weight_J*weight_E);
1565 amrex::Gpu::Atomic::AddNoRet(&Syy_arr(k_J, 0, 0, comp_yy), fpyy*weight_J*weight_E);
1566 }
1567 amrex::Gpu::Atomic::AddNoRet(&Sxy_arr(k_J, 0, 0, comp_yy), fpxy*weight_J*weight_E);
1568 amrex::Gpu::Atomic::AddNoRet(&Syx_arr(k_J, 0, 0, comp_yy), fpyx*weight_J*weight_E);
1569 }
1570 for (int kE=0; kE<=depos_order-1; kE++) {
1571 const amrex::Real weight_E = weight_cell[ms][kE];
1572 const int comp_yz = depos_order - k + kE + SegShift;
1573 amrex::Gpu::Atomic::AddNoRet(&Sxz_arr(k_J, 0, 0, comp_yz), fpxz*weight_J*weight_E);
1574 amrex::Gpu::Atomic::AddNoRet(&Syz_arr(k_J, 0, 0, comp_yz), fpyz*weight_J*weight_E);
1575 }
1576 }
1577
1578 }
1579
1580 // Deposit Szx, Szy, and Szz for this segment
1581 for (int k=0; k<=depos_order-1; k++) {
1582
1583 const int k_J = lo.x + k0_cell[ns] + k;
1584 const amrex::Real weight_J = weight_cell[ns][k];
1585 for (int ms=0; ms<num_segments; ms++) {
1586 const int SegShift = max_crossings + SegNum[ms] - SegNum[ns];
1587 for (int kE=0; kE<=depos_order-1; kE++) {
1588 const int comp_zz = depos_order - 1 - k + kE + SegShift;
1589 if (comp_zz > width_zz) { break; } // Reduced deposit for diagonal mass matrices
1590 const amrex::Real weight_E = weight_cell[ms][kE];
1591 amrex::Gpu::Atomic::AddNoRet(&Szz_arr(k_J, 0, 0, comp_zz), fpzz*weight_J*weight_E);
1592 }
1593 for (int kE=0; kE<=depos_order; kE++) {
1594 const amrex::Real weight_E = weight_node[ms][kE];
1595 const int comp_zy = depos_order-1 - k + kE + SegShift;
1596 amrex::Gpu::Atomic::AddNoRet(&Szx_arr(k_J, 0, 0, comp_zy), fpzx*weight_J*weight_E);
1597 amrex::Gpu::Atomic::AddNoRet(&Szy_arr(k_J, 0, 0, comp_zy), fpzy*weight_J*weight_E);
1598 }
1599 }
1600
1601 }
1602
1603 }
1604
1605 }
1606
1607#endif
1608}
1609
1638template <int depos_order, bool full_mass_matrices>
1639void doVillasenorSigmaDeposition ([[maybe_unused]] const amrex::ParticleReal* xp_n_data,
1640 [[maybe_unused]] const amrex::ParticleReal* yp_n_data,
1641 [[maybe_unused]] const amrex::ParticleReal* zp_n_data,
1642 const GetParticlePosition<PIdx>& GetPosition,
1643 [[maybe_unused]] const int* nsuborbits,
1644 const amrex::ParticleReal* wp,
1645 const amrex::ParticleReal* uxp_n,
1646 const amrex::ParticleReal* uyp_n,
1647 const amrex::ParticleReal* uzp_n,
1648 const amrex::ParticleReal* uxp_nph,
1649 const amrex::ParticleReal* uyp_nph,
1650 const amrex::ParticleReal* uzp_nph,
1651 const int max_crossings,
1652 amrex::Array4<amrex::Real> const& Sxx_arr,
1653 amrex::Array4<amrex::Real> const& Sxy_arr,
1654 amrex::Array4<amrex::Real> const& Sxz_arr,
1655 amrex::Array4<amrex::Real> const& Syx_arr,
1656 amrex::Array4<amrex::Real> const& Syy_arr,
1657 amrex::Array4<amrex::Real> const& Syz_arr,
1658 amrex::Array4<amrex::Real> const& Szx_arr,
1659 amrex::Array4<amrex::Real> const& Szy_arr,
1660 amrex::Array4<amrex::Real> const& Szz_arr,
1661 const amrex::Array4<amrex::Real const>& Bx_arr,
1662 const amrex::Array4<amrex::Real const>& By_arr,
1663 const amrex::Array4<amrex::Real const>& Bz_arr,
1664 const amrex::IndexType Bx_type,
1665 const amrex::IndexType By_type,
1666 const amrex::IndexType Bz_type,
1667 const long np_to_deposit,
1668 const amrex::Real dt,
1669 const amrex::XDim3& dinv,
1670 const amrex::XDim3& xyzmin,
1671 const amrex::GpuArray<amrex::GpuArray<double,2>, AMREX_SPACEDIM> & domain_double,
1672 const amrex::GpuArray<amrex::GpuArray<bool,2>, AMREX_SPACEDIM> & do_cropping,
1673 const amrex::Dim3 lo,
1674 const amrex::Real qs,
1675 const amrex::Real ms)
1676{
1677 using namespace amrex::literals;
1678
1679 const amrex::Real invvol = dinv.x*dinv.y*dinv.z;
1680
1681 // Loop over particles and deposit mass matrices
1682 amrex::ParallelFor(
1683 np_to_deposit,
1684 [=] AMREX_GPU_DEVICE (long const ip) {
1685
1686 // Skip mass matrix deposition for particles with suborbits.
1687 if (nsuborbits && nsuborbits[ip] > 1) { return; }
1688
1689 amrex::ParticleReal xp_nph, yp_nph, zp_nph;
1690 GetPosition(ip, xp_nph, yp_nph, zp_nph);
1691
1692 // Compute magnetic field on particle
1693 amrex::ParticleReal Bxp = 0.0;
1694 amrex::ParticleReal Byp = 0.0;
1695 amrex::ParticleReal Bzp = 0.0;
1696 const int depos_order_perp = 1;
1697 const int depos_order_para = 1;
1698 doDirectGatherVectorField<depos_order_perp,depos_order_para>(
1699 xp_nph, yp_nph, zp_nph,
1700 Bxp, Byp, Bzp,
1701 Bx_arr, By_arr, Bz_arr,
1702 Bx_type, By_type, Bz_type,
1703 dinv, xyzmin, lo, /*n_rz_azimuthal_modes=*/0 );
1704
1705 // Compute inverse Lorentz factor, the average of gamma at time levels n and n+1
1706 const amrex::ParticleReal gaminv = GetImplicitGammaInverse(uxp_n[ip], uyp_n[ip], uzp_n[ip],
1707 uxp_nph[ip], uyp_nph[ip], uzp_nph[ip]);
1708
1709 // Compute current density kernels to deposit
1710 const amrex::Real wq_invvol = qs*wp[ip]*invvol;
1711 const amrex::Real rhop = wq_invvol*gaminv;
1712
1713 // Set the Mass Matrices kernels
1714 amrex::ParticleReal fpxx, fpxy, fpxz;
1715 amrex::ParticleReal fpyx, fpyy, fpyz;
1716 amrex::ParticleReal fpzx, fpzy, fpzz;
1717 setMassMatricesKernels( qs, ms, dt, rhop,
1718 uxp_nph[ip], uyp_nph[ip], uzp_nph[ip],
1719 Bxp, Byp, Bzp,
1720 fpxx, fpxy, fpxz,
1721 fpyx, fpyy, fpyz,
1722 fpzx, fpzy, fpzz );
1723
1724 amrex::ParticleReal const xp_n = (xp_n_data ? xp_n_data[ip] : 0._prt);
1725 amrex::ParticleReal const yp_n = (yp_n_data ? yp_n_data[ip] : 0._prt);
1726 amrex::ParticleReal const zp_n = (zp_n_data ? zp_n_data[ip] : 0._prt);
1727
1728 // Compute position at time n + 1
1729 amrex::ParticleReal const xp_np1 = 2._prt*xp_nph - xp_n;
1730 amrex::ParticleReal const yp_np1 = 2._prt*yp_nph - yp_n;
1731 amrex::ParticleReal const zp_np1 = 2._prt*zp_nph - zp_n;
1732
1733 // Pass dummy arrays for Jx, Jy, Jz (which will not be used)
1734 amrex::Array4<amrex::Real> const dummy_Jx{};
1735 amrex::Array4<amrex::Real> const dummy_Jy{};
1736 amrex::Array4<amrex::Real> const dummy_Jz{};
1737
1738 //NOLINTNEXTLINE(readability-suspicious-call-argument)
1739 doVillasenorJandSigmaDepositionKernel<depos_order,full_mass_matrices,/*deposit_J=*/false>(
1740 xp_n, yp_n, zp_n,
1741 xp_np1, yp_np1, zp_np1,
1742 wq_invvol,
1743 uxp_nph[ip], uyp_nph[ip], uzp_nph[ip],
1744 gaminv,
1745 fpxx, fpxy, fpxz,
1746 fpyx, fpyy, fpyz,
1747 fpzx, fpzy, fpzz,
1748 dummy_Jx, dummy_Jy, dummy_Jz,
1749 max_crossings,
1750 Sxx_arr, Sxy_arr, Sxz_arr,
1751 Syx_arr, Syy_arr, Syz_arr,
1752 Szx_arr, Szy_arr, Szz_arr,
1753 dt, dinv, xyzmin, domain_double, do_cropping, lo );
1754
1755 });
1756}
1757
1758#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
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
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 upx, const amrex::ParticleReal upy, const amrex::ParticleReal upz, const amrex::ParticleReal Bpx, const amrex::ParticleReal Bpy, const amrex::ParticleReal Bpz, 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
doDirectSigmaDeposition
void doDirectSigmaDeposition(const GetParticlePosition< PIdx > &GetPosition, const int *nsuborbits, 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 &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::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 mass matrices for thread thread_num
Definition MassMatricesDeposition.H:522
doVillasenorSigmaDeposition
void doVillasenorSigmaDeposition(const amrex::ParticleReal *xp_n_data, const amrex::ParticleReal *yp_n_data, const amrex::ParticleReal *zp_n_data, const GetParticlePosition< PIdx > &GetPosition, const int *nsuborbits, 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, 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 mass matrices for thread thread_num.
Definition MassMatricesDeposition.H:1639
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:112
RigidAdvanceMode::vz
@ vz
Definition RigidInjectedParticleContainer.H:27
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::IntVectND< 3 >::TheZeroVector
__host__ static __device__ constexpr IntVectND< dim > TheZeroVector() noexcept
amrex::Real
amrex_real Real
amrex::ParticleReal
amrex_particle_real ParticleReal
amrex::Array4
ArrayND< T, 4, true > Array4
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:255
ablastr::constant::SI::inv_c2_v
constexpr auto inv_c2_v
inverse of the square of the vacuum speed of light [s^2/m^2] (variable template)
Definition constant.H:146
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::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