diff options
Diffstat (limited to 'Source')
| -rw-r--r-- | Source/FieldSolver/WarpXPushFieldsEM.cpp | 14 | ||||
| -rw-r--r-- | Source/FieldSolver/WarpX_FDTD.H | 154 |
2 files changed, 137 insertions, 31 deletions
diff --git a/Source/FieldSolver/WarpXPushFieldsEM.cpp b/Source/FieldSolver/WarpXPushFieldsEM.cpp index 1df05bc0f..dfb9f0211 100644 --- a/Source/FieldSolver/WarpXPushFieldsEM.cpp +++ b/Source/FieldSolver/WarpXPushFieldsEM.cpp @@ -175,20 +175,21 @@ WarpX::EvolveB (int lev, PatchType patch_type, amrex::Real a_dt) warpx_push_bz_nodal(j,k,l,Bzfab,Exfab,Eyfab,dtsdx,dtsdy); }); } else if (WarpX::maxwell_fdtd_solver_id == 0) { + const long nmodes = n_rz_azimuthal_modes; amrex::ParallelFor(tbx, [=] AMREX_GPU_DEVICE (int j, int k, int l) { - warpx_push_bx_yee(j,k,l,Bxfab,Eyfab,Ezfab,dtsdy,dtsdz); + warpx_push_bx_yee(j,k,l,Bxfab,Eyfab,Ezfab,dtsdx,dtsdy,dtsdz,dxinv,xmin,nmodes); }); amrex::ParallelFor(tby, [=] AMREX_GPU_DEVICE (int j, int k, int l) { - warpx_push_by_yee(j,k,l,Byfab,Exfab,Ezfab,dtsdx,dtsdz); + warpx_push_by_yee(j,k,l,Byfab,Exfab,Ezfab,dtsdx,dtsdz,nmodes); }); amrex::ParallelFor(tbz, [=] AMREX_GPU_DEVICE (int j, int k, int l) { - warpx_push_bz_yee(j,k,l,Bzfab,Exfab,Eyfab,dtsdx,dtsdy,dxinv,xmin); + warpx_push_bz_yee(j,k,l,Bzfab,Exfab,Eyfab,dtsdx,dtsdy,dxinv,xmin,nmodes); }); } else if (WarpX::maxwell_fdtd_solver_id == 1) { Real betaxy, betaxz, betayx, betayz, betazx, betazy; @@ -372,20 +373,21 @@ WarpX::EvolveE (int lev, PatchType patch_type, amrex::Real a_dt) warpx_push_ez_nodal(j,k,l,Ezfab,Bxfab,Byfab,jzfab,mu_c2_dt,dtsdx_c2,dtsdy_c2); }); } else { + const long nmodes = n_rz_azimuthal_modes; amrex::ParallelFor(tex, [=] AMREX_GPU_DEVICE (int j, int k, int l) { - warpx_push_ex_yee(j,k,l,Exfab,Byfab,Bzfab,jxfab,mu_c2_dt,dtsdy_c2,dtsdz_c2); + warpx_push_ex_yee(j,k,l,Exfab,Byfab,Bzfab,jxfab,mu_c2_dt,dtsdx_c2,dtsdy_c2,dtsdz_c2,dxinv,xmin,nmodes); }); amrex::ParallelFor(tey, [=] AMREX_GPU_DEVICE (int j, int k, int l) { - warpx_push_ey_yee(j,k,l,Eyfab,Bxfab,Bzfab,jyfab,mu_c2_dt,dtsdx_c2,dtsdz_c2,xmin); + warpx_push_ey_yee(j,k,l,Eyfab,Bxfab,Bzfab,jyfab,Exfab,mu_c2_dt,dtsdx_c2,dtsdz_c2,xmin,nmodes); }); amrex::ParallelFor(tez, [=] AMREX_GPU_DEVICE (int j, int k, int l) { - warpx_push_ez_yee(j,k,l,Ezfab,Bxfab,Byfab,jzfab,mu_c2_dt,dtsdx_c2,dtsdy_c2,dxinv,xmin); + warpx_push_ez_yee(j,k,l,Ezfab,Bxfab,Byfab,jzfab,mu_c2_dt,dtsdx_c2,dtsdy_c2,dxinv,xmin,nmodes); }); } diff --git a/Source/FieldSolver/WarpX_FDTD.H b/Source/FieldSolver/WarpX_FDTD.H index e2d40f4c2..4a5139d50 100644 --- a/Source/FieldSolver/WarpX_FDTD.H +++ b/Source/FieldSolver/WarpX_FDTD.H @@ -8,36 +8,75 @@ using namespace amrex; AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE 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) + Real dtsdx, Real dtsdy, Real dtsdz, Real dxinv, Real rmin, const long nmodes) { #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)); -#elif (defined WARPX_DIM_2D) || (defined WARPX_DIM_RZ) - // Note that the 2D Cartesian and RZ versions are the same +#elif (defined WARPX_DIM_2D) Bx(i,j,0) += + dtsdz * (Ey(i,j+1,0) - Ey(i,j,0)); +#if (defined WARPX_DIM_RZ) + if (i != 0 || rmin != 0.) { + Bx(i,j,0,0) += + dtsdz * (Ey(i,j+1,0,0) - Ey(i,j,0,0)); + } else { + Bx(i,j,0,0) = 0.; + } + for (int imode=1 ; imode < nmodes ; imode++) { + if (i == 0 && rmin == 0) { + if (imode == 1) { + // For the mode m = 1, the bulk equation diverges on axis + // (due to the 1/r terms). The following expressions regularize + // these divergences by assuming, on axis : + // Ez/r = 0/r + dEz/dr + Bx(i,j,2*imode) = Bx(i,j,0,2*imode) - imode*dtsdx*Ez(i+1,j,0,2*imode+1) & + + dtsdz*(Ey(i,j+1,0,2*imode) - Ey(i,j,0,2*imode)); + Bx(i,j,2*imode+1) = Bx(i,j,0,2*imode+1) + imode*dtsdx*Ez(i+1,j,0,2*imode) & + + dtsdz*(Ey(i,j+1,0,2*imode+1) - Ey(i,j,0,2*imode+1)); + } else { + Bx(i,j,0,2*imode) = 0.; + Bx(i,j,0,2*imode+1) = 0.; + } + } else { + // Br(i,j,m) = Br(i,j,m) + I*m*dt*Ez(i,j,m)/r + dtsdz*(Et(i,j+1,m) - Et(i,j,m)) + const Real r = rmin*dxinv + i; + Bx(i,j,2*imode) = Bx(i,j,0,2*imode) - imode*dtsdx*Ez(i,j,0,2*imode+1)/r & + + dtsdz*(Ey(i,j+1,0,2*imode) - Ey(i,j,0,2*imode)); + Bx(i,j,2*imode+1) = Bx(i,j,0,2*imode+1) + imode*dtsdx*Ez(i,j,0,2*imode)/r & + + dtsdz*(Ey(i,j+1,0,2*imode+1) - Ey(i,j,0,2*imode+1)); + } + } +#endif #endif } AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE 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) + Real dtsdx, Real dtsdz, const long nmodes) { #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)); #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)); + // Note that the 2D Cartesian and RZ mode 0 are the same + By(i,j,0,0) += + dtsdx * (Ez(i+1,j ,0,0) - Ez(i,j,0,0)) + - dtsdz * (Ex(i ,j+1,0,0) - Ex(i,j,0,0)); +#if (defined WARPX_DIM_RZ) + for (int imode=1 ; imode < nmodes ; imode++) { + // Bt(i,j,m) = Bt(i,j,m) + dtsdr*(Ez(i+1,j,m) - Ez(i,j,m)) - dtsdz*(Er(i,j+1,m) - Er(i,j,m)) + By(i,j,0,2*imode) += + dtsdx*(Ez(i+1,j ,0,2*imode) - Ez(i,j,0,2*imode)) + - dtsdz*(Ex(i ,j+1,0,2*imode) - Ex(i,j,0,2*imode)); + By(i,j,0,2*imode+1) += + dtsdx*(Ez(i+1,j ,0,2*imode+1) - Ez(i,j,0,2*imode+1)) + - dtsdz*(Ex(i ,j+1,0,2*imode+1) - Ex(i,j,0,2*imode+1)); + } +#endif #endif } AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE 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) + Real dtsdx, Real dtsdy, Real dxinv, Real rmin, const long nmodes) { #if defined WARPX_DIM_3D Bz(i,j,k) += - dtsdx * (Ey(i+1,j ,k) - Ey(i,j,k)) @@ -45,46 +84,94 @@ void warpx_push_bz_yee(int i, int j, int k, Array4<Real> const& Bz, #elif defined WARPX_DIM_2D Bz(i,j,0) += - dtsdx * (Ey(i+1,j,0) - Ey(i,j,0)); #elif defined WARPX_DIM_RZ + const Real r = rmin*dxinv + i + 0.5; const Real ru = 1. + 0.5/(rmin*dxinv + i + 0.5); const 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)); + Bz(i,j,0,0) += - dtsdx*(ru*Ey(i+1,j,0,0) - rd*Ey(i,j,0,0)); + for (int imode=1 ; imode < nmodes ; imode++) { + // Bz(i,j,m) = Bz(i,j,m) - dtsdr*(ru*Et(i+1,j,m) - rd*Et(i,j,m)) - I*m*dt*Er(i,j,m)/r + Bz(i,j,0,2*imode) += - dtsdx*(ru*Ey(i+1,j,0,2*imode) - rd*Ey(i,j,0,2*imode)) + imode*dtsdx*Ex(i,j,2*imode+1)/r; + Bz(i,j,0,2*imode+1) += - dtsdx*(ru*Ey(i+1,j,0,2*imode+1) - rd*Ey(i,j,0,2*imode+1)) - imode*dtsdx*Ex(i,j,2*imode)/r; + } #endif } AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE 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) + Real mu_c2_dt, Real dtsdx_c2, Real dtsdy_c2, Real dtsdz_c2, Real dxinv, Real rmin, const long nmodes) { #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); #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); + // Note that the 2D Cartesian and RZ mode 0 are the same + Ex(i,j,0,0) += - dtsdz_c2 * (By(i,j,0,0) - By(i,j-1,0,0)) + - mu_c2_dt * Jx(i,j,0,0); +#if (defined WARPX_DIM_RZ) + const Real r = rmin*dxinv+ i + 0.5; + for (int imode=1 ; imode < nmodes ; imode++) { + // Er(i,j,m) = Er(i,j,m) - I*m*dt*Bz(i,j,m)/r - dtsdz*(Bt(i,j,m) - Bt(i,j-1,m)) - mudt*Jr(i,j,m) + Ex(i,j,0,2*imode) += - dtsdz_c2*(By(i,j,0,2*imode) - By(i,j-1,0,2*imode)) + imode*dtsdx_c2*Bz(i,j,0,2*imode+1)/r + - mu_c2_dt*Jx(i,j,0,2*imode); + Ex(i,j,0,2*imode+1) += - dtsdz_c2*(By(i,j,0,2*imode+1) - By(i,j-1,0,2*imode+1)) - imode*dtsdx_c2*Bz(i,j,0,2*imode)/r + - mu_c2_dt*Jx(i,j,0,2*imode+1); + } +#endif #endif } AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE 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) + Array4<Real> const& Ex, + Real mu_c2_dt, Real dtsdx_c2, Real dtsdz_c2, Real rmin, const long nmodes) { #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); -#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_DIM_RZ - if (i != 0 || rmin != 0.) -#endif - { +#elif (defined WARPX_DIM_2D) Ey(i,j,0) += - dtsdx_c2 * (Bz(i,j,0) - Bz(i-1,j,0)) + dtsdz_c2 * (Bx(i,j,0) - Bx(i,j-1,0)) - mu_c2_dt * Jy(i,j,0); +#elif (defined WARPX_DIM_RZ) + if (i != 0 || rmin != 0.) { + Ey(i,j,0,0) += - dtsdx_c2 * (Bz(i,j,0,0) - Bz(i-1,j,0,0)) + + dtsdz_c2 * (Bx(i,j,0,0) - Bx(i,j-1,0,0)) + - mu_c2_dt * Jy(i,j,0,0); + } else { + Ey(i,j,0,0) = 0.; + } + for (int imode=1 ; imode < nmodes ; imode++) { + if (i == 0 && rmin == 0) { + if (imode == 1) { + // The bulk equation could in principle be used here since it does not diverge + // on axis. However, it typically gives poor results e.g. for the propagation + // of a laser pulse (the field is spuriously reduced on axis). For this reason + // a modified on-axis condition is used here: we use the fact that + // Etheta(r=0,m=1) should equal -iEr(r=0,m=1), for the fields Er and Et to be + // independent of theta at r=0. Now with linear interpolation: + // Er(r=0,m=1) = 0.5*[Er(r=dr/2,m=1) + Er(r=-dr/2,m=1)] + // And using the rule applying for the guards cells + // Er(r=-dr/2,m=1) = Er(r=dr/2,m=1). Thus: Et(i,j,m) = -i*Er(i,j,m) + Ey(i,j,0,2*imode) = Ex(i,j,0,2*imode+1); + Ey(i,j,0,2*imode+1) = -Ex(i,0,j,2*imode); + } else { + // Etheta should remain 0 on axis, for modes different than m=1 + Ey(i,j,0,2*imode) = 0.; + Ey(i,j,0,2*imode+1) = 0.; + } + } else { + // Et(i,j,m) = Et(i,j,m) - dtsdr*(Bz(i,j,m) - Bz(i-1,j,m)) + dtsdz*(Br(i,j,m) - Br(i,j-1,m)) - mudt*Jt(i,j,m) + Ey(i,j,0,2*imode) += - dtsdx_c2 * (Bz(i,j,0,2*imode) - Bz(i-1,j,0,2*imode)) + + dtsdz_c2 * (Bx(i,j,0,2*imode) - Bx(i,j-1,0,2*imode)) + - mu_c2_dt * Jy(i,j,0,2*imode); + Ey(i,j,0,2*imode+1) += - dtsdx_c2 * (Bz(i,j,0,2*imode+1) - Bz(i-1,j,0,2*imode+1)) + + dtsdz_c2 * (Bx(i,j,0,2*imode+1) - Bx(i,j-1,0,2*imode+1)) + - mu_c2_dt * Jy(i,j,0,2*imode+1); + } } #endif } @@ -92,7 +179,7 @@ void warpx_push_ey_yee(int i, int j, int k, Array4<Real> const& Ey, AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE 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) + Real mu_c2_dt, Real dtsdx_c2, Real dtsdy_c2, Real dxinv, Real rmin, const long nmodes) { #if defined WARPX_DIM_3D Ez(i,j,k) += + dtsdx_c2 * (By(i,j,k) - By(i-1,j ,k)) @@ -105,11 +192,28 @@ void warpx_push_ez_yee(int i, int j, int k, Array4<Real> const& Ez, if (i != 0 || rmin != 0.) { const Real ru = 1. + 0.5/(rmin*dxinv + i); const Real rd = 1. - 0.5/(rmin*dxinv + i); - Ez(i,j,0) += + dtsdx_c2 * (ru*By(i,j,0) - rd*By(i-1,j,0)) - - mu_c2_dt * Jz(i,j,0); + Ez(i,j,0,0) += + dtsdx_c2 * (ru*By(i,j,0,0) - rd*By(i-1,j,0,0)) + - mu_c2_dt * Jz(i,j,0,0); } else { - Ez(i,j,0) += + 4.*dtsdx_c2 * By(i,j,0) - - mu_c2_dt * Jz(i,j,0); + Ez(i,j,0,0) += + 4.*dtsdx_c2 * By(i,j,0,0) + - mu_c2_dt * Jz(i,j,0,0); + } + for (int imode=1 ; imode < nmodes ; imode++) { + if (i == 0 && rmin == 0) { + Ez(i,j,0,2*imode) = 0.; + Ez(i,j,0,2*imode+1) = 0.; + } else { + const Real r = rmin*dxinv + i + 0.5; + const Real ru = 1. + 0.5/(rmin*dxinv + i); + const Real rd = 1. - 0.5/(rmin*dxinv + i); + // Ez(i,j,m) = Ez(i,j,m) + dtsdr*(ru*Bt(i,j,m) - rd*Bt(i-1,j,m)) + I*m*dt*Br(i,j,m)/r - mudt*Jz(i,j,m) + Ez(i,j,0,2*imode) += + dtsdx_c2 * (ru*By(i,j,0,2*imode) - rd*By(i-1,j,0,2*imode)) + - imode*dtsdx_c2*Bx(i,j,0,2*imode+1)/r + - mu_c2_dt * Jz(i,j,0,2*imode); + Ez(i,j,0,2*imode+1) += + dtsdx_c2 * (ru*By(i,j,0,2*imode+1) - rd*By(i-1,j,0,2*imode+1)) + + imode*dtsdx_c2*Bx(i,j,0,2*imode)/r + - mu_c2_dt * Jz(i,j,0,2*imode+1); + } } #endif } |