31 void calc_M(amrex::Real arr [], amrex::Real ex, amrex::Real ey, amrex::Real ez,
32 amrex::Real bx, amrex::Real by, amrex::Real bz,
33 amrex::Real
xi_c2, amrex::Real c2)
35 using namespace amrex;
36 const Real ee = ex*ex+ey*ey+ez*ez;
37 const Real bb_c2 = c2*(bx*bx+by*by+bz*bz);
38 const Real eb = ex*bx+ey*by+ez*bz;
39 arr[0] = -2._rt*
xi_c2*( 2._rt*bx*(ee-bb_c2) - 7._rt*ex*eb );
40 arr[1] = -2._rt*
xi_c2*( 2._rt*by*(ee-bb_c2) - 7._rt*ey*eb );
41 arr[2] = -2._rt*
xi_c2*( 2._rt*bz*(ee-bb_c2) - 7._rt*ez*eb );
88 const amrex::Real
dz,
const amrex::Real
dt,
const amrex::Real
xi_c2)
91 using namespace amrex;
96 constexpr amrex::Real c2i = 1._rt/c2;
97 const amrex::Real dxi = 1._rt/
dx;
98 const amrex::Real dzi = 1._rt/
dz;
100 #if defined(WARPX_DIM_3D)
101 const amrex::Real dyi = 1._rt/dy;
105 amrex::Real Mpx [3] = {0._rt,0._rt,0._rt};
106 amrex::Real Mnx [3] = {0._rt,0._rt,0._rt};
107 amrex::Real Mpy [3] = {0._rt,0._rt,0._rt};
108 amrex::Real Mny [3] = {0._rt,0._rt,0._rt};
109 amrex::Real Mpz [3] = {0._rt,0._rt,0._rt};
110 amrex::Real Mnz [3] = {0._rt,0._rt,0._rt};
114 calc_M(Mpx, tmpEx(j+1,k,l), tmpEy(j+1,k,l), tmpEz(j+1,k,l),
115 Bx(j+1,k,l), By(j+1,k,l), Bz(j+1,k,l),
xi_c2, c2);
116 calc_M(Mnx, tmpEx(j-1,k,l), tmpEy(j-1,k,l), tmpEz(j-1,k,l),
117 Bx(j-1,k,l), By(j-1,k,l), Bz(j-1,k,l),
xi_c2, c2);
118 calc_M(Mpy, tmpEx(j,k+1,l), tmpEy(j,k+1,l), tmpEz(j,k+1,l),
119 Bx(j,k+1,l), By(j,k+1,l), Bz(j,k+1,l),
xi_c2, c2);
120 calc_M(Mny, tmpEx(j,k-1,l), tmpEy(j,k-1,l), tmpEz(j,k-1,l),
121 Bx(j,k-1,l), By(j,k-1,l), Bz(j,k-1,l),
xi_c2, c2);
122 calc_M(Mpz, tmpEx(j,k,l+1), tmpEy(j,k,l+1), tmpEz(j,k,l+1),
123 Bx(j,k,l+1), By(j,k,l+1), Bz(j,k,l+1),
xi_c2, c2);
124 calc_M(Mnz, tmpEx(j,k,l-1), tmpEy(j,k,l-1), tmpEz(j,k,l-1),
125 Bx(j,k,l-1), By(j,k,l-1), Bz(j,k,l-1),
xi_c2, c2);
129 const amrex::Real VxM[3] = {
130 0.5_rt*( (Mpy[2]-Mny[2])*dyi - (Mpz[1]-Mnz[1])*dzi ),
131 0.5_rt*( (Mpz[0]-Mnz[0])*dzi - (Mpx[2]-Mnx[2])*dxi ),
132 0.5_rt*( (Mpx[1]-Mnx[1])*dxi - (Mpy[0]-Mny[0])*dyi ),
135 const amrex::Real VxE[3] = {
136 0.5_rt*( (tmpEz(j,k+1,l)-tmpEz(j,k-1,l) )*dyi - (tmpEy(j,k,l+1)-tmpEy(j,k,l-1) )*dzi ),
137 0.5_rt*( (tmpEx(j,k,l+1)-tmpEx(j,k,l-1) )*dzi - (tmpEz(j+1,k,l)-tmpEz(j-1,k,l) )*dxi ),
138 0.5_rt*( (tmpEy(j+1,k,l)-tmpEy(j-1,k,l) )*dxi - (tmpEx(j,k+1,l)-tmpEx(j,k-1,l) )*dyi ),
141 const amrex::Real VxB[3] = {
142 0.5_rt*( (Bz(j,k+1,l)-Bz(j,k-1,l) )*dyi - (By(j,k,l+1)-By(j,k,l-1) )*dzi ),
143 0.5_rt*( (Bx(j,k,l+1)-Bx(j,k,l-1) )*dzi - (Bz(j+1,k,l)-Bz(j-1,k,l) )*dxi ),
144 0.5_rt*( (By(j+1,k,l)-By(j-1,k,l) )*dxi - (Bx(j,k+1,l)-Bx(j,k-1,l) )*dyi ),
149 const amrex::Real ex = tmpEx(j,k,l);
150 const amrex::Real ey = tmpEy(j,k,l);
151 const amrex::Real ez = tmpEz(j,k,l);
152 const amrex::Real bx = Bx(j,k,l);
153 const amrex::Real by = By(j,k,l);
154 const amrex::Real bz = Bz(j,k,l);
158 const amrex::Real ee = ex*ex + ey*ey + ez*ez;
159 const amrex::Real bb = bx*bx + by*by + bz*bz;
160 const amrex::Real eb = ex*bx + ey*by + ez*bz;
161 const amrex::Real EVxE = ex*VxE[0] + ey*VxE[1] + ez*VxE[2];
162 const amrex::Real BVxE = bx*VxE[0] + by*VxE[1] + bz*VxE[2];
163 const amrex::Real EVxB = ex*VxB[0] + ey*VxB[1] + ez*VxB[2];
164 const amrex::Real BVxB = bx*VxB[0] + by*VxB[1] + bz*VxB[2];
165 const amrex::Real Emu0J = ex*mu0jx + ey*mu0jy + ez*mu0jz;
166 const amrex::Real Bmu0J = bx*mu0jx + by*mu0jy + bz*mu0jz;
170 const amrex::Real Alpha[3] = {
171 2._rt*
xi_c2*( -7._rt*bx*EVxE - 7._rt*VxE[0]*eb + 4._rt*ex*BVxE ) + VxM[0],
172 2._rt*
xi_c2*( -7._rt*by*EVxE - 7._rt*VxE[1]*eb + 4._rt*ey*BVxE ) + VxM[1],
173 2._rt*
xi_c2*( -7._rt*bz*EVxE - 7._rt*VxE[2]*eb + 4._rt*ez*BVxE ) + VxM[2]
176 const amrex::Real Omega[3] = {
177 Alpha[0] + 2._rt*
xi_c2*( 4._rt*ex*(EVxB + Emu0J) + 2._rt*(VxB[0] + mu0jx)*( ee - c2*bb )
178 + 7._rt*c2*bx*(BVxB +Bmu0J) ),
179 Alpha[1] + 2._rt*
xi_c2*( 4._rt*ey*(EVxB + Emu0J) + 2._rt*(VxB[1] + mu0jy)*( ee - c2*bb )
180 + 7._rt*c2*by*(BVxB + Bmu0J) ),
181 Alpha[2] + 2._rt*
xi_c2*( 4._rt*ez*(EVxB + Emu0J) + 2._rt*(VxB[2] + mu0jz)*( ee - c2*bb )
182 + 7._rt*c2*bz*(BVxB + Bmu0J) )
187 const amrex::Real a00 =
beta +
xi_c2*( 8._rt*c2i*ex*ex + 14._rt*bx*bx );
189 const amrex::Real a11 =
beta +
xi_c2*( 8._rt*c2i*ey*ey + 14._rt*by*by );
191 const amrex::Real a22 =
beta +
xi_c2*( 8._rt*c2i*ez*ez + 14._rt*bz*bz );
193 const amrex::Real a01 =
xi_c2*( 2._rt*c2i*ex*ey + 14*bx*by );
195 const amrex::Real a02 =
xi_c2*( 2._rt*c2i*ex*ez + 14._rt*bx*bz );
197 const amrex::Real a12 =
xi_c2*( 2._rt*c2i*ez*ey + 14._rt*bz*by );
199 const amrex::Real detA = a00*( a11*a22 - a12*a12 ) - a01*( a01*a22 - a02*a12 )+a02*( a01*a12 - a02*a11 );
203 const amrex::Real invAx[3] = {a22*a11 - a12*a12, a12*a02 - a22*a01, a12*a01 - a11*a02};
205 const amrex::Real invAy[3] = {a02*a12 - a22*a01, a00*a22 - a02*a02, a01*a02 - a12*a00};
207 const amrex::Real invAz[3] = {a12*a01 - a02*a11, a02*a01 - a12*a00, a11*a00 - a01*a01};
211 const amrex::Real dEx = (-1._rt/detA)*(invAx[0]*Omega[0] +
215 const amrex::Real dEy = (-1._rt/detA)*(invAy[0]*Omega[0] +
219 const amrex::Real dEz = (-1._rt/detA)*(invAz[0]*Omega[0] +
225 Ex(j,k,l) = Ex(j,k,l) + 0.5_rt*
dt*dEx;
227 Ey(j,k,l) = Ey(j,k,l) + 0.5_rt*
dt*dEy;
229 Ez(j,k,l) = Ez(j,k,l) + 0.5_rt*
dt*dEz;
237 amrex::Real Mpx [3] = {0._rt,0._rt,0._rt};
238 amrex::Real Mnx [3] = {0._rt,0._rt,0._rt};
239 amrex::Real Mpz [3] = {0._rt,0._rt,0._rt};
240 amrex::Real Mnz [3] = {0._rt,0._rt,0._rt};
244 calc_M(Mpx, tmpEx(j+1,k,0), tmpEy(j+1,k,0), tmpEz(j+1,k,0),
245 Bx(j+1,k,0), By(j+1,k,0), Bz(j+1,k,0),
xi_c2, c2);
246 calc_M(Mnx, tmpEx(j-1,k,0), tmpEy(j-1,k,0), tmpEz(j-1,k,0),
247 Bx(j-1,k,0), By(j-1,k,0), Bz(j-1,k,0),
xi_c2, c2);
248 calc_M(Mpz, tmpEx(j,k+1,0), tmpEy(j,k+1,0), tmpEz(j,k+1,0),
249 Bx(j,k+1,0), By(j,k+1,0), Bz(j,k+1,0),
xi_c2, c2);
250 calc_M(Mnz, tmpEx(j,k-1,0), tmpEy(j,k-1,0), tmpEz(j,k-1,0),
251 Bx(j,k-1,0), By(j,k-1,0), Bz(j,k-1,0),
xi_c2, c2);
255 const amrex::Real VxM[3] = {
256 -0.5_rt*(Mpz[1]-Mnz[1])*dzi,
257 0.5_rt*( (Mpz[0]-Mnz[0])*dzi - (Mpx[2]-Mnx[2])*dxi ),
258 0.5_rt*(Mpx[1]-Mnx[1])*dxi,
261 const amrex::Real VxE[3] = {
262 -0.5_rt*(tmpEy(j,k+1,0)-tmpEy(j,k-1,0) )*dzi,
263 0.5_rt*( (tmpEx(j,k+1,0)-tmpEx(j,k-1,0) )*dzi - (tmpEz(j+1,k,0)-tmpEz(j-1,k,0) )*dxi ),
264 0.5_rt*(tmpEy(j+1,k,0)-tmpEy(j-1,k,0) )*dxi,
267 const amrex::Real VxB[3] = {
268 -0.5_rt*(By(j,k+1,0)-By(j,k-1,0) )*dzi,
269 0.5_rt*( (Bx(j,k+1,0)-Bx(j,k-1,0) )*dzi - (Bz(j+1,k,0)-Bz(j-1,k,0) )*dxi ),
270 0.5_rt*(By(j+1,k,0)-By(j-1,k,0) )*dxi,
275 const amrex::Real ex = tmpEx(j,k,0);
276 const amrex::Real ey = tmpEy(j,k,0);
277 const amrex::Real ez = tmpEz(j,k,0);
278 const amrex::Real bx = Bx(j,k,0);
279 const amrex::Real by = By(j,k,0);
280 const amrex::Real bz = Bz(j,k,0);
284 const amrex::Real ee = ex*ex + ey*ey + ez*ez;
285 const amrex::Real bb = bx*bx + by*by + bz*bz;
286 const amrex::Real eb = ex*bx + ey*by + ez*bz;
287 const amrex::Real EVxE = ex*VxE[0] + ey*VxE[1] + ez*VxE[2];
288 const amrex::Real BVxE = bx*VxE[0] + by*VxE[1] + bz*VxE[2];
289 const amrex::Real EVxB = ex*VxB[0] + ey*VxB[1] + ez*VxB[2];
290 const amrex::Real BVxB = bx*VxB[0] + by*VxB[1] + bz*VxB[2];
291 const amrex::Real Emu0J = ex*mu0jx + ey*mu0jy + ez*mu0jz;
292 const amrex::Real Bmu0J = bx*mu0jx + by*mu0jy + bz*mu0jz;
296 const amrex::Real Alpha[3] = {
297 2._rt*
xi_c2*( -7._rt*bx*EVxE - 7._rt*VxE[0]*eb + 4._rt*ex*BVxE ) + VxM[0],
298 2._rt*
xi_c2*( -7._rt*by*EVxE - 7._rt*VxE[1]*eb + 4._rt*ey*BVxE ) + VxM[1],
299 2._rt*
xi_c2*( -7._rt*bz*EVxE - 7._rt*VxE[2]*eb + 4._rt*ez*BVxE ) + VxM[2]
302 const amrex::Real Omega[3] = {
303 Alpha[0] + 2._rt*
xi_c2*( 4._rt*ex*(EVxB + Emu0J) + 2._rt*(VxB[0] + mu0jx)*( ee - c2*bb )
304 + 7._rt*c2*bx*(BVxB + Bmu0J) ),
305 Alpha[1] + 2._rt*
xi_c2*( 4._rt*ey*(EVxB + Emu0J) + 2._rt*(VxB[1] + mu0jy)*( ee - c2*bb )
306 + 7._rt*c2*by*(BVxB + Bmu0J) ),
307 Alpha[2] + 2._rt*
xi_c2*( 4._rt*ez*(EVxB +Emu0J) + 2._rt*(VxB[2] + mu0jz)*( ee - c2*bb )
308 + 7._rt*c2*bz*(BVxB + Bmu0J) )
313 const amrex::Real a00 =
beta +
xi_c2*( 8._rt*c2i*ex*ex + 14._rt*bx*bx );
315 const amrex::Real a11 =
beta +
xi_c2*( 8._rt*c2i*ey*ey + 14._rt*by*by );
317 const amrex::Real a22 =
beta +
xi_c2*( 8._rt*c2i*ez*ez + 14._rt*bz*bz );
319 const amrex::Real a01 =
xi_c2*( 2._rt*c2i*ex*ey + 14._rt*bx*by );
321 const amrex::Real a02 =
xi_c2*( 2._rt*c2i*ex*ez + 14._rt*bx*bz );
323 const amrex::Real a12 =
xi_c2*( 2._rt*c2i*ez*ey + 14._rt*bz*by );
325 const amrex::Real detA = a00*( a11*a22 - a12*a12 ) - a01*( a01*a22 - a02*a12 ) + a02*( a01*a12 - a02*a11 );
329 const amrex::Real invAx[3] = {a22*a11 - a12*a12, a12*a02 - a22*a01, a12*a01 - a11*a02};
331 const amrex::Real invAy[3] = {a02*a12 - a22*a01, a00*a22 - a02*a02, a01*a02 - a12*a00};
333 const amrex::Real invAz[3] = {a12*a01 - a02*a11, a02*a01 - a12*a00, a11*a00 - a01*a01};
337 const amrex::Real dEx = (-1._rt/detA)*(invAx[0]*Omega[0] +
341 const amrex::Real dEy = (-1._rt/detA)*(invAy[0]*Omega[0] +
345 const amrex::Real dEz = (-1._rt/detA)*(invAz[0]*Omega[0] +
351 Ex(j,k,0) = Ex(j,k,0) + 0.5_rt*
dt*dEx;
353 Ey(j,k,0) = Ey(j,k,0) + 0.5_rt*
dt*dEy;
355 Ez(j,k,0) = Ez(j,k,0) + 0.5_rt*
dt*dEz;
#define AMREX_GPU_HOST_DEVICE
AMREX_GPU_HOST_DEVICE AMREX_INLINE void warpx_hybrid_QED_push(int j, int k, int l, amrex::Array4< amrex::Real > const &Ex, amrex::Array4< amrex::Real > const &Ey, amrex::Array4< amrex::Real > const &Ez, amrex::Array4< amrex::Real > const &Bx, amrex::Array4< amrex::Real > const &By, amrex::Array4< amrex::Real const > const &Bz, amrex::Array4< amrex::Real > const &tmpEx, amrex::Array4< amrex::Real > const &tmpEy, amrex::Array4< amrex::Real > const &tmpEz, amrex::Array4< amrex::Real > const &Jx, amrex::Array4< amrex::Real > const &Jy, amrex::Array4< amrex::Real > const &Jz, const amrex::Real dx, const amrex::Real dy, const amrex::Real dz, const amrex::Real dt, const amrex::Real xi_c2)
Definition: WarpX_QED_K.H:80
AMREX_GPU_HOST_DEVICE AMREX_INLINE void calc_M(amrex::Real arr[], amrex::Real ex, amrex::Real ey, amrex::Real ez, amrex::Real bx, amrex::Real by, amrex::Real bz, amrex::Real xi_c2, amrex::Real c2)
Definition: WarpX_QED_K.H:31
static constexpr auto xi_c2
xi times c2 = xi*c*c. This should be usable for single precision instead of xi; very close to smalles...
Definition: constant.H:67
static constexpr auto c
vacuum speed of light [m/s]
Definition: constant.H:44
static constexpr auto ep0
vacuum permittivity: dielectric permittivity of vacuum [F/m]
Definition: constant.H:46
static constexpr auto mu0
vacuum permeability: magnetic permeability of vacuum = 4.0e-7 * pi [H/m]
Definition: constant.H:48
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void ignore_unused(const Ts &...)
int dz
Definition: compute_domain.py:36
float dt
Definition: stencil.py:442
tuple dx
lab frame
Definition: stencil.py:429
beta
Definition: stencil.py:434