diff options
Diffstat (limited to 'Source')
| -rw-r--r-- | Source/BoundaryConditions/PMLComponent.H | 3 | ||||
| -rw-r--r-- | Source/BoundaryConditions/WarpX_PML_kernels.H | 484 | ||||
| -rw-r--r-- | Source/FieldSolver/FiniteDifferenceSolver/EvolveBPML.cpp | 131 | ||||
| -rw-r--r-- | Source/FieldSolver/FiniteDifferenceSolver/EvolveEPML.cpp | 206 | ||||
| -rw-r--r-- | Source/FieldSolver/FiniteDifferenceSolver/EvolveFPML.cpp | 113 | ||||
| -rw-r--r-- | Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceSolver.H | 35 | ||||
| -rw-r--r-- | Source/FieldSolver/FiniteDifferenceSolver/Make.package | 3 | ||||
| -rw-r--r-- | Source/FieldSolver/WarpXPushFieldsEM.cpp | 268 | ||||
| -rw-r--r-- | Source/FieldSolver/WarpX_FDTD.H | 58 | ||||
| -rw-r--r-- | Source/WarpX.cpp | 11 |
10 files changed, 525 insertions, 787 deletions
diff --git a/Source/BoundaryConditions/PMLComponent.H b/Source/BoundaryConditions/PMLComponent.H index 4e43cf34e..6fa4f8af7 100644 --- a/Source/BoundaryConditions/PMLComponent.H +++ b/Source/BoundaryConditions/PMLComponent.H @@ -15,7 +15,8 @@ struct PMLComp { enum { xy=0, xz=1, xx=2, yz=0, yx=1, yy=2, - zx=0, zy=1, zz=2 }; + zx=0, zy=1, zz=2, + x=0, y=1, z=2 }; // Used for the PML components of F }; #endif diff --git a/Source/BoundaryConditions/WarpX_PML_kernels.H b/Source/BoundaryConditions/WarpX_PML_kernels.H index 8a573c4b9..75db5ba1d 100644 --- a/Source/BoundaryConditions/WarpX_PML_kernels.H +++ b/Source/BoundaryConditions/WarpX_PML_kernels.H @@ -12,490 +12,6 @@ using namespace amrex; - -// PML BX YEE -AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE -void warpx_push_pml_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) - Bx(i,j,k,0) += - dtsdy * (Ez(i,j+1,k,0) + Ez(i,j+1,k,1) + Ez(i,j+1,k,2) - - Ez(i,j,k,0) - Ez(i,j,k,1) - Ez(i,j,k,2)) ; - Bx(i,j,k,1) += dtsdz * (Ey(i,j,k+1,0) + Ey(i,j,k+1,1) + Ey(i,j,k+1,2) - - Ey(i,j,k,0) - Ey(i,j,k,1) - Ey(i,j,k,2)) ; -#else - Bx(i,j,k,1) += dtsdz * (Ey(i,j+1,k,0) + Ey(i,j+1,k,1) + Ey(i,j+1,k,2) - - Ey(i,j,k,0) - Ey(i,j,k,1) - Ey(i,j,k,2)) ; -#endif - -} - -//PML BY YEE -AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE -void warpx_push_pml_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) - By(i,j,k,0) += -dtsdz * (Ex(i,j,k+1,0) + Ex(i,j,k+1,1) + Ex(i,j,k+1,2) - - Ex(i,j,k,0) - Ex(i,j,k,1) - Ex(i,j,k,2)); -#else - By(i,j,k,0) += -dtsdz * (Ex(i,j+1,k,0) + Ex(i,j+1,k,1) + Ex(i,j+1,k,2) - - Ex(i,j,k,0) - Ex(i,j,k,1) - Ex(i,j,k,2)); -#endif - By(i,j,k,1) += dtsdx * (Ez(i+1,j,k,0) + Ez(i+1,j,k,1) + Ez(i+1,j,k,2) - - Ez(i,j,k,0) - Ez(i,j,k,1) - Ez(i,j,k,2)); -} - - -//PML BZ YEE -AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE -void warpx_push_pml_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) -{ - Bz(i,j,k,0) += -dtsdx * (Ey(i+1,j,k,0) + Ey(i+1,j,k,1) + Ey(i+1,j,k,2) - -Ey(i,j,k,0) - Ey(i,j,k,1) - Ey(i,j,k,2) ); -#if (AMREX_SPACEDIM == 3) - Bz(i,j,k,1) += dtsdy * (Ex(i,j+1,k,0) + Ex(i,j+1,k,1) + Ex(i,j+1,k,2) - -Ex(i,j,k,0) - Ex(i,j,k,1) - Ex(i,j,k,2) ); -#endif -} - - -//PML BX CKC -AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE -void warpx_push_pml_bx_ckc(int i, int j, int k, Array4<Real> const&Bx, - Array4<Real const> const& Ey, - Array4<Real const> const& Ez, - Real betaxy, Real betaxz, Real betayx, Real betayz, - Real betazx, Real betazy, Real gammax, Real gammay, - Real gammaz, Real alphax, Real alphay, Real alphaz) -{ -#if (AMREX_SPACEDIM == 3) - Bx(i,j,k,0) += -alphay*(Ez(i ,j+1,k ,0) + Ez(i ,j+1,k ,1) + Ez(i ,j+1,k ,2) - - Ez(i ,j ,k ,0) - Ez(i ,j ,k ,1) - Ez(i ,j ,k ,2)) - -betayx*(Ez(i+1,j+1,k ,0) + Ez(i+1,j+1,k ,1) + Ez(i+1,j+1,k ,2) - - Ez(i+1,j ,k ,0) - Ez(i+1,j ,k ,1) - Ez(i+1,j ,k ,2) - + Ez(i-1,j+1,k ,0) + Ez(i-1,j+1,k ,1) + Ez(i-1,j+1,k ,2) - - Ez(i-1,j ,k ,0) - Ez(i-1,j ,k ,1) - Ez(i-1,j ,k ,2)) - -betayz*(Ez(i ,j+1,k+1,0) + Ez(i ,j+1,k+1,1) + Ez(i ,j+1,k+1,2) - - Ez(i ,j ,k+1,0) - Ez(i ,j ,k+1,1) - Ez(i ,j ,k+1,2) - + Ez(i ,j+1,k-1,0) + Ez(i ,j+1,k-1,1) + Ez(i ,j+1,k-1,2) - - Ez(i ,j ,k-1,0) - Ez(i ,j ,k-1,1) - Ez(i ,j ,k-1,2)) - -gammay*(Ez(i+1,j+1,k+1,0) + Ez(i+1,j+1,k+1,1) + Ez(i+1,j+1,k+1,2) - - Ez(i+1,j ,k+1,0) - Ez(i+1,j ,k+1,1) - Ez(i+1,j ,k+1,2) - + Ez(i-1,j+1,k+1,0) + Ez(i-1,j+1,k+1,1) + Ez(i-1,j+1,k+1,2) - - Ez(i-1,j ,k+1,0) - Ez(i-1,j ,k+1,1) - Ez(i-1,j ,k+1,2) - + Ez(i+1,j+1,k-1,0) + Ez(i+1,j+1,k-1,1) + Ez(i+1,j+1,k-1,2) - - Ez(i+1,j ,k-1,0) - Ez(i+1,j ,k-1,1) - Ez(i+1,j ,k-1,2) - + Ez(i-1,j+1,k-1,0) + Ez(i-1,j+1,k-1,1) + Ez(i-1,j+1,k-1,2) - - Ez(i-1,j ,k-1,0) - Ez(i-1,j ,k-1,1) - Ez(i-1,j ,k-1,2)); - - Bx(i,j,k,1) += alphaz*(Ey(i ,j ,k+1,0) + Ey(i ,j ,k+1,1) + Ey(i ,j ,k+1,2) - - Ey(i ,j ,k ,0) - Ey(i ,j ,k ,1) - Ey(i ,j ,k ,2)) - +betazx*(Ey(i+1,j ,k+1,0) + Ey(i+1,j ,k+1,1) + Ey(i+1,j ,k+1,2) - - Ey(i+1,j ,k ,0) - Ey(i+1,j ,k ,1) - Ey(i+1,j ,k ,2) - + Ey(i-1,j ,k+1,0) + Ey(i-1,j ,k+1,1) + Ey(i-1,j ,k+1,2) - - Ey(i-1,j ,k ,0) - Ey(i-1,j ,k ,1) - Ey(i-1,j ,k ,2)) - +betazy*(Ey(i ,j+1,k+1,0) + Ey(i ,j+1,k+1,1) + Ey(i ,j+1,k+1,2) - - Ey(i ,j+1,k ,0) - Ey(i ,j+1,k ,1) - Ey(i ,j+1,k ,2) - + Ey(i, j-1,k+1,0) + Ey(i ,j-1,k+1,1) + Ey(i ,j-1,k+1,2) - - Ey(i ,j-1,k ,0) - Ey(i ,j-1,k ,1) - Ey(i ,j-1,k ,2)) - +gammaz*(Ey(i+1,j+1,k+1,0) + Ey(i+1,j+1,k+1,1) + Ey(i+1,j+1,k+1,2) - - Ey(i+1,j+1,k ,0) - Ey(i+1,j+1,k ,1) - Ey(i+1,j+1,k ,2) - + Ey(i-1,j+1,k+1,0) + Ey(i-1,j+1,k+1,1) + Ey(i-1,j+1,k+1,2) - - Ey(i-1,j+1,k ,0) - Ey(i-1,j+1,k ,1) - Ey(i-1,j+1,k ,2) - + Ey(i+1,j-1,k+1,0) + Ey(i+1,j-1,k+1,1) + Ey(i+1,j-1,k+1,2) - - Ey(i+1,j-1,k ,0) - Ey(i+1,j-1,k ,1) - Ey(i+1,j-1,k ,2) - + Ey(i-1,j-1,k+1,0) + Ey(i-1,j-1,k+1,1) + Ey(i-1,j-1,k+1,2) - - Ey(i-1,j-1,k ,0) - Ey(i-1,j-1,k ,1) - Ey(i-1,j-1,k ,2)); - -#else - - Bx(i,j,k,1) += alphaz*(Ey(i ,j+1,k ,0) + Ey(i ,j+1,k ,1) + Ey(i ,j+1,k ,2) - - Ey(i ,j ,k ,0) - Ey(i ,j ,k ,1) - Ey(i ,j ,k ,2)) - +betazx*(Ey(i+1,j+1,k ,0) + Ey(i+1,j+1,k ,1) + Ey(i+1,j+1,k ,2) - - Ey(i+1,j ,k ,0) - Ey(i+1,j ,k ,1) - Ey(i+1,j ,k ,2) - + Ey(i-1,j+1,k ,0) + Ey(i-1,j+1,k ,1) + Ey(i-1,j+1,k ,2) - - Ey(i-1,j ,k ,0) - Ey(i-1,j ,k ,1) - Ey(i-1,j ,k ,2)); - - -#endif -} - - -AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE -void warpx_push_pml_by_ckc(int i, int j, int k, Array4<Real> const&By, - Array4<Real const> const& Ex, - Array4<Real const> const& Ez, - Real betaxy, Real betaxz, Real betayx, Real betayz, - Real betazx, Real betazy, Real gammax, Real gammay, - Real gammaz, Real alphax, Real alphay, Real alphaz) -{ -#if (AMREX_SPACEDIM==3) - By(i,j,k,0) += -alphaz*(Ex(i ,j ,k+1,0) + Ex(i ,j ,k+1,1) + Ex(i ,j ,k+1,2) - - Ex(i ,j ,k ,0) - Ex(i ,j ,k ,1) - Ex(i ,j ,k ,2)) - -betazx*(Ex(i+1,j ,k+1,0) + Ex(i+1,j ,k+1,1) + Ex(i+1,j ,k+1,2) - - Ex(i+1,j ,k ,0) - Ex(i+1,j ,k ,1) - Ex(i+1,j ,k ,2) - + Ex(i-1,j ,k+1,0) + Ex(i-1,j ,k+1,1) + Ex(i-1,j ,k+1,2) - - Ex(i-1,j ,k ,0) - Ex(i-1,j ,k ,1) - Ex(i-1,j ,k ,2)) - -betazy*(Ex(i ,j+1,k+1,0) + Ex(i ,j+1,k+1,1) + Ex(i ,j+1,k+1,2) - - Ex(i ,j+1,k ,0) - Ex(i ,j+1,k ,1) - Ex(i ,j+1,k ,2) - + Ex(i ,j-1,k+1,0) + Ex(i ,j-1,k+1,1) + Ex(i ,j-1,k+1,2) - - Ex(i ,j-1,k ,0) - Ex(i ,j-1,k ,1) - Ex(i ,j-1,k ,2)) - -gammaz*(Ex(i+1,j+1,k+1,0) + Ex(i+1,j+1,k+1,1) + Ex(i+1,j+1,k+1,2) - - Ex(i+1,j+1,k ,0) - Ex(i+1,j+1,k ,1) - Ex(i+1,j+1,k ,2) - + Ex(i-1,j+1,k+1,0) + Ex(i-1,j+1,k+1,1) + Ex(i-1,j+1,k+1,2) - - Ex(i-1,j+1,k ,0) - Ex(i-1,j+1,k ,1) - Ex(i-1,j+1,k ,2) - + Ex(i+1,j-1,k+1,0) + Ex(i+1,j-1,k+1,1) + Ex(i+1,j-1,k+1,2) - - Ex(i+1,j-1,k ,0) - Ex(i+1,j-1,k ,1) - Ex(i+1,j-1,k ,2) - + Ex(i-1,j-1,k+1,0) + Ex(i-1,j-1,k+1,1) + Ex(i-1,j-1,k+1,2) - - Ex(i-1,j-1,k ,0) - Ex(i-1,j-1,k ,1) - Ex(i-1,j-1,k ,2)); - - By(i,j,k,1) += alphax*(Ez(i+1,j ,k ,0) + Ez(i+1,j ,k ,1) + Ez(i+1,j ,k ,2) - - Ez(i ,j ,k ,0) - Ez(i ,j ,k ,1) - Ez(i ,j ,k ,2)) - +betaxy*(Ez(i+1,j+1,k ,0) + Ez(i+1,j+1,k ,1) + Ez(i+1,j+1,k ,2) - - Ez(i ,j+1,k ,0) - Ez(i ,j+1,k ,1) - Ez(i ,j+1,k ,2) - + Ez(i+1,j-1,k ,0) + Ez(i+1,j-1,k ,1) + Ez(i+1,j-1,k ,2) - - Ez(i ,j-1,k ,0) - Ez(i ,j-1,k ,1) - Ez(i ,j-1,k ,2)) - +betaxz*(Ez(i+1,j ,k+1,0) + Ez(i+1,j ,k+1,1) + Ez(i+1,j ,k+1,2) - - Ez(i ,j ,k+1,0) - Ez(i ,j ,k+1,1) - Ez(i ,j ,k+1,2) - + Ez(i+1,j ,k-1,0) + Ez(i+1,j ,k-1,1) + Ez(i+1,j ,k-1,2) - - Ez(i ,j ,k-1,0) - Ez(i ,j ,k-1,1) - Ez(i ,j ,k-1,2)) - +gammax*(Ez(i+1,j+1,k+1,0) + Ez(i+1,j+1,k+1,1) + Ez(i+1,j+1,k+1,2) - - Ez(i ,j+1,k+1,0) - Ez(i ,j+1,k+1,1) - Ez(i ,j+1,k+1,2) - + Ez(i+1,j-1,k+1,0) + Ez(i+1,j-1,k+1,1) + Ez(i+1,j-1,k+1,2) - - Ez(i ,j-1,k+1,0) - Ez(i ,j-1,k+1,1) - Ez(i ,j-1,k+1,2) - + Ez(i+1,j+1,k-1,0) + Ez(i+1,j+1,k-1,1) + Ez(i+1,j+1,k-1,2) - - Ez(i ,j+1,k-1,0) - Ez(i ,j+1,k-1,1) - Ez(i ,j+1,k-1,2) - + Ez(i+1,j-1,k-1,0) + Ez(i+1,j-1,k-1,1) + Ez(i+1,j-1,k-1,2) - - Ez(i ,j-1,k-1,0) - Ez(i ,j-1,k-1,1) - Ez(i ,j-1,k-1,2)); -#else - By(i,j,k,0) += -alphaz*(Ex(i ,j+1,k ,0) + Ex(i ,j+1,k ,1) + Ex(i ,j+1,k ,2) - - Ex(i ,j ,k ,0) - Ex(i ,j ,k ,1) - Ex(i ,j ,k ,2)) - -betazx*(Ex(i+1,j+1,k ,0) + Ex(i+1,j+1,k ,1) + Ex(i+1,j+1,k ,2) - - Ex(i+1,j ,k ,0) - Ex(i+1,j ,k ,1) - Ex(i+1,j ,k ,2) - + Ex(i-1,j+1,k ,0) + Ex(i-1,j+1,k ,1) + Ex(i-1,j+1,k ,2) - - Ex(i-1,j ,k ,0) - Ex(i-1,j ,k ,1) - Ex(i-1,j ,k ,2)); - - By(i,j,k,1) += alphax*(Ez(i+1,j ,k ,0) + Ez(i+1,j ,k ,1) + Ez(i+1,j ,k ,2) - - Ez(i ,j ,k ,0) - Ez(i ,j ,k ,1) - Ez(i ,j ,k ,2)) - +betaxz*(Ez(i+1,j+1,k ,0) + Ez(i+1,j+1,k ,1) + Ez(i+1,j+1,k ,2) - - Ez(i ,j+1,k ,0) - Ez(i ,j+1,k ,1) - Ez(i ,j+1,k ,2) - + Ez(i+1,j-1,k ,0) + Ez(i+1,j-1,k ,1) + Ez(i+1,j-1,k ,2) - - Ez(i ,j-1,k ,0) - Ez(i ,j-1,k ,1) - Ez(i ,j-1,k ,2)); - - -#endif - -} - - -AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE -void warpx_push_pml_bz_ckc(int i, int j, int k, Array4<Real> const&Bz, - Array4<Real const> const& Ex, - Array4<Real const> const& Ey, - Real betaxy, Real betaxz, Real betayx, Real betayz, - Real betazx, Real betazy, Real gammax, Real gammay, - Real gammaz, Real alphax, Real alphay, Real alphaz) -{ -#if (AMREX_SPACEDIM==3) - Bz(i,j,k,0) += -alphax*(Ey(i+1,j ,k ,0) + Ey(i+1,j ,k ,1) + Ey(i+1,j ,k ,2) - - Ey(i ,j ,k ,0) - Ey(i ,j ,k ,1) - Ey(i ,j ,k ,2)) - -betaxy*(Ey(i+1,j+1,k ,0) + Ey(i+1,j+1,k ,1) + Ey(i+1,j+1,k ,2) - - Ey(i ,j+1,k ,0) - Ey(i ,j+1,k ,1) - Ey(i ,j+1,k ,2) - + Ey(i+1,j-1,k ,0) + Ey(i+1,j-1,k ,1) + Ey(i+1,j-1,k ,2) - - Ey(i ,j-1,k ,0) - Ey(i ,j-1,k ,1) - Ey(i ,j-1,k ,2)) - -betaxz*(Ey(i+1,j ,k+1,0) + Ey(i+1,j ,k+1,1) + Ey(i+1,j ,k+1,2) - - Ey(i ,j ,k+1,0) - Ey(i ,j ,k+1,1) - Ey(i ,j ,k+1,2) - + Ey(i+1,j ,k-1,0) + Ey(i+1,j ,k-1,1) + Ey(i+1,j ,k-1,2) - - Ey(i ,j ,k-1,0) - Ey(i ,j ,k-1,1) - Ey(i ,j ,k-1,2)) - -gammax*(Ey(i+1,j+1,k+1,0) + Ey(i+1,j+1,k+1,1) + Ey(i+1,j+1,k+1,2) - - Ey(i ,j+1,k+1,0) - Ey(i ,j+1,k+1,1) - Ey(i ,j+1,k+1,2) - + Ey(i+1,j-1,k+1,0) + Ey(i+1,j-1,k+1,1) + Ey(i+1,j-1,k+1,2) - - Ey(i ,j-1,k+1,0) - Ey(i ,j-1,k+1,1) - Ey(i ,j-1,k+1,2) - + Ey(i+1,j+1,k-1,0) + Ey(i+1,j+1,k-1,1) + Ey(i+1,j+1,k-1,2) - - Ey(i ,j+1,k-1,0) - Ey(i ,j+1,k-1,1) - Ey(i ,j+1,k-1,2) - + Ey(i+1,j-1,k-1,0) + Ey(i+1,j-1,k-1,1) + Ey(i+1,j-1,k-1,2) - - Ey(i ,j-1,k-1,0) - Ey(i ,j-1,k-1,1) - Ey(i ,j-1,k-1,2)); - - Bz(i,j,k,1) += alphay*(Ex(i ,j+1,k ,0) + Ex(i ,j+1,k ,1) + Ex(i ,j+1,k ,2) - - Ex(i ,j ,k ,0) - Ex(i ,j ,k ,1) - Ex(i ,j ,k ,2)) - +betayx*(Ex(i+1,j+1,k ,0) + Ex(i+1,j+1,k ,1) + Ex(i+1,j+1,k ,2) - - Ex(i+1,j ,k ,0) - Ex(i+1,j ,k ,1) - Ex(i+1,j ,k ,2) - + Ex(i-1,j+1,k ,0) + Ex(i-1,j+1,k ,1) + Ex(i-1,j+1,k ,2) - - Ex(i-1,j ,k ,0) - Ex(i-1,j ,k ,1) - Ex(i-1,j ,k ,2)) - +betayz*(Ex(i ,j+1,k+1,0) + Ex(i ,j+1,k+1,1) + Ex(i ,j+1,k+1,2) - - Ex(i ,j ,k+1,0) - Ex(i ,j ,k+1,1) - Ex(i ,j ,k+1,2) - + Ex(i ,j+1,k-1,0) + Ex(i ,j+1,k-1,1) + Ex(i ,j+1,k-1,2) - - Ex(i ,j ,k-1,0) - Ex(i ,j ,k-1,1) - Ex(i ,j ,k-1,2)) - +gammay*(Ex(i+1,j+1,k+1,0) + Ex(i+1,j+1,k+1,1) + Ex(i+1,j+1,k+1,2) - - Ex(i+1,j ,k+1,0) - Ex(i+1,j ,k+1,1) - Ex(i+1,j ,k+1,2) - + Ex(i-1,j+1,k+1,0) + Ex(i-1,j+1,k+1,1) + Ex(i-1,j+1,k+1,2) - - Ex(i-1,j ,k+1,0) - Ex(i-1,j ,k+1,1) - Ex(i-1,j ,k+1,2) - + Ex(i+1,j+1,k-1,0) + Ex(i+1,j+1,k-1,1) + Ex(i+1,j+1,k-1,2) - - Ex(i+1,j ,k-1,0) - Ex(i+1,j ,k-1,1) - Ex(i+1,j ,k-1,2) - + Ex(i-1,j+1,k-1,0) + Ex(i-1,j+1,k-1,1) + Ex(i-1,j+1,k-1,2) - - Ex(i-1,j ,k-1,0) - Ex(i-1,j ,k-1,1) - Ex(i-1,j ,k-1,2)); - - -#else - Bz(i,j,k,0) += -alphax*(Ey(i+1,j ,k ,0) + Ey(i+1,j ,k ,1) + Ey(i+1,j ,k ,2) - - Ey(i ,j ,k ,0) - Ey(i ,j ,k ,1) - Ey(i ,j ,k ,2)) - -betaxz*(Ey(i+1,j+1,k ,0) + Ey(i+1,j+1,k ,1) + Ey(i+1,j+1,k ,2) - - Ey(i ,j+1,k ,0) - Ey(i ,j+1,k ,1) - Ey(i ,j+1,k ,2) - + Ey(i+1,j-1,k ,0) + Ey(i+1,j-1,k ,1) + Ey(i+1,j-1,k ,2) - - Ey(i ,j-1,k ,0) - Ey(i ,j-1,k ,1) - Ey(i ,j-1,k ,2)); -#endif - -} - - - - - -//PML EX YEE -AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE -void warpx_push_pml_ex_yee (int i, int j, int k, Array4<Real> const& Ex, - Array4<Real const> const& By, - Array4<Real const> const& Bz, - Real dtsdy_c2, Real dtsdz_c2) -{ -#if (AMREX_SPACEDIM == 3) - Ex(i,j,k,0) += dtsdy_c2 * ( Bz(i,j ,k,0) + Bz(i,j ,k,1) - - Bz(i,j-1,k,0) - Bz(i,j-1,k,1) ); - Ex(i,j,k,1) += -dtsdz_c2 * (By(i,j,k ,0) + By(i,j,k ,1) - - By(i,j,k-1,0) - By(i,j,k-1,1) ); - -#else - Ex(i,j,k,1) += -dtsdz_c2 * (By(i,j ,k,0) + By(i,j ,k,1) - - By(i,j-1,k,0) - By(i,j-1,k,1) ); -#endif -} - -//PML EY YEE -AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE -void warpx_push_pml_ey_yee (int i, int j, int k, Array4<Real> const& Ey, - Array4<Real const> const& Bx, - Array4<Real const> const& Bz, - Real dtsdx_c2, Real dtsdz_c2) -{ -#if (AMREX_SPACEDIM == 3) - Ey(i,j,k,0) += dtsdz_c2 * (Bx(i,j,k ,0) + Bx(i,j,k ,1) - - Bx(i,j,k-1,0) - Bx(i,j,k-1,1)); -#else - Ey(i,j,k,0) += dtsdz_c2 * (Bx(i,j ,k,0) + Bx(i,j ,k,1) - - Bx(i,j-1,k,0) - Bx(i,j-1,k,1)); -#endif - Ey(i,j,k,1) += -dtsdx_c2 * (Bz(i ,j,k,0) + Bz(i ,j,k,1) - - Bz(i-1,j,k,0) - Bz(i-1,j,k,1) ); -} - -//PML EZ YEE -AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE -void warpx_push_pml_ez_yee (int i, int j, int k, Array4<Real> const& Ez, - Array4<Real const> const& Bx, - Array4<Real const> const& By, - Real dtsdx_c2, Real dtsdy_c2) -{ - Ez(i,j,k,0) += dtsdx_c2 * (By(i, j,k,0) + By(i, j,k,1) - - By(i-1,j,k,0) - By(i-1,j,k,1) ); -#if (AMREX_SPACEDIM==3) - Ez(i,j,k,1) += - dtsdy_c2 * (Bx(i,j, k,0) + Bx(i,j, k,1) - - Bx(i,j-1,k,0) - Bx(i,j-1,k,1) ); -#endif -} - - -// PML EX F YEE -AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE -void warpx_push_pml_ex_f_yee(int i, int j, int k, - Array4<Real> const& Exfab, - Array4<Real const> const& F_fab, - Real dtsdx_c2) -{ - Exfab(i,j,k,2) += dtsdx_c2 * (F_fab(i+1,j,k,0) - F_fab(i,j,k,0) - + F_fab(i+1,j,k,1) - F_fab(i,j,k,1) - + F_fab(i+1,j,k,2) - F_fab(i,j,k,2)); -} - -// PML EY F YEE -AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE -void warpx_push_pml_ey_f_yee(int i, int j, int k, - Array4<Real> const& Eyfab, - Array4<Real const> const& F_fab, - Real dtsdy_c2) -{ -#if (AMREX_SPACEDIM == 3) - Eyfab(i,j,k,2) += dtsdy_c2 * (F_fab(i,j+1,k,0) - F_fab(i,j,k,0) - + F_fab(i,j+1,k,1) - F_fab(i,j,k,1) - + F_fab(i,j+1,k,2) - F_fab(i,j,k,2)); -#endif -} - -// PML EZ F YEE -AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE -void warpx_push_pml_ez_f_yee(int i, int j, int k, - Array4<Real> const& Ezfab, - Array4<Real const> const& F_fab, - Real dtsdz_c2) -{ -#if (AMREX_SPACEDIM == 3) - Ezfab(i,j,k,2) += dtsdz_c2 * (F_fab(i,j,k+1,0) - F_fab(i,j,k,0) - + F_fab(i,j,k+1,1) - F_fab(i,j,k,1) - + F_fab(i,j,k+1,2) - F_fab(i,j,k,2)); -#else - Ezfab(i,j,k,2) += dtsdz_c2 * (F_fab(i,j+1,k,0) - F_fab(i,j,k,0) - + F_fab(i,j+1,k,1) - F_fab(i,j,k,1) - + F_fab(i,j+1,k,2) - F_fab(i,j,k,2)); -#endif -} - - -AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE -void warpx_push_pml_ex_f_ckc(int i, int j, int k, - Array4<Real> const& Exfab, - Array4<Real const> const& F, - Real alphax, Real betaxy, Real betaxz,Real gammax) -{ -#if (AMREX_SPACEDIM==3) - Exfab(i,j,k,2) += alphax*( - F(i+1,j ,k ,0) + F(i+1,j ,k ,1) + F(i+1,j ,k ,2) - - F(i ,j ,k ,0) - F(i ,j ,k ,1) - F(i ,j ,k ,2)) - +betaxy*( - F(i+1,j+1,k ,0) + F(i+1,j+1,k ,1) + F(i+1,j+1,k ,2) - - F(i ,j+1,k ,0) - F(i ,j+1,k ,1) - F(i ,j+1,k ,2) - + F(i+1,j-1,k ,0) + F(i+1,j-1,k ,1) + F(i+1,j-1,k ,2) - - F(i ,j-1,k ,0) - F(i ,j-1,k ,1) - F(i ,j-1,k ,2)) - +betaxz*( - F(i+1,j ,k+1,0) + F(i+1,j ,k+1,1) + F(i+1,j ,k+1,2) - - F(i ,j ,k+1,0) - F(i ,j ,k+1,1) - F(i ,j ,k+1,2) - + F(i+1,j ,k-1,0) + F(i+1,j ,k-1,1) + F(i+1,j ,k-1,2) - - F(i ,j ,k-1,0) - F(i ,j ,k-1,1) - F(i ,j ,k-1,2)) - +gammax*( - F(i+1,j+1,k+1,0) + F(i+1,j+1,k+1,1) + F(i+1,j+1,k+1,2) - - F(i ,j+1,k+1,0) - F(i ,j+1,k+1,1) - F(i ,j+1,k+1,2) - + F(i+1,j-1,k+1,0) + F(i+1,j-1,k+1,1) + F(i+1,j-1,k+1,2) - - F(i ,j-1,k+1,0) - F(i ,j-1,k+1,1) - F(i ,j-1,k+1,2) - + F(i+1,j+1,k-1,0) + F(i+1,j+1,k-1,1) + F(i+1,j+1,k-1,2) - - F(i ,j+1,k-1,0) - F(i ,j+1,k-1,1) - F(i ,j+1,k-1,2) - + F(i+1,j-1,k-1,0) + F(i+1,j-1,k-1,1) + F(i+1,j-1,k-1,2) - - F(i ,j-1,k-1,0) - F(i ,j-1,k-1,1) - F(i ,j-1,k-1,2)); -#else - Exfab(i,j,k,2) += alphax*( - F(i+1,j ,k ,0) + F(i+1,j ,k ,1) + F(i+1,j ,k ,2) - - F(i ,j ,k ,0) - F(i ,j ,k ,1) - F(i ,j ,k ,2)) - +betaxz*( - F(i+1,j+1,k ,0) + F(i+1,j+1,k ,1) + F(i+1,j+1,k ,2) - - F(i ,j+1,k ,0) - F(i ,j+1,k ,1) - F(i ,j+1,k ,2) - + F(i+1,j-1,k ,0) + F(i+1,j-1,k ,1) + F(i+1,j-1,k ,2) - - F(i ,j-1,k ,0) - F(i ,j-1,k ,1) - F(i ,j-1,k ,2)); -#endif - -} - -AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE -void warpx_push_pml_ey_f_ckc(int i, int j, int k, - Array4<Real> const& Eyfab, - Array4<Real const> const& F, - Real alphay, Real betayx, Real betayz,Real gammay) -{ -#if (AMREX_SPACEDIM==3) - Eyfab(i,j,k,2) += alphay*( - F(i ,j+1,k ,0) + F(i ,j+1,k ,1) + F(i ,j+1,k ,2) - - F(i ,j ,k ,0) - F(i ,j ,k ,1) - F(i ,j ,k ,2)) - +betayx*( - F(i+1,j+1,k ,0) + F(i+1,j+1,k ,1) + F(i+1,j+1,k ,2) - - F(i+1,j ,k ,0) - F(i+1,j ,k ,1) - F(i+1,j ,k ,2) - + F(i-1,j+1,k ,0) + F(i-1,j+1,k ,1) + F(i-1,j+1,k ,2) - - F(i-1,j ,k ,0) - F(i-1,j ,k ,1) - F(i-1,j ,k ,2)) - +betayz*( - F(i ,j+1,k+1,0) + F(i ,j+1,k+1,1) + F(i ,j+1,k+1,2) - - F(i ,j ,k+1,0) - F(i ,j ,k+1,1) - F(i ,j ,k+1,2) - + F(i ,j+1,k-1,0) + F(i ,j+1,k-1,1) + F(i ,j+1,k-1,2) - - F(i ,j ,k-1,0) - F(i ,j ,k-1,1) - F(i ,j ,k-1,2)) - +gammay*( - F(i+1,j+1,k+1,0) + F(i+1,j+1,k+1,1) + F(i+1,j+1,k+1,2) - - F(i+1,j ,k+1,0) - F(i+1,j ,k+1,1) - F(i+1,j ,k+1,2) - + F(i-1,j+1,k+1,0) + F(i-1,j+1,k+1,1) + F(i-1,j+1,k+1,2) - - F(i-1,j ,k+1,0) - F(i-1,j ,k+1,1) - F(i-1,j ,k+1,2) - + F(i+1,j+1,k-1,0) + F(i+1,j+1,k-1,1) + F(i+1,j+1,k-1,2) - - F(i+1,j ,k-1,0) - F(i+1,j ,k-1,1) - F(i+1,j ,k-1,2) - + F(i-1,j+1,k-1,0) + F(i-1,j+1,k-1,1) + F(i-1,j+1,k-1,2) - - F(i-1,j ,k-1,0) - F(i-1,j ,k-1,1) - F(i-1,j ,k-1,2)); -#endif -} - -AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE -void warpx_push_pml_ez_f_ckc(int i, int j, int k, - Array4<Real> const& Ezfab, - Array4<Real const> const& F, - Real alphaz, Real betazx, Real betazy,Real gammaz) -{ -#if (AMREX_SPACEDIM==3) - Ezfab(i,j,k,2) += alphaz*( - F(i ,j ,k+1,0) + F(i ,j ,k+1,1) + F(i ,j ,k+1,2) - - F(i ,j ,k ,0) - F(i ,j ,k ,1) - F(i ,j ,k ,2)) - +betazx*( - F(i+1,j ,k+1,0) + F(i+1,j ,k+1,1) + F(i+1,j ,k+1,2) - - F(i+1,j ,k ,0) - F(i+1,j ,k ,1) - F(i+1,j ,k ,2) - + F(i-1,j ,k+1,0) + F(i-1,j ,k+1,1) + F(i-1,j ,k+1,2) - - F(i-1,j ,k ,0) - F(i-1,j ,k ,1) - F(i-1,j ,k ,2)) - +betazy*( - F(i ,j+1,k+1,0) + F(i ,j+1,k+1,1) + F(i ,j+1,k+1,2) - - F(i ,j+1,k ,0) - F(i ,j+1,k ,1) - F(i ,j+1,k ,2) - + F(i ,j-1,k+1,0) + F(i ,j-1,k+1,1) + F(i ,j-1,k+1,2) - - F(i ,j-1,k ,0) - F(i ,j-1,k ,1) - F(i ,j-1,k ,2)) - +gammaz*( - F(i+1,j+1,k+1,0) + F(i+1,j+1,k+1,1) + F(i+1,j+1,k+1,2) - - F(i+1,j+1,k ,0) - F(i+1,j+1,k ,1) - F(i+1,j+1,k ,2) - + F(i-1,j+1,k+1,0) + F(i-1,j+1,k+1,1) + F(i-1,j+1,k+1,2) - - F(i-1,j+1,k ,0) - F(i-1,j+1,k ,1) - F(i-1,j+1,k ,2) - + F(i+1,j-1,k+1,0) + F(i+1,j-1,k+1,1) + F(i+1,j-1,k+1,2) - - F(i+1,j-1,k ,0) - F(i+1,j-1,k ,1) - F(i+1,j-1,k ,2) - + F(i-1,j-1,k+1,0) + F(i-1,j-1,k+1,1) + F(i-1,j-1,k+1,2) - - F(i-1,j-1,k ,0) - F(i-1,j-1,k ,1) - F(i-1,j-1,k ,2)); -#else - Ezfab(i,j,k,2) += alphaz*( - F(i ,j+1,k ,0) + F(i ,j+1,k ,1) + F(i ,j+1,k ,2) - - F(i ,j ,k ,0) - F(i ,j ,k ,1) - F(i ,j ,k ,2)) - +betazx*( - F(i+1,j+1,k ,0) + F(i+1,j+1,k ,1) + F(i+1,j+1,k ,2) - - F(i+1,j ,k ,0) - F(i+1,j ,k ,1) - F(i+1,j ,k ,2) - + F(i-1,j+1,k ,0) + F(i-1,j+1,k ,1) + F(i-1,j+1,k ,2) - - F(i-1,j ,k ,0) - F(i-1,j ,k ,1) - F(i-1,j ,k ,2)); -#endif - -} - - -AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE -void warpx_push_pml_F(int i, int j, int k, Array4<Real> const& F_fab, - Array4<Real const> const& Exfab, - Array4<Real const> const& Eyfab, - Array4<Real const> const& Ezfab, - Real dtsdx, Real dtsdy, Real dtsdz) -{ - F_fab(i,j,k,0) += dtsdx * (Exfab(i,j,k,0) - Exfab(i-1,j,k,0) - + Exfab(i,j,k,1) - Exfab(i-1,j,k,1) - + Exfab(i,j,k,2) - Exfab(i-1,j,k,2)); -#if (AMREX_SPACEDIM==3) - F_fab(i,j,k,1) += dtsdy * (Eyfab(i,j,k,0) - Eyfab(i,j-1,k,0) - + Eyfab(i,j,k,1) - Eyfab(i,j-1,k,1) - + Eyfab(i,j,k,2) - Eyfab(i,j-1,k,2)); - F_fab(i,j,k,2) += dtsdz * (Ezfab(i,j,k,0) - Ezfab(i,j,k-1,0) - + Ezfab(i,j,k,1) - Ezfab(i,j,k-1,1) - + Ezfab(i,j,k,2) - Ezfab(i,j,k-1,2)); -#else - F_fab(i,j,k,2) += dtsdz * (Ezfab(i,j,k,0) - Ezfab(i,j-1,k,0) - + Ezfab(i,j,k,1) - Ezfab(i,j-1,k,1) - + Ezfab(i,j,k,2) - Ezfab(i,j-1,k,2)); -#endif -} - - AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void warpx_damp_pml_ex (int i, int j, int k, Array4<Real> const& Ex, const Real* const sigma_fac_y, diff --git a/Source/FieldSolver/FiniteDifferenceSolver/EvolveBPML.cpp b/Source/FieldSolver/FiniteDifferenceSolver/EvolveBPML.cpp new file mode 100644 index 000000000..6282892ab --- /dev/null +++ b/Source/FieldSolver/FiniteDifferenceSolver/EvolveBPML.cpp @@ -0,0 +1,131 @@ +/* Copyright 2020 Remi Lehe + * + * This file is part of WarpX. + * + * License: BSD-3-Clause-LBNL + */ + +#include "Utils/WarpXAlgorithmSelection.H" +#include "FieldSolver/FiniteDifferenceSolver/FiniteDifferenceSolver.H" +#ifdef WARPX_DIM_RZ +# include "FieldSolver/FiniteDifferenceSolver/FiniteDifferenceAlgorithms/CylindricalYeeAlgorithm.H" +#else +# include "FieldSolver/FiniteDifferenceSolver/FiniteDifferenceAlgorithms/CartesianYeeAlgorithm.H" +# include "FieldSolver/FiniteDifferenceSolver/FiniteDifferenceAlgorithms/CartesianCKCAlgorithm.H" +# include "FieldSolver/FiniteDifferenceSolver/FiniteDifferenceAlgorithms/CartesianNodalAlgorithm.H" +#endif +#include "BoundaryConditions/PMLComponent.H" +#include <AMReX_Gpu.H> + +using namespace amrex; + +/** + * \brief Update the B field, over one timestep + */ +void FiniteDifferenceSolver::EvolveBPML ( + std::array< amrex::MultiFab*, 3 > Bfield, + std::array< amrex::MultiFab*, 3 > const Efield, + amrex::Real const dt ) { + + // Select algorithm (The choice of algorithm is a runtime option, + // but we compile code for each algorithm, using templates) +#ifdef WARPX_DIM_RZ + amrex::Abort("PML are not implemented in cylindrical geometry."); +#else + if (m_do_nodal) { + + EvolveBPMLCartesian <CartesianNodalAlgorithm> ( Bfield, Efield, dt ); + + } else if (m_fdtd_algo == MaxwellSolverAlgo::Yee) { + + EvolveBPMLCartesian <CartesianYeeAlgorithm> ( Bfield, Efield, dt ); + + } else if (m_fdtd_algo == MaxwellSolverAlgo::CKC) { + + EvolveBPMLCartesian <CartesianCKCAlgorithm> ( Bfield, Efield, dt ); + + } else { + amrex::Abort("Unknown algorithm"); + } +#endif +} + + +#ifndef WARPX_DIM_RZ + +template<typename T_Algo> +void FiniteDifferenceSolver::EvolveBPMLCartesian ( + std::array< amrex::MultiFab*, 3 > Bfield, + std::array< amrex::MultiFab*, 3 > const Efield, + amrex::Real const dt ) { + + // Loop through the grids, and over the tiles within each grid +#ifdef _OPENMP +#pragma omp parallel if (amrex::Gpu::notInLaunchRegion()) +#endif + for ( MFIter mfi(*Bfield[0], TilingIfNotGPU()); mfi.isValid(); ++mfi ) { + + // Extract field data for this grid/tile + Array4<Real> const& Bx = Bfield[0]->array(mfi); + Array4<Real> const& By = Bfield[1]->array(mfi); + Array4<Real> const& Bz = Bfield[2]->array(mfi); + Array4<Real> const& Ex = Efield[0]->array(mfi); + Array4<Real> const& Ey = Efield[1]->array(mfi); + Array4<Real> const& Ez = Efield[2]->array(mfi); + + // Extract stencil coefficients + Real const * const AMREX_RESTRICT coefs_x = m_stencil_coefs_x.dataPtr(); + int const n_coefs_x = m_stencil_coefs_x.size(); + Real const * const AMREX_RESTRICT coefs_y = m_stencil_coefs_y.dataPtr(); + int const n_coefs_y = m_stencil_coefs_y.size(); + Real const * const AMREX_RESTRICT coefs_z = m_stencil_coefs_z.dataPtr(); + int const n_coefs_z = m_stencil_coefs_z.size(); + + // Extract tileboxes for which to loop + Box const& tbx = mfi.tilebox(Bfield[0]->ixType().ixType()); + Box const& tby = mfi.tilebox(Bfield[1]->ixType().ixType()); + Box const& tbz = mfi.tilebox(Bfield[2]->ixType().ixType()); + + // Loop over the cells and update the fields + amrex::ParallelFor(tbx, tby, tbz, + + [=] AMREX_GPU_DEVICE (int i, int j, int k){ + Bx(i, j, k, PMLComp::xz) += dt * ( + T_Algo::UpwardDz(Ey, coefs_z, n_coefs_z, i, j, k, PMLComp::yx) + + T_Algo::UpwardDz(Ey, coefs_z, n_coefs_z, i, j, k, PMLComp::yy) + + T_Algo::UpwardDz(Ey, coefs_z, n_coefs_z, i, j, k, PMLComp::yz) ); + Bx(i, j, k, PMLComp::xy) -= dt * ( + T_Algo::UpwardDy(Ez, coefs_y, n_coefs_y, i, j, k, PMLComp::zx) + + T_Algo::UpwardDy(Ez, coefs_y, n_coefs_y, i, j, k, PMLComp::zy) + + T_Algo::UpwardDy(Ez, coefs_y, n_coefs_y, i, j, k, PMLComp::zz) ); + }, + + [=] AMREX_GPU_DEVICE (int i, int j, int k){ + By(i, j, k, PMLComp::yx) += dt * ( + T_Algo::UpwardDx(Ez, coefs_x, n_coefs_x, i, j, k, PMLComp::zx) + + T_Algo::UpwardDx(Ez, coefs_x, n_coefs_x, i, j, k, PMLComp::zy) + + T_Algo::UpwardDx(Ez, coefs_x, n_coefs_x, i, j, k, PMLComp::zz) ); + By(i, j, k, PMLComp::yz) -= dt * ( + T_Algo::UpwardDz(Ex, coefs_z, n_coefs_z, i, j, k, PMLComp::xx) + + T_Algo::UpwardDz(Ex, coefs_z, n_coefs_z, i, j, k, PMLComp::xy) + + T_Algo::UpwardDz(Ex, coefs_z, n_coefs_z, i, j, k, PMLComp::xz) ); + }, + + [=] AMREX_GPU_DEVICE (int i, int j, int k){ + Bz(i, j, k, PMLComp::zy) += dt * ( + T_Algo::UpwardDy(Ex, coefs_y, n_coefs_y, i, j, k, PMLComp::xx) + + T_Algo::UpwardDy(Ex, coefs_y, n_coefs_y, i, j, k, PMLComp::xy) + + T_Algo::UpwardDy(Ex, coefs_y, n_coefs_y, i, j, k, PMLComp::xz) ); + Bz(i, j, k, PMLComp::zx) -= dt * ( + T_Algo::UpwardDx(Ey, coefs_x, n_coefs_x, i, j, k, PMLComp::yx) + + T_Algo::UpwardDx(Ey, coefs_x, n_coefs_x, i, j, k, PMLComp::yy) + + T_Algo::UpwardDx(Ey, coefs_x, n_coefs_x, i, j, k, PMLComp::yz) ); + } + + ); + + } + +} + +#endif // corresponds to ifndef WARPX_DIM_RZ diff --git a/Source/FieldSolver/FiniteDifferenceSolver/EvolveEPML.cpp b/Source/FieldSolver/FiniteDifferenceSolver/EvolveEPML.cpp new file mode 100644 index 000000000..34144ff0d --- /dev/null +++ b/Source/FieldSolver/FiniteDifferenceSolver/EvolveEPML.cpp @@ -0,0 +1,206 @@ +/* Copyright 2020 Remi Lehe + * + * This file is part of WarpX. + * + * License: BSD-3-Clause-LBNL + */ + +#include "Utils/WarpXAlgorithmSelection.H" +#include "FieldSolver/FiniteDifferenceSolver/FiniteDifferenceSolver.H" +#ifdef WARPX_DIM_RZ +# include "FieldSolver/FiniteDifferenceSolver/FiniteDifferenceAlgorithms/CylindricalYeeAlgorithm.H" +#else +# include "FieldSolver/FiniteDifferenceSolver/FiniteDifferenceAlgorithms/CartesianYeeAlgorithm.H" +# include "FieldSolver/FiniteDifferenceSolver/FiniteDifferenceAlgorithms/CartesianCKCAlgorithm.H" +# include "FieldSolver/FiniteDifferenceSolver/FiniteDifferenceAlgorithms/CartesianNodalAlgorithm.H" +#endif +#include "BoundaryConditions/PML.H" +#include "BoundaryConditions/PML_current.H" +#include "BoundaryConditions/PMLComponent.H" +#include "Utils/WarpXConst.H" +#include <AMReX_Gpu.H> + +using namespace amrex; + +/** + * \brief Update the E field, over one timestep + */ +void FiniteDifferenceSolver::EvolveEPML ( + std::array< amrex::MultiFab*, 3 > Efield, + std::array< amrex::MultiFab*, 3 > const Bfield, + std::array< amrex::MultiFab*, 3 > const Jfield, + amrex::MultiFab* const Ffield, + MultiSigmaBox const& sigba, + amrex::Real const dt, bool pml_has_particles ) { + + // Select algorithm (The choice of algorithm is a runtime option, + // but we compile code for each algorithm, using templates) +#ifdef WARPX_DIM_RZ + amrex::Abort("PML are not implemented in cylindrical geometry."); +#else + if (m_do_nodal) { + + EvolveEPMLCartesian <CartesianNodalAlgorithm> ( + Efield, Bfield, Jfield, Ffield, sigba, dt, pml_has_particles ); + + } else if (m_fdtd_algo == MaxwellSolverAlgo::Yee) { + + EvolveEPMLCartesian <CartesianYeeAlgorithm> ( + Efield, Bfield, Jfield, Ffield, sigba, dt, pml_has_particles ); + + } else if (m_fdtd_algo == MaxwellSolverAlgo::CKC) { + + EvolveEPMLCartesian <CartesianCKCAlgorithm> ( + Efield, Bfield, Jfield, Ffield, sigba, dt, pml_has_particles ); + + } else { + amrex::Abort("Unknown algorithm"); + } +#endif +} + + +#ifndef WARPX_DIM_RZ + +template<typename T_Algo> +void FiniteDifferenceSolver::EvolveEPMLCartesian ( + std::array< amrex::MultiFab*, 3 > Efield, + std::array< amrex::MultiFab*, 3 > const Bfield, + std::array< amrex::MultiFab*, 3 > const Jfield, + amrex::MultiFab* const Ffield, + MultiSigmaBox const& sigba, + amrex::Real const dt, bool pml_has_particles ) { + + Real constexpr c2 = PhysConst::c * PhysConst::c; + + // Loop through the grids, and over the tiles within each grid +#ifdef _OPENMP +#pragma omp parallel if (amrex::Gpu::notInLaunchRegion()) +#endif + for ( MFIter mfi(*Efield[0], TilingIfNotGPU()); mfi.isValid(); ++mfi ) { + + // Extract field data for this grid/tile + Array4<Real> const& Ex = Efield[0]->array(mfi); + Array4<Real> const& Ey = Efield[1]->array(mfi); + Array4<Real> const& Ez = Efield[2]->array(mfi); + Array4<Real> const& Bx = Bfield[0]->array(mfi); + Array4<Real> const& By = Bfield[1]->array(mfi); + Array4<Real> const& Bz = Bfield[2]->array(mfi); + + // Extract stencil coefficients + Real const * const AMREX_RESTRICT coefs_x = m_stencil_coefs_x.dataPtr(); + int const n_coefs_x = m_stencil_coefs_x.size(); + Real const * const AMREX_RESTRICT coefs_y = m_stencil_coefs_y.dataPtr(); + int const n_coefs_y = m_stencil_coefs_y.size(); + Real const * const AMREX_RESTRICT coefs_z = m_stencil_coefs_z.dataPtr(); + int const n_coefs_z = m_stencil_coefs_z.size(); + + // Extract tileboxes for which to loop + Box const& tex = mfi.tilebox(Efield[0]->ixType().ixType()); + Box const& tey = mfi.tilebox(Efield[1]->ixType().ixType()); + Box const& tez = mfi.tilebox(Efield[2]->ixType().ixType()); + + // Loop over the cells and update the fields + amrex::ParallelFor(tex, tey, tez, + + [=] AMREX_GPU_DEVICE (int i, int j, int k){ + Ex(i, j, k, PMLComp::xz) -= c2 * dt * ( + T_Algo::DownwardDz(By, coefs_z, n_coefs_z, i, j, k, PMLComp::yx) + + T_Algo::DownwardDz(By, coefs_z, n_coefs_z, i, j, k, PMLComp::yz) ); + Ex(i, j, k, PMLComp::xy) += c2 * dt * ( + T_Algo::DownwardDy(Bz, coefs_y, n_coefs_y, i, j, k, PMLComp::zx) + + T_Algo::DownwardDy(Bz, coefs_y, n_coefs_y, i, j, k, PMLComp::zy) ); + }, + + [=] AMREX_GPU_DEVICE (int i, int j, int k){ + Ey(i, j, k, PMLComp::yx) -= c2 * dt * ( + T_Algo::DownwardDx(Bz, coefs_x, n_coefs_x, i, j, k, PMLComp::zx) + + T_Algo::DownwardDx(Bz, coefs_x, n_coefs_x, i, j, k, PMLComp::zy) ); + Ey(i, j, k, PMLComp::yz) += c2 * dt * ( + T_Algo::DownwardDz(Bx, coefs_z, n_coefs_z, i, j, k, PMLComp::xy) + + T_Algo::DownwardDz(Bx, coefs_z, n_coefs_z, i, j, k, PMLComp::xz) ); + }, + + [=] AMREX_GPU_DEVICE (int i, int j, int k){ + Ez(i, j, k, PMLComp::zy) -= c2 * dt * ( + T_Algo::DownwardDy(Bx, coefs_y, n_coefs_y, i, j, k, PMLComp::xy) + + T_Algo::DownwardDy(Bx, coefs_y, n_coefs_y, i, j, k, PMLComp::xz) ); + Ez(i, j, k, PMLComp::zx) += c2 * dt * ( + T_Algo::DownwardDx(By, coefs_x, n_coefs_x, i, j, k, PMLComp::yx) + + T_Algo::DownwardDx(By, coefs_x, n_coefs_x, i, j, k, PMLComp::yz) ); + } + + ); + + // If F is not a null pointer, further update E using the grad(F) term + // (hyperbolic correction for errors in charge conservation) + if (Ffield) { + // Extract field data for this grid/tile + Array4<Real> const& F = Ffield->array(mfi); + + // Loop over the cells and update the fields + amrex::ParallelFor(tex, tey, tez, + + [=] AMREX_GPU_DEVICE (int i, int j, int k){ + Ex(i, j, k, PMLComp::xx) += c2 * dt * ( + T_Algo::UpwardDx(F, coefs_x, n_coefs_x, i, j, k, PMLComp::x) + + T_Algo::UpwardDx(F, coefs_x, n_coefs_x, i, j, k, PMLComp::y) + + T_Algo::UpwardDx(F, coefs_x, n_coefs_x, i, j, k, PMLComp::z) ); + }, + [=] AMREX_GPU_DEVICE (int i, int j, int k){ + Ey(i, j, k, PMLComp::yy) += c2 * dt * ( + T_Algo::UpwardDy(F, coefs_y, n_coefs_y, i, j, k, PMLComp::x) + + T_Algo::UpwardDy(F, coefs_y, n_coefs_y, i, j, k, PMLComp::y) + + T_Algo::UpwardDy(F, coefs_y, n_coefs_y, i, j, k, PMLComp::z) ); + }, + [=] AMREX_GPU_DEVICE (int i, int j, int k){ + Ez(i, j, k, PMLComp::zz) += c2 * dt * ( + T_Algo::UpwardDz(F, coefs_z, n_coefs_z, i, j, k, PMLComp::x) + + T_Algo::UpwardDz(F, coefs_z, n_coefs_z, i, j, k, PMLComp::y) + + T_Algo::UpwardDz(F, coefs_z, n_coefs_z, i, j, k, PMLComp::z) ); + } + ); + } + + // Update the E field in the PML, using the current + // deposited by the particles in the PML + if (pml_has_particles) { + + // Extract field data for this grid/tile + Array4<Real> const& Jx = Jfield[0]->array(mfi); + Array4<Real> const& Jy = Jfield[1]->array(mfi); + Array4<Real> const& Jz = Jfield[2]->array(mfi); + const Real* sigmaj_x = sigba[mfi].sigma[0].data(); + const Real* sigmaj_y = sigba[mfi].sigma[1].data(); + const Real* sigmaj_z = sigba[mfi].sigma[2].data(); + int const x_lo = sigba[mfi].sigma[0].lo(); +#if (AMREX_SPACEDIM == 3) + int const y_lo = sigba[mfi].sigma[1].lo(); + int const z_lo = sigba[mfi].sigma[2].lo(); +#else + int const y_lo = 0; + int const z_lo = sigba[mfi].sigma[1].lo(); +#endif + const Real mu_c2_dt = (PhysConst::mu0*PhysConst::c*PhysConst::c) * dt; + + amrex::ParallelFor( tex, tey, tez, + [=] AMREX_GPU_DEVICE (int i, int j, int k) { + push_ex_pml_current(i, j, k, Ex, Jx, + sigmaj_y, sigmaj_z, y_lo, z_lo, mu_c2_dt); + }, + [=] AMREX_GPU_DEVICE (int i, int j, int k) { + push_ey_pml_current(i, j, k, Ey, Jy, + sigmaj_x, sigmaj_z, x_lo, z_lo, mu_c2_dt); + }, + [=] AMREX_GPU_DEVICE (int i, int j, int k) { + push_ez_pml_current(i, j, k, Ez, Jz, + sigmaj_x, sigmaj_y, x_lo, y_lo, mu_c2_dt); + } + ); + } + + } + +} + +#endif // corresponds to ifndef WARPX_DIM_RZ diff --git a/Source/FieldSolver/FiniteDifferenceSolver/EvolveFPML.cpp b/Source/FieldSolver/FiniteDifferenceSolver/EvolveFPML.cpp new file mode 100644 index 000000000..b36c50c35 --- /dev/null +++ b/Source/FieldSolver/FiniteDifferenceSolver/EvolveFPML.cpp @@ -0,0 +1,113 @@ +/* Copyright 2020 Remi Lehe + * + * This file is part of WarpX. + * + * License: BSD-3-Clause-LBNL + */ + +#include "Utils/WarpXAlgorithmSelection.H" +#include "FieldSolver/FiniteDifferenceSolver/FiniteDifferenceSolver.H" +#ifdef WARPX_DIM_RZ +# include "FieldSolver/FiniteDifferenceSolver/FiniteDifferenceAlgorithms/CylindricalYeeAlgorithm.H" +#else +# include "FieldSolver/FiniteDifferenceSolver/FiniteDifferenceAlgorithms/CartesianYeeAlgorithm.H" +# include "FieldSolver/FiniteDifferenceSolver/FiniteDifferenceAlgorithms/CartesianCKCAlgorithm.H" +# include "FieldSolver/FiniteDifferenceSolver/FiniteDifferenceAlgorithms/CartesianNodalAlgorithm.H" +#endif +#include "BoundaryConditions/PMLComponent.H" +#include <AMReX_Gpu.H> + +using namespace amrex; + +/** + * \brief Update the E field, over one timestep + */ +void FiniteDifferenceSolver::EvolveFPML ( + amrex::MultiFab* Ffield, + std::array< amrex::MultiFab*, 3 > const Efield, + amrex::Real const dt ) { + + // Select algorithm (The choice of algorithm is a runtime option, + // but we compile code for each algorithm, using templates) +#ifdef WARPX_DIM_RZ + amrex::Abort("PML are not implemented in cylindrical geometry."); +#else + if (m_do_nodal) { + + EvolveFPMLCartesian <CartesianNodalAlgorithm> ( Ffield, Efield, dt ); + + } else if (m_fdtd_algo == MaxwellSolverAlgo::Yee) { + + EvolveFPMLCartesian <CartesianYeeAlgorithm> ( Ffield, Efield, dt ); + + } else if (m_fdtd_algo == MaxwellSolverAlgo::CKC) { + + EvolveFPMLCartesian <CartesianCKCAlgorithm> ( Ffield, Efield, dt ); + + } else { + amrex::Abort("Unknown algorithm"); + } +#endif +} + + +#ifndef WARPX_DIM_RZ + +template<typename T_Algo> +void FiniteDifferenceSolver::EvolveFPMLCartesian ( + amrex::MultiFab* Ffield, + std::array< amrex::MultiFab*, 3 > const Efield, + amrex::Real const dt ) { + + // Loop through the grids, and over the tiles within each grid +#ifdef _OPENMP +#pragma omp parallel if (amrex::Gpu::notInLaunchRegion()) +#endif + for ( MFIter mfi(*Ffield, TilingIfNotGPU()); mfi.isValid(); ++mfi ) { + + // Extract field data for this grid/tile + Array4<Real> const& F = Ffield->array(mfi); + Array4<Real> const& Ex = Efield[0]->array(mfi); + Array4<Real> const& Ey = Efield[1]->array(mfi); + Array4<Real> const& Ez = Efield[2]->array(mfi); + + // Extract stencil coefficients + Real const * const AMREX_RESTRICT coefs_x = m_stencil_coefs_x.dataPtr(); + int const n_coefs_x = m_stencil_coefs_x.size(); + Real const * const AMREX_RESTRICT coefs_y = m_stencil_coefs_y.dataPtr(); + int const n_coefs_y = m_stencil_coefs_y.size(); + Real const * const AMREX_RESTRICT coefs_z = m_stencil_coefs_z.dataPtr(); + int const n_coefs_z = m_stencil_coefs_z.size(); + + // Extract tileboxes for which to loop + Box const& tf = mfi.tilebox(Ffield->ixType().ixType()); + + // Loop over the cells and update the fields + amrex::ParallelFor(tf, + + [=] AMREX_GPU_DEVICE (int i, int j, int k){ + + F(i, j, k, PMLComp::x) += dt * ( + T_Algo::DownwardDx(Ex, coefs_x, n_coefs_x, i, j, k, PMLComp::xx) + + T_Algo::DownwardDx(Ex, coefs_x, n_coefs_x, i, j, k, PMLComp::xy) + + T_Algo::DownwardDx(Ex, coefs_x, n_coefs_x, i, j, k, PMLComp::xz) ); + + F(i, j, k, PMLComp::y) += dt * ( + T_Algo::DownwardDy(Ey, coefs_y, n_coefs_y, i, j, k, PMLComp::yx) + + T_Algo::DownwardDy(Ey, coefs_y, n_coefs_y, i, j, k, PMLComp::yy) + + T_Algo::DownwardDy(Ey, coefs_y, n_coefs_y, i, j, k, PMLComp::yz) ); + + F(i, j, k, PMLComp::z) += dt * ( + T_Algo::DownwardDz(Ez, coefs_z, n_coefs_z, i, j, k, PMLComp::zx) + + T_Algo::DownwardDz(Ez, coefs_z, n_coefs_z, i, j, k, PMLComp::zy) + + T_Algo::DownwardDz(Ez, coefs_z, n_coefs_z, i, j, k, PMLComp::zz) ); + + } + + ); + + } + +} + +#endif // corresponds to ifndef WARPX_DIM_RZ diff --git a/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceSolver.H b/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceSolver.H index b23e3bbed..e60ceed3e 100644 --- a/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceSolver.H +++ b/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceSolver.H @@ -9,6 +9,7 @@ #define WARPX_FINITE_DIFFERENCE_SOLVER_H_ #include <AMReX_MultiFab.H> +#include "BoundaryConditions/PML.H" /** * \brief Top-level class for the electromagnetic finite-difference solver @@ -53,6 +54,20 @@ class FiniteDifferenceSolver void ComputeDivE ( const std::array<std::unique_ptr<amrex::MultiFab>,3>& Efield, amrex::MultiFab& divE ); + void EvolveBPML ( std::array< amrex::MultiFab*, 3 > Bfield, + std::array< amrex::MultiFab*, 3 > const Efield, + amrex::Real const dt ); + + void EvolveEPML ( std::array< amrex::MultiFab*, 3 > Efield, + std::array< amrex::MultiFab*, 3 > const Bfield, + std::array< amrex::MultiFab*, 3 > const Jfield, + amrex::MultiFab* const Ffield, + MultiSigmaBox const& sigba, + amrex::Real const dt, bool pml_has_particles ); + + void EvolveFPML ( amrex::MultiFab* Ffield, + std::array< amrex::MultiFab*, 3 > const Efield, + amrex::Real const dt ); private: @@ -130,6 +145,26 @@ class FiniteDifferenceSolver const std::array<std::unique_ptr<amrex::MultiFab>,3>& Efield, amrex::MultiFab& divE ); + template< typename T_Algo > + void EvolveBPMLCartesian ( + std::array< amrex::MultiFab*, 3 > Bfield, + std::array< amrex::MultiFab*, 3 > const Efield, + amrex::Real const dt ); + + template< typename T_Algo > + void EvolveEPMLCartesian ( + std::array< amrex::MultiFab*, 3 > Efield, + std::array< amrex::MultiFab*, 3 > const Bfield, + std::array< amrex::MultiFab*, 3 > const Jfield, + amrex::MultiFab* const Ffield, + MultiSigmaBox const& sigba, + amrex::Real const dt, bool pml_has_particles ); + + template< typename T_Algo > + void EvolveFPMLCartesian ( amrex::MultiFab* Ffield, + std::array< amrex::MultiFab*, 3 > const Efield, + amrex::Real const dt ); + #endif }; diff --git a/Source/FieldSolver/FiniteDifferenceSolver/Make.package b/Source/FieldSolver/FiniteDifferenceSolver/Make.package index 9d95e04cc..b267463bd 100644 --- a/Source/FieldSolver/FiniteDifferenceSolver/Make.package +++ b/Source/FieldSolver/FiniteDifferenceSolver/Make.package @@ -3,5 +3,8 @@ CEXE_sources += EvolveB.cpp CEXE_sources += EvolveE.cpp CEXE_sources += EvolveF.cpp CEXE_sources += ComputeDivE.cpp +CEXE_sources += EvolveBPML.cpp +CEXE_sources += EvolveEPML.cpp +CEXE_sources += EvolveFPML.cpp VPATH_LOCATIONS += $(WARPX_HOME)/Source/FieldSolver/FiniteDifferenceSolver diff --git a/Source/FieldSolver/WarpXPushFieldsEM.cpp b/Source/FieldSolver/WarpXPushFieldsEM.cpp index f63102b34..67fb40c01 100644 --- a/Source/FieldSolver/WarpXPushFieldsEM.cpp +++ b/Source/FieldSolver/WarpXPushFieldsEM.cpp @@ -140,81 +140,24 @@ void WarpX::EvolveB (int lev, PatchType patch_type, amrex::Real a_dt) { + // Evolve B field in regular cells if (patch_type == PatchType::fine) { m_fdtd_solver_fp[lev]->EvolveB( Bfield_fp[lev], Efield_fp[lev], a_dt ); } else { m_fdtd_solver_cp[lev]->EvolveB( Bfield_cp[lev], Efield_cp[lev], a_dt ); } - const int patch_level = (patch_type == PatchType::fine) ? lev : lev-1; - const std::array<Real,3>& dx = WarpX::CellSize(patch_level); - const Real dtsdx = a_dt/dx[0], dtsdy = a_dt/dx[1], dtsdz = a_dt/dx[2]; - - if (do_pml && pml[lev]->ok()) - { - const auto& pml_B = (patch_type == PatchType::fine) ? pml[lev]->GetB_fp() : pml[lev]->GetB_cp(); - const auto& pml_E = (patch_type == PatchType::fine) ? pml[lev]->GetE_fp() : pml[lev]->GetE_cp(); - -#ifdef _OPENMP -#pragma omp parallel if (Gpu::notInLaunchRegion()) -#endif - for ( MFIter mfi(*pml_B[0], TilingIfNotGPU()); mfi.isValid(); ++mfi ) - { - const Box& tbx = mfi.tilebox( pml_B[0]->ixType().toIntVect() ); - const Box& tby = mfi.tilebox( pml_B[1]->ixType().toIntVect() ); - const Box& tbz = mfi.tilebox( pml_B[2]->ixType().toIntVect() ); - auto const& pml_Bxfab = pml_B[0]->array(mfi); - auto const& pml_Byfab = pml_B[1]->array(mfi); - auto const& pml_Bzfab = pml_B[2]->array(mfi); - auto const& pml_Exfab = pml_E[0]->array(mfi); - auto const& pml_Eyfab = pml_E[1]->array(mfi); - auto const& pml_Ezfab = pml_E[2]->array(mfi); - if (WarpX::maxwell_fdtd_solver_id == 0) { - amrex::ParallelFor(tbx, tby, tbz, - [=] AMREX_GPU_DEVICE (int i, int j, int k) { - warpx_push_pml_bx_yee(i,j,k,pml_Bxfab,pml_Eyfab,pml_Ezfab, - dtsdy,dtsdz); - }, - [=] AMREX_GPU_DEVICE (int i, int j, int k) { - warpx_push_pml_by_yee(i,j,k,pml_Byfab,pml_Exfab,pml_Ezfab, - dtsdx,dtsdz); - }, - [=] AMREX_GPU_DEVICE (int i, int j, int k) { - warpx_push_pml_bz_yee(i,j,k,pml_Bzfab,pml_Exfab,pml_Eyfab, - dtsdx,dtsdy); - }); - } else if (WarpX::maxwell_fdtd_solver_id == 1) { - Real betaxy, betaxz, betayx, betayz, betazx, betazy; - Real gammax, gammay, gammaz; - Real alphax, alphay, alphaz; - warpx_calculate_ckc_coefficients(dtsdx, dtsdy, dtsdz, - betaxy, betaxz, betayx, betayz, - betazx, betazy, gammax, gammay, - gammaz, alphax, alphay, alphaz); - - amrex::ParallelFor(tbx, tby, tbz, - [=] AMREX_GPU_DEVICE (int i, int j, int k) { - warpx_push_pml_bx_ckc(i,j,k,pml_Bxfab,pml_Eyfab,pml_Ezfab, - betaxy, betaxz, betayx, betayz, - betazx, betazy, gammax, gammay, - gammaz, alphax, alphay, alphaz); - }, - [=] AMREX_GPU_DEVICE (int i, int j, int k) { - warpx_push_pml_by_ckc(i,j,k,pml_Byfab,pml_Exfab,pml_Ezfab, - betaxy, betaxz, betayx, betayz, - betazx, betazy, gammax, gammay, - gammaz, alphax, alphay, alphaz); - }, - [=] AMREX_GPU_DEVICE (int i, int j, int k) { - warpx_push_pml_bz_ckc(i,j,k,pml_Bzfab,pml_Exfab,pml_Eyfab, - betaxy, betaxz, betayx, betayz, - betazx, betazy, gammax, gammay, - gammaz, alphax, alphay, alphaz); - }); - - } + // Evolve B field in PML cells + if (do_pml && pml[lev]->ok()) { + if (patch_type == PatchType::fine) { + m_fdtd_solver_fp[lev]->EvolveBPML( + pml[lev]->GetB_fp(), pml[lev]->GetE_fp(), a_dt ); + } else { + m_fdtd_solver_cp[lev]->EvolveBPML( + pml[lev]->GetB_cp(), pml[lev]->GetE_cp(), a_dt ); } } + } void @@ -240,7 +183,7 @@ WarpX::EvolveE (int lev, amrex::Real a_dt) void WarpX::EvolveE (int lev, PatchType patch_type, amrex::Real a_dt) { - + // Evolve E field in regular cells if (patch_type == PatchType::fine) { m_fdtd_solver_fp[lev]->EvolveE( Efield_fp[lev], Bfield_fp[lev], current_fp[lev], F_fp[lev], a_dt ); @@ -249,142 +192,20 @@ WarpX::EvolveE (int lev, PatchType patch_type, amrex::Real a_dt) current_cp[lev], F_cp[lev], a_dt ); } - const Real mu_c2_dt = (PhysConst::mu0*PhysConst::c*PhysConst::c) * a_dt; - const Real c2dt = (PhysConst::c*PhysConst::c) * a_dt; - - const int patch_level = (patch_type == PatchType::fine) ? lev : lev-1; - const std::array<Real,3>& dx = WarpX::CellSize(patch_level); - const Real dtsdx_c2 = c2dt/dx[0], dtsdy_c2 = c2dt/dx[1], dtsdz_c2 = c2dt/dx[2]; - - MultiFab* F; - if (patch_type == PatchType::fine) - { - F = F_fp[lev].get(); - } - else if (patch_type == PatchType::coarse) - { - F = F_cp[lev].get(); - } - - if (do_pml && pml[lev]->ok()) - { - const auto& pml_B = (patch_type == PatchType::fine) ? pml[lev]->GetB_fp() : pml[lev]->GetB_cp(); - const auto& pml_E = (patch_type == PatchType::fine) ? pml[lev]->GetE_fp() : pml[lev]->GetE_cp(); - const auto& pml_j = (patch_type == PatchType::fine) ? pml[lev]->Getj_fp() : pml[lev]->Getj_cp(); - const auto& pml_F = (patch_type == PatchType::fine) ? pml[lev]->GetF_fp() : pml[lev]->GetF_cp(); - const auto& sigba = (patch_type == PatchType::fine) ? pml[lev]->GetMultiSigmaBox_fp() - : pml[lev]->GetMultiSigmaBox_cp(); -#ifdef _OPENMP -#pragma omp parallel if (Gpu::notInLaunchRegion()) -#endif - for ( MFIter mfi(*pml_E[0], TilingIfNotGPU()); mfi.isValid(); ++mfi ) - { - const Box& tex = mfi.tilebox( pml_E[0]->ixType().toIntVect() ); - const Box& tey = mfi.tilebox( pml_E[1]->ixType().toIntVect() ); - const Box& tez = mfi.tilebox( pml_E[2]->ixType().toIntVect() ); - - auto const& pml_Exfab = pml_E[0]->array(mfi); - auto const& pml_Eyfab = pml_E[1]->array(mfi); - auto const& pml_Ezfab = pml_E[2]->array(mfi); - auto const& pml_Bxfab = pml_B[0]->array(mfi); - auto const& pml_Byfab = pml_B[1]->array(mfi); - auto const& pml_Bzfab = pml_B[2]->array(mfi); - - amrex::ParallelFor(tex, tey, tez, - [=] AMREX_GPU_DEVICE (int i, int j, int k) { - warpx_push_pml_ex_yee(i,j,k,pml_Exfab,pml_Byfab,pml_Bzfab, - dtsdy_c2,dtsdz_c2); - }, - [=] AMREX_GPU_DEVICE (int i, int j, int k) { - warpx_push_pml_ey_yee(i,j,k,pml_Eyfab,pml_Bxfab,pml_Bzfab, - dtsdx_c2,dtsdz_c2); - }, - [=] AMREX_GPU_DEVICE (int i, int j, int k) { - warpx_push_pml_ez_yee(i,j,k,pml_Ezfab,pml_Bxfab,pml_Byfab, - dtsdx_c2,dtsdy_c2); - }); - - if (pml_has_particles) { - // Update the E field in the PML, using the current - // deposited by the particles in the PML - auto const& pml_jxfab = pml_j[0]->array(mfi); - auto const& pml_jyfab = pml_j[1]->array(mfi); - auto const& pml_jzfab = pml_j[2]->array(mfi); - const Real* sigmaj_x = sigba[mfi].sigma[0].data(); - const Real* sigmaj_y = sigba[mfi].sigma[1].data(); - const Real* sigmaj_z = sigba[mfi].sigma[2].data(); - - int const x_lo = sigba[mfi].sigma[0].lo(); -#if (AMREX_SPACEDIM == 3) - int const y_lo = sigba[mfi].sigma[1].lo(); - int const z_lo = sigba[mfi].sigma[2].lo(); -#else - int const y_lo = 0; - int const z_lo = sigba[mfi].sigma[1].lo(); -#endif - amrex::ParallelFor( tex, tey, tez, - [=] AMREX_GPU_DEVICE (int i, int j, int k) { - push_ex_pml_current(i,j,k, - pml_Exfab, pml_jxfab, sigmaj_y, sigmaj_z, - y_lo, z_lo, mu_c2_dt); - }, - [=] AMREX_GPU_DEVICE (int i, int j, int k) { - push_ey_pml_current(i,j,k, - pml_Eyfab, pml_jyfab, sigmaj_x, sigmaj_z, - x_lo, z_lo, mu_c2_dt); - }, - [=] AMREX_GPU_DEVICE (int i, int j, int k) { - push_ez_pml_current(i,j,k, - pml_Ezfab, pml_jzfab, sigmaj_x, sigmaj_y, - x_lo, y_lo, mu_c2_dt); - } - ); - } - - - if (pml_F) - { - - auto const& pml_F_fab = pml_F->array(mfi); - - if (WarpX::maxwell_fdtd_solver_id == 0) { - - amrex::ParallelFor(tex, tey, tez, - [=] AMREX_GPU_DEVICE (int i, int j, int k) { - warpx_push_pml_ex_f_yee(i,j,k,pml_Exfab,pml_F_fab,dtsdx_c2); - }, - [=] AMREX_GPU_DEVICE (int i, int j, int k) { - warpx_push_pml_ey_f_yee(i,j,k,pml_Eyfab,pml_F_fab,dtsdy_c2); - }, - [=] AMREX_GPU_DEVICE (int i, int j, int k) { - warpx_push_pml_ez_f_yee(i,j,k,pml_Ezfab,pml_F_fab,dtsdz_c2); - }); - - } else if (WarpX::maxwell_fdtd_solver_id == 1) { - - Real betaxy, betaxz, betayx, betayz, betazx, betazy; - Real gammax, gammay, gammaz; - Real alphax, alphay, alphaz; - warpx_calculate_ckc_coefficients(dtsdx_c2, dtsdy_c2, dtsdz_c2, - betaxy, betaxz, betayx, betayz, - betazx, betazy, gammax, gammay, - gammaz, alphax, alphay, alphaz); - amrex::ParallelFor(tex, tey, tez, - [=] AMREX_GPU_DEVICE (int i, int j, int k) { - warpx_push_pml_ex_f_ckc(i,j,k,pml_Exfab,pml_F_fab, - alphax,betaxy,betaxz,gammax); - }, - [=] AMREX_GPU_DEVICE (int i, int j, int k) { - warpx_push_pml_ey_f_ckc(i,j,k,pml_Eyfab,pml_F_fab, - alphay,betayx,betayz,gammay); - }, - [=] AMREX_GPU_DEVICE (int i, int j, int k) { - warpx_push_pml_ez_f_ckc(i,j,k,pml_Ezfab,pml_F_fab, - alphaz,betazx,betazy,gammaz); - }); - - } - } + // Evolve E field in PML cells + if (do_pml && pml[lev]->ok()) { + if (patch_type == PatchType::fine) { + m_fdtd_solver_fp[lev]->EvolveEPML( + pml[lev]->GetE_fp(), pml[lev]->GetB_fp(), + pml[lev]->Getj_fp(), pml[lev]->GetF_fp(), + pml[lev]->GetMultiSigmaBox_fp(), + a_dt, pml_has_particles ); + } else { + m_fdtd_solver_cp[lev]->EvolveEPML( + pml[lev]->GetE_cp(), pml[lev]->GetB_cp(), + pml[lev]->Getj_cp(), pml[lev]->GetF_cp(), + pml[lev]->GetMultiSigmaBox_cp(), + a_dt, pml_has_particles ); } } } @@ -418,6 +239,7 @@ WarpX::EvolveF (int lev, PatchType patch_type, amrex::Real a_dt, DtType a_dt_typ const int rhocomp = (a_dt_type == DtType::FirstHalf) ? 0 : 1; + // Evolve F field in regular cells if (patch_type == PatchType::fine) { m_fdtd_solver_fp[lev]->EvolveF( F_fp[lev], Efield_fp[lev], rho_fp[lev], rhocomp, a_dt ); @@ -426,37 +248,17 @@ WarpX::EvolveF (int lev, PatchType patch_type, amrex::Real a_dt, DtType a_dt_typ rho_cp[lev], rhocomp, a_dt ); } - const int patch_level = (patch_type == PatchType::fine) ? lev : lev-1; - const auto& dx = WarpX::CellSize(patch_level); - const std::array<Real,3> dtsdx {a_dt/dx[0], a_dt/dx[1], a_dt/dx[2]}; - - if (do_pml && pml[lev]->ok()) - { - const auto& pml_F = (patch_type == PatchType::fine) ? pml[lev]->GetF_fp() : pml[lev]->GetF_cp(); - const auto& pml_E = (patch_type == PatchType::fine) ? pml[lev]->GetE_fp() : pml[lev]->GetE_cp(); - -#ifdef _OPENMP -#pragma omp parallel if (Gpu::notInLaunchRegion()) -#endif - for ( MFIter mfi(*pml_F, TilingIfNotGPU()); mfi.isValid(); ++mfi ) - { - const Box& bx = mfi.tilebox(); - - auto const& pml_F_fab = pml_F->array(mfi); - auto const& pml_Exfab = pml_E[0]->array(mfi); - auto const& pml_Eyfab = pml_E[1]->array(mfi); - auto const& pml_Ezfab = pml_E[2]->array(mfi); - - amrex::ParallelFor(bx, - [=] AMREX_GPU_DEVICE (int i, int j, int k) - { - warpx_push_pml_F(i, j, k, pml_F_fab, pml_Exfab, - pml_Eyfab, pml_Ezfab, - dtsdx[0], dtsdx[1], dtsdx[2]); - }); - + // Evolve F field in PML cells + if (do_pml && pml[lev]->ok()) { + if (patch_type == PatchType::fine) { + m_fdtd_solver_fp[lev]->EvolveFPML( + pml[lev]->GetF_fp(), pml[lev]->GetE_fp(), a_dt ); + } else { + m_fdtd_solver_cp[lev]->EvolveFPML( + pml[lev]->GetF_cp(), pml[lev]->GetE_cp(), a_dt ); } } + } #ifdef WARPX_DIM_RZ diff --git a/Source/FieldSolver/WarpX_FDTD.H b/Source/FieldSolver/WarpX_FDTD.H index 847323a0f..a49da2e08 100644 --- a/Source/FieldSolver/WarpX_FDTD.H +++ b/Source/FieldSolver/WarpX_FDTD.H @@ -9,64 +9,6 @@ #include <AMReX_FArrayBox.H> -AMREX_FORCE_INLINE -void warpx_calculate_ckc_coefficients(amrex::Real dtsdx, amrex::Real dtsdy, amrex::Real dtsdz, - amrex::Real &betaxy, amrex::Real &betaxz, amrex::Real &betayx, - amrex::Real &betayz, amrex::Real &betazx, amrex::Real &betazy, - amrex::Real &gammax, amrex::Real &gammay, amrex::Real &gammaz, - amrex::Real &alphax, amrex::Real &alphay, amrex::Real &alphaz) -{ - using namespace amrex; - - // Cole-Karkkainen-Cowan push - // computes coefficients according to Cowan - PRST-AB 16, 041303 (2013) -#if defined WARPX_DIM_3D - const Real delta = std::max( { dtsdx,dtsdy,dtsdz } ); - const Real rx = (dtsdx/delta)*(dtsdx/delta); - const Real ry = (dtsdy/delta)*(dtsdy/delta); - const Real rz = (dtsdz/delta)*(dtsdz/delta); - const Real beta = 0.125*(1. - rx*ry*rz/(ry*rz + rz*rx + rx*ry)); - betaxy = ry*beta; - betaxz = rz*beta; - betayx = rx*beta; - betayz = rz*beta; - betazx = rx*beta; - betazy = ry*beta; - gammax = ry*rz*(0.0625 - 0.125*ry*rz/(ry*rz + rz*rx + rx*ry)); - gammay = rx*rz*(0.0625 - 0.125*rx*rz/(ry*rz + rz*rx + rx*ry)); - gammaz = rx*ry*(0.0625 - 0.125*rx*ry/(ry*rz + rz*rx + rx*ry)); - alphax = 1. - 2.*betaxy - 2.*betaxz - 4.*gammax; - alphay = 1. - 2.*betayx - 2.*betayz - 4.*gammay; - alphaz = 1. - 2.*betazx - 2.*betazy - 4.*gammaz; - - betaxy *= dtsdx; - betaxz *= dtsdx; - betayx *= dtsdy; - betayz *= dtsdy; - betazx *= dtsdz; - betazy *= dtsdz; - alphax *= dtsdx; - alphay *= dtsdy; - alphaz *= dtsdz; - gammax *= dtsdx; - gammay *= dtsdy; - gammaz *= dtsdz; -#elif defined WARPX_DIM_XZ - const Real delta = std::max(dtsdx,dtsdz); - const Real rx = (dtsdx/delta)*(dtsdx/delta); - const Real rz = (dtsdz/delta)*(dtsdz/delta); - betaxz = 0.125*rz; - betazx = 0.125*rx; - alphax = 1. - 2.*betaxz; - alphaz = 1. - 2.*betazx; - - betaxz *= dtsdx; - betazx *= dtsdz; - alphax *= dtsdx; - alphaz *= dtsdz; -#endif -} - AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void warpx_computedivb(int i, int j, int k, int dcomp, amrex::Array4<amrex::Real> const& divB, diff --git a/Source/WarpX.cpp b/Source/WarpX.cpp index 28a72624b..82ccabb2d 100644 --- a/Source/WarpX.cpp +++ b/Source/WarpX.cpp @@ -282,17 +282,6 @@ WarpX::WarpX () // Sanity checks. Must be done after calling the MultiParticleContainer // constructor, as it reads additional parameters // (e.g., use_fdtd_nci_corr) - -#ifndef WARPX_USE_PSATD - AMREX_ALWAYS_ASSERT_WITH_MESSAGE( - not ( do_pml && do_nodal ), - "PML + do_nodal for finite-difference not implemented" - ); -#endif - AMREX_ALWAYS_ASSERT_WITH_MESSAGE( - not ( do_dive_cleaning && do_nodal ), - "divE cleaning + do_nodal not implemented" - ); #ifdef WARPX_USE_PSATD AMREX_ALWAYS_ASSERT(use_fdtd_nci_corr == 0); AMREX_ALWAYS_ASSERT(do_subcycling == 0); |