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
172
173
174
175
176
|
#include <limits>
#include <WarpX.H>
#include <WarpXConst.H>
#include <WarpX_f.H>
#include <ParticleContainer.H>
#include <ParticleIterator.H>
LaserContainer::LaserContainer (AmrCore* amr_core, int ispecies)
: WarpXParticleContainer(amr_core, ispecies)
{
charge = PhysConst::q_e; // note that q_e is defined to be positive.
mass = std::numeric_limits<Real>::max();
}
void
LaserContainer::AllocData ()
{
// have to resize here, not in the constructor because GDB was not
// ready in constructor.
m_particles.resize(GDB().finestLevel()+1);
}
void
LaserContainer::InitData ()
{
BoxLib::Abort("LaserContainer::InitData: Not implemented");
}
void
LaserContainer::Evolve (int lev,
const MultiFab&, const MultiFab&, const MultiFab&,
const MultiFab&, const MultiFab&, const MultiFab&,
MultiFab& jx, MultiFab& jy, MultiFab& jz, Real dt)
{
BL_PROFILE("Laser::Evolve()");
BL_PROFILE_VAR_NS("Laser::Evolve::Copy", blp_copy);
BL_PROFILE_VAR_NS("PICSAR::LaserParticlePush", blp_pxr_pp);
BL_PROFILE_VAR_NS("PICSAR::LaserCurrentDepo", blp_pxr_cd);
const Geometry& gm = GDB().Geom(lev);
const BoxArray& ba = jx.boxArray();
#if (BL_SPACEDIM == 3)
const Real* dx = gm.CellSize();
#elif (BL_SPACEDIM == 2)
Real dx[3] = { gm.CellSize(0), std::numeric_limits<Real>::quiet_NaN(), gm.CellSize(1) };
#endif
#if (BL_SPACEDIM == 3)
long ngx_j = jx.nGrow();
long ngy_j = ngx_j;
long ngz_j = ngx_j;
#elif (BL_SPACEDIM == 2)
long ngx_j = jx.nGrow();;
long ngy_j = 0;
long ngz_j = ngx_j;
#endif
BL_ASSERT(OnSameGrids(lev,jx));
{
Array<Real> xp, yp, zp, wp, uxp, uyp, uzp, giv;
PartIterInfo info {lev, do_tiling, tile_size};
for (PartIter pti(*this, info); pti.isValid(); ++pti)
{
const int gid = pti.index();
const Box& tbx = pti.tilebox();
const Box& vbx = pti.validbox();
const long np = pti.numParticles();
// Data on the grid
FArrayBox& jxfab = jx[gid];
FArrayBox& jyfab = jy[gid];
FArrayBox& jzfab = jz[gid];
xp.resize(np);
yp.resize(np);
zp.resize(np);
wp.resize(np);
uxp.resize(np);
uyp.resize(np);
uzp.resize(np);
giv.resize(np);
//
// copy data from particle container to temp arrays
//
BL_PROFILE_VAR_START(blp_copy);
pti.foreach([&](int i, ParticleType& p) {
#if (BL_SPACEDIM == 3)
xp[i] = p.m_pos[0];
yp[i] = p.m_pos[1];
zp[i] = p.m_pos[2];
#elif (BL_SPACEDIM == 2)
xp[i] = p.m_pos[0];
yp[i] = std::numeric_limits<Real>::quiet_NaN();
zp[i] = p.m_pos[1];
#endif
wp[i] = p.m_data[PIdx::w];
uxp[i] = p.m_data[PIdx::ux];
uyp[i] = p.m_data[PIdx::uy];
uzp[i] = p.m_data[PIdx::uz];
});
BL_PROFILE_VAR_STOP(blp_copy);
const Box& box = BoxLib::enclosedCells(ba[gid]);
BL_ASSERT(box == vbx);
#if (BL_SPACEDIM == 3)
long nx = box.length(0);
long ny = box.length(1);
long nz = box.length(2);
#elif (BL_SPACEDIM == 2)
long nx = box.length(0);
long ny = 0;
long nz = box.length(1);
#endif
RealBox grid_box = RealBox( box, gm.CellSize(), gm.ProbLo() );
#if (BL_SPACEDIM == 3)
const Real* xyzmin = grid_box.lo();
#elif (BL_SPACEDIM == 2)
Real xyzmin[3] = { grid_box.lo(0), std::numeric_limits<Real>::quiet_NaN(), grid_box.lo(1) };
#endif
//
// Particle Push
//
BL_PROFILE_VAR_START(blp_pxr_pp);
warpx_laser_pusher(&np, xp.data(), yp.data(), zp.data(),
uxp.data(), uyp.data(), uzp.data(), giv.data(),
&this->charge, &dt,
&WarpX::laser_pusher_algo);
BL_PROFILE_VAR_STOP(blp_pxr_pp);
//
// Current Deposition
// xxxxx this part needs to be thread safe if we have OpenMP over tiles
//
long lvect = 8;
BL_PROFILE_VAR_START(blp_pxr_cd);
warpx_current_deposition(jxfab.dataPtr(), jyfab.dataPtr(), jzfab.dataPtr(),
&np, xp.data(), yp.data(), zp.data(),
uxp.data(), uyp.data(), uzp.data(),
giv.data(), wp.data(), &this->charge,
&xyzmin[0], &xyzmin[1], &xyzmin[2],
&dt, &dx[0], &dx[1], &dx[2], &nx, &ny, &nz,
&ngx_j, &ngy_j, &ngz_j,
&WarpX::nox,&WarpX::noy,&WarpX::noz,
&lvect,&WarpX::current_deposition_algo);
BL_PROFILE_VAR_STOP(blp_pxr_cd);
//
// copy particle data back
//
BL_PROFILE_VAR_START(blp_copy);
pti.foreach([&](int i, ParticleType& p) {
BL_ASSERT(p.m_id > 0);
#if (BL_SPACEDIM == 3)
p.m_pos[0] = xp[i];
p.m_pos[1] = yp[i];
p.m_pos[2] = zp[i];
#elif (BL_SPACEDIM == 2)
p.m_pos[0] = xp[i];
p.m_pos[1] = zp[i];
#endif
p.m_data[PIdx::ux] = uxp[i];
p.m_data[PIdx::uy] = uyp[i];
p.m_data[PIdx::uz] = uzp[i];
});
BL_PROFILE_VAR_STOP(blp_copy);
}
}
}
|
