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
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
|
/* Copyright 2019-2020 David Grote
*
* This file is part of WarpX.
*
* License: BSD-3-Clause-LBNL
*/
#include "WarpX.H"
#include "PsatdAlgorithmRZ.H"
#include "Utils/WarpXConst.H"
#include "Utils/WarpXProfilerWrapper.H"
#include <cmath>
using amrex::operator""_rt;
/* \brief Initialize coefficients for the update equation */
PsatdAlgorithmRZ::PsatdAlgorithmRZ (SpectralKSpaceRZ const & spectral_kspace,
amrex::DistributionMapping const & dm,
int const n_rz_azimuthal_modes, int const norder_z,
bool const nodal, amrex::Real const dt,
bool const update_with_rho)
// Initialize members of base class
: SpectralBaseAlgorithmRZ(spectral_kspace, dm,
norder_z, nodal),
m_dt(dt),
m_update_with_rho(update_with_rho)
{
// Allocate the arrays of coefficients
amrex::BoxArray const & ba = spectral_kspace.spectralspace_ba;
C_coef = SpectralRealCoefficients(ba, dm, n_rz_azimuthal_modes, 0);
S_ck_coef = SpectralRealCoefficients(ba, dm, n_rz_azimuthal_modes, 0);
X1_coef = SpectralRealCoefficients(ba, dm, n_rz_azimuthal_modes, 0);
X2_coef = SpectralRealCoefficients(ba, dm, n_rz_azimuthal_modes, 0);
X3_coef = SpectralRealCoefficients(ba, dm, n_rz_azimuthal_modes, 0);
coefficients_initialized = false;
}
/* Advance the E and B field in spectral space (stored in `f`)
* over one time step */
void
PsatdAlgorithmRZ::pushSpectralFields(SpectralFieldDataRZ & f)
{
bool const update_with_rho = m_update_with_rho;
if (not coefficients_initialized) {
// This is called from here since it needs the kr values
// which can be obtained from the SpectralFieldDataRZ
InitializeSpectralCoefficients(f);
coefficients_initialized = true;
}
// Loop over boxes
for (amrex::MFIter mfi(f.fields); mfi.isValid(); ++mfi){
amrex::Box const & bx = f.fields[mfi].box();
// Extract arrays for the fields to be updated
amrex::Array4<Complex> const& fields = f.fields[mfi].array();
// Extract arrays for the coefficients
amrex::Array4<const amrex::Real> const& C_arr = C_coef[mfi].array();
amrex::Array4<const amrex::Real> const& S_ck_arr = S_ck_coef[mfi].array();
amrex::Array4<const amrex::Real> const& X1_arr = X1_coef[mfi].array();
amrex::Array4<const amrex::Real> const& X2_arr = X2_coef[mfi].array();
amrex::Array4<const amrex::Real> const& X3_arr = X3_coef[mfi].array();
// Extract pointers for the k vectors
auto const & kr_modes = f.getKrArray(mfi);
amrex::Real const* kr_arr = kr_modes.dataPtr();
amrex::Real const* modified_kz_arr = modified_kz_vec[mfi].dataPtr();
int const nr = bx.length(0);
amrex::Real const dt = m_dt;
// Loop over indices within one box
// Note that k = 0
int const modes = f.n_rz_azimuthal_modes;
amrex::ParallelFor(bx, modes,
[=] AMREX_GPU_DEVICE(int i, int j, int k, int mode) noexcept
{
// All of the fields of each mode are grouped together
using Idx = SpectralFieldIndex;
auto const Ep_m = Idx::Ex + Idx::n_fields*mode;
auto const Em_m = Idx::Ey + Idx::n_fields*mode;
auto const Ez_m = Idx::Ez + Idx::n_fields*mode;
auto const Bp_m = Idx::Bx + Idx::n_fields*mode;
auto const Bm_m = Idx::By + Idx::n_fields*mode;
auto const Bz_m = Idx::Bz + Idx::n_fields*mode;
auto const Jp_m = Idx::Jx + Idx::n_fields*mode;
auto const Jm_m = Idx::Jy + Idx::n_fields*mode;
auto const Jz_m = Idx::Jz + Idx::n_fields*mode;
auto const rho_old_m = Idx::rho_old + Idx::n_fields*mode;
auto const rho_new_m = Idx::rho_new + Idx::n_fields*mode;
// Record old values of the fields to be updated
Complex const Ep_old = fields(i,j,k,Ep_m);
Complex const Em_old = fields(i,j,k,Em_m);
Complex const Ez_old = fields(i,j,k,Ez_m);
Complex const Bp_old = fields(i,j,k,Bp_m);
Complex const Bm_old = fields(i,j,k,Bm_m);
Complex const Bz_old = fields(i,j,k,Bz_m);
// Shortcut for the values of J and rho
Complex const Jp = fields(i,j,k,Jp_m);
Complex const Jm = fields(i,j,k,Jm_m);
Complex const Jz = fields(i,j,k,Jz_m);
Complex const rho_old = fields(i,j,k,rho_old_m);
Complex const rho_new = fields(i,j,k,rho_new_m);
// k vector values, and coefficients
// The k values for each mode are grouped together
int const ir = i + nr*mode;
amrex::Real const kr = kr_arr[ir];
amrex::Real const kz = modified_kz_arr[j];
constexpr amrex::Real c2 = PhysConst::c*PhysConst::c;
constexpr amrex::Real inv_ep0 = 1._rt/PhysConst::ep0;
Complex const I = Complex{0._rt,1._rt};
amrex::Real const C = C_arr(i,j,k,mode);
amrex::Real const S_ck = S_ck_arr(i,j,k,mode);
amrex::Real const X1 = X1_arr(i,j,k,mode);
amrex::Real const X2 = X2_arr(i,j,k,mode);
amrex::Real const X3 = X3_arr(i,j,k,mode);
Complex rho_diff;
if (update_with_rho) {
rho_diff = X2*rho_new - X3*rho_old;
} else {
Complex const divE = kr*(Ep_old - Em_old) + I*kz*Ez_old;
Complex const divJ = kr*(Jp - Jm) + I*kz*Jz;
rho_diff = (X2 - X3)*PhysConst::ep0*divE - X2*dt*divJ;
}
// Update E (see WarpX online documentation: theory section)
fields(i,j,k,Ep_m) = C*Ep_old
+ S_ck*(-c2*I*kr/2._rt*Bz_old + c2*kz*Bp_old - inv_ep0*Jp)
+ 0.5_rt*kr*rho_diff;
fields(i,j,k,Em_m) = C*Em_old
+ S_ck*(-c2*I*kr/2._rt*Bz_old - c2*kz*Bm_old - inv_ep0*Jm)
- 0.5_rt*kr*rho_diff;
fields(i,j,k,Ez_m) = C*Ez_old
+ S_ck*(c2*I*kr*Bp_old + c2*I*kr*Bm_old - inv_ep0*Jz)
- I*kz*rho_diff;
// Update B (see WarpX online documentation: theory section)
fields(i,j,k,Bp_m) = C*Bp_old
- S_ck*(-I*kr/2._rt*Ez_old + kz*Ep_old)
+ X1*(-I*kr/2._rt*Jz + kz*Jp);
fields(i,j,k,Bm_m) = C*Bm_old
- S_ck*(-I*kr/2._rt*Ez_old - kz*Em_old)
+ X1*(-I*kr/2._rt*Jz - kz*Jm);
fields(i,j,k,Bz_m) = C*Bz_old
- S_ck*I*(kr*Ep_old + kr*Em_old)
+ X1*I*(kr*Jp + kr*Jm);
});
}
}
void PsatdAlgorithmRZ::InitializeSpectralCoefficients (SpectralFieldDataRZ const & f)
{
// Fill them with the right values:
// Loop over boxes and allocate the corresponding coefficients
// for each box owned by the local MPI proc
for (amrex::MFIter mfi(f.fields); mfi.isValid(); ++mfi){
amrex::Box const & bx = f.fields[mfi].box();
// Extract pointers for the k vectors
amrex::Real const* const modified_kz = modified_kz_vec[mfi].dataPtr();
// Extract arrays for the coefficients
amrex::Array4<amrex::Real> const& C = C_coef[mfi].array();
amrex::Array4<amrex::Real> const& S_ck = S_ck_coef[mfi].array();
amrex::Array4<amrex::Real> const& X1 = X1_coef[mfi].array();
amrex::Array4<amrex::Real> const& X2 = X2_coef[mfi].array();
amrex::Array4<amrex::Real> const& X3 = X3_coef[mfi].array();
auto const & kr_modes = f.getKrArray(mfi);
amrex::Real const* kr_arr = kr_modes.dataPtr();
int const nr = bx.length(0);
amrex::Real const dt = m_dt;
// Loop over indices within one box
int const modes = f.n_rz_azimuthal_modes;
amrex::ParallelFor(bx, modes,
[=] AMREX_GPU_DEVICE(int i, int j, int k, int mode) noexcept
{
// Calculate norm of vector
int const ir = i + nr*mode;
amrex::Real const kr = kr_arr[ir];
amrex::Real const kz = modified_kz[j];
amrex::Real const k_norm = std::sqrt(kr*kr + kz*kz);
// Calculate coefficients
constexpr amrex::Real c = PhysConst::c;
constexpr amrex::Real ep0 = PhysConst::ep0;
if (k_norm != 0){
C(i,j,k,mode) = std::cos(c*k_norm*dt);
S_ck(i,j,k,mode) = std::sin(c*k_norm*dt)/(c*k_norm);
X1(i,j,k,mode) = (1._rt - C(i,j,k,mode))/(ep0 * c*c * k_norm*k_norm);
X2(i,j,k,mode) = (1._rt - S_ck(i,j,k,mode)/dt)/(ep0 * k_norm*k_norm);
X3(i,j,k,mode) = (C(i,j,k,mode) - S_ck(i,j,k,mode)/dt)/(ep0 * k_norm*k_norm);
} else { // Handle k_norm = 0, by using the analytical limit
C(i,j,k,mode) = 1._rt;
S_ck(i,j,k,mode) = dt;
X1(i,j,k,mode) = 0.5_rt * dt*dt / ep0;
X2(i,j,k,mode) = c*c * dt*dt / (6._rt*ep0);
X3(i,j,k,mode) = - c*c * dt*dt / (3._rt*ep0);
}
});
}
}
void
PsatdAlgorithmRZ::CurrentCorrection (const int lev,
SpectralFieldDataRZ& field_data,
std::array<std::unique_ptr<amrex::MultiFab>,3>& current,
const std::unique_ptr<amrex::MultiFab>& rho)
{
// Profiling
WARPX_PROFILE( "PsatdAlgorithmRZ::CurrentCorrection" );
using Idx = SpectralFieldIndex;
// Forward Fourier transform of J and rho
field_data.ForwardTransform( lev,
*current[0], Idx::Jx,
*current[1], Idx::Jy);
field_data.ForwardTransform( lev, *current[2], Idx::Jz, 0);
field_data.ForwardTransform( lev, *rho, Idx::rho_old, 0 );
field_data.ForwardTransform( lev, *rho, Idx::rho_new, 1 );
// Loop over boxes
for (amrex::MFIter mfi(field_data.fields); mfi.isValid(); ++mfi){
amrex::Box const & bx = field_data.fields[mfi].box();
// Extract arrays for the fields to be updated
amrex::Array4<Complex> fields = field_data.fields[mfi].array();
// Extract pointers for the k vectors
auto const & kr_modes = field_data.getKrArray(mfi);
amrex::Real const* kr_arr = kr_modes.dataPtr();
amrex::Real const* modified_kz_arr = modified_kz_vec[mfi].dataPtr();
int const nr = bx.length(0);
// Local copy of member variables before GPU loop
amrex::Real const dt = m_dt;
// Loop over indices within one box
int const modes = field_data.n_rz_azimuthal_modes;
ParallelFor(bx, modes,
[=] AMREX_GPU_DEVICE(int i, int j, int k, int mode) noexcept
{
// All of the fields of each mode are grouped together
auto const Jp_m = Idx::Jx + Idx::n_fields*mode;
auto const Jm_m = Idx::Jy + Idx::n_fields*mode;
auto const Jz_m = Idx::Jz + Idx::n_fields*mode;
auto const rho_old_m = Idx::rho_old + Idx::n_fields*mode;
auto const rho_new_m = Idx::rho_new + Idx::n_fields*mode;
// Shortcuts for the values of J and rho
Complex const Jp = fields(i,j,k,Jp_m);
Complex const Jm = fields(i,j,k,Jm_m);
Complex const Jz = fields(i,j,k,Jz_m);
Complex const rho_old = fields(i,j,k,rho_old_m);
Complex const rho_new = fields(i,j,k,rho_new_m);
// k vector values, and coefficients
// The k values for each mode are grouped together
int const ir = i + nr*mode;
amrex::Real const kr = kr_arr[ir];
amrex::Real const kz = modified_kz_arr[j];
amrex::Real const k_norm2 = kr*kr + kz*kz;
constexpr Complex I = Complex{0._rt,1._rt};
// Correct J
if ( k_norm2 != 0._rt )
{
Complex const F = - ((rho_new - rho_old)/dt + I*kz*Jz + kr*(Jp - Jm))/k_norm2;
fields(i,j,k,Jp_m) += +0.5_rt*kr*F;
fields(i,j,k,Jm_m) += -0.5_rt*kr*F;
fields(i,j,k,Jz_m) += -I*kz*F;
}
});
}
// Backward Fourier transform of J
field_data.BackwardTransform( lev,
*current[0], Idx::Jx,
*current[1], Idx::Jy);
field_data.BackwardTransform( lev,
*current[2], Idx::Jz, 0 );
}
void
PsatdAlgorithmRZ::VayDeposition (const int lev /**/,
SpectralFieldDataRZ& /*field_data*/,
std::array<std::unique_ptr<amrex::MultiFab>,3>& /*current*/)
{
amrex::Abort("Vay deposition not implemented in RZ geometry");
}
|
