1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
|
/* Copyright 2019 Axel Huebl, David Grote, Maxence Thevenet
* Weiqun Zhang
*
* This file is part of WarpX.
*
* License: BSD-3-Clause-LBNL
*/
#ifndef CHARGEDEPOSITION_H_
#define CHARGEDEPOSITION_H_
#include "Particles/Pusher/GetAndSetPosition.H"
#include "Particles/ShapeFactors.H"
#include <AMReX.H>
/* \brief Charge Deposition for thread thread_num
* /param GetPosition : A functor for returning the particle position.
* \param wp : Pointer to array of particle weights.
* \param ion_lev : Pointer to array of particle ionization level. This is
required to have the charge of each macroparticle
since q is a scalar. For non-ionizable species,
ion_lev is a null pointer.
* \param rho_fab : FArrayBox of charge density, either full array or tile.
* \param np_to_depose : Number of particles for which current is deposited.
* \param dx : 3D cell size
* \param xyzmin : Physical lower bounds of domain.
* \param lo : Index lower bounds of domain.
* /param q : species charge.
*/
template <int depos_order>
void doChargeDepositionShapeN(const GetParticlePosition& GetPosition,
const amrex::ParticleReal * const wp,
const int * const ion_lev,
amrex::FArrayBox& rho_fab,
const long np_to_depose,
const std::array<amrex::Real,3>& dx,
const std::array<amrex::Real, 3> xyzmin,
const amrex::Dim3 lo,
const amrex::Real q,
const long n_rz_azimuthal_modes)
{
using namespace amrex;
// Whether ion_lev is a null pointer (do_ionization=0) or a real pointer
// (do_ionization=1)
const bool do_ionization = ion_lev;
const amrex::Real dxi = 1.0/dx[0];
const amrex::Real dzi = 1.0/dx[2];
#if (AMREX_SPACEDIM == 2)
const amrex::Real invvol = dxi*dzi;
#elif (defined WARPX_DIM_3D)
const amrex::Real dyi = 1.0/dx[1];
const amrex::Real invvol = dxi*dyi*dzi;
#endif
const amrex::Real xmin = xyzmin[0];
#if (defined WARPX_DIM_3D)
const amrex::Real ymin = xyzmin[1];
#endif
const amrex::Real zmin = xyzmin[2];
amrex::Array4<amrex::Real> const& rho_arr = rho_fab.array();
amrex::IntVect const rho_type = rho_fab.box().type();
constexpr int zdir = (AMREX_SPACEDIM - 1);
constexpr int NODE = amrex::IndexType::NODE;
constexpr int CELL = amrex::IndexType::CELL;
// Loop over particles and deposit into rho_fab
amrex::ParallelFor(
np_to_depose,
[=] AMREX_GPU_DEVICE (long ip) {
// --- Get particle quantities
amrex::Real wq = q*wp[ip]*invvol;
if (do_ionization){
wq *= ion_lev[ip];
}
amrex::ParticleReal xp, yp, zp;
GetPosition(ip, xp, yp, zp);
// --- Compute shape factors
// x direction
// Get particle position in grid coordinates
#if (defined WARPX_DIM_RZ)
const amrex::Real rp = std::sqrt(xp*xp + yp*yp);
amrex::Real costheta;
amrex::Real sintheta;
if (rp > 0.) {
costheta = xp/rp;
sintheta = yp/rp;
} else {
costheta = 1.;
sintheta = 0.;
}
const Complex xy0 = Complex{costheta, sintheta};
const amrex::Real x = (rp - xmin)*dxi;
#else
const amrex::Real x = (xp - xmin)*dxi;
#endif
// Compute shape factor along x
// i: leftmost grid point that the particle touches
amrex::Real sx[depos_order + 1];
int i = 0;
Compute_shape_factor< depos_order > const compute_shape_factor;
if (rho_type[0] == NODE) {
i = compute_shape_factor(sx, x);
} else if (rho_type[0] == CELL) {
i = compute_shape_factor(sx, x - 0.5_rt);
}
#if (defined WARPX_DIM_3D)
// y direction
const amrex::Real y = (yp - ymin)*dyi;
amrex::Real sy[depos_order + 1];
int j = 0;
if (rho_type[1] == NODE) {
j = compute_shape_factor(sy, y);
} else if (rho_type[1] == CELL) {
j = compute_shape_factor(sy, y - 0.5_rt);
}
#endif
// z direction
const amrex::Real z = (zp - zmin)*dzi;
amrex::Real sz[depos_order + 1];
int k = 0;
if (rho_type[zdir] == NODE) {
k = compute_shape_factor(sz, z);
} else if (rho_type[zdir] == CELL) {
k = compute_shape_factor(sz, z - 0.5_rt);
}
// Deposit charge into rho_arr
#if (defined WARPX_DIM_XZ) || (defined WARPX_DIM_RZ)
for (int iz=0; iz<=depos_order; iz++){
for (int ix=0; ix<=depos_order; ix++){
amrex::Gpu::Atomic::Add(
&rho_arr(lo.x+i+ix, lo.y+k+iz, 0, 0),
sx[ix]*sz[iz]*wq);
#if (defined WARPX_DIM_RZ)
Complex xy = xy0; // Throughout the following loop, xy takes the value e^{i m theta}
for (int imode=1 ; imode < n_rz_azimuthal_modes ; imode++) {
// The factor 2 on the weighting comes from the normalization of the modes
amrex::Gpu::Atomic::Add( &rho_arr(lo.x+i+ix, lo.y+k+iz, 0, 2*imode-1), 2._rt*sx[ix]*sz[iz]*wq*xy.real());
amrex::Gpu::Atomic::Add( &rho_arr(lo.x+i+ix, lo.y+k+iz, 0, 2*imode ), 2._rt*sx[ix]*sz[iz]*wq*xy.imag());
xy = xy*xy0;
}
#endif
}
}
#elif (defined WARPX_DIM_3D)
for (int iz=0; iz<=depos_order; iz++){
for (int iy=0; iy<=depos_order; iy++){
for (int ix=0; ix<=depos_order; ix++){
amrex::Gpu::Atomic::Add(
&rho_arr(lo.x+i+ix, lo.y+j+iy, lo.z+k+iz),
sx[ix]*sy[iy]*sz[iz]*wq);
}
}
}
#endif
}
);
#ifndef WARPX_DIM_RZ
amrex::ignore_unused(n_rz_azimuthal_modes);
#endif
}
#endif // CHARGEDEPOSITION_H_
|
