diff options
Diffstat (limited to 'Source/Particles/interpolate_cic.F90')
| -rw-r--r-- | Source/Particles/interpolate_cic.F90 | 371 |
1 files changed, 371 insertions, 0 deletions
diff --git a/Source/Particles/interpolate_cic.F90 b/Source/Particles/interpolate_cic.F90 new file mode 100644 index 000000000..bc2c4f37e --- /dev/null +++ b/Source/Particles/interpolate_cic.F90 @@ -0,0 +1,371 @@ +module warpx_ES_interpolate_cic + + use iso_c_binding + use amrex_fort_module, only : amrex_real + + implicit none + +contains + +! This routine interpolates the electric field to the particle positions +! using cloud-in-cell interpolation. The electric fields are assumed to be +! node-centered. +! +! Arguments: +! particles : a pointer to the particle array-of-structs +! ns : the stride length of particle struct (the size of the struct in number of reals) +! np : the number of particles +! Ex_p : the electric field in the x-direction at the particle positions (output) +! Ey_p : the electric field in the y-direction at the particle positions (output) +! Ez_p : the electric field in the z-direction at the particle positions (output) +! Ex, Ey, Ez: Fabs conting the electric field on the mesh +! lo : a pointer to the lo corner of this valid box, in index space +! hi : a pointer to the hi corner of this valid box, in index space +! plo : the real position of the left-hand corner of the problem domain +! dx : the mesh spacing +! ng : the number of ghost cells for the E-field +! + subroutine warpx_interpolate_cic_3d(particles, ns, np, & + Ex_p, Ey_p, Ez_p, & + Ex, Ey, Ez, & + lo, hi, plo, dx, ng) & + bind(c,name='warpx_interpolate_cic_3d') + integer, value, intent(in) :: ns, np + real(amrex_real), intent(in) :: particles(ns,np) + real(amrex_real), intent(inout) :: Ex_p(np), Ey_p(np), Ez_p(np) + integer, intent(in) :: ng + integer, intent(in) :: lo(3) + integer, intent(in) :: hi(3) + real(amrex_real), intent(in) :: Ex(lo(1)-ng:hi(1)+ng, lo(2)-ng:hi(2)+ng, lo(3)-ng:hi(3)+ng) + real(amrex_real), intent(in) :: Ey(lo(1)-ng:hi(1)+ng, lo(2)-ng:hi(2)+ng, lo(3)-ng:hi(3)+ng) + real(amrex_real), intent(in) :: Ez(lo(1)-ng:hi(1)+ng, lo(2)-ng:hi(2)+ng, lo(3)-ng:hi(3)+ng) + real(amrex_real), intent(in) :: plo(3) + real(amrex_real), intent(in) :: dx(3) + + integer i, j, k, n + real(amrex_real) wx_lo, wy_lo, wz_lo, wx_hi, wy_hi, wz_hi + real(amrex_real) lx, ly, lz + real(amrex_real) inv_dx(3) + inv_dx = 1.0d0/dx + + do n = 1, np + lx = (particles(1, n) - plo(1))*inv_dx(1) + ly = (particles(2, n) - plo(2))*inv_dx(2) + lz = (particles(3, n) - plo(3))*inv_dx(3) + + i = floor(lx) + j = floor(ly) + k = floor(lz) + + wx_hi = lx - i + wy_hi = ly - j + wz_hi = lz - k + + wx_lo = 1.0d0 - wx_hi + wy_lo = 1.0d0 - wy_hi + wz_lo = 1.0d0 - wz_hi + + Ex_p(n) = wx_lo*wy_lo*wz_lo*Ex(i, j, k ) + & + wx_lo*wy_lo*wz_hi*Ex(i, j, k+1) + & + wx_lo*wy_hi*wz_lo*Ex(i, j+1, k ) + & + wx_lo*wy_hi*wz_hi*Ex(i, j+1, k+1) + & + wx_hi*wy_lo*wz_lo*Ex(i+1, j, k ) + & + wx_hi*wy_lo*wz_hi*Ex(i+1, j, k+1) + & + wx_hi*wy_hi*wz_lo*Ex(i+1, j+1, k ) + & + wx_hi*wy_hi*wz_hi*Ex(i+1, j+1, k+1) + + Ey_p(n) = wx_lo*wy_lo*wz_lo*Ey(i, j, k ) + & + wx_lo*wy_lo*wz_hi*Ey(i, j, k+1) + & + wx_lo*wy_hi*wz_lo*Ey(i, j+1, k ) + & + wx_lo*wy_hi*wz_hi*Ey(i, j+1, k+1) + & + wx_hi*wy_lo*wz_lo*Ey(i+1, j, k ) + & + wx_hi*wy_lo*wz_hi*Ey(i+1, j, k+1) + & + wx_hi*wy_hi*wz_lo*Ey(i+1, j+1, k ) + & + wx_hi*wy_hi*wz_hi*Ey(i+1, j+1, k+1) + + Ez_p(n) = wx_lo*wy_lo*wz_lo*Ez(i, j, k ) + & + wx_lo*wy_lo*wz_hi*Ez(i, j, k+1) + & + wx_lo*wy_hi*wz_lo*Ez(i, j+1, k ) + & + wx_lo*wy_hi*wz_hi*Ez(i, j+1, k+1) + & + wx_hi*wy_lo*wz_lo*Ez(i+1, j, k ) + & + wx_hi*wy_lo*wz_hi*Ez(i+1, j, k+1) + & + wx_hi*wy_hi*wz_lo*Ez(i+1, j+1, k ) + & + wx_hi*wy_hi*wz_hi*Ez(i+1, j+1, k+1) + + end do + + end subroutine warpx_interpolate_cic_3d + + + subroutine warpx_interpolate_cic_2d(particles, ns, np, & + Ex_p, Ey_p, & + Ex, Ey, & + lo, hi, plo, dx, ng) & + bind(c,name='warpx_interpolate_cic_2d') + integer, value, intent(in) :: ns, np + real(amrex_real), intent(in) :: particles(ns,np) + real(amrex_real), intent(inout) :: Ex_p(np), Ey_p(np) + integer, intent(in) :: ng + integer, intent(in) :: lo(2) + integer, intent(in) :: hi(2) + real(amrex_real), intent(in) :: Ex(lo(1)-ng:hi(1)+ng, lo(2)-ng:hi(2)+ng) + real(amrex_real), intent(in) :: Ey(lo(1)-ng:hi(1)+ng, lo(2)-ng:hi(2)+ng) + real(amrex_real), intent(in) :: plo(2) + real(amrex_real), intent(in) :: dx(2) + + integer i, j, n + real(amrex_real) wx_lo, wy_lo, wx_hi, wy_hi + real(amrex_real) lx, ly + real(amrex_real) inv_dx(2) + inv_dx = 1.0d0/dx + + do n = 1, np + lx = (particles(1, n) - plo(1))*inv_dx(1) + ly = (particles(2, n) - plo(2))*inv_dx(2) + + i = floor(lx) + j = floor(ly) + + wx_hi = lx - i + wy_hi = ly - j + + wx_lo = 1.0d0 - wx_hi + wy_lo = 1.0d0 - wy_hi + + Ex_p(n) = wx_lo*wy_lo*Ex(i, j ) + & + wx_lo*wy_hi*Ex(i, j+1) + & + wx_hi*wy_lo*Ex(i+1, j ) + & + wx_hi*wy_hi*Ex(i+1, j+1) + + Ey_p(n) = wx_lo*wy_lo*Ey(i, j ) + & + wx_lo*wy_hi*Ey(i, j+1) + & + wx_hi*wy_lo*Ey(i+1, j ) + & + wx_hi*wy_hi*Ey(i+1, j+1) + + end do + + end subroutine warpx_interpolate_cic_2d + + + subroutine warpx_interpolate_cic_two_levels_3d(particles, ns, np, & + Ex_p, Ey_p, Ez_p, & + Ex, Ey, Ez, & + lo, hi, dx, & + cEx, cEy, cEz, & + mask, & + clo, chi, cdx, & + plo, ng, lev) & + bind(c,name='warpx_interpolate_cic_two_levels_3d') + integer, value, intent(in) :: ns, np + real(amrex_real), intent(in) :: particles(ns,np) + real(amrex_real), intent(inout) :: Ex_p(np), Ey_p(np), Ez_p(np) + integer, intent(in) :: ng, lev + integer, intent(in) :: lo(3), hi(3) + integer, intent(in) :: clo(3), chi(3) + real(amrex_real), intent(in) :: Ex(lo(1)-ng:hi(1)+ng, lo(2)-ng:hi(2)+ng, lo(3)-ng:hi(3)+ng) + real(amrex_real), intent(in) :: Ey(lo(1)-ng:hi(1)+ng, lo(2)-ng:hi(2)+ng, lo(3)-ng:hi(3)+ng) + real(amrex_real), intent(in) :: Ez(lo(1)-ng:hi(1)+ng, lo(2)-ng:hi(2)+ng, lo(3)-ng:hi(3)+ng) + real(amrex_real), intent(in) :: cEx(clo(1)-ng:chi(1)+ng, clo(2)-ng:chi(2)+ng, clo(3)-ng:chi(3)+ng) + real(amrex_real), intent(in) :: cEy(clo(1)-ng:chi(1)+ng, clo(2)-ng:chi(2)+ng, clo(3)-ng:chi(3)+ng) + real(amrex_real), intent(in) :: cEz(clo(1)-ng:chi(1)+ng, clo(2)-ng:chi(2)+ng, clo(3)-ng:chi(3)+ng) + integer(c_int), intent(in) :: mask (lo(1):hi(1),lo(2):hi(2),lo(3):hi(3)) + real(amrex_real), intent(in) :: plo(3) + real(amrex_real), intent(in) :: dx(3), cdx(3) + + integer i, j, k, n + real(amrex_real) wx_lo, wy_lo, wz_lo, wx_hi, wy_hi, wz_hi + real(amrex_real) lx, ly, lz + real(amrex_real) inv_dx(3), inv_cdx(3) + inv_dx = 1.0d0/dx + inv_cdx = 1.0d0/cdx + + do n = 1, np + + lx = (particles(1, n) - plo(1))*inv_dx(1) + ly = (particles(2, n) - plo(2))*inv_dx(2) + lz = (particles(3, n) - plo(3))*inv_dx(3) + + i = floor(lx) + j = floor(ly) + k = floor(lz) + +! use the coarse E if near the level boundary + if (lev .eq. 1 .and. mask(i,j,k) .eq. 1) then + + lx = (particles(1, n) - plo(1))*inv_cdx(1) + ly = (particles(2, n) - plo(2))*inv_cdx(2) + lz = (particles(3, n) - plo(3))*inv_cdx(3) + + i = floor(lx) + j = floor(ly) + k = floor(lz) + + wx_hi = lx - i + wy_hi = ly - j + wz_hi = lz - k + + wx_lo = 1.0d0 - wx_hi + wy_lo = 1.0d0 - wy_hi + wz_lo = 1.0d0 - wz_hi + + Ex_p(n) = wx_lo*wy_lo*wz_lo*cEx(i, j, k ) + & + wx_lo*wy_lo*wz_hi*cEx(i, j, k+1) + & + wx_lo*wy_hi*wz_lo*cEx(i, j+1, k ) + & + wx_lo*wy_hi*wz_hi*cEx(i, j+1, k+1) + & + wx_hi*wy_lo*wz_lo*cEx(i+1, j, k ) + & + wx_hi*wy_lo*wz_hi*cEx(i+1, j, k+1) + & + wx_hi*wy_hi*wz_lo*cEx(i+1, j+1, k ) + & + wx_hi*wy_hi*wz_hi*cEx(i+1, j+1, k+1) + + Ey_p(n) = wx_lo*wy_lo*wz_lo*cEy(i, j, k ) + & + wx_lo*wy_lo*wz_hi*cEy(i, j, k+1) + & + wx_lo*wy_hi*wz_lo*cEy(i, j+1, k ) + & + wx_lo*wy_hi*wz_hi*cEy(i, j+1, k+1) + & + wx_hi*wy_lo*wz_lo*cEy(i+1, j, k ) + & + wx_hi*wy_lo*wz_hi*cEy(i+1, j, k+1) + & + wx_hi*wy_hi*wz_lo*cEy(i+1, j+1, k ) + & + wx_hi*wy_hi*wz_hi*cEy(i+1, j+1, k+1) + + Ez_p(n) = wx_lo*wy_lo*wz_lo*cEz(i, j, k ) + & + wx_lo*wy_lo*wz_hi*cEz(i, j, k+1) + & + wx_lo*wy_hi*wz_lo*cEz(i, j+1, k ) + & + wx_lo*wy_hi*wz_hi*cEz(i, j+1, k+1) + & + wx_hi*wy_lo*wz_lo*cEz(i+1, j, k ) + & + wx_hi*wy_lo*wz_hi*cEz(i+1, j, k+1) + & + wx_hi*wy_hi*wz_lo*cEz(i+1, j+1, k ) + & + wx_hi*wy_hi*wz_hi*cEz(i+1, j+1, k+1) + +! otherwise use the fine + else + + wx_hi = lx - i + wy_hi = ly - j + wz_hi = lz - k + + wx_lo = 1.0d0 - wx_hi + wy_lo = 1.0d0 - wy_hi + wz_lo = 1.0d0 - wz_hi + + Ex_p(n) = wx_lo*wy_lo*wz_lo*Ex(i, j, k ) + & + wx_lo*wy_lo*wz_hi*Ex(i, j, k+1) + & + wx_lo*wy_hi*wz_lo*Ex(i, j+1, k ) + & + wx_lo*wy_hi*wz_hi*Ex(i, j+1, k+1) + & + wx_hi*wy_lo*wz_lo*Ex(i+1, j, k ) + & + wx_hi*wy_lo*wz_hi*Ex(i+1, j, k+1) + & + wx_hi*wy_hi*wz_lo*Ex(i+1, j+1, k ) + & + wx_hi*wy_hi*wz_hi*Ex(i+1, j+1, k+1) + + Ey_p(n) = wx_lo*wy_lo*wz_lo*Ey(i, j, k ) + & + wx_lo*wy_lo*wz_hi*Ey(i, j, k+1) + & + wx_lo*wy_hi*wz_lo*Ey(i, j+1, k ) + & + wx_lo*wy_hi*wz_hi*Ey(i, j+1, k+1) + & + wx_hi*wy_lo*wz_lo*Ey(i+1, j, k ) + & + wx_hi*wy_lo*wz_hi*Ey(i+1, j, k+1) + & + wx_hi*wy_hi*wz_lo*Ey(i+1, j+1, k ) + & + wx_hi*wy_hi*wz_hi*Ey(i+1, j+1, k+1) + + Ez_p(n) = wx_lo*wy_lo*wz_lo*Ez(i, j, k ) + & + wx_lo*wy_lo*wz_hi*Ez(i, j, k+1) + & + wx_lo*wy_hi*wz_lo*Ez(i, j+1, k ) + & + wx_lo*wy_hi*wz_hi*Ez(i, j+1, k+1) + & + wx_hi*wy_lo*wz_lo*Ez(i+1, j, k ) + & + wx_hi*wy_lo*wz_hi*Ez(i+1, j, k+1) + & + wx_hi*wy_hi*wz_lo*Ez(i+1, j+1, k ) + & + wx_hi*wy_hi*wz_hi*Ez(i+1, j+1, k+1) + + end if + + end do + + end subroutine warpx_interpolate_cic_two_levels_3d + + + subroutine warpx_interpolate_cic_two_levels_2d(particles, ns, np, & + Ex_p, Ey_p, & + Ex, Ey, & + lo, hi, dx, & + cEx, cEy, & + mask, & + clo, chi, cdx, & + plo, ng, lev) & + bind(c,name='warpx_interpolate_cic_two_levels_2d') + integer, value, intent(in) :: ns, np + real(amrex_real), intent(in) :: particles(ns,np) + real(amrex_real), intent(inout) :: Ex_p(np), Ey_p(np) + integer, intent(in) :: ng, lev + integer, intent(in) :: lo(2), hi(2) + integer, intent(in) :: clo(2), chi(2) + real(amrex_real), intent(in) :: Ex(lo(1)-ng:hi(1)+ng, lo(2)-ng:hi(2)+ng) + real(amrex_real), intent(in) :: Ey(lo(1)-ng:hi(1)+ng, lo(2)-ng:hi(2)+ng) + real(amrex_real), intent(in) :: cEx(clo(1)-ng:chi(1)+ng, clo(2)-ng:chi(2)+ng) + real(amrex_real), intent(in) :: cEy(clo(1)-ng:chi(1)+ng, clo(2)-ng:chi(2)+ng) + integer(c_int), intent(in) :: mask (lo(1):hi(1),lo(2):hi(2)) + real(amrex_real), intent(in) :: plo(2) + real(amrex_real), intent(in) :: dx(2), cdx(2) + + integer i, j, n + real(amrex_real) wx_lo, wy_lo, wx_hi, wy_hi + real(amrex_real) lx, ly + real(amrex_real) inv_dx(2), inv_cdx(2) + inv_dx = 1.0d0/dx + inv_cdx = 1.0d0/cdx + + do n = 1, np + + lx = (particles(1, n) - plo(1))*inv_dx(1) + ly = (particles(2, n) - plo(2))*inv_dx(2) + + i = floor(lx) + j = floor(ly) + +! use the coarse E if near the level boundary + if (lev .eq. 1 .and. mask(i,j) .eq. 1) then + + lx = (particles(1, n) - plo(1))*inv_cdx(1) + ly = (particles(2, n) - plo(2))*inv_cdx(2) + + i = floor(lx) + j = floor(ly) + + wx_hi = lx - i + wy_hi = ly - j + + wx_lo = 1.0d0 - wx_hi + wy_lo = 1.0d0 - wy_hi + + Ex_p(n) = wx_lo*wy_lo*cEx(i, j ) + & + wx_lo*wy_hi*cEx(i, j+1) + & + wx_hi*wy_lo*cEx(i+1, j ) + & + wx_hi*wy_hi*cEx(i+1, j+1) + + Ey_p(n) = wx_lo*wy_lo*cEy(i, j ) + & + wx_lo*wy_hi*cEy(i, j+1) + & + wx_hi*wy_lo*cEy(i+1, j ) + & + wx_hi*wy_hi*cEy(i+1, j+1) + +! otherwise use the fine + else + + wx_hi = lx - i + wy_hi = ly - j + + wx_lo = 1.0d0 - wx_hi + wy_lo = 1.0d0 - wy_hi + + Ex_p(n) = wx_lo*wy_lo*Ex(i, j ) + & + wx_lo*wy_hi*Ex(i, j+1) + & + wx_hi*wy_lo*Ex(i+1, j ) + & + wx_hi*wy_hi*Ex(i+1, j+1) + + Ey_p(n) = wx_lo*wy_lo*Ey(i, j ) + & + wx_lo*wy_hi*Ey(i, j+1) + & + wx_hi*wy_lo*Ey(i+1, j ) + & + wx_hi*wy_hi*Ey(i+1, j+1) + + end if + + end do + + end subroutine warpx_interpolate_cic_two_levels_2d + +end module warpx_ES_interpolate_cic |