WarpX
WarpX_QED_K.H
Go to the documentation of this file.
1 /* Copyright 2019-2020 Glenn Richardson, Maxence Thevenet
2  *
3  * This file is part of WarpX.
4  *
5  * License: BSD-3-Clause-LBNL
6  */
7 
8 #ifndef WARPX_QED_K_h
9 #define WARPX_QED_K_h
10 
11 #include "Utils/WarpXConst.H"
12 
13 #include <AMReX.H>
14 #include <AMReX_FArrayBox.H>
15 
16 #include <cmath>
17 
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)
34 {
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 );
42 }
43 
44 
81  int j, int k, int l, amrex::Array4<amrex::Real> const& Ex,
85  amrex::Array4<amrex::Real> const& tmpEy, amrex::Array4<amrex::Real> const& tmpEz,
87  amrex::Array4<amrex::Real> const& Jz, const amrex::Real dx, const amrex::Real dy,
88  const amrex::Real dz, const amrex::Real dt, const amrex::Real xi_c2)
89 {
90 
91 using namespace amrex;
92 
93 // Defining constants to be used in the calculations
94 
95 constexpr amrex::Real c2 = PhysConst::c * PhysConst::c;
96 constexpr amrex::Real c2i = 1._rt/c2;
97 const amrex::Real dxi = 1._rt/dx;
98 const amrex::Real dzi = 1._rt/dz;
99 
100 #if defined(WARPX_DIM_3D)
101 const amrex::Real dyi = 1._rt/dy;
102 
103  // Picking out points for stencil to be used in curl function of M
104 
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};
111 
112  // Calculating the M-field at the chosen stencil points
113 
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);
126 
127  // Calculating necessary curls
128 
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 ),
133  };
134 
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 ),
139  };
140 
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 ),
145  };
146 
147  // Defining comapct values for QED corrections
148 
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);
155  const amrex::Real mu0jx = PhysConst::mu0*Jx(j,k,l);
156  const amrex::Real mu0jy = PhysConst::mu0*Jy(j,k,l);
157  const amrex::Real mu0jz = PhysConst::mu0*Jz(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;
167 
168  const amrex::Real beta = 4._rt*xi_c2*( c2i*ee - bb ) + PhysConst::ep0;
169 
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]
174  };
175 
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) )
183  };
184 
185  // Calculating matrix values for the QED correction algorithm
186 
187  const amrex::Real a00 = beta + xi_c2*( 8._rt*c2i*ex*ex + 14._rt*bx*bx );
188 
189  const amrex::Real a11 = beta + xi_c2*( 8._rt*c2i*ey*ey + 14._rt*by*by );
190 
191  const amrex::Real a22 = beta + xi_c2*( 8._rt*c2i*ez*ez + 14._rt*bz*bz );
192 
193  const amrex::Real a01 = xi_c2*( 2._rt*c2i*ex*ey + 14*bx*by );
194 
195  const amrex::Real a02 = xi_c2*( 2._rt*c2i*ex*ez + 14._rt*bx*bz );
196 
197  const amrex::Real a12 = xi_c2*( 2._rt*c2i*ez*ey + 14._rt*bz*by );
198 
199  const amrex::Real detA = a00*( a11*a22 - a12*a12 ) - a01*( a01*a22 - a02*a12 )+a02*( a01*a12 - a02*a11 );
200 
201  // Calculating the rows of the inverse matrix using the general 3x3 inverse form
202 
203  const amrex::Real invAx[3] = {a22*a11 - a12*a12, a12*a02 - a22*a01, a12*a01 - a11*a02};
204 
205  const amrex::Real invAy[3] = {a02*a12 - a22*a01, a00*a22 - a02*a02, a01*a02 - a12*a00};
206 
207  const amrex::Real invAz[3] = {a12*a01 - a02*a11, a02*a01 - a12*a00, a11*a00 - a01*a01};
208 
209  // Calculating the final QED corrections by mutliplying the Omega vector with the inverse matrix
210 
211  const amrex::Real dEx = (-1._rt/detA)*(invAx[0]*Omega[0] +
212  invAx[1]*Omega[1] +
213  invAx[2]*Omega[2]);
214 
215  const amrex::Real dEy = (-1._rt/detA)*(invAy[0]*Omega[0] +
216  invAy[1]*Omega[1] +
217  invAy[2]*Omega[2]);
218 
219  const amrex::Real dEz = (-1._rt/detA)*(invAz[0]*Omega[0] +
220  invAz[1]*Omega[1] +
221  invAz[2]*Omega[2]);
222 
223  // Adding the QED corrections to the original fields
224 
225  Ex(j,k,l) = Ex(j,k,l) + 0.5_rt*dt*dEx;
226 
227  Ey(j,k,l) = Ey(j,k,l) + 0.5_rt*dt*dEy;
228 
229  Ez(j,k,l) = Ez(j,k,l) + 0.5_rt*dt*dEz;
230 
231 
232 // 2D case - follows naturally from 3D case
233 #else
234 
235  // Picking out points for stencil to be used in curl function of M
236 
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};
241 
242  // Calculating the M-field at the chosen stencil points
243 
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);
252 
253  // Calculating necessary curls
254 
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,
259  };
260 
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,
265  };
266 
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,
271  };
272 
273  // Defining comapct values for QED corrections
274 
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);
281  const amrex::Real mu0jx = PhysConst::mu0*Jx(j,k,0);
282  const amrex::Real mu0jy = PhysConst::mu0*Jy(j,k,0);
283  const amrex::Real mu0jz = PhysConst::mu0*Jz(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;
293 
294  const amrex::Real beta = 4._rt*xi_c2*( c2i*ee - bb ) + PhysConst::ep0;
295 
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]
300  };
301 
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) )
309  };
310 
311  // Calculating matrix values for the QED correction algorithm
312 
313  const amrex::Real a00 = beta + xi_c2*( 8._rt*c2i*ex*ex + 14._rt*bx*bx );
314 
315  const amrex::Real a11 = beta + xi_c2*( 8._rt*c2i*ey*ey + 14._rt*by*by );
316 
317  const amrex::Real a22 = beta + xi_c2*( 8._rt*c2i*ez*ez + 14._rt*bz*bz );
318 
319  const amrex::Real a01 = xi_c2*( 2._rt*c2i*ex*ey + 14._rt*bx*by );
320 
321  const amrex::Real a02 = xi_c2*( 2._rt*c2i*ex*ez + 14._rt*bx*bz );
322 
323  const amrex::Real a12 = xi_c2*( 2._rt*c2i*ez*ey + 14._rt*bz*by );
324 
325  const amrex::Real detA = a00*( a11*a22 - a12*a12 ) - a01*( a01*a22 - a02*a12 ) + a02*( a01*a12 - a02*a11 );
326 
327  // Calculating inverse matrix values using the general 3x3 form
328 
329  const amrex::Real invAx[3] = {a22*a11 - a12*a12, a12*a02 - a22*a01, a12*a01 - a11*a02};
330 
331  const amrex::Real invAy[3] = {a02*a12 - a22*a01, a00*a22 - a02*a02, a01*a02 - a12*a00};
332 
333  const amrex::Real invAz[3] = {a12*a01 - a02*a11, a02*a01 - a12*a00, a11*a00 - a01*a01};
334 
335  // Calculating the final QED corrections by mutliplying the Omega vector with the inverse matrix
336 
337  const amrex::Real dEx = (-1._rt/detA)*(invAx[0]*Omega[0] +
338  invAx[1]*Omega[1] +
339  invAx[2]*Omega[2]);
340 
341  const amrex::Real dEy = (-1._rt/detA)*(invAy[0]*Omega[0] +
342  invAy[1]*Omega[1] +
343  invAy[2]*Omega[2]);
344 
345  const amrex::Real dEz = (-1._rt/detA)*(invAz[0]*Omega[0] +
346  invAz[1]*Omega[1] +
347  invAz[2]*Omega[2]);
348 
349  // Adding the QED corrections to the original fields
350 
351  Ex(j,k,0) = Ex(j,k,0) + 0.5_rt*dt*dEx;
352 
353  Ey(j,k,0) = Ey(j,k,0) + 0.5_rt*dt*dEy;
354 
355  Ez(j,k,0) = Ez(j,k,0) + 0.5_rt*dt*dEz;
356 
357  amrex::ignore_unused(l, dy);
358 #endif
359 
360 }
361 
362 #endif //WARPX_QED_K_h
#define AMREX_INLINE
#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