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
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
|
/* Copyright 2020 Remi Lehe
*
* This file is part of WarpX.
*
* License: BSD-3-Clause-LBNL
*/
#ifndef WARPX_FINITE_DIFFERENCE_SOLVER_H_
#define WARPX_FINITE_DIFFERENCE_SOLVER_H_
#include "EmbeddedBoundary/WarpXFaceInfoBox_fwd.H"
#include "FiniteDifferenceSolver_fwd.H"
#include "BoundaryConditions/PML_fwd.H"
#include "Evolve/WarpXDtType.H"
#include "HybridPICModel/HybridPICModel_fwd.H"
#include "MacroscopicProperties/MacroscopicProperties_fwd.H"
#include <AMReX_GpuContainers.H>
#include <AMReX_REAL.H>
#include <AMReX_BaseFwd.H>
#include <array>
#include <memory>
/**
* \brief Top-level class for the electromagnetic finite-difference solver
*
* Stores the coefficients of the finite-difference stencils,
* and has member functions to update fields over one time step.
*/
class FiniteDifferenceSolver
{
public:
// Constructor
/** \brief Initialize the finite-difference Maxwell solver (for a given refinement level)
*
* This function initializes the stencil coefficients for the chosen finite-difference algorithm
*
* \param fdtd_algo Identifies the chosen algorithm, as defined in WarpXAlgorithmSelection.H
* \param cell_size Cell size along each dimension, for the chosen refinement level
* \param grid_type Whether the solver is applied to a collocated or staggered grid
*/
FiniteDifferenceSolver (
int fdtd_algo,
std::array<amrex::Real,3> cell_size,
short grid_type );
void EvolveB ( std::array< std::unique_ptr<amrex::MultiFab>, 3 >& Bfield,
std::array< std::unique_ptr<amrex::MultiFab>, 3 > const& Efield,
std::unique_ptr<amrex::MultiFab> const& Gfield,
std::array< std::unique_ptr<amrex::MultiFab>, 3 > const& face_areas,
std::array< std::unique_ptr<amrex::MultiFab>, 3 > const& area_mod,
std::array< std::unique_ptr<amrex::MultiFab>, 3 >& ECTRhofield,
std::array< std::unique_ptr<amrex::MultiFab>, 3 >& Venl,
std::array< std::unique_ptr<amrex::iMultiFab>, 3 >& flag_info_cell,
std::array< std::unique_ptr<amrex::LayoutData<FaceInfoBox> >, 3 >& borrowing,
int lev, amrex::Real dt );
void EvolveE ( std::array< std::unique_ptr<amrex::MultiFab>, 3 >& Efield,
std::array< std::unique_ptr<amrex::MultiFab>, 3 > const& Bfield,
std::array< std::unique_ptr<amrex::MultiFab>, 3 > const& Jfield,
std::array< std::unique_ptr<amrex::MultiFab>, 3 > const& edge_lengths,
std::array< std::unique_ptr<amrex::MultiFab>, 3 > const& face_areas,
std::array< std::unique_ptr<amrex::MultiFab>, 3 >& ECTRhofield,
std::unique_ptr<amrex::MultiFab> const& Ffield,
int lev, amrex::Real dt );
void EvolveF ( std::unique_ptr<amrex::MultiFab>& Ffield,
std::array< std::unique_ptr<amrex::MultiFab>, 3 > const& Efield,
std::unique_ptr<amrex::MultiFab> const& rhofield,
int rhocomp,
amrex::Real dt );
void EvolveG (std::unique_ptr<amrex::MultiFab>& Gfield,
std::array<std::unique_ptr<amrex::MultiFab>,3> const& Bfield,
amrex::Real dt);
void EvolveECTRho ( std::array< std::unique_ptr<amrex::MultiFab>, 3 > const& Efield,
std::array< std::unique_ptr<amrex::MultiFab>, 3 > const& edge_lengths,
std::array< std::unique_ptr<amrex::MultiFab>, 3 > const& face_areas,
std::array< std::unique_ptr<amrex::MultiFab>, 3 >& ECTRhofield,
int lev );
void ApplySilverMuellerBoundary(
std::array< std::unique_ptr<amrex::MultiFab>, 3 >& Efield,
std::array< std::unique_ptr<amrex::MultiFab>, 3 >& Bfield,
amrex::Box domain_box,
amrex::Real dt,
amrex::Vector<int> field_boundary_lo,
amrex::Vector<int> field_boundary_hi);
void ComputeDivE ( const std::array<std::unique_ptr<amrex::MultiFab>,3>& Efield,
amrex::MultiFab& divE );
/**
* \brief Macroscopic E-update for non-vacuum medium using the user-selected
* finite-difference algorithm and macroscopic sigma-method defined in
* WarpXAlgorithmSelection.H
*
* \param[out] Efield vector of electric field MultiFabs updated at a given level
* \param[in] Bfield vector of magnetic field MultiFabs at a given level
* \param[in] Jfield vector of current density MultiFabs at a given level
* \param[in] edge_lengths length of edges along embedded boundaries
* \param[in] dt timestep of the simulation
* \param[in] macroscopic_properties contains user-defined properties of the medium.
*/
void MacroscopicEvolveE ( std::array< std::unique_ptr<amrex::MultiFab>, 3>& Efield,
std::array< std::unique_ptr<amrex::MultiFab>, 3> const& Bfield,
std::array< std::unique_ptr<amrex::MultiFab>, 3 > const& Jfield,
std::array< std::unique_ptr<amrex::MultiFab>, 3 > const& edge_lengths,
amrex::Real dt,
std::unique_ptr<MacroscopicProperties> const& macroscopic_properties);
void EvolveBPML ( std::array< amrex::MultiFab*, 3 > Bfield,
std::array< amrex::MultiFab*, 3 > Efield,
amrex::Real dt,
bool dive_cleaning);
void EvolveEPML ( std::array< amrex::MultiFab*, 3 > Efield,
std::array< amrex::MultiFab*, 3 > Bfield,
std::array< amrex::MultiFab*, 3 > Jfield,
std::array< amrex::MultiFab*, 3 > edge_lengths,
amrex::MultiFab* Ffield,
MultiSigmaBox const& sigba,
amrex::Real dt, bool pml_has_particles );
void EvolveFPML ( amrex::MultiFab* Ffield,
std::array< amrex::MultiFab*, 3 > Efield,
amrex::Real dt );
/**
* \brief E-update in the hybrid PIC algorithm as described in
* Winske et al. (2003) Eq. 10.
* https://link.springer.com/chapter/10.1007/3-540-36530-3_8
*
* \param[out] Efield vector of electric field MultiFabs updated at a given level
* \param[in] Jfield vector of total current MultiFabs at a given level
* \param[in] Jifield vector of ion current density MultiFabs at a given level
* \param[in] Bfield vector of magnetic field MultiFabs at a given level
* \param[in] rhofield scalar ion charge density Multifab at a given level
* \param[in] Pefield scalar electron pressure MultiFab at a given level
* \param[in] edge_lengths length of edges along embedded boundaries
* \param[in] lev level number for the calculation
* \param[in] hybrid_pic_model instance of the hybrid-PIC model
* \param[in] include_resistivity_term boolean flag for whether to include resistivity
*/
void HybridPICSolveE ( std::array< std::unique_ptr<amrex::MultiFab>, 3>& Efield,
std::array< std::unique_ptr<amrex::MultiFab>, 3>& Jfield,
std::array< std::unique_ptr<amrex::MultiFab>, 3 > const& Jifield,
std::array< std::unique_ptr<amrex::MultiFab>, 3> const& Bfield,
std::unique_ptr<amrex::MultiFab> const& rhofield,
std::unique_ptr<amrex::MultiFab> const& Pefield,
std::array< std::unique_ptr<amrex::MultiFab>, 3 > const& edge_lengths,
int lev, HybridPICModel const* hybrid_pic_model,
bool include_resistivity_term );
/**
* \brief Calculation of total current using Ampere's law (without
* displacement current): J = (curl x B) / mu0.
*
* \param[out] Jfield vector of current MultiFabs at a given level
* \param[in] Bfield vector of magnetic field MultiFabs at a given level
* \param[in] edge_lengths length of edges along embedded boundaries
* \param[in] lev level number for the calculation
*/
void CalculateCurrentAmpere (
std::array< std::unique_ptr<amrex::MultiFab>, 3>& Jfield,
std::array< std::unique_ptr<amrex::MultiFab>, 3> const& Bfield,
std::array< std::unique_ptr<amrex::MultiFab>, 3 > const& edge_lengths,
int lev );
private:
int m_fdtd_algo;
short m_grid_type;
#ifdef WARPX_DIM_RZ
amrex::Real m_dr, m_rmin;
int m_nmodes;
// host-only
amrex::Vector<amrex::Real> m_h_stencil_coefs_r, m_h_stencil_coefs_z;
// device copy after init
amrex::Gpu::DeviceVector<amrex::Real> m_stencil_coefs_r;
amrex::Gpu::DeviceVector<amrex::Real> m_stencil_coefs_z;
#else
// host-only
amrex::Vector<amrex::Real> m_h_stencil_coefs_x, m_h_stencil_coefs_y, m_h_stencil_coefs_z;
// device copy after init
amrex::Gpu::DeviceVector<amrex::Real> m_stencil_coefs_x;
amrex::Gpu::DeviceVector<amrex::Real> m_stencil_coefs_y;
amrex::Gpu::DeviceVector<amrex::Real> m_stencil_coefs_z;
#endif
public:
// The member functions below contain extended __device__ lambda.
// In order to compile with nvcc, they need to be public.
#ifdef WARPX_DIM_RZ
template< typename T_Algo >
void EvolveBCylindrical (
std::array< std::unique_ptr<amrex::MultiFab>, 3 >& Bfield,
std::array< std::unique_ptr<amrex::MultiFab>, 3 > const& Efield,
int lev,
amrex::Real dt );
template< typename T_Algo >
void EvolveECylindrical (
std::array< std::unique_ptr<amrex::MultiFab>, 3 >& Efield,
std::array< std::unique_ptr<amrex::MultiFab>, 3 > const& Bfield,
std::array< std::unique_ptr<amrex::MultiFab>, 3 > const& Jfield,
std::unique_ptr<amrex::MultiFab> const& Ffield,
int lev,
amrex::Real dt );
template< typename T_Algo >
void EvolveFCylindrical (
std::unique_ptr<amrex::MultiFab>& Ffield,
std::array< std::unique_ptr<amrex::MultiFab>, 3 > const& Efield,
std::unique_ptr<amrex::MultiFab> const& rhofield,
int rhocomp,
amrex::Real dt );
template< typename T_Algo >
void ComputeDivECylindrical (
const std::array<std::unique_ptr<amrex::MultiFab>,3>& Efield,
amrex::MultiFab& divE );
template<typename T_Algo>
void HybridPICSolveECylindrical (
std::array< std::unique_ptr<amrex::MultiFab>, 3>& Efield,
std::array< std::unique_ptr<amrex::MultiFab>, 3> const& Jfield,
std::array< std::unique_ptr<amrex::MultiFab>, 3> const& Jifield,
std::array< std::unique_ptr<amrex::MultiFab>, 3> const& Bfield,
std::unique_ptr<amrex::MultiFab> const& rhofield,
std::unique_ptr<amrex::MultiFab> const& Pefield,
std::array< std::unique_ptr<amrex::MultiFab>, 3 > const& edge_lengths,
int lev, HybridPICModel const* hybrid_pic_model,
bool include_resistivity_term );
template<typename T_Algo>
void CalculateCurrentAmpereCylindrical (
std::array< std::unique_ptr<amrex::MultiFab>, 3 >& Jfield,
std::array< std::unique_ptr<amrex::MultiFab>, 3 > const& Bfield,
std::array< std::unique_ptr<amrex::MultiFab>, 3 > const& edge_lengths,
int lev
);
#else
template< typename T_Algo >
void EvolveBCartesian (
std::array< std::unique_ptr<amrex::MultiFab>, 3 >& Bfield,
std::array< std::unique_ptr<amrex::MultiFab>, 3 > const& Efield,
std::unique_ptr<amrex::MultiFab> const& Gfield,
int lev, amrex::Real dt );
template< typename T_Algo >
void EvolveECartesian (
std::array< std::unique_ptr<amrex::MultiFab>, 3 >& Efield,
std::array< std::unique_ptr<amrex::MultiFab>, 3 > const& Bfield,
std::array< std::unique_ptr<amrex::MultiFab>, 3 > const& Jfield,
std::array< std::unique_ptr<amrex::MultiFab>, 3 > const& edge_lengths,
std::unique_ptr<amrex::MultiFab> const& Ffield,
int lev, amrex::Real dt );
template< typename T_Algo >
void EvolveFCartesian (
std::unique_ptr<amrex::MultiFab>& Ffield,
std::array< std::unique_ptr<amrex::MultiFab>, 3 > const& Efield,
std::unique_ptr<amrex::MultiFab> const& rhofield,
int rhocomp,
amrex::Real dt );
template< typename T_Algo >
void EvolveGCartesian (
std::unique_ptr<amrex::MultiFab>& Gfield,
std::array<std::unique_ptr<amrex::MultiFab>,3> const& Bfield,
amrex::Real dt);
void EvolveRhoCartesianECT (
std::array< std::unique_ptr<amrex::MultiFab>, 3 > const& Efield,
std::array< std::unique_ptr<amrex::MultiFab>, 3 > const& edge_lengths,
std::array< std::unique_ptr<amrex::MultiFab>, 3 > const& face_areas,
std::array< std::unique_ptr<amrex::MultiFab>, 3 >& ECTRhofield, int lev);
void EvolveBCartesianECT (
std::array< std::unique_ptr<amrex::MultiFab>, 3 >& Bfield,
std::array< std::unique_ptr<amrex::MultiFab>, 3 > const& face_areas,
std::array< std::unique_ptr<amrex::MultiFab>, 3 > const& area_mod,
std::array< std::unique_ptr<amrex::MultiFab>, 3 >& ECTRhofield,
std::array< std::unique_ptr<amrex::MultiFab>, 3 >& Venl,
std::array< std::unique_ptr<amrex::iMultiFab>, 3 >& flag_info_cell,
std::array< std::unique_ptr<amrex::LayoutData<FaceInfoBox> >, 3 >& borrowing,
int lev, amrex::Real dt
);
template< typename T_Algo >
void ComputeDivECartesian (
const std::array<std::unique_ptr<amrex::MultiFab>,3>& Efield,
amrex::MultiFab& divE );
template< typename T_Algo, typename T_MacroAlgo >
void MacroscopicEvolveECartesian (
std::array< std::unique_ptr< amrex::MultiFab>, 3>& Efield,
std::array< std::unique_ptr< amrex::MultiFab>, 3> const& Bfield,
std::array< std::unique_ptr< amrex::MultiFab>, 3> const& Jfield,
std::array< std::unique_ptr<amrex::MultiFab>, 3 > const& edge_lengths,
amrex::Real dt,
std::unique_ptr<MacroscopicProperties> const& macroscopic_properties);
template< typename T_Algo >
void EvolveBPMLCartesian (
std::array< amrex::MultiFab*, 3 > Bfield,
std::array< amrex::MultiFab*, 3 > Efield,
amrex::Real dt,
bool dive_cleaning);
template< typename T_Algo >
void EvolveEPMLCartesian (
std::array< amrex::MultiFab*, 3 > Efield,
std::array< amrex::MultiFab*, 3 > Bfield,
std::array< amrex::MultiFab*, 3 > Jfield,
std::array< amrex::MultiFab*, 3 > edge_lengths,
amrex::MultiFab* Ffield,
MultiSigmaBox const& sigba,
amrex::Real dt, bool pml_has_particles );
template< typename T_Algo >
void EvolveFPMLCartesian ( amrex::MultiFab* Ffield,
std::array< amrex::MultiFab*, 3 > Efield,
amrex::Real dt );
template<typename T_Algo>
void HybridPICSolveECartesian (
std::array< std::unique_ptr<amrex::MultiFab>, 3>& Efield,
std::array< std::unique_ptr<amrex::MultiFab>, 3> const& Jfield,
std::array< std::unique_ptr<amrex::MultiFab>, 3> const& Jifield,
std::array< std::unique_ptr<amrex::MultiFab>, 3> const& Bfield,
std::unique_ptr<amrex::MultiFab> const& rhofield,
std::unique_ptr<amrex::MultiFab> const& Pefield,
std::array< std::unique_ptr<amrex::MultiFab>, 3 > const& edge_lengths,
int lev, HybridPICModel const* hybrid_pic_model,
bool include_resistivity_term );
template<typename T_Algo>
void CalculateCurrentAmpereCartesian (
std::array< std::unique_ptr<amrex::MultiFab>, 3 >& Jfield,
std::array< std::unique_ptr<amrex::MultiFab>, 3 > const& Bfield,
std::array< std::unique_ptr<amrex::MultiFab>, 3 > const& edge_lengths,
int lev
);
#endif
};
#endif // WARPX_FINITE_DIFFERENCE_SOLVER_H_
|