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
|
#include "Utils/WarpXAlgorithmSelection.H"
#include "FiniteDifferenceSolver.H"
#ifdef WARPX_DIM_RZ
// currently works only for 3D
#else
# include "FiniteDifferenceAlgorithms/CartesianYeeAlgorithm.H"
# include "FiniteDifferenceAlgorithms/CartesianCKCAlgorithm.H"
#endif
#include "Utils/WarpXConst.H"
#include <AMReX_Gpu.H>
#include <WarpX.H>
using namespace amrex;
void FiniteDifferenceSolver::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,
amrex::Real const dt, std::unique_ptr<MacroscopicProperties> const& macroscopic_properties ) {
// Select algorithm (The choice of algorithm is a runtime option,
// but we compile code for each algorithm, using templates)
#ifdef WARPX_DIM_RZ
amrex::Abort("currently macro E-push does not work for RZ");
#else
if (m_do_nodal) {
amrex::Abort(" macro E-push does not work for nodal ");
} else if (m_fdtd_algo == MaxwellSolverAlgo::Yee) {
MacroscopicEvolveECartesian <CartesianYeeAlgorithm> ( Efield, Bfield, Jfield, dt,
macroscopic_properties );
} else if (m_fdtd_algo == MaxwellSolverAlgo::CKC) {
// Note : EvolveE is the same for CKC and Yee.
// In the templated Yee and CKC calls, the core operations for EvolveE is the same.
MacroscopicEvolveECartesian <CartesianCKCAlgorithm> ( Efield, Bfield, Jfield, dt,
macroscopic_properties );
} else {
amrex::Abort("Unknown algorithm");
}
#endif
}
#ifndef WARPX_DIM_RZ
template<typename T_Algo>
void FiniteDifferenceSolver::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,
amrex::Real const dt, std::unique_ptr<MacroscopicProperties> const& macroscopic_properties ) {
const int ¯oscopic_solver_algo = WarpX::macroscopic_solver_algo;
Real sigma = macroscopic_properties->sigma();
Real const mu = macroscopic_properties->mu();
Real const epsilon = macroscopic_properties->epsilon();
Real alpha = 0._rt;
Real beta = 0._rt;
Real fac1 = 0._rt;
Real inv_fac = 0._rt;
if (macroscopic_solver_algo == MacroscopicSolverAlgo::BackwardEuler) {
fac1 = sigma * dt / epsilon;
inv_fac = 1._rt / ( 1._rt + fac1);
alpha = inv_fac;
beta = dt * inv_fac / epsilon;
} else if (macroscopic_solver_algo == MacroscopicSolverAlgo::LaxWendroff) {
fac1 = 0.5_rt * sigma * dt / epsilon;
inv_fac = 1._rt / ( 1._rt + fac1);
alpha = (1.0_rt - fac1) * inv_fac;
beta = dt * inv_fac / epsilon;
}
// Loop through the grids, and over the tiles within each grid
#ifdef _OPENMP
#pragma omp parallel if (amrex::Gpu::notInLaunchRegion())
#endif
for ( MFIter mfi(*Efield[0], TilingIfNotGPU()); mfi.isValid(); ++mfi ) {
// Extract field data for this grid/tile
Array4<Real> const& Ex = Efield[0]->array(mfi);
Array4<Real> const& Ey = Efield[1]->array(mfi);
Array4<Real> const& Ez = Efield[2]->array(mfi);
Array4<Real> const& Bx = Bfield[0]->array(mfi);
Array4<Real> const& By = Bfield[1]->array(mfi);
Array4<Real> const& Bz = Bfield[2]->array(mfi);
// Extract stencil coefficients
Real const * const AMREX_RESTRICT coefs_x = m_stencil_coefs_x.dataPtr();
int const n_coefs_x = m_stencil_coefs_x.size();
Real const * const AMREX_RESTRICT coefs_y = m_stencil_coefs_y.dataPtr();
int const n_coefs_y = m_stencil_coefs_y.size();
Real const * const AMREX_RESTRICT coefs_z = m_stencil_coefs_z.dataPtr();
int const n_coefs_z = m_stencil_coefs_z.size();
// Extract tileboxes for which to loop
Box const& tex = mfi.tilebox(Efield[0]->ixType().toIntVect());
Box const& tey = mfi.tilebox(Efield[1]->ixType().toIntVect());
Box const& tez = mfi.tilebox(Efield[2]->ixType().toIntVect());
// Loop over the cells and update the fields
amrex::ParallelFor(tex, tey, tez,
[=] AMREX_GPU_DEVICE (int i, int j, int k){
Ex(i, j, k) = alpha * Ex(i, j, k) + (beta/mu)
* ( - T_Algo::DownwardDz(By, coefs_z, n_coefs_z, i, j, k)
+ T_Algo::DownwardDy(Bz, coefs_y, n_coefs_y, i, j, k));
},
[=] AMREX_GPU_DEVICE (int i, int j, int k){
Ey(i, j, k) = alpha * Ey(i, j, k) + (beta/mu)
* ( - T_Algo::DownwardDx(Bz, coefs_x, n_coefs_x, i, j, k)
+ T_Algo::DownwardDz(Bx, coefs_z, n_coefs_z, i, j, k));
},
[=] AMREX_GPU_DEVICE (int i, int j, int k){
Ez(i, j, k) = alpha * Ez(i, j, k) + (beta/mu)
* ( - T_Algo::DownwardDy(Bx, coefs_y, n_coefs_y, i, j, k)
+ T_Algo::DownwardDx(By, coefs_x, n_coefs_x, i, j, k));
}
);
// update E using J, if source currents are specified.
if (Jfield[0]) {
Array4<Real> const& jx = Jfield[0]->array(mfi);
Array4<Real> const& jy = Jfield[1]->array(mfi);
Array4<Real> const& jz = Jfield[2]->array(mfi);
amrex::ParallelFor(tex, tey, tez,
[=] AMREX_GPU_DEVICE (int i, int j, int k) {
Ex(i, j, k) += -beta * jx(i, j, k);
},
[=] AMREX_GPU_DEVICE (int i, int j, int k) {
Ey(i, j, k) += -beta * jy(i, j, k);
},
[=] AMREX_GPU_DEVICE (int i, int j, int k) {
Ez(i, j, k) += -beta * jz(i, j, k);
}
);
}
}
}
#endif // corresponds to ifndef WARPX_DIM_RZ
|
