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