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
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 int col_base = depos_order + aa - iz;
283 int Nc;
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
520template <int depos_order, bool full_mass_matrices>
521void doDirectSigmaDeposition (const GetParticlePosition<PIdx>& GetPosition,
522 const amrex::ParticleReal* wp,
523 const amrex::ParticleReal* uxp_n,
524 const amrex::ParticleReal* uyp_n,
525 const amrex::ParticleReal* uzp_n,
526 const amrex::ParticleReal* uxp_nph,
527 const amrex::ParticleReal* uyp_nph,
528 const amrex::ParticleReal* uzp_nph,
529 amrex::Array4<amrex::Real> const& Sxx_arr,
530 amrex::Array4<amrex::Real> const& Sxy_arr,
531 amrex::Array4<amrex::Real> const& Sxz_arr,
532 amrex::Array4<amrex::Real> const& Syx_arr,
533 amrex::Array4<amrex::Real> const& Syy_arr,
534 amrex::Array4<amrex::Real> const& Syz_arr,
535 amrex::Array4<amrex::Real> const& Szx_arr,
536 amrex::Array4<amrex::Real> const& Szy_arr,
537 amrex::Array4<amrex::Real> const& Szz_arr,
538 const amrex::IntVect& jx_type,
539 const amrex::IntVect& jy_type,
540 const amrex::IntVect& jz_type,
541 const amrex::Array4<amrex::Real const>& Bx_arr,
542 const amrex::Array4<amrex::Real const>& By_arr,
543 const amrex::Array4<amrex::Real const>& Bz_arr,
544 const amrex::IndexType Bx_type,
545 const amrex::IndexType By_type,
546 const amrex::IndexType Bz_type,
547 const long np_to_deposit,
548 const amrex::Real dt,
549 const amrex::XDim3& dinv,
550 const amrex::XDim3& xyzmin,
551 const amrex::Dim3 lo,
552 const amrex::Real qs,
553 const amrex::Real ms)
554{
555 using namespace amrex::literals;
556
557 const amrex::Real invvol = dinv.x*dinv.y*dinv.z;
558
559 // Loop over particles and deposit mass matrices
560 amrex::ParallelFor(
561 np_to_deposit,
562 [=] AMREX_GPU_DEVICE (long ip) {
563 amrex::ParticleReal xp_nph, yp_nph, zp_nph;
564 GetPosition(ip, xp_nph, yp_nph, zp_nph);
565
566 // Compute magnetic field on particle
567 amrex::ParticleReal Bxp = 0.0;
568 amrex::ParticleReal Byp = 0.0;
569 amrex::ParticleReal Bzp = 0.0;
570 const int depos_order_perp = 1;
571 const int depos_order_para = 1;
572 const int n_rz_azimuthal_modes = 0;
573 doDirectGatherVectorField<depos_order_perp,depos_order_para>(
574 xp_nph, yp_nph, zp_nph,
575 Bxp, Byp, Bzp,
576 Bx_arr, By_arr, Bz_arr,
577 Bx_type, By_type, Bz_type,
578 dinv, xyzmin, lo, n_rz_azimuthal_modes );
579
580 // Compute inverse Lorentz factor, the average of gamma at time levels n and n+1
581 const amrex::ParticleReal gaminv = GetImplicitGammaInverse(uxp_n[ip], uyp_n[ip], uzp_n[ip],
582 uxp_nph[ip], uyp_nph[ip], uzp_nph[ip]);
583
584 // Compute current density kernels to deposit
585 const amrex::Real rhop = qs*wp[ip]*invvol*gaminv;
586 const amrex::Real wqx = rhop*uxp_nph[ip];
587 const amrex::Real wqy = rhop*uyp_nph[ip];
588 const amrex::Real wqz = rhop*uzp_nph[ip];
589
590 // Set the Mass Matrices kernels
591 amrex::ParticleReal fpxx, fpxy, fpxz;
592 amrex::ParticleReal fpyx, fpyy, fpyz;
593 amrex::ParticleReal fpzx, fpzy, fpzz;
594 setMassMatricesKernels( qs, ms, dt, rhop,
595 uxp_nph[ip], uyp_nph[ip], uzp_nph[ip],
596 Bxp, Byp, Bzp,
597 fpxx, fpxy, fpxz,
598 fpyx, fpyy, fpyz,
599 fpzx, fpzy, fpzz );
600
601 // Pass dummy arrays for Jx, Jy, Jz (which will not be used)
602 amrex::Array4<amrex::Real> const dummy_Jx{};
603 amrex::Array4<amrex::Real> const dummy_Jy{};
604 amrex::Array4<amrex::Real> const dummy_Jz{};
605
606 doDirectJandSigmaDepositionKernel<depos_order,full_mass_matrices,/*deposit_J=*/false>(
607 xp_nph, yp_nph, zp_nph,
608 wqx, wqy, wqz,
609 fpxx, fpxy, fpxz,
610 fpyx, fpyy, fpyz,
611 fpzx, fpzy, fpzz,
612 dummy_Jx, dummy_Jy, dummy_Jz,
613 Sxx_arr, Sxy_arr, Sxz_arr,
614 Syx_arr, Syy_arr, Syz_arr,
615 Szx_arr, Szy_arr, Szz_arr,
616 jx_type, jy_type, jz_type,
617 dinv, xyzmin, lo );
618
619 }
620 );
621}
622
648template <int depos_order, bool full_mass_matrices, bool deposit_J>
649AMREX_GPU_HOST_DEVICE AMREX_INLINE
650void doVillasenorJandSigmaDepositionKernel ( [[maybe_unused]] const amrex::ParticleReal xp_old,
651 [[maybe_unused]] const amrex::ParticleReal yp_old,
652 [[maybe_unused]] const amrex::ParticleReal zp_old,
653 [[maybe_unused]] const amrex::ParticleReal xp_new,
654 [[maybe_unused]] const amrex::ParticleReal yp_new,
655 [[maybe_unused]] const amrex::ParticleReal zp_new,
656 const amrex::ParticleReal wq_invvol,
657 [[maybe_unused]] const amrex::ParticleReal uxp_mid,
658 [[maybe_unused]] const amrex::ParticleReal uyp_mid,
659 [[maybe_unused]] const amrex::ParticleReal uzp_mid,
660 [[maybe_unused]] const amrex::ParticleReal gaminv,
661 const amrex::ParticleReal fpxx,
662 [[maybe_unused]] const amrex::ParticleReal fpxy,
663 [[maybe_unused]] const amrex::ParticleReal fpxz,
664 [[maybe_unused]] const amrex::ParticleReal fpyx,
665 const amrex::ParticleReal fpyy,
666 [[maybe_unused]] const amrex::ParticleReal fpyz,
667 [[maybe_unused]] const amrex::ParticleReal fpzx,
668 [[maybe_unused]] const amrex::ParticleReal fpzy,
669 const amrex::ParticleReal fpzz,
670 amrex::Array4<amrex::Real> const& Jx_arr,
671 amrex::Array4<amrex::Real> const& Jy_arr,
672 amrex::Array4<amrex::Real> const& Jz_arr,
673 [[maybe_unused]] int max_crossings,
674 amrex::Array4<amrex::Real> const& Sxx_arr,
675 [[maybe_unused]] amrex::Array4<amrex::Real> const& Sxy_arr,
676 [[maybe_unused]] amrex::Array4<amrex::Real> const& Sxz_arr,
677 [[maybe_unused]] amrex::Array4<amrex::Real> const& Syx_arr,
678 amrex::Array4<amrex::Real> const& Syy_arr,
679 [[maybe_unused]] amrex::Array4<amrex::Real> const& Syz_arr,
680 [[maybe_unused]] amrex::Array4<amrex::Real> const& Szx_arr,
681 [[maybe_unused]] amrex::Array4<amrex::Real> const& Szy_arr,
682 amrex::Array4<amrex::Real> const& Szz_arr,
683 const amrex::Real dt,
684 const amrex::XDim3& dinv,
685 const amrex::XDim3& xyzmin,
686 const amrex::GpuArray<amrex::GpuArray<double,2>, AMREX_SPACEDIM> & domain_double,
687 const amrex::GpuArray<amrex::GpuArray<bool,2>, AMREX_SPACEDIM> & do_cropping,
688 const amrex::Dim3 lo )
689{
690
691 using namespace amrex::literals;
692
693#if (AMREX_SPACEDIM > 1)
694 amrex::Real constexpr one_third = 1.0_rt / 3.0_rt;
695 amrex::Real constexpr one_sixth = 1.0_rt / 6.0_rt;
696#endif
697
698 // computes current and old position in grid units
699#if defined(WARPX_DIM_RZ) || defined(WARPX_DIM_RCYLINDER)
700 amrex::Real const xp_mid = (xp_new + xp_old)*0.5_rt;
701 amrex::Real const yp_mid = (yp_new + yp_old)*0.5_rt;
702 amrex::Real const rp_new = std::sqrt(xp_new*xp_new + yp_new*yp_new);
703 amrex::Real const rp_old = std::sqrt(xp_old*xp_old + yp_old*yp_old);
704 amrex::Real const rp_mid = (rp_new + rp_old)/2._rt;
705 amrex::Real const costheta_mid = (rp_mid > 0._rt ? xp_mid/rp_mid : 1._rt);
706 amrex::Real const sintheta_mid = (rp_mid > 0._rt ? yp_mid/rp_mid : 0._rt);
707
708 // Keep these double to avoid bug in single precision
709 double x_new = (rp_new - xyzmin.x)*dinv.x;
710 double const x_old = (rp_old - xyzmin.x)*dinv.x;
711 amrex::Real const vx = (rp_new - rp_old)/dt;
712 amrex::Real const vy = (-uxp_mid*sintheta_mid + uyp_mid*costheta_mid)*gaminv;
713#if defined(WARPX_DIM_RCYLINDER)
714 amrex::Real const vz = uzp_mid*gaminv;
715#endif
716#elif defined(WARPX_DIM_RSPHERE)
717 amrex::Real const xp_mid = (xp_new + xp_old)*0.5_rt;
718 amrex::Real const yp_mid = (yp_new + yp_old)*0.5_rt;
719 amrex::Real const zp_mid = (zp_new + zp_old)*0.5_rt;
720 amrex::Real const rpxy_new = std::sqrt(xp_new*xp_new + yp_new*yp_new);
721 amrex::Real const rp_new = std::sqrt(xp_new*xp_new + yp_new*yp_new + zp_new*zp_new);
722 amrex::Real const rpxy_old = std::sqrt(xp_old*xp_old + yp_old*yp_old);
723 amrex::Real const rp_old = std::sqrt(xp_old*xp_old + yp_old*yp_old + zp_old*zp_old);
724 amrex::Real const rpxy_mid = (rpxy_new + rpxy_old)*0.5_rt;
725 amrex::Real const rp_mid = (rp_new + rp_old)*0.5_rt;
726 amrex::Real const costheta_mid = (rpxy_mid > 0._rt ? xp_mid/rpxy_mid : 1._rt);
727 amrex::Real const sintheta_mid = (rpxy_mid > 0._rt ? yp_mid/rpxy_mid : 0._rt);
728 amrex::Real const cosphi_mid = (rp_mid > 0._rt ? rpxy_mid/rp_mid : 1._rt);
729 amrex::Real const sinphi_mid = (rp_mid > 0._rt ? zp_mid/rp_mid : 0._rt);
730
731 // Keep these double to avoid bug in single precision
732 double x_new = (rp_new - xyzmin.x)*dinv.x;
733 double const x_old = (rp_old - xyzmin.x)*dinv.x;
734 amrex::Real const vx = (rp_new - rp_old)/dt;
735 amrex::Real const vy = (-uxp_mid*sintheta_mid + uyp_mid*costheta_mid)*gaminv;
736 amrex::Real const vz = (-uxp_mid*costheta_mid*cosphi_mid - uyp_mid*sintheta_mid*cosphi_mid + uzp_mid*sinphi_mid)*gaminv;
737#elif defined(WARPX_DIM_XZ)
738 // Keep these double to avoid bug in single precision
739 double x_new = (xp_new - xyzmin.x)*dinv.x;
740 double const x_old = (xp_old - xyzmin.x)*dinv.x;
741 amrex::Real const vx = (xp_new - xp_old)/dt;
742 amrex::Real const vy = uyp_mid*gaminv;
743#elif defined(WARPX_DIM_1D_Z)
744 amrex::Real const vx = uxp_mid*gaminv;
745 amrex::Real const vy = uyp_mid*gaminv;
746#elif defined(WARPX_DIM_3D)
747 // Keep these double to avoid bug in single precision
748 double x_new = (xp_new - xyzmin.x)*dinv.x;
749 double const x_old = (xp_old - xyzmin.x)*dinv.x;
750 double y_new = (yp_new - xyzmin.y)*dinv.y;
751 double const y_old = (yp_old - xyzmin.y)*dinv.y;
752 amrex::Real const vx = (xp_new - xp_old)/dt;
753 amrex::Real const vy = (yp_new - yp_old)/dt;
754#endif
755
756#if !defined(WARPX_DIM_RCYLINDER) && !defined(WARPX_DIM_RSPHERE)
757 // Keep these double to avoid bug in single precision
758 double z_new = (zp_new - xyzmin.z)*dinv.z;
759 double const z_old = (zp_old - xyzmin.z)*dinv.z;
760 amrex::Real const vz = (zp_new - zp_old)/dt;
761#endif
762
763 // Define velocity kernels to deposit
764 amrex::Real const wqx = wq_invvol*vx;
765 amrex::Real const wqy = wq_invvol*vy;
766 amrex::Real const wqz = wq_invvol*vz;
767
768 // Compute total change in particle position (always do before cropping)
769#if !defined(WARPX_DIM_1D_Z)
770 const double dxp = x_new - x_old;
771#endif
772#if defined(WARPX_DIM_3D)
773 const double dyp = y_new - y_old;
774#endif
775#if !defined(WARPX_DIM_RCYLINDER) && !defined(WARPX_DIM_RSPHERE)
776 const double dzp = z_new - z_old;
777#endif
778
779 // Crop particle orbits at absorbing domain boundaries
780#if defined(WARPX_DIM_3D)
781 ParticleUtils::crop_at_boundary(x_old, y_old, z_old, x_new, y_new, z_new,
782 domain_double[0][0], domain_double[0][1], do_cropping[0]);
783 ParticleUtils::crop_at_boundary(y_old, z_old, x_old, y_new, z_new, x_new,
784 domain_double[1][0], domain_double[1][1], do_cropping[1]);
785 ParticleUtils::crop_at_boundary(z_old, x_old, y_old, z_new, x_new, y_new,
786 domain_double[2][0], domain_double[2][1], do_cropping[2]);
787#elif defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ)
788 ParticleUtils::crop_at_boundary(x_old, z_old, x_new, z_new,
789 domain_double[0][0], domain_double[0][1], do_cropping[0]);
790 ParticleUtils::crop_at_boundary(z_old, x_old, z_new, x_new,
791 domain_double[1][0], domain_double[1][1], do_cropping[1]);
792#elif defined(WARPX_DIM_1D_Z)
793 ParticleUtils::crop_at_boundary(z_new, domain_double[0][0], domain_double[0][1], do_cropping[0]);
794#elif defined(WARPX_DIM_RCYLINDER) || defined(WARPX_DIM_RSPHERE)
795 ParticleUtils::crop_at_boundary(x_new, domain_double[0][0], domain_double[0][1], do_cropping[0]);
796#endif
797
798 // 1) Determine the number of segments.
799 // 2) Loop over segments and deposit current.
800
801 // cell crossings are defined at cell edges if depos_order is odd
802 // cell crossings are defined at cell centers if depos_order is even
803
804 int num_segments = 1;
805 double shift = 0.0;
806 if ( (depos_order % 2) == 0 ) { shift = 0.5; }
807
808#if defined(WARPX_DIM_3D)
809
810 // compute cell crossings in X-direction
811 const auto i_old = static_cast<int>(x_old-shift);
812 const auto i_new = static_cast<int>(x_new-shift);
813 const int cell_crossings_x = std::abs(i_new-i_old);
814 num_segments += cell_crossings_x;
815
816 // compute cell crossings in Y-direction
817 const auto j_old = static_cast<int>(y_old-shift);
818 const auto j_new = static_cast<int>(y_new-shift);
819 const int cell_crossings_y = std::abs(j_new-j_old);
820 num_segments += cell_crossings_y;
821
822 // compute cell crossings in Z-direction
823 const auto k_old = static_cast<int>(z_old-shift);
824 const auto k_new = static_cast<int>(z_new-shift);
825 const int cell_crossings_z = std::abs(k_new-k_old);
826 num_segments += cell_crossings_z;
827
828 // Compute initial particle cell locations in each direction
829 // used to find the position at cell crossings.
830 // Keep these double to avoid bug in single precision
831 const auto dirX_sign = static_cast<double>(dxp < 0. ? -1. : 1.);
832 const auto dirY_sign = static_cast<double>(dyp < 0. ? -1. : 1.);
833 const auto dirZ_sign = static_cast<double>(dzp < 0. ? -1. : 1.);
834 double Xcell = 0., Ycell = 0., Zcell = 0.;
835 if (num_segments > 1) {
836 Xcell = static_cast<double>(i_old) + shift + 0.5*(1.-dirX_sign);
837 Ycell = static_cast<double>(j_old) + shift + 0.5*(1.-dirY_sign);
838 Zcell = static_cast<double>(k_old) + shift + 0.5*(1.-dirZ_sign);
839 }
840
841 // loop over the number of segments and deposit
842 const Compute_shape_factor< depos_order-1 > compute_shape_factor_cell;
843 const Compute_shape_factor_pair< depos_order > compute_shape_factors_node;
844 double dxp_seg, dyp_seg, dzp_seg;
845 double x0_new, y0_new, z0_new;
846 double x0_old = x_old;
847 double y0_old = y_old;
848 double z0_old = z_old;
849
850 for (int ns=0; ns<num_segments; ns++) {
851
852 if (ns == num_segments-1) { // final segment
853
854 x0_new = x_new;
855 y0_new = y_new;
856 z0_new = z_new;
857 dxp_seg = x0_new - x0_old;
858 dyp_seg = y0_new - y0_old;
859 dzp_seg = z0_new - z0_old;
860
861 }
862 else {
863
864 x0_new = Xcell + dirX_sign;
865 y0_new = Ycell + dirY_sign;
866 z0_new = Zcell + dirZ_sign;
867 dxp_seg = x0_new - x0_old;
868 dyp_seg = y0_new - y0_old;
869 dzp_seg = z0_new - z0_old;
870
871 if ( (dyp == 0. || std::abs(dxp_seg) < std::abs(dxp/dyp*dyp_seg))
872 && (dzp == 0. || std::abs(dxp_seg) < std::abs(dxp/dzp*dzp_seg)) ) {
873 Xcell = x0_new;
874 dyp_seg = dyp/dxp*dxp_seg;
875 dzp_seg = dzp/dxp*dxp_seg;
876 y0_new = y0_old + dyp_seg;
877 z0_new = z0_old + dzp_seg;
878 }
879 else if (dzp == 0. || std::abs(dyp_seg) < std::abs(dyp/dzp*dzp_seg)) {
880 Ycell = y0_new;
881 dxp_seg = dxp/dyp*dyp_seg;
882 dzp_seg = dzp/dyp*dyp_seg;
883 x0_new = x0_old + dxp_seg;
884 z0_new = z0_old + dzp_seg;
885 }
886 else {
887 Zcell = z0_new;
888 dxp_seg = dxp/dzp*dzp_seg;
889 dyp_seg = dyp/dzp*dzp_seg;
890 x0_new = x0_old + dxp_seg;
891 y0_new = y0_old + dyp_seg;
892 }
893
894 }
895
896 // Compute the segment factors (each equal to dt_seg/dt for nonzero dxp, dyp, or dzp)
897 const auto seg_factor_x = static_cast<amrex::Real>(dxp == 0. ? 1._rt : dxp_seg/dxp);
898 const auto seg_factor_y = static_cast<amrex::Real>(dyp == 0. ? 1._rt : dyp_seg/dyp);
899 const auto seg_factor_z = static_cast<amrex::Real>(dzp == 0. ? 1._rt : dzp_seg/dzp);
900
901 // Compute cell-based weights using the average segment position
902 // Keep these double to avoid bug in single precision
903 double sx_cell[depos_order] = {0.};
904 double sy_cell[depos_order] = {0.};
905 double sz_cell[depos_order] = {0.};
906 double const x0_bar = (x0_new + x0_old)/2.0;
907 double const y0_bar = (y0_new + y0_old)/2.0;
908 double const z0_bar = (z0_new + z0_old)/2.0;
909 const int i0_cell = compute_shape_factor_cell( sx_cell, x0_bar-0.5 );
910 const int j0_cell = compute_shape_factor_cell( sy_cell, y0_bar-0.5 );
911 const int k0_cell = compute_shape_factor_cell( sz_cell, z0_bar-0.5 );
912
913 if constexpr (depos_order >= 3) { // higher-order correction to the cell-based weights
914 const Compute_shape_factor_pair<depos_order-1> compute_shape_factors_cell;
915 double sx_old_cell[depos_order] = {0.};
916 double sx_new_cell[depos_order] = {0.};
917 double sy_old_cell[depos_order] = {0.};
918 double sy_new_cell[depos_order] = {0.};
919 double sz_old_cell[depos_order] = {0.};
920 double sz_new_cell[depos_order] = {0.};
921 const int i0_cell_2 = compute_shape_factors_cell( sx_old_cell, sx_new_cell, x0_old-0.5, x0_new-0.5 );
922 const int j0_cell_2 = compute_shape_factors_cell( sy_old_cell, sy_new_cell, y0_old-0.5, y0_new-0.5 );
923 const int k0_cell_2 = compute_shape_factors_cell( sz_old_cell, sz_new_cell, z0_old-0.5, z0_new-0.5 );
924 amrex::ignore_unused(i0_cell_2, j0_cell_2, k0_cell_2);
925 for (int m=0; m<depos_order; m++) {
926 sx_cell[m] = (4.0*sx_cell[m] + sx_old_cell[m] + sx_new_cell[m])/6.0;
927 sy_cell[m] = (4.0*sy_cell[m] + sy_old_cell[m] + sy_new_cell[m])/6.0;
928 sz_cell[m] = (4.0*sz_cell[m] + sz_old_cell[m] + sz_new_cell[m])/6.0;
929 }
930 }
931
932 // Compute node-based weights using the old and new segment positions
933 // Keep these double to avoid bug in single precision
934 double sx_old_node[depos_order+1] = {0.};
935 double sx_new_node[depos_order+1] = {0.};
936 double sy_old_node[depos_order+1] = {0.};
937 double sy_new_node[depos_order+1] = {0.};
938 double sz_old_node[depos_order+1] = {0.};
939 double sz_new_node[depos_order+1] = {0.};
940 const int i0_node = compute_shape_factors_node( sx_old_node, sx_new_node, x0_old, x0_new );
941 const int j0_node = compute_shape_factors_node( sy_old_node, sy_new_node, y0_old, y0_new );
942 const int k0_node = compute_shape_factors_node( sz_old_node, sz_new_node, z0_old, z0_new );
943
944 // deposit Jx and Sxx for this segment
945 amrex::Real weight;
946 for (int i=0; i<=depos_order-1; i++) {
947 for (int j=0; j<=depos_order; j++) {
948 for (int k=0; k<=depos_order; k++) {
949 weight = sx_cell[i]*( sy_old_node[j]*sz_old_node[k]*one_third
950 + sy_old_node[j]*sz_new_node[k]*one_sixth
951 + sy_new_node[j]*sz_old_node[k]*one_sixth
952 + sy_new_node[j]*sz_new_node[k]*one_third )*seg_factor_x;
953 if constexpr (deposit_J) {
954 amrex::Gpu::Atomic::AddNoRet( &Jx_arr(lo.x+i0_cell+i, lo.y+j0_node+j, lo.z+k0_node+k), wqx*weight);
955 }
956 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);
957 }
958 }
959 }
960
961 // deposit Jy and Syy or this segment
962 for (int i=0; i<=depos_order; i++) {
963 for (int j=0; j<=depos_order-1; j++) {
964 for (int k=0; k<=depos_order; k++) {
965 weight = sy_cell[j]*( sx_old_node[i]*sz_old_node[k]*one_third
966 + sx_old_node[i]*sz_new_node[k]*one_sixth
967 + sx_new_node[i]*sz_old_node[k]*one_sixth
968 + sx_new_node[i]*sz_new_node[k]*one_third )*seg_factor_y;
969 if constexpr (deposit_J) {
970 amrex::Gpu::Atomic::AddNoRet( &Jy_arr(lo.x+i0_node+i, lo.y+j0_cell+j, lo.z+k0_node+k), wqy*weight);
971 }
972 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);
973 }
974 }
975 }
976
977 // deposit Jz and Sz for this segment
978 for (int i=0; i<=depos_order; i++) {
979 for (int j=0; j<=depos_order; j++) {
980 for (int k=0; k<=depos_order-1; k++) {
981 weight = sz_cell[k]*( sx_old_node[i]*sy_old_node[j]*one_third
982 + sx_old_node[i]*sy_new_node[j]*one_sixth
983 + sx_new_node[i]*sy_old_node[j]*one_sixth
984 + sx_new_node[i]*sy_new_node[j]*one_third )*seg_factor_z;
985 if constexpr (deposit_J) {
986 amrex::Gpu::Atomic::AddNoRet( &Jz_arr(lo.x+i0_node+i, lo.y+j0_node+j, lo.z+k0_cell+k), wqz*weight);
987 }
988 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);
989 }
990 }
991 }
992
993 // update old segment values
994 if (ns < num_segments-1) {
995 x0_old = x0_new;
996 y0_old = y0_new;
997 z0_old = z0_new;
998 }
999
1000 } // end loop over segments
1001
1002#elif defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ)
1003
1004 // compute cell crossings in X-direction
1005 const auto i_old = static_cast<int>(x_old-shift);
1006 const auto i_new = static_cast<int>(x_new-shift);
1007 const int cell_crossings_x = std::abs(i_new-i_old);
1008 num_segments += cell_crossings_x;
1009
1010 // compute cell crossings in Z-direction
1011 const auto k_old = static_cast<int>(z_old-shift);
1012 const auto k_new = static_cast<int>(z_new-shift);
1013 const int cell_crossings_z = std::abs(k_new-k_old);
1014 num_segments += cell_crossings_z;
1015
1016 // Compute initial particle cell locations in each direction
1017 // used to find the position at cell crossings.
1018 // Keep these double to avoid bug in single precision
1019 const auto dirX_sign = static_cast<double>(dxp < 0. ? -1. : 1.);
1020 const auto dirZ_sign = static_cast<double>(dzp < 0. ? -1. : 1.);
1021 double Xcell = 0., Zcell = 0.;
1022 if (num_segments > 1) {
1023 Xcell = static_cast<double>(i_old) + shift + 0.5*(1.-dirX_sign);
1024 Zcell = static_cast<double>(k_old) + shift + 0.5*(1.-dirZ_sign);
1025 }
1026
1027 // loop over the number of segments and deposit
1028 const Compute_shape_factor< depos_order-1 > compute_shape_factor_cell;
1029 const Compute_shape_factor_pair< depos_order > compute_shape_factors_node;
1030 double dxp_seg, dzp_seg;
1031 double x0_new, z0_new;
1032 double x0_old = x_old;
1033 double z0_old = z_old;
1034
1035 constexpr int num_segments_max = 1 + 4*AMREX_SPACEDIM;
1036 AMREX_ALWAYS_ASSERT_WITH_MESSAGE( num_segments <= num_segments_max,
1037 "Error: num_segments must be less than or equal to 1 + 4*AMREX_SPACEDIM.");
1038
1039 // Save the start index and interpolation weights for each segment
1040 int i0_cell[num_segments_max];
1041 int i0_node[num_segments_max];
1042 int k0_cell[num_segments_max];
1043 int k0_node[num_segments_max];
1044 amrex::Real weight_cellX_nodeZ[num_segments_max][depos_order][depos_order+1];
1045 amrex::Real weight_nodeX_cellZ[num_segments_max][depos_order+1][depos_order];
1046 amrex::Real weight_nodeX_nodeZ[num_segments_max][depos_order+1][depos_order+1];
1047
1048 const auto i_mid = static_cast<int>(0.5*(x_new+x_old)-shift);
1049 const auto k_mid = static_cast<int>(0.5*(z_new+z_old)-shift);
1050 int SegNumX[num_segments_max];
1051 int SegNumZ[num_segments_max];
1052
1053 for (int ns=0; ns<num_segments; ns++) {
1054
1055 if (ns == num_segments-1) { // final segment
1056
1057 x0_new = x_new;
1058 z0_new = z_new;
1059 dxp_seg = x0_new - x0_old;
1060 dzp_seg = z0_new - z0_old;
1061
1062 }
1063 else {
1064
1065 x0_new = Xcell + dirX_sign;
1066 z0_new = Zcell + dirZ_sign;
1067 dxp_seg = x0_new - x0_old;
1068 dzp_seg = z0_new - z0_old;
1069
1070 if (dzp == 0. || std::abs(dxp_seg) < std::abs(dxp/dzp*dzp_seg)) {
1071 Xcell = x0_new;
1072 dzp_seg = dzp/dxp*dxp_seg;
1073 z0_new = z0_old + dzp_seg;
1074 }
1075 else {
1076 Zcell = z0_new;
1077 dxp_seg = dxp/dzp*dzp_seg;
1078 x0_new = x0_old + dxp_seg;
1079 }
1080
1081 }
1082
1083 // Compute the segment factors (each equal to dt_seg/dt for nonzero dxp, or dzp)
1084 const auto seg_factor_x = static_cast<amrex::Real>(dxp == 0. ? 1._rt : dxp_seg/dxp);
1085 const auto seg_factor_z = static_cast<amrex::Real>(dzp == 0. ? 1._rt : dzp_seg/dzp);
1086
1087 // Compute cell-based weights using the average segment position
1088 // Keep these double to avoid bug in single precision
1089 double sx_cell[depos_order] = {0.};
1090 double sz_cell[depos_order] = {0.};
1091 double const x0_bar = (x0_new + x0_old)/2.0;
1092 double const z0_bar = (z0_new + z0_old)/2.0;
1093 i0_cell[ns] = compute_shape_factor_cell( sx_cell, x0_bar-0.5 );
1094 k0_cell[ns] = compute_shape_factor_cell( sz_cell, z0_bar-0.5 );
1095
1096 // Set the segment number for the mass matrix component calc
1097 if constexpr (full_mass_matrices) {
1098 const auto i0_mid = static_cast<int>(x0_bar-shift);
1099 const auto k0_mid = static_cast<int>(z0_bar-shift);
1100 SegNumX[ns] = 1 + i0_mid - i_mid;
1101 SegNumZ[ns] = 1 + k0_mid - k_mid;
1102 }
1103
1104 if constexpr (depos_order >= 3) { // higher-order correction to the cell-based weights
1105 const Compute_shape_factor_pair<depos_order-1> compute_shape_factors_cell;
1106 double sx_old_cell[depos_order] = {0.};
1107 double sx_new_cell[depos_order] = {0.};
1108 double sz_old_cell[depos_order] = {0.};
1109 double sz_new_cell[depos_order] = {0.};
1110 const int i0_cell_2 = compute_shape_factors_cell( sx_old_cell, sx_new_cell, x0_old-0.5, x0_new-0.5 );
1111 const int k0_cell_2 = compute_shape_factors_cell( sz_old_cell, sz_new_cell, z0_old-0.5, z0_new-0.5 );
1112 amrex::ignore_unused(i0_cell_2, k0_cell_2);
1113 for (int m=0; m<depos_order; m++) {
1114 sx_cell[m] = (4.0*sx_cell[m] + sx_old_cell[m] + sx_new_cell[m])/6.0;
1115 sz_cell[m] = (4.0*sz_cell[m] + sz_old_cell[m] + sz_new_cell[m])/6.0;
1116 }
1117 }
1118
1119 // Compute node-based weights using the old and new segment positions
1120 // Keep these double to avoid bug in single precision
1121 double sx_old_node[depos_order+1] = {0.};
1122 double sx_new_node[depos_order+1] = {0.};
1123 double sz_old_node[depos_order+1] = {0.};
1124 double sz_new_node[depos_order+1] = {0.};
1125 i0_node[ns] = compute_shape_factors_node( sx_old_node, sx_new_node, x0_old, x0_new );
1126 k0_node[ns] = compute_shape_factors_node( sz_old_node, sz_new_node, z0_old, z0_new );
1127
1128 // deposit Jx and Sx for this segment
1129 amrex::Real weight;
1130 for (int i=0; i<=depos_order-1; i++) {
1131 for (int k=0; k<=depos_order; k++) {
1132 const int i_J = lo.x + i0_cell[ns] + i;
1133 const int k_J = lo.y + k0_node[ns] + k;
1134 weight = sx_cell[i]*(sz_old_node[k] + sz_new_node[k])/2.0_rt*seg_factor_x;
1135 if constexpr (deposit_J) {
1136 amrex::Gpu::Atomic::AddNoRet(&Jx_arr(i_J, k_J, 0, 0), wqx*weight);
1137 }
1138 if constexpr (full_mass_matrices) { weight_cellX_nodeZ[ns][i][k] = weight; }
1139 else {
1140 amrex::Gpu::Atomic::AddNoRet(&Sxx_arr(i_J, k_J, 0, 0), fpxx*weight*weight);
1141 }
1142 }
1143 }
1144
1145 // deposit out-of-plane Jy and Sy for this segment
1146 const auto seg_factor_y = std::min(seg_factor_x,seg_factor_z);
1147 for (int i=0; i<=depos_order; i++) {
1148 for (int k=0; k<=depos_order; k++) {
1149 const int i_J = lo.x + i0_node[ns] + i;
1150 const int k_J = lo.y + k0_node[ns] + k;
1151 weight = ( sx_old_node[i]*sz_old_node[k]*one_third
1152 + sx_old_node[i]*sz_new_node[k]*one_sixth
1153 + sx_new_node[i]*sz_old_node[k]*one_sixth
1154 + sx_new_node[i]*sz_new_node[k]*one_third )*seg_factor_y;
1155 if constexpr (deposit_J) {
1156 amrex::Gpu::Atomic::AddNoRet(&Jy_arr(i_J, k_J, 0, 0), wqy*weight);
1157 }
1158 if constexpr (full_mass_matrices) { weight_nodeX_nodeZ[ns][i][k] = weight; }
1159 else {
1160 amrex::Gpu::Atomic::AddNoRet(&Syy_arr(i_J, k_J, 0, 0), fpyy*weight*weight);
1161 }
1162 }
1163 }
1164
1165 // deposit Jz and Szz for this segment
1166 for (int i=0; i<=depos_order; i++) {
1167 for (int k=0; k<=depos_order-1; k++) {
1168 const int i_J = lo.x + i0_node[ns] + i;
1169 const int k_J = lo.y + k0_cell[ns] + k;
1170 weight = sz_cell[k]*(sx_old_node[i] + sx_new_node[i])/2.0_rt*seg_factor_z;
1171 if constexpr (deposit_J) {
1172 amrex::Gpu::Atomic::AddNoRet(&Jz_arr(i_J, k_J, 0, 0), wqz*weight);
1173 }
1174 if constexpr (full_mass_matrices) { weight_nodeX_cellZ[ns][i][k] = weight; }
1175 else {
1176 amrex::Gpu::Atomic::AddNoRet(&Szz_arr(i_J, k_J, 0, 0), fpzz*weight*weight);
1177 }
1178 }
1179 }
1180
1181 // update old segment values
1182 if (ns < num_segments-1) {
1183 x0_old = x0_new;
1184 z0_old = z0_new;
1185 }
1186
1187 } // end loop over segments
1188
1189 if constexpr (full_mass_matrices) {
1190
1191 const int Ncomp_base = 2*depos_order + 2*max_crossings;
1192 const int Ncomp_xx0 = Ncomp_base - 1;
1193 const int Ncomp_xz0 = Ncomp_base;
1194 const int Ncomp_zx0 = Ncomp_base;
1195 const int Ncomp_zz0 = 1 + Ncomp_base;
1196 const int Ncomp_xy0 = Ncomp_base;
1197 const int Ncomp_yx0 = Ncomp_base;
1198 const int Ncomp_yz0 = 1 + Ncomp_base;
1199 const int Ncomp_yy0 = 1 + Ncomp_base;
1200 const int Ncomp_zy0 = 1 + Ncomp_base;
1201
1202 const int width_xx1 = depos_order + max_crossings; // (Ncomp_xx[1] - 1)/2
1203 const int width_zz1 = width_xx1 - 1; // (Ncomp_zz[1] - 1)/2
1204 const int width_yy1 = width_xx1; // (Ncomp_yy[1] - 1)/2
1205
1206 // Loop over segments and deposit full mass matrices
1207 for (int ns=0; ns<num_segments; ns++) {
1208
1209 // Deposit Sxx, Sxz, and Sxy for this segment
1210 for (int i = 0; i <= depos_order - 1; i++) {
1211 for (int k = 0; k <= depos_order; k++) {
1212 const int i_J = lo.x + i0_cell[ns] + i;
1213 const int k_J = lo.y + k0_node[ns] + k;
1214 const amrex::Real weight_J = weight_cellX_nodeZ[ns][i][k];
1215 for (int ms=0; ms<num_segments; ms++) {
1216 const int SegShiftX = max_crossings + SegNumX[ms] - SegNumX[ns];
1217 const int SegShiftZ = max_crossings + SegNumZ[ms] - SegNumZ[ns];
1218 // Deposit Sxx
1219 for (int kE = 0; kE <= depos_order; kE++) {
1220 const int row_xx = depos_order - k + kE + SegShiftZ;
1221 const int above_diag = (row_xx > width_xx1) ? 1 : 0;
1222 for (int iE = 0; iE <= depos_order - 1; iE++) {
1223 const int col_xx = depos_order - 1 - i + iE + SegShiftX;
1224 if (col_xx > Ncomp_xx0 - row_xx - above_diag) { break; } // Reduced deposit for diagonal mass matrices
1225 const int comp_xx = col_xx + Ncomp_xx0*row_xx;
1226 const amrex::Real weight_E = weight_cellX_nodeZ[ms][iE][kE];
1227 amrex::Gpu::Atomic::AddNoRet(&Sxx_arr(i_J, k_J, 0, comp_xx), fpxx*weight_J*weight_E);
1228 }
1229 }
1230 // Deposit Sxz
1231 for (int iE = 0; iE <= depos_order; iE++) {
1232 for (int kE = 0; kE <= depos_order - 1; kE++) {
1233 const amrex::Real weight_E = weight_nodeX_cellZ[ms][iE][kE];
1234 const int comp_xz = depos_order - 1 - i + iE + SegShiftX
1235 + Ncomp_xz0*(depos_order - k + kE + SegShiftZ);
1236 amrex::Gpu::Atomic::AddNoRet(&Sxz_arr(i_J, k_J, 0, comp_xz), fpxz*weight_J*weight_E);
1237 }
1238 }
1239 // Deposit Sxy
1240 for (int iE = 0; iE <= depos_order; iE++) {
1241 for (int kE = 0; kE <= depos_order; kE++) {
1242 const amrex::Real weight_E = weight_nodeX_nodeZ[ms][iE][kE];
1243 const int comp_xy = depos_order - 1 - i + iE + SegShiftX
1244 + Ncomp_xy0*(depos_order - k + kE + SegShiftZ);
1245 amrex::Gpu::Atomic::AddNoRet(&Sxy_arr(i_J, k_J, 0, comp_xy), fpxy*weight_J*weight_E);
1246 }
1247 }
1248 }
1249 }
1250 }
1251
1252 // Deposit Szx, Szz, and Szy for this segment
1253 for (int i=0; i<=depos_order; i++) {
1254 for (int k=0; k<=depos_order-1; k++) {
1255 const int i_J = lo.x + i0_node[ns] + i;
1256 const int k_J = lo.y + k0_cell[ns] + k;
1257 const amrex::Real weight_J = weight_nodeX_cellZ[ns][i][k];
1258 for (int ms=0; ms<num_segments; ms++) {
1259 const int SegShiftX = max_crossings + SegNumX[ms] - SegNumX[ns];
1260 const int SegShiftZ = max_crossings + SegNumZ[ms] - SegNumZ[ns];
1261 // Deposit Szx
1262 for (int iE=0; iE<=depos_order-1; iE++) {
1263 for (int kE=0; kE<=depos_order; kE++) {
1264 const amrex::Real weight_E = weight_cellX_nodeZ[ms][iE][kE];
1265 const int comp_zx = depos_order - i + iE + SegShiftX
1266 + Ncomp_zx0*(depos_order-1 - k + kE + SegShiftZ);
1267 amrex::Gpu::Atomic::AddNoRet( &Szx_arr(i_J, k_J, 0, comp_zx), fpzx*weight_J*weight_E);
1268 }
1269 }
1270 // Deposit Szz
1271 for (int kE=0; kE<=depos_order-1; kE++) {
1272 const int row_zz = depos_order - 1 - k + kE + SegShiftZ;
1273 const int above_diag = (row_zz > width_zz1) ? 1 : 0;
1274 for (int iE=0; iE<=depos_order; iE++) {
1275 const int col_zz = depos_order - i + iE + SegShiftX;
1276 if (col_zz > Ncomp_zz0 - 2 - row_zz - above_diag) { break; } // Reduced deposit for diagonal mass matrices
1277 const int comp_zz = col_zz + Ncomp_zz0*row_zz;
1278 const amrex::Real weight_E = weight_nodeX_cellZ[ms][iE][kE];
1279 amrex::Gpu::Atomic::AddNoRet( &Szz_arr(i_J, k_J, 0, comp_zz), fpzz*weight_J*weight_E);
1280 }
1281 }
1282 // Deposit Szy
1283 for (int iE=0; iE<=depos_order; iE++) {
1284 for (int kE=0; kE<=depos_order; kE++) {
1285 const amrex::Real weight_E = weight_nodeX_nodeZ[ms][iE][kE];
1286 const int comp_zy = depos_order - i + iE + SegShiftX
1287 + Ncomp_zy0*(depos_order-1 - k + kE + SegShiftZ);
1288 amrex::Gpu::Atomic::AddNoRet( &Szy_arr(i_J, k_J, 0, comp_zy), fpzy*weight_J*weight_E);
1289 }
1290 }
1291 }
1292 }
1293 }
1294
1295 // Deposit Syx, Syz, and Syy for this segment
1296 for (int i=0; i<=depos_order; i++) {
1297 for (int k=0; k<=depos_order; k++) {
1298 const int i_J = lo.x + i0_node[ns] + i;
1299 const int k_J = lo.y + k0_node[ns] + k;
1300 const amrex::Real weight_J = weight_nodeX_nodeZ[ns][i][k];
1301 for (int ms=0; ms<num_segments; ms++) {
1302 const int SegShiftX = max_crossings + SegNumX[ms] - SegNumX[ns];
1303 const int SegShiftZ = max_crossings + SegNumZ[ms] - SegNumZ[ns];
1304 // Deposit Syx
1305 for (int iE=0; iE<=depos_order-1; iE++) {
1306 for (int kE=0; kE<=depos_order; kE++) {
1307 const amrex::Real weight_E = weight_cellX_nodeZ[ms][iE][kE];
1308 const int comp_yx = depos_order - i + iE + SegShiftX
1309 + Ncomp_yx0*(depos_order - k + kE + SegShiftZ);
1310 amrex::Gpu::Atomic::AddNoRet( &Syx_arr(i_J, k_J, 0, comp_yx), fpyx*weight_J*weight_E);
1311 }
1312 }
1313 // Deposit Syz
1314 for (int iE=0; iE<=depos_order; iE++) {
1315 for (int kE=0; kE<=depos_order-1; kE++) {
1316 const amrex::Real weight_E = weight_nodeX_cellZ[ms][iE][kE];
1317 const int comp_yz = depos_order - i + iE + SegShiftX
1318 + Ncomp_yz0*(depos_order - k + kE + SegShiftZ);
1319 amrex::Gpu::Atomic::AddNoRet( &Syz_arr(i_J, k_J, 0, comp_yz), fpyz*weight_J*weight_E);
1320 }
1321 }
1322 // Deposit Syy
1323 for (int kE=0; kE<=depos_order; kE++) {
1324 const int row_yy = depos_order - k + kE + SegShiftZ;
1325 const int above_diag = (row_yy > width_yy1) ? 1 : 0;
1326 for (int iE=0; iE<=depos_order; iE++) {
1327 const int col_yy = depos_order - i + iE + SegShiftX;
1328 if (col_yy > Ncomp_yy0 - 1 - row_yy - above_diag) { break; } // Reduced deposit for diagonal mass matrices
1329 const int comp_yy = col_yy + Ncomp_yy0*row_yy;
1330 const amrex::Real weight_E = weight_nodeX_nodeZ[ms][iE][kE];
1331 amrex::Gpu::Atomic::AddNoRet( &Syy_arr(i_J, k_J, 0, comp_yy), fpyy*weight_J*weight_E);
1332 }
1333 }
1334
1335 }
1336 }
1337 }
1338
1339 }
1340
1341 }
1342
1343#elif defined(WARPX_DIM_RCYLINDER) || defined(WARPX_DIM_RSPHERE)
1344
1345 // compute cell crossings in X-direction
1346 const auto i_old = static_cast<int>(x_old-shift);
1347 const auto i_new = static_cast<int>(x_new-shift);
1348 const int cell_crossings_x = std::abs(i_new-i_old);
1349 num_segments += cell_crossings_x;
1350
1351 // Compute the initial cell location used to find the cell crossings.
1352 // Keep these double to avoid bug in single precision
1353 const auto dirX_sign = static_cast<double>(dxp < 0. ? -1. : 1.);
1354 double Xcell = static_cast<double>(i_old) + shift + 0.5*(1.-dirX_sign);
1355
1356 // loop over the number of segments and deposit
1357 const Compute_shape_factor< depos_order-1 > compute_shape_factor_cell;
1358 const Compute_shape_factor_pair< depos_order > compute_shape_factors_node;
1359 double dxp_seg;
1360 double x0_new;
1361 double x0_old = x_old;
1362
1363 for (int ns=0; ns<num_segments; ns++) {
1364
1365 if (ns == num_segments-1) { // final segment
1366 x0_new = x_new;
1367 dxp_seg = x0_new - x0_old;
1368 }
1369 else {
1370 Xcell = Xcell + dirX_sign;
1371 x0_new = Xcell;
1372 dxp_seg = x0_new - x0_old;
1373 }
1374
1375 // Compute the segment factor (equal to dt_seg/dt for nonzero dxp)
1376 const auto seg_factor = static_cast<amrex::Real>(dxp == 0. ? 1._rt : dxp_seg/dxp);
1377
1378 // Compute cell-based weights using the average segment position
1379 // Keep these double to avoid bug in single precision
1380 double sx_cell[depos_order] = {0.};
1381 double const x0_bar = (x0_new + x0_old)/2.0;
1382 const int i0_cell = compute_shape_factor_cell( sx_cell, x0_bar-0.5 );
1383
1384 if constexpr (depos_order >= 3) { // higher-order correction to the cell-based weights
1385 const Compute_shape_factor_pair<depos_order-1> compute_shape_factors_cell;
1386 double sx_old_cell[depos_order] = {0.};
1387 double sx_new_cell[depos_order] = {0.};
1388 const int i0_cell_2 = compute_shape_factors_cell( sx_old_cell, sx_new_cell, x0_old-0.5, x0_new-0.5 );
1389 amrex::ignore_unused(i0_cell_2);
1390 for (int m=0; m<depos_order; m++) {
1391 sx_cell[m] = (4.0*sx_cell[m] + sx_old_cell[m] + sx_new_cell[m])/6.0;
1392 }
1393 }
1394
1395 // Compute node-based weights using the old and new segment positions
1396 // Keep these double to avoid bug in single precision
1397 double sx_old_node[depos_order+1] = {0.};
1398 double sx_new_node[depos_order+1] = {0.};
1399 const int i0_node = compute_shape_factors_node( sx_old_node, sx_new_node, x0_old, x0_new );
1400
1401 // deposit out-of-plane Jy, Jz, Syy, and Szz for this segment
1402 for (int i=0; i<=depos_order; i++) {
1403 const amrex::Real weight = 0.5_rt*(sx_old_node[i] + sx_new_node[i])*seg_factor;
1404 if constexpr (deposit_J) {
1405 amrex::Gpu::Atomic::AddNoRet( &Jy_arr(lo.x+i0_node+i, 0, 0), wqy*weight);
1406 amrex::Gpu::Atomic::AddNoRet( &Jz_arr(lo.x+i0_node+i, 0, 0), wqz*weight);
1407 }
1408 //
1409 amrex::Gpu::Atomic::AddNoRet( &Syy_arr(lo.x+i0_node+i, 0, 0, 0), fpyy*weight*weight);
1410 amrex::Gpu::Atomic::AddNoRet( &Szz_arr(lo.x+i0_node+i, 0, 0, 0), fpzz*weight*weight);
1411 }
1412
1413 // deposit Jx and Sxx for this segment
1414 for (int i=0; i<=depos_order-1; i++) {
1415 const amrex::Real weight = sx_cell[i]*seg_factor;
1416 if constexpr (deposit_J) {
1417 amrex::Gpu::Atomic::AddNoRet( &Jx_arr(lo.x+i0_cell+i, 0, 0), wqx*weight);
1418 }
1419 //
1420 amrex::Gpu::Atomic::AddNoRet( &Sxx_arr(lo.x+i0_cell+i, 0, 0, 0), fpxx*weight*weight);
1421 }
1422
1423 // update old segment values
1424 if (ns < num_segments-1) {
1425 x0_old = x0_new;
1426 }
1427
1428 }
1429
1430#elif defined(WARPX_DIM_1D_Z)
1431
1432 // compute cell crossings in Z-direction
1433 const auto k_old = static_cast<int>(z_old-shift);
1434 const auto k_new = static_cast<int>(z_new-shift);
1435 const int cell_crossings_z = std::abs(k_new-k_old);
1436 num_segments += cell_crossings_z;
1437
1438 // Compute initial particle cell location used to find cell crossings.
1439 // Keep these double to avoid bug in single precision
1440 const auto dirZ_sign = static_cast<double>(dzp < 0. ? -1. : 1.);
1441 double Zcell = static_cast<double>(k_old) + shift + 0.5*(1.-dirZ_sign);
1442
1443 // loop over the number of segments and deposit
1444 const Compute_shape_factor< depos_order-1 > compute_shape_factor_cell;
1445 const Compute_shape_factor_pair< depos_order > compute_shape_factors_node;
1446 double dzp_seg;
1447 double z0_new;
1448 double z0_old = z_old;
1449
1450 constexpr int num_segments_max = 1 + 4*AMREX_SPACEDIM;
1451 AMREX_ALWAYS_ASSERT_WITH_MESSAGE( num_segments <= num_segments_max,
1452 "Error: num_segments must be less than or equal to 1 + 4*AMREX_SPACEDIM.");
1453
1454 // Save the start index and interpolation weights for each segment
1455 int k0_cell[num_segments_max];
1456 int k0_node[num_segments_max];
1457 amrex::Real weight_cell[num_segments_max][depos_order];
1458 amrex::Real weight_node[num_segments_max][depos_order+1];
1459
1460 const auto k_mid = static_cast<int>(0.5*(z_new+z_old)-shift);
1461 int SegNum[num_segments_max];
1462
1463 for (int ns=0; ns<num_segments; ns++) {
1464
1465 if (ns == num_segments-1) { // final segment
1466 z0_new = z_new;
1467 dzp_seg = z0_new - z0_old;
1468 }
1469 else {
1470 Zcell = Zcell + dirZ_sign;
1471 z0_new = Zcell;
1472 dzp_seg = z0_new - z0_old;
1473 }
1474
1475 // Compute the segment factor (equal to dt_seg/dt for nonzero dzp)
1476 const auto seg_factor = static_cast<amrex::Real>(dzp == 0. ? 1._rt : dzp_seg/dzp);
1477
1478 // Compute cell-based weights using the average segment position
1479 // Keep these double to avoid bug in single precision
1480 double sz_cell[depos_order] = {0.};
1481 double const z0_bar = (z0_new + z0_old)/2.0;
1482 k0_cell[ns] = compute_shape_factor_cell( sz_cell, z0_bar-0.5 );
1483
1484 // Set the segment number for the mass matrix component calc
1485 if constexpr (full_mass_matrices) {
1486 const auto k0_mid = static_cast<int>(z0_bar-shift);
1487 SegNum[ns] = 1 + k0_mid - k_mid;
1488 }
1489
1490 if constexpr (depos_order >= 3) { // higher-order correction to the cell-based weights
1491 const Compute_shape_factor_pair<depos_order-1> compute_shape_factors_cell;
1492 double sz_old_cell[depos_order] = {0.};
1493 double sz_new_cell[depos_order] = {0.};
1494 const int k0_cell_2 = compute_shape_factors_cell( sz_old_cell, sz_new_cell, z0_old-0.5, z0_new-0.5 );
1495 amrex::ignore_unused(k0_cell_2);
1496 for (int m=0; m<depos_order; m++) {
1497 sz_cell[m] = (4.0*sz_cell[m] + sz_old_cell[m] + sz_new_cell[m])/6.0;
1498 }
1499 }
1500
1501 // Compute node-based weights using the old and new segment positions
1502 // Keep these double to avoid bug in single precision
1503 double sz_old_node[depos_order+1] = {0.};
1504 double sz_new_node[depos_order+1] = {0.};
1505 k0_node[ns] = compute_shape_factors_node( sz_old_node, sz_new_node, z0_old, z0_new );
1506
1507 // deposit out-of-plane Jx, Jy, Sx, and Sy for this segment
1508 for (int k=0; k<=depos_order; k++) {
1509 const amrex::Real weight = 0.5_rt*(sz_old_node[k] + sz_new_node[k])*seg_factor;
1510 const int k_J = lo.x + k0_node[ns] + k;
1511 if constexpr (deposit_J) {
1512 amrex::Gpu::Atomic::AddNoRet(&Jx_arr(k_J, 0, 0), wqx*weight);
1513 amrex::Gpu::Atomic::AddNoRet(&Jy_arr(k_J, 0, 0), wqy*weight);
1514 }
1515 if constexpr (full_mass_matrices) { weight_node[ns][k] = weight; }
1516 else {
1517 amrex::Gpu::Atomic::AddNoRet(&Sxx_arr(k_J, 0, 0, 0), fpxx*weight*weight);
1518 amrex::Gpu::Atomic::AddNoRet(&Syy_arr(k_J, 0, 0, 0), fpyy*weight*weight);
1519 }
1520 }
1521
1522 // deposit Jz and Szz for this segment
1523 for (int k=0; k<=depos_order-1; k++) {
1524 const amrex::Real weight = sz_cell[k]*seg_factor;
1525 const int k_J = lo.x + k0_cell[ns] + k;
1526 if constexpr (deposit_J) {
1527 amrex::Gpu::Atomic::AddNoRet(&Jz_arr(k_J, 0, 0), wqz*weight);
1528 }
1529 if constexpr (full_mass_matrices) { weight_cell[ns][k] = weight; }
1530 else {
1531 amrex::Gpu::Atomic::AddNoRet(&Szz_arr(k_J, 0, 0, 0), fpzz*weight*weight);
1532 }
1533 }
1534
1535 // update old segment values
1536 if (ns < num_segments-1) {
1537 z0_old = z0_new;
1538 }
1539
1540 }
1541
1542 if constexpr (full_mass_matrices) {
1543
1544 const int width_zz = depos_order - 1 + max_crossings;
1545 const int width_yy = depos_order + max_crossings;
1546
1547 // Loop over segments and deposit full mass matrices
1548 for (int ns=0; ns<num_segments; ns++) {
1549
1550 // Deposit Sxx, Sxy, Sxz, Syx, Syy, and Syz for this segment
1551 for (int k=0; k<=depos_order; k++) {
1552
1553 const int k_J = lo.x + k0_node[ns] + k;
1554 const amrex::Real weight_J = weight_node[ns][k];
1555 for (int ms=0; ms<num_segments; ms++) {
1556 const int SegShift = max_crossings + SegNum[ms] - SegNum[ns];
1557 for (int kE=0; kE<=depos_order; kE++) {
1558 const amrex::Real weight_E = weight_node[ms][kE];
1559 const int comp_yy = depos_order - k + kE + SegShift;
1560 if (comp_yy <= width_yy) { // Reduced deposit for diagonal mass matrices
1561 amrex::Gpu::Atomic::AddNoRet(&Sxx_arr(k_J, 0, 0, comp_yy), fpxx*weight_J*weight_E);
1562 amrex::Gpu::Atomic::AddNoRet(&Syy_arr(k_J, 0, 0, comp_yy), fpyy*weight_J*weight_E);
1563 }
1564 amrex::Gpu::Atomic::AddNoRet(&Sxy_arr(k_J, 0, 0, comp_yy), fpxy*weight_J*weight_E);
1565 amrex::Gpu::Atomic::AddNoRet(&Syx_arr(k_J, 0, 0, comp_yy), fpyx*weight_J*weight_E);
1566 }
1567 for (int kE=0; kE<=depos_order-1; kE++) {
1568 const amrex::Real weight_E = weight_cell[ms][kE];
1569 const int comp_yz = depos_order - k + kE + SegShift;
1570 amrex::Gpu::Atomic::AddNoRet(&Sxz_arr(k_J, 0, 0, comp_yz), fpxz*weight_J*weight_E);
1571 amrex::Gpu::Atomic::AddNoRet(&Syz_arr(k_J, 0, 0, comp_yz), fpyz*weight_J*weight_E);
1572 }
1573 }
1574
1575 }
1576
1577 // Deposit Szx, Szy, and Szz for this segment
1578 for (int k=0; k<=depos_order-1; k++) {
1579
1580 const int k_J = lo.x + k0_cell[ns] + k;
1581 const amrex::Real weight_J = weight_cell[ns][k];
1582 for (int ms=0; ms<num_segments; ms++) {
1583 const int SegShift = max_crossings + SegNum[ms] - SegNum[ns];
1584 for (int kE=0; kE<=depos_order-1; kE++) {
1585 const int comp_zz = depos_order - 1 - k + kE + SegShift;
1586 if (comp_zz > width_zz) { break; } // Reduced deposit for diagonal mass matrices
1587 const amrex::Real weight_E = weight_cell[ms][kE];
1588 amrex::Gpu::Atomic::AddNoRet(&Szz_arr(k_J, 0, 0, comp_zz), fpzz*weight_J*weight_E);
1589 }
1590 for (int kE=0; kE<=depos_order; kE++) {
1591 const amrex::Real weight_E = weight_node[ms][kE];
1592 const int comp_zy = depos_order-1 - k + kE + SegShift;
1593 amrex::Gpu::Atomic::AddNoRet(&Szx_arr(k_J, 0, 0, comp_zy), fpzx*weight_J*weight_E);
1594 amrex::Gpu::Atomic::AddNoRet(&Szy_arr(k_J, 0, 0, comp_zy), fpzy*weight_J*weight_E);
1595 }
1596 }
1597
1598 }
1599
1600 }
1601
1602 }
1603
1604#endif
1605}
1606
1635template <int depos_order, bool full_mass_matrices>
1636void doVillasenorSigmaDeposition ([[maybe_unused]] const amrex::ParticleReal* xp_n_data,
1637 [[maybe_unused]] const amrex::ParticleReal* yp_n_data,
1638 [[maybe_unused]] const amrex::ParticleReal* zp_n_data,
1639 const GetParticlePosition<PIdx>& GetPosition,
1640 [[maybe_unused]] const int* nsuborbits,
1641 const amrex::ParticleReal* wp,
1642 const amrex::ParticleReal* uxp_n,
1643 const amrex::ParticleReal* uyp_n,
1644 const amrex::ParticleReal* uzp_n,
1645 const amrex::ParticleReal* uxp_nph,
1646 const amrex::ParticleReal* uyp_nph,
1647 const amrex::ParticleReal* uzp_nph,
1648 const int max_crossings,
1649 amrex::Array4<amrex::Real> const& Sxx_arr,
1650 amrex::Array4<amrex::Real> const& Sxy_arr,
1651 amrex::Array4<amrex::Real> const& Sxz_arr,
1652 amrex::Array4<amrex::Real> const& Syx_arr,
1653 amrex::Array4<amrex::Real> const& Syy_arr,
1654 amrex::Array4<amrex::Real> const& Syz_arr,
1655 amrex::Array4<amrex::Real> const& Szx_arr,
1656 amrex::Array4<amrex::Real> const& Szy_arr,
1657 amrex::Array4<amrex::Real> const& Szz_arr,
1658 const amrex::Array4<amrex::Real const>& Bx_arr,
1659 const amrex::Array4<amrex::Real const>& By_arr,
1660 const amrex::Array4<amrex::Real const>& Bz_arr,
1661 const amrex::IndexType Bx_type,
1662 const amrex::IndexType By_type,
1663 const amrex::IndexType Bz_type,
1664 const long np_to_deposit,
1665 const amrex::Real dt,
1666 const amrex::XDim3& dinv,
1667 const amrex::XDim3& xyzmin,
1668 const amrex::GpuArray<amrex::GpuArray<double,2>, AMREX_SPACEDIM> & domain_double,
1669 const amrex::GpuArray<amrex::GpuArray<bool,2>, AMREX_SPACEDIM> & do_cropping,
1670 const amrex::Dim3 lo,
1671 const amrex::Real qs,
1672 const amrex::Real ms)
1673{
1674 using namespace amrex::literals;
1675
1676 const amrex::Real invvol = dinv.x*dinv.y*dinv.z;
1677
1678 // Loop over particles and deposit mass matrices
1679 amrex::ParallelFor(
1680 np_to_deposit,
1681 [=] AMREX_GPU_DEVICE (long const ip) {
1682
1683 // Skip mass matrix deposition for particles with suborbits.
1684 if (nsuborbits && nsuborbits[ip] > 1) { return; }
1685
1686 // Skip particles with zero weight.
1687 if (wp[ip] == 0.) { 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 const int n_rz_azimuthal_modes = 0;
1699 doDirectGatherVectorField<depos_order_perp,depos_order_para>(
1700 xp_nph, yp_nph, zp_nph,
1701 Bxp, Byp, Bzp,
1702 Bx_arr, By_arr, Bz_arr,
1703 Bx_type, By_type, Bz_type,
1704 dinv, xyzmin, lo, n_rz_azimuthal_modes );
1705
1706 // Compute inverse Lorentz factor, the average of gamma at time levels n and n+1
1707 const amrex::ParticleReal gaminv = GetImplicitGammaInverse(uxp_n[ip], uyp_n[ip], uzp_n[ip],
1708 uxp_nph[ip], uyp_nph[ip], uzp_nph[ip]);
1709
1710 // Compute current density kernels to deposit
1711 const amrex::Real wq_invvol = qs*wp[ip]*invvol;
1712 const amrex::Real rhop = wq_invvol*gaminv;
1713
1714 // Set the Mass Matrices kernels
1715 amrex::ParticleReal fpxx, fpxy, fpxz;
1716 amrex::ParticleReal fpyx, fpyy, fpyz;
1717 amrex::ParticleReal fpzx, fpzy, fpzz;
1718 setMassMatricesKernels( qs, ms, dt, rhop,
1719 uxp_nph[ip], uyp_nph[ip], uzp_nph[ip],
1720 Bxp, Byp, Bzp,
1721 fpxx, fpxy, fpxz,
1722 fpyx, fpyy, fpyz,
1723 fpzx, fpzy, fpzz );
1724
1725 amrex::ParticleReal const xp_n = (xp_n_data ? xp_n_data[ip] : 0._prt);
1726 amrex::ParticleReal const yp_n = (yp_n_data ? yp_n_data[ip] : 0._prt);
1727 amrex::ParticleReal const zp_n = (zp_n_data ? zp_n_data[ip] : 0._prt);
1728
1729 // Compute position at time n + 1
1730 amrex::ParticleReal const xp_np1 = 2._prt*xp_nph - xp_n;
1731 amrex::ParticleReal const yp_np1 = 2._prt*yp_nph - yp_n;
1732 amrex::ParticleReal const zp_np1 = 2._prt*zp_nph - zp_n;
1733
1734 // Pass dummy arrays for Jx, Jy, Jz (which will not be used)
1735 amrex::Array4<amrex::Real> const dummy_Jx{};
1736 amrex::Array4<amrex::Real> const dummy_Jy{};
1737 amrex::Array4<amrex::Real> const dummy_Jz{};
1738
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:650
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
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:1636
doDirectSigmaDeposition
void doDirectSigmaDeposition(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 &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:521
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::c
constexpr auto c
vacuum speed of light [m/s]
Definition constant.H:153
amrex::Gpu::Atomic::AddNoRet
__host__ __device__ AMREX_FORCE_INLINE void AddNoRet(T *sum, T value) noexcept
amrex::ignore_unused
__host__ __device__ void ignore_unused(const Ts &...)
amrex::ParallelFor
std::enable_if_t< std::is_integral_v< T > > ParallelFor(TypeList< CTOs... > ctos, std::array< int, sizeof...(CTOs)> const &runtime_options, T N, F &&f)
amrex::shift
__host__ __device__ BoxND< dim > shift(const BoxND< dim > &b, int dir, int nzones) noexcept
amrex::IndexType
IndexTypeND< 3 > IndexType
amrex::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