diff options
Diffstat (limited to 'Source')
| -rw-r--r-- | Source/FieldSolver/WarpX_FDTD.H | 91 |
1 files changed, 41 insertions, 50 deletions
diff --git a/Source/FieldSolver/WarpX_FDTD.H b/Source/FieldSolver/WarpX_FDTD.H index e2a9d0abb..ee13875f6 100644 --- a/Source/FieldSolver/WarpX_FDTD.H +++ b/Source/FieldSolver/WarpX_FDTD.H @@ -10,10 +10,10 @@ void warpx_push_bx_yee(int i, int j, int k, Array4<Real> const& Bx, Array4<Real const> const& Ey, Array4<Real const> const& Ez, Real dtsdy, Real dtsdz) { -#if (AMREX_SPACEDIM == 3) +#if defined WARPX_DIM_3D Bx(i,j,k) += - dtsdy * (Ez(i,j+1,k ) - Ez(i,j,k)) + dtsdz * (Ey(i,j ,k+1) - Ey(i,j,k)); -#else +#elif (defined WARPX_DIM_2D) || (defined WARPX_DIM_RZ) // Note that the 2D Cartesian and RZ versions are the same Bx(i,j,0) += + dtsdz * (Ey(i,j+1,0) - Ey(i,j,0)); #endif @@ -24,10 +24,10 @@ void warpx_push_by_yee(int i, int j, int k, Array4<Real> const& By, Array4<Real const> const& Ex, Array4<Real const> const& Ez, Real dtsdx, Real dtsdz) { -#if (AMREX_SPACEDIM == 3) +#if defined WARPX_DIM_3D By(i,j,k) += + dtsdx * (Ez(i+1,j,k ) - Ez(i,j,k)) - dtsdz * (Ex(i ,j,k+1) - Ex(i,j,k)); -#else +#elif (defined WARPX_DIM_2D) || (defined WARPX_DIM_RZ) // Note that the 2D Cartesian and RZ versions are the same By(i,j,0) += + dtsdx * (Ez(i+1,j ,0) - Ez(i,j,0)) - dtsdz * (Ex(i ,j+1,0) - Ex(i,j,0)); @@ -39,18 +39,16 @@ void warpx_push_bz_yee(int i, int j, int k, Array4<Real> const& Bz, Array4<Real const> const& Ex, Array4<Real const> const& Ey, Real dtsdx, Real dtsdy, Real dxinv, Real rmin) { -#if (AMREX_SPACEDIM == 3) +#if defined WARPX_DIM_3D Bz(i,j,k) += - dtsdx * (Ey(i+1,j ,k) - Ey(i,j,k)) + dtsdy * (Ex(i ,j+1,k) - Ex(i,j,k)); -#else -#ifndef WARPX_RZ +#elif defined WARPX_DIM_2D Bz(i,j,0) += - dtsdx * (Ey(i+1,j,0) - Ey(i,j,0)); -#else +#elif defined WARPX_DIM_RZ Real ru = 1. + 0.5/(rmin*dxinv + i + 0.5); Real rd = 1. - 0.5/(rmin*dxinv + i + 0.5); Bz(i,j,0) += - dtsdx * (ru*Ey(i+1,j,0) - rd*Ey(i,j,0)); #endif -#endif } AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE @@ -58,11 +56,11 @@ void warpx_push_ex_yee(int i, int j, int k, Array4<Real> const& Ex, Array4<Real const> const& By, Array4<Real const> const& Bz, Array4<Real const> const& Jx, Real mu_c2_dt, Real dtsdy_c2, Real dtsdz_c2) { -#if (AMREX_SPACEDIM == 3) +#if defined WARPX_DIM_3D Ex(i,j,k) += + dtsdy_c2 * (Bz(i,j,k) - Bz(i,j-1,k )) - dtsdz_c2 * (By(i,j,k) - By(i,j ,k-1)) - mu_c2_dt * Jx(i,j,k); -#else +#elif (defined WARPX_DIM_2D) || (defined WARPX_DIM_RZ) // Note that the 2D Cartesian and RZ versions are the same Ex(i,j,0) += - dtsdz_c2 * (By(i,j,0) - By(i,j-1,0)) - mu_c2_dt * Jx(i,j,0); @@ -74,13 +72,13 @@ void warpx_push_ey_yee(int i, int j, int k, Array4<Real> const& Ey, Array4<Real const> const& Bx, Array4<Real const> const& Bz, Array4<Real const> const& Jy, Real mu_c2_dt, Real dtsdx_c2, Real dtsdz_c2, Real rmin) { -#if (AMREX_SPACEDIM == 3) +#if defined WARPX_DIM_3D Ey(i,j,k) += - dtsdx_c2 * (Bz(i,j,k) - Bz(i-1,j,k)) + dtsdz_c2 * (Bx(i,j,k) - Bx(i,j,k-1)) - mu_c2_dt * Jy(i,j,k); -#else +#elif (defined WARPX_DIM_2D) || (defined WARPX_DIM_RZ) // 2D Cartesian and RZ are the same, except that the axis is skipped with RZ -#ifdef WARPX_RZ +#ifdef WARPX_DIM_RZ if (i != 0 || rmin != 0.) #endif { @@ -96,15 +94,14 @@ void warpx_push_ez_yee(int i, int j, int k, Array4<Real> const& Ez, Array4<Real const> const& Bx, Array4<Real const> const& By, Array4<Real const> const& Jz, Real mu_c2_dt, Real dtsdx_c2, Real dtsdy_c2, Real dxinv, Real rmin) { -#if (AMREX_SPACEDIM == 3) +#if defined WARPX_DIM_3D Ez(i,j,k) += + dtsdx_c2 * (By(i,j,k) - By(i-1,j ,k)) - dtsdy_c2 * (Bx(i,j,k) - Bx(i ,j-1,k)) - mu_c2_dt * Jz(i,j,k); -#else -#ifndef WARPX_RZ +#elif defined WARPX_DIM_2D Ez(i,j,0) += + dtsdx_c2 * (By(i,j,0) - By(i-1,j,0)) - mu_c2_dt * Jz(i,j,0); -#else +#elif defined WARPX_DIM_RZ if (i != 0 || rmin != 0.) { Real ru = 1. + 0.5/(rmin*dxinv + i); Real rd = 1. - 0.5/(rmin*dxinv + i); @@ -115,16 +112,15 @@ void warpx_push_ez_yee(int i, int j, int k, Array4<Real> const& Ez, - mu_c2_dt * Jz(i,j,0); } #endif -#endif } AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void warpx_push_ex_f_yee(int j, int k, int l, Array4<Real> const& Ex, Array4<Real const> const& F, Real dtsdx_c2) { -#if (AMREX_SPACEDIM == 3) +#if defined WARPX_DIM_3D Ex(j,k,l) += + dtsdx_c2 * (F(j+1,k,l) - F(j,k,l)); -#else +#elif (defined WARPX_DIM_2D) || (defined WARPX_DIM_RZ) Ex(j,k,0) += + dtsdx_c2 * (F(j+1,k,0) - F(j ,k,0)); #endif } @@ -133,7 +129,7 @@ AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void warpx_push_ey_f_yee(int j, int k, int l, Array4<Real> const& Ey, Array4<Real const> const& F, Real dtsdy_c2) { -#if (AMREX_SPACEDIM == 3) +#if defined WARPX_DIM_3D Ey(j,k,l) += + dtsdy_c2 * (F(j,k+1,l) - F(j,k,l)); #endif } @@ -142,9 +138,9 @@ AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void warpx_push_ez_f_yee(int j, int k, int l, Array4<Real> const& Ez, Array4<Real const> const& F, Real dtsdz_c2) { -#if (AMREX_SPACEDIM == 3) +#if defined WARPX_DIM_3D Ez(j,k,l) += + dtsdz_c2 * (F(j,k,l+1) - F(j,k,l)); -#else +#elif (defined WARPX_DIM_2D) || (defined WARPX_DIM_RZ) Ez(j,k,0) += + dtsdz_c2 * (F(j,k+1,0) - F(j,k,0)); #endif } @@ -156,7 +152,7 @@ static void warpx_calculate_ckc_coefficients(Real dtsdx, Real dtsdy, Real dtsdz, { // Cole-Karkkainen-Cowan push // computes coefficients according to Cowan - PRST-AB 16, 041303 (2013) -#if (AMREX_SPACEDIM == 3) +#if defined WARPX_DIM_3D Real delta = std::max( { dtsdx,dtsdy,dtsdz } ); Real rx = (dtsdx/delta)*(dtsdx/delta); Real ry = (dtsdy/delta)*(dtsdy/delta); @@ -187,7 +183,7 @@ static void warpx_calculate_ckc_coefficients(Real dtsdx, Real dtsdy, Real dtsdz, gammax *= dtsdx; gammay *= dtsdy; gammaz *= dtsdz; -#else +#elif defined WARPX_DIM_2D Real delta = std::max(dtsdx,dtsdz); Real rx = (dtsdx/delta)*(dtsdx/delta); Real rz = (dtsdz/delta)*(dtsdz/delta); @@ -210,7 +206,7 @@ void warpx_push_bx_ckc(int j, int k, int l, Array4<Real> const& Bx, Real gammax, Real gammay, Real gammaz, Real alphax, Real alphay, Real alphaz) { -#if (AMREX_SPACEDIM == 3) +#if defined WARPX_DIM_3D Bx(j,k,l) += - alphay * (Ez(j ,k+1,l ) - Ez(j, k ,l )) - betayx * (Ez(j+1,k+1,l ) - Ez(j+1,k ,l ) + Ez(j-1,k+1,l ) - Ez(j-1,k ,l )) @@ -229,7 +225,7 @@ void warpx_push_bx_ckc(int j, int k, int l, Array4<Real> const& Bx, + Ey(j-1,k+1,l+1) - Ey(j-1,k+1,l ) + Ey(j+1,k-1,l+1) - Ey(j+1,k-1,l ) + Ey(j-1,k-1,l+1) - Ey(j-1,k-1,l )); -#else +#elif defined WARPX_DIM_2D Bx(j,k,0) += + alphaz * (Ey(j ,k+1,0) - Ey(j, k,0)) + betazx * (Ey(j+1,k+1,0) - Ey(j+1,k,0) + Ey(j-1,k+1,0) - Ey(j-1,k,0)); @@ -243,7 +239,7 @@ void warpx_push_by_ckc(int j, int k, int l, Array4<Real> const& By, Real gammax, Real gammay, Real gammaz, Real alphax, Real alphay, Real alphaz) { -#if (AMREX_SPACEDIM == 3) +#if defined WARPX_DIM_3D By(j,k,l) += + alphax * (Ez(j+1,k ,l ) - Ez(j, k, l )) + betaxy * (Ez(j+1,k+1,l ) - Ez(j ,k+1,l ) + Ez(j+1,k-1,l ) - Ez(j ,k-1,l )) @@ -262,7 +258,7 @@ void warpx_push_by_ckc(int j, int k, int l, Array4<Real> const& By, + Ex(j-1,k+1,l+1) - Ex(j-1,k+1,l ) + Ex(j+1,k-1,l+1) - Ex(j+1,k-1,l ) + Ex(j-1,k-1,l+1) - Ex(j-1,k-1,l )); -#else +#elif defined WARPX_DIM_2D By(j,k,0) += + alphax * (Ez(j+1,k ,0) - Ez(j,k ,0)) + betaxz * (Ez(j+1,k+1,0) - Ez(j,k+1,0) + Ez(j+1,k-1,0) - Ez(j,k-1,0)) @@ -279,7 +275,7 @@ void warpx_push_bz_ckc(int j, int k, int l, Array4<Real> const& Bz, Real gammax, Real gammay, Real gammaz, Real alphax, Real alphay, Real alphaz) { -#if (AMREX_SPACEDIM == 3) +#if defined WARPX_DIM_3D Bz(j,k,l) += - alphax * (Ey(j+1,k ,l ) - Ey(j ,k ,l )) - betaxy * (Ey(j+1,k+1,l ) - Ey(j ,k+1,l ) + Ey(j+1,k-1,l ) - Ey(j ,k-1,l )) @@ -298,8 +294,7 @@ void warpx_push_bz_ckc(int j, int k, int l, Array4<Real> const& Bz, + Ex(j-1,k+1,l+1) - Ex(j-1,k ,l+1) + Ex(j+1,k+1,l-1) - Ex(j+1,k ,l-1) + Ex(j-1,k+1,l-1) - Ex(j-1,k ,l-1)); - -#else +#elif defined WARPX_DIM_2D Bz(j,k,0) += - alphax * (Ey(j+1,k ,0) - Ey(j,k ,0)) - betaxz * (Ey(j+1,k+1,0) - Ey(j,k+1,0) + Ey(j+1,k-1,0) - Ey(j,k-1,0)); @@ -313,7 +308,7 @@ void warpx_push_ex_f_ckc(int j, int k, int l, Array4<Real> const& Ex, Real gammax, Real gammay, Real gammaz, Real alphax, Real alphay, Real alphaz) { -#if (AMREX_SPACEDIM == 3) +#if defined WARPX_DIM_3D Ex(j,k,l) += + alphax * (F(j+1,k ,l ) - F(j, k, l )) + betaxy * (F(j+1,k+1,l ) - F(j,k+1,l ) + F(j+1,k-1,l ) - F(j,k-1,l )) @@ -323,7 +318,7 @@ void warpx_push_ex_f_ckc(int j, int k, int l, Array4<Real> const& Ex, + F(j+1,k-1,l+1) - F(j,k-1,l+1) + F(j+1,k+1,l-1) - F(j,k+1,l-1) + F(j+1,k-1,l-1) - F(j,k-1,l-1)); -#else +#elif defined WARPX_DIM_2D Ex(j,k,0) += + alphax * (F(j+1,k ,0) - F(j,k ,0)) + betaxz * (F(j+1,k+1,0) - F(j,k+1,0) + F(j+1,k-1,0) - F(j,k-1,0)); @@ -337,7 +332,7 @@ void warpx_push_ey_f_ckc(int j, int k, int l, Array4<Real> const& Ey, Real gammax, Real gammay, Real gammaz, Real alphax, Real alphay, Real alphaz) { -#if (AMREX_SPACEDIM == 3) +#if defined WARPX_DIM_3D Ey(j,k,l) += + alphay * (F(j ,k+1,l ) - F(j ,k,l )) + betayx * (F(j+1,k+1,l ) - F(j+1,k,l ) + F(j-1,k+1,l ) - F(j-1,k,l )) @@ -357,7 +352,7 @@ void warpx_push_ez_f_ckc(int j, int k, int l, Array4<Real> const& Ez, Real gammax, Real gammay, Real gammaz, Real alphax, Real alphay, Real alphaz) { -#if (AMREX_SPACEDIM == 3) +#if defined WARPX_DIM_3D Ez(j,k,l) += + alphaz * (F(j ,k ,l+1) - F(j, k, l)) + betazx * (F(j+1,k ,l+1) - F(j+1,k ,l) + F(j-1,k ,l+1) - F(j-1,k ,l)) @@ -367,7 +362,7 @@ void warpx_push_ez_f_ckc(int j, int k, int l, Array4<Real> const& Ez, + F(j-1,k+1,l+1) - F(j-1,k+1,l) + F(j+1,k-1,l+1) - F(j+1,k-1,l) + F(j-1,k-1,l+1) - F(j-1,k-1,l)); -#else +#elif defined WARPX_DIM_2D Ez(j,k,0) += + alphaz * (F(j ,k+1,0) - F(j ,k,0)) + betazx * (F(j+1,k+1,0) - F(j+1,k,0) + F(j-1,k+1,0) - F(j-1,k,0)); @@ -378,46 +373,43 @@ AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void warpx_computedivb(int i, int j, int k, int dcomp, Array4<Real> const& divB, Array4<Real const> const& Bx, Array4<Real const> const& By, Array4<Real const> const& Bz, Real dxinv, Real dyinv, Real dzinv -#ifdef WARPX_RZ +#ifdef WARPX_DIM_RZ ,const Real rmin #endif ) { -#if (AMREX_SPACEDIM == 3) +#if defined WARPX_DIM_3D divB(i,j,k,dcomp) = (Bx(i+1,j ,k ) - Bx(i,j,k))*dxinv + (By(i ,j+1,k ) - By(i,j,k))*dyinv + (Bz(i ,j ,k+1) - Bz(i,j,k))*dzinv; -#else -#ifndef WARPX_RZ +#elif defined WARPX_DIM_2D divB(i,j,0,dcomp) = (Bx(i+1,j ,0) - Bx(i,j,0))*dxinv + (Bz(i ,j+1,0) - Bz(i,j,0))*dzinv; -#else +#elif defined WARPX_DIM_RZ Real ru = 1. + 0.5/(rmin*dxinv + i + 0.5); Real rd = 1. - 0.5/(rmin*dxinv + i + 0.5); divB(i,j,0,dcomp) = (ru*Bx(i+1,j,0) - rd*Bx(i,j,0))*dxinv + (Bz(i,j+1,0) - Bz(i,j,0))*dzinv; #endif -#endif } AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void warpx_computedive(int i, int j, int k, int dcomp, Array4<Real> const& divE, Array4<Real const> const& Ex, Array4<Real const> const& Ey, Array4<Real const> const& Ez, Real dxinv, Real dyinv, Real dzinv -#ifdef WARPX_RZ +#ifdef WARPX_DIM_RZ ,const Real rmin #endif ) { -#if (AMREX_SPACEDIM == 3) +#if defined WARPX_DIM_3D divE(i,j,k,dcomp) = (Ex(i,j,k) - Ex(i-1,j,k))*dxinv + (Ey(i,j,k) - Ey(i,j-1,k))*dyinv + (Ez(i,j,k) - Ez(i,j,k-1))*dzinv; -#else -#ifndef WARPX_RZ +#elif defined WARPX_DIM_2D divE(i,j,0,dcomp) = (Ex(i,j,0) - Ex(i-1,j,0))*dxinv + (Ez(i,j,0) - Ez(i,j-1,0))*dzinv; -#else +#elif defined WARPX_DIM_RZ if (i == 0 && rmin == 0.) { // the bulk equation diverges on axis // (due to the 1/r terms). The following expressions regularize @@ -431,7 +423,6 @@ void warpx_computedive(int i, int j, int k, int dcomp, Array4<Real> const& divE, + (Ez(i,j,0) - Ez(i,j-1,0))*dzinv; } #endif -#endif } #endif |