diff options
Diffstat (limited to 'Source/FieldSolver/WarpXPushFieldsEM.cpp')
-rw-r--r-- | Source/FieldSolver/WarpXPushFieldsEM.cpp | 91 |
1 files changed, 77 insertions, 14 deletions
diff --git a/Source/FieldSolver/WarpXPushFieldsEM.cpp b/Source/FieldSolver/WarpXPushFieldsEM.cpp index 11d598b98..b3d8e9d76 100644 --- a/Source/FieldSolver/WarpXPushFieldsEM.cpp +++ b/Source/FieldSolver/WarpXPushFieldsEM.cpp @@ -175,18 +175,19 @@ 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, tby, tbz, [=] 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_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_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; @@ -366,18 +367,19 @@ 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, tey, tez, [=] 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_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_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); }); } @@ -647,6 +649,8 @@ WarpX::ApplyInverseVolumeScalingToCurrentDensity (MultiFab* Jx, MultiFab* Jy, Mu const Real rmin = xyzmin[0]; const int irmin = lo.x; + const long nmodes = n_rz_azimuthal_modes; + // Rescale current in r-z mode since the inverse volume factor was not // included in the current deposition. amrex::ParallelFor(tbr, tbt, tbz, @@ -656,13 +660,29 @@ WarpX::ApplyInverseVolumeScalingToCurrentDensity (MultiFab* Jx, MultiFab* Jy, Mu // to the cells above the axis. // Note that Jr(i==0) is at 1/2 dr. if (rmin == 0. && 0 <= i && i < ngJ) { - Jr_arr(i,j,0) -= Jr_arr(-1-i,j,0); + Jr_arr(i,j,0,0) -= Jr_arr(-1-i,j,0,0); } // Apply the inverse volume scaling // Since Jr is not cell centered in r, no need for distinction // between on axis and off-axis factors const amrex::Real r = std::abs(rmin + (i - irmin + 0.5)*dr); - Jr_arr(i,j,0) /= (2.*MathConst::pi*r); + Jr_arr(i,j,0,0) /= (2.*MathConst::pi*r); + + for (int imode=1 ; imode < nmodes ; imode++) { + const Real ifact = ( (imode%2) == 0 ? +1. : -1.); + // Wrap the current density deposited in the guard cells around + // to the cells above the axis. + // Note that Jr(i==0) is at 1/2 dr. + if (rmin == 0. && 0 <= i && i < ngJ) { + Jr_arr(i,j,0,2*imode-1) -= ifact*Jr_arr(-1-i,j,0,2*imode-1); + Jr_arr(i,j,0,2*imode) -= ifact*Jr_arr(-1-i,j,0,2*imode); + } + // Apply the inverse volume scaling + // Since Jr is not cell centered in r, no need for distinction + // between on axis and off-axis factors + Jr_arr(i,j,0,2*imode-1) /= (2.*MathConst::pi*r); + Jr_arr(i,j,0,2*imode) /= (2.*MathConst::pi*r); + } }, [=] AMREX_GPU_DEVICE (int i, int j, int k) { @@ -670,16 +690,37 @@ WarpX::ApplyInverseVolumeScalingToCurrentDensity (MultiFab* Jx, MultiFab* Jy, Mu // to the cells above the axis. // Jt is located on the boundary if (rmin == 0. && 0 < i && i <= ngJ) { - Jt_arr(i,j,0) += Jt_arr(-i,j,0); + Jt_arr(i,j,0,0) += Jt_arr(-i,j,0,0); } // Apply the inverse volume scaling // Jt is forced to zero on axis. const amrex::Real r = std::abs(rmin + (i - irmin)*dr); if (r == 0.) { - Jt_arr(i,j,0) = 0.; + Jt_arr(i,j,0,0) = 0.; } else { - Jt_arr(i,j,0) /= (2.*MathConst::pi*r); + Jt_arr(i,j,0,0) /= (2.*MathConst::pi*r); + } + + for (int imode=1 ; imode < nmodes ; imode++) { + const Real ifact = ( (imode%2) == 0 ? +1. : -1.); + // Wrap the current density deposited in the guard cells around + // to the cells above the axis. + // Jt is located on the boundary + if (rmin == 0. && 0 < i && i <= ngJ) { + Jt_arr(i,j,0,2*imode-1) += ifact*Jt_arr(-i,j,0,2*imode-1); + Jt_arr(i,j,0,2*imode) += ifact*Jt_arr(-i,j,0,2*imode); + } + + // Apply the inverse volume scaling + // Jt is forced to zero on axis. + if (r == 0.) { + Jt_arr(i,j,0,2*imode-1) = 0.; + Jt_arr(i,j,0,2*imode) = 0.; + } else { + Jt_arr(i,j,0,2*imode-1) /= (2.*MathConst::pi*r); + Jt_arr(i,j,0,2*imode) /= (2.*MathConst::pi*r); + } } }, [=] AMREX_GPU_DEVICE (int i, int j, int k) @@ -688,17 +729,39 @@ WarpX::ApplyInverseVolumeScalingToCurrentDensity (MultiFab* Jx, MultiFab* Jy, Mu // to the cells above the axis. // Jz is located on the boundary if (rmin == 0. && 0 < i && i <= ngJ) { - Jz_arr(i,j,0) += Jz_arr(-i,j,0); + Jz_arr(i,j,0,0) += Jz_arr(-i,j,0,0); } // Apply the inverse volume scaling const amrex::Real r = std::abs(rmin + (i - irmin)*dr); if (r == 0.) { // Verboncoeur JCP 164, 421-427 (2001) : corrected volume on axis - Jz_arr(i,j,0) /= (MathConst::pi*dr/3.); + Jz_arr(i,j,0,0) /= (MathConst::pi*dr/3.); } else { - Jz_arr(i,j,0) /= (2.*MathConst::pi*r); + Jz_arr(i,j,0,0) /= (2.*MathConst::pi*r); } + + for (int imode=1 ; imode < nmodes ; imode++) { + const Real ifact = ( (imode%2) == 0 ? +1. : -1.); + // Wrap the current density deposited in the guard cells around + // to the cells above the axis. + // Jz is located on the boundary + if (rmin == 0. && 0 < i && i <= ngJ) { + Jz_arr(i,j,0,2*imode-1) += ifact*Jz_arr(-i,j,0,2*imode-1); + Jz_arr(i,j,0,2*imode) += ifact*Jz_arr(-i,j,0,2*imode); + } + + // Apply the inverse volume scaling + if (r == 0.) { + // Verboncoeur JCP 164, 421-427 (2001) : corrected volume on axis + Jz_arr(i,j,0,2*imode-1) /= (MathConst::pi*dr/3.); + Jz_arr(i,j,0,2*imode) /= (MathConst::pi*dr/3.); + } else { + Jz_arr(i,j,0,2*imode-1) /= (2.*MathConst::pi*r); + Jz_arr(i,j,0,2*imode) /= (2.*MathConst::pi*r); + } + } + }); } } |