Fast Fourier Transform periodic interpolation method for superposition sums in a periodic unit cell

TL;DR

This work addresses the challenge of efficiently evaluating periodic sums that arise from non-uniform sources in infinite periodic unit cells. It introduces the FFT-PIM approach, which splits the periodic Green's function into a near-zone part with few images and a far-zone part computed on a sparse grid with interpolation, while the near-zone contribution is handled via an FFT-based method with a corrective step. The method achieves $O(N\log N)$ time and $O(N)$ memory and applies to 1D, 2D, and 3D periodicities for static and dynamic problems, with or without phase shifts, including the neutral-source (NPSP) special case. Numerical results on CPU and GPU validate high accuracy and substantial speedups, and the authors release an open-source PUFF package for the far-zone computations, enabling scalable evaluation of periodic sums in electromagnetics, acoustics, micromagnetics, and density functional theory contexts.

Abstract

We propose a Fast Fourier Transform based Periodic Interpolation Method (FFT-PIM), a flexible and computationally efficient approach for computing the scalar potential given by a superposition sum in a unit cell of an infinitely periodic array. Under the same umbrella, FFT-PIM allows computing the potential for 1D, 2D, and 3D periodicities for dynamic and static problems, including problems with and without a periodic phase shift. The computational complexity of the FFT-PIM is of $O(N \log N)$ for $N$ spatially coinciding sources and observer points. The FFT-PIM uses rapidly converging series representations of the Green's function serving as a kernel in the superposition sum. Based on these representations, the FFT-PIM splits the potential into its near-zone component, which includes a small number of images surrounding the unit cell of interest, and far-zone component, which includes the rest of an infinite number of images. The far-zone component is evaluated by projecting the non-uniform sources onto a sparse uniform grid, performing superposition sums on this sparse grid, and interpolating the potential from the uniform grid to the non-uniform observation points. The near-zone component is evaluated using an FFT-based method, which is adapted to efficiently handle non-uniform source-observer distributions within the periodic unit cell. The FFT-PIM can be used for a broad range of applications, such as periodic problems involving integral equations in computational electromagnetic and acoustic, micromagnetic solvers, and density functional theory solvers.

Figures (10)

• Figure 1: Illustration for a 2D periodic problem consisting of an infinite 2D array of cubes. The central red cube is the zeroth unit cell. Surrounding it along the $x$ and $y$ axes are the $1^{st}$ order image cubes in green and the $2^{nd}$ order image cubes in blue.
• Figure 2: Convergence of PGFs, $L_x,L_y,L_z=1$, $x=y=z=0.5$. (a) The relative error of the sum of the first $m$ terms from Eq. \ref{['eq:pgf']}, with $k_0=-1-j$ and $k_{x0}=1-j, k_{y0}=1+j, k_{z0}=-1+j$, (b) The relative error of the sum of the first $m$ terms from Eq. \ref{['eq:lgf']}. Exponential convergence is achieved in both cases.
• Figure 3: Illustration of the source (black) and observer (red) grids for the far-zone PSP component. Black and red circles are source/observer grid points on the Cartesian lattices. Black and red dots are source/observer points with arbitrary distribution. Black arrows denote projections from source points to source grid points. Green arrow is the grid-to-grid interaction, performed by convolution. Red arrows represent interpolations from observer grid points to observer points.
• Figure 4: Illustration for the error correction step 4 for the near-zone PSP component. The black dot is a source point and red dot is an observer point. The green region is the error-correction region $\Omega^{ER}$ . The black circles are grid points. For the observers in $\Omega^{ER}$, PSPs are calculated through direct calculations, this is done by subtracting the grid interaction inside the green region (dotted arrows) and adding the direct interactions (dashed purple arrow). The shaded green region and grey dot is the image of the right boundary in the $i_d=-1$ image cell, which is close to the left boundary of the unit cell within $\Omega^{ER}$.
• Figure 5: Potential on the $x$ axis of a coaxial structure . (a) Illustration of coaxial structure unit cell ($L=1$) and potential of the non-periodic unit cell. Inner radius $r_1=1$ with negative line charge density $\rho_1=-1$, outer radius $r_2=2$ with positive line charge density $\rho_2=1/2$. (b) PSP with a 1D periodicity along the $x$ axis with $L_x$ = 1. The blue line is the magnitude of PSP for the NPSP case with $k_0=k_{x0}=0$. The red, green and black lines are real part, imaginary part, and magnitude of PSP for the dynamic case with $k_0=1$, $k_{x0}=1-j$.