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
|
#include <AMReX_ParallelDescriptor.H>
#include <AMReX_MLMG.H>
#include <AMReX_MLNodeLaplacian.H>
#include <WarpX.H>
using namespace amrex;
void
WarpX::InitSpaceChargeField (WarpXParticleContainer& pc)
{
#ifdef WARPX_DIM_RZ
amrex::Abort("The initialization of space-charge field has not yet been implemented in RZ geometry.");
#endif
// Allocate fields for charge and potential
const int num_levels = max_level + 1;
Vector<std::unique_ptr<MultiFab> > rho(num_levels);
Vector<std::unique_ptr<MultiFab> > phi(num_levels);
const int ng = WarpX::nox;
for (int lev = 0; lev <= max_level; lev++) {
BoxArray nba = boxArray(lev);
nba.surroundingNodes();
rho[lev].reset(new MultiFab(nba, dmap[lev], 1, ng)); // Make ng big enough/use rho from sim
phi[lev].reset(new MultiFab(nba, dmap[lev], 1, 0));
phi[lev]->setVal(0.);
}
// Deposit particle charge density (source of Poisson solver)
bool const local = false;
bool const reset = true;
bool const do_rz_volume_scaling = true;
pc.DepositCharge(rho, local, reset, do_rz_volume_scaling);
// Compute the potential phi, by solving the Poisson equation
computePhi( rho, phi );
// Compute the corresponding electric field, from the potential phi
computeE( Efield_fp, phi );
}
/* Compute the potential `phi` by solving the Poisson equation with `rho` as
a source. This uses the amrex solver.
\param[in] rho The charge density a given species
\param[out] phi The potential to be computed by this function
*/
void
WarpX::computePhi(const amrex::Vector<std::unique_ptr<amrex::MultiFab> >& rho,
amrex::Vector<std::unique_ptr<amrex::MultiFab> >& phi) const
{
// Define the boundary conditions
Array<LinOpBCType,AMREX_SPACEDIM> lobc, hibc;
for (int idim=0; idim<AMREX_SPACEDIM; idim++){
if ( Geom(0).isPeriodic(idim) ) {
lobc[idim] = LinOpBCType::Periodic;
hibc[idim] = LinOpBCType::Periodic;
} else {
// Use Dirichlet boundary condition by default.
// Ideally, we would often want open boundary conditions here.
lobc[idim] = LinOpBCType::Dirichlet;
hibc[idim] = LinOpBCType::Dirichlet;
}
}
// Define the linear operator (Poisson operator)
MLNodeLaplacian linop( Geom(), boxArray(), DistributionMap() );
linop.setDomainBC( lobc, hibc );
for (int lev = 0; lev <= max_level; lev++) {
BoxArray cba = boxArray(lev);
cba.enclosedCells();
MultiFab sigma(cba, DistributionMap(lev), 1, 0);
sigma.setVal(-PhysConst::ep0);
linop.setSigma(lev, sigma);
}
// Solve the Poisson equation
MLMG mlmg(linop);
mlmg.setVerbose(1);
const Real reltol = 1.e-11;
mlmg.solve( GetVecOfPtrs(phi), GetVecOfConstPtrs(rho), reltol, 0.0);
}
/* \bried Compute the electric field that corresponds to `phi`, and
add it to the set of MultiFab `E`.
\param[inout] E Electric field on the grid
\param[in] phi The potential from which to compute the electric field
*/
void
WarpX::computeE(amrex::Vector<std::array<std::unique_ptr<amrex::MultiFab>, 3> >& E,
const amrex::Vector<std::unique_ptr<amrex::MultiFab> >& phi) const
{
for (int lev = 0; lev <= max_level; lev++) {
const Real* dx = Geom(lev).CellSize();
#ifdef _OPENMP
#pragma omp parallel if (Gpu::notInLaunchRegion())
#endif
for ( MFIter mfi(*phi[lev], TilingIfNotGPU()); mfi.isValid(); ++mfi )
{
const Real inv_dx = 1./dx[0];
#if (AMREX_SPACEDIM == 3)
const Real inv_dy = 1./dx[1];
const Real inv_dz = 1./dx[2];
#else
const Real inv_dz = 1./dx[1];
#endif
const Box& tbx = mfi.tilebox(Ex_nodal_flag);
const Box& tby = mfi.tilebox(Ey_nodal_flag);
const Box& tbz = mfi.tilebox(Ez_nodal_flag);
const auto& phi_arr = phi[lev]->array(mfi);
const auto& Ex_arr = (*E[lev][0])[mfi].array();
const auto& Ey_arr = (*E[lev][1])[mfi].array();
const auto& Ez_arr = (*E[lev][2])[mfi].array();
#if (AMREX_SPACEDIM == 3)
amrex::ParallelFor( tbx, tby, tbz,
[=] AMREX_GPU_DEVICE (int i, int j, int k) {
Ex_arr(i,j,k) += -inv_dx*( phi_arr(i+1,j,k) - phi_arr(i,j,k) );
},
[=] AMREX_GPU_DEVICE (int i, int j, int k) {
Ey_arr(i,j,k) += -inv_dy*( phi_arr(i,j+1,k) - phi_arr(i,j,k) );
},
[=] AMREX_GPU_DEVICE (int i, int j, int k) {
Ez_arr(i,j,k) += -inv_dz*( phi_arr(i,j,k+1) - phi_arr(i,j,k) );
}
);
#else
amrex::ParallelFor( tbx, tbz,
[=] AMREX_GPU_DEVICE (int i, int j, int k) {
Ex_arr(i,j,k) += -inv_dx*( phi_arr(i+1,j,k) - phi_arr(i,j,k) );
},
[=] AMREX_GPU_DEVICE (int i, int j, int k) {
Ez_arr(i,j,k) += -inv_dz*( phi_arr(i,j+1,k) - phi_arr(i,j,k) );
}
);
#endif
}
}
}
|
