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
|
//
// Each problem must have its own version of MyParticleContainer::InitData()
// to initialize the particle data on this level
//
#include <cmath>
#include <BLProfiler.H>
#include <ParticleContainer.H>
#include <WarpXConst.H>
void
MyParticleContainer::InitData()
{
BL_PROFILE("MyPC::InitData()");
charge = -PhysConst::q_e;
mass = PhysConst::m_e;
m_particles.resize(m_gdb->finestLevel()+1);
const int lev = 0;
const Geometry& geom = m_gdb->Geom(lev);
const Real* dx = geom.CellSize();
Real weight, ux, uy, uz;
Real particle_xmin, particle_xmax, particle_ymin, particle_ymax, particle_zmin, particle_zmax;
int n_part_per_cell;
{
ParmParse pp("langmuirwave");
n_part_per_cell = 1;
pp.query("num_particles_per_cell", n_part_per_cell);
weight = 1.e25;
pp.query("n_e", weight);
weight *= dx[0]*dx[1]*dx[2]/n_part_per_cell;
pp.query("particle_xmin", particle_xmin);
pp.query("particle_xmax", particle_xmax);
pp.query("particle_ymin", particle_ymin);
pp.query("particle_ymax", particle_ymax);
pp.query("particle_zmin", particle_zmin);
pp.query("particle_zmax", particle_zmax);
ux = 0.;
uy = 0.;
uz = 0.;
pp.query("ux", ux);
pp.query("uy", uy);
pp.query("uz", uz);
Real gamma = 1./std::sqrt(1.0 - ux*ux - uy*uy - uz*uz);
ux *= PhysConst::c*gamma;
uy *= PhysConst::c*gamma;
uz *= PhysConst::c*gamma;
}
const BoxArray& ba = m_gdb->ParticleBoxArray(lev);
const DistributionMapping& dm = m_gdb->ParticleDistributionMap(lev);
MultiFab dummy_mf(ba, 1, 0, dm, Fab_noallocate);
for (MFIter mfi(dummy_mf,false); mfi.isValid(); ++mfi)
{
int gid = mfi.index();
Box grid = ba[gid];
RealBox grid_box = RealBox(grid,dx,geom.ProbLo());
int nx = grid.length(0), ny = grid.length(1), nz = grid.length(2);
for (int k = 0; k < nz; k++) {
for (int j = 0; j < ny; j++) {
for (int i = 0; i < nx; i++) {
for (int i_part=0; i_part<n_part_per_cell;i_part++) {
Real particle_shift = (0.5+i_part)/n_part_per_cell;
Real x = grid_box.lo(0) + (i + particle_shift)*dx[0];
Real y = grid_box.lo(1) + (j + particle_shift)*dx[1];
Real z = grid_box.lo(2) + (k + particle_shift)*dx[2];
if (x >= particle_xmax || x < particle_xmin ||
y >= particle_ymax || y < particle_ymin ||
z >= particle_zmax || z < particle_zmin ) continue;
ParticleType p;
p.m_id = ParticleBase::NextID();
p.m_cpu = ParallelDescriptor::MyProc();
p.m_lev = lev;
p.m_grid = gid;
p.m_pos[0] = x;
p.m_pos[1] = y;
p.m_pos[2] = z;
for (int i = 0; i < BL_SPACEDIM; i++) {
BL_ASSERT(p.m_pos[i] < grid_box.hi(i));
}
p.m_data[PIdx::w] = weight;
for (int i = 1; i < PIdx::nattribs; i++) {
p.m_data[i] = 0;
}
p.m_data[PIdx::ux] = ux;
p.m_data[PIdx::uy] = uy;
p.m_data[PIdx::uz] = uz;
if (!ParticleBase::Where(p,m_gdb)) // this will set m_cell
{
BoxLib::Abort("invalid particle");
}
BL_ASSERT(p.m_lev >= 0 && p.m_lev <= m_gdb->finestLevel());
//
// Add it to the appropriate PBox at the appropriate level.
//
m_particles[p.m_lev][p.m_grid].push_back(p);
}
}
}
}
}
//
// We still need to redistribute in order to define each particle's cell, grid and level, but this
// shouldn't require any inter-node communication because the particles should already be in the right grid.
//
Redistribute(true);
}
|
