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