3.1.1 FFT Discrete convolution method

This method is mostly used as soon as a spatial discretization based on Finite Differences (FD) [52] can be defined over a structured mesh. This occurs for example, if one considers a magnetic body with rectangular geometry and subdivides it into a collection of square rectangular prisms with edges $d_x,d_y,d_z$ parallel to the coordinate axes. To start our discussion we recall the fact that the solution of magnetostatic problem (3.1)-(3.2) can be obtained in terms of the scalar potential $\varphi$ such that ${\mathbf{H}_\text{m}}=-\nabla\varphi$, solution of the following boundary value problem:

$$\left\{\begin{aligned} &\nabla^2\varphi=\nabla\cdot\textbf{M}&... ... \Omega}=-\textbf{n}\cdot\textbf{M} \end{aligned} \right. \quad.$$

The boundary value problem (3.3) admits the following solution [13]:

$$\varphi(\textbf{r})=-\frac{1}{4\pi}\int_\Omega \frac{\nabla'\cdo... ...)\cdot\textbf{n}}{\vert\textbf{r}-\textbf{r}'\vert} dS_{\textbf{r}'} \quad.$$ (3.4)

By using the divergence theorem, the surface integral can be rewritten as volume integral. Equation (3.4) can be put in the compact form:

$$\varphi(\textbf{r})=\frac{1}{4\pi}\int_\Omega \left[ -\frac{\nab... ...')}{\vert\textbf{r}-\textbf{r}'\vert}\right) \right] dV_{\textbf{r}'} \quad.$$ (3.5)

By using the fact that $\nabla\cdot(f\textbf{v})=f\nabla\cdot\textbf{v}+\textbf{v}\cdot\nabla f$, one ends up with:

$$\varphi(\textbf{r})=\frac{1}{4\pi}\int_\Omega \nabla'\left( \fra... ...textbf{r}'\vert}\right)\cdot\textbf{M}(\textbf{r}') dV_{\textbf{r}'} \quad.$$ (3.6)

Thus, the magnetostatic field ${\mathbf{H}_\text{m}}$ is given by

$${\mathbf{H}_\text{m}}(\textbf{r}) =-\nabla\varphi=-\frac{1}{4\pi}\nabla\int_\Omega \nabla'\left( \... ...f{r}-\textbf{r}'\vert}\right)\cdot\textbf{M}(\textbf{r}') dV_{\textbf{r}'}=$$

$$=-\int_\Omega \mathcal{N}(\textbf{r}-\textbf{r}')\cdot\textbf{M}(\textbf{r}') dV_{\textbf{r}'} \quad,$$ (3.7)

where $\mathcal{N}(\textbf{r}-\textbf{r}')$ is the demagnetizing tensor. The product $-\mathcal{N}(\textbf{r}-\textbf{r}')\cdot\textbf{M}(\textbf{r}') dV_{\textbf{r}'}$ gives the magnetostatic field produced at location $\textbf{r}$ by an elementary magnetic moment situated at location $\textbf{r}'$. The expression (3.7) remains formally unchanged if one assumes suitable discretization over $N$ elementary cells. For instance, if we subdivide the magnetic body into $N=n_x n_y n_z$ rectangular square prisms with edges parallel to the coordinate axes, with $n_x,n_y,n_z$ cells along the $x,y,z$ axis respectively, each cell can be uniquely determined by means of three indexes $i,j,k$. As far as magnetization within the cells is concerned, there are two approaches proposed in literature. The first is often referred to as constant volume charges^3.1 method, that is, $\nabla\cdot\textbf{M}$ is supposed to be constant within each cell. The second approach supposes the magnetization $\textbf{M}$ uniform within each cell. A comparison between these two approach has been performed in Ref. [53] with respect to the solution of the standard problem no. 2 (see $\mu-$mag website [79] for details). Thus, apart from this particular choice, one can rewrite Eq. (3.7) as a discrete convolution:

$${\mathbf{H}_\text{m}}_{;i,j,k}=-\sum_{i,j,k\neq i',j',k'} N_{i-i',j-j',k-k'} \cdot \textbf{M}_{i',j',k'} d_x d_y d_z \quad,$$ (3.8)

where $N_{p,q,r}$ is the demagnetizing tensor associated to the prism cell $(p,q,r)$:

$$N_{p,q,r}=\left(\begin{array}{ccc} N_{xx} & N_{xy} &N_{xz} N_{yx} &N_{yy} &N_{yz} N_{zx} &N_{zy} &N_{zz} \end{array} \right) \quad.$$ (3.9)

The discrete convolution (3.8) can be computed by means of the Discrete Fourier Transform (DFT), which can be implemented very efficiently with the well-known algorithm referred to as Fast Fourier Transform [60]. In fact, the time-domain convolution can be changed into a scalar product in frequency space using the Fourier transform. To properly take care of the finite size effect of the system and, therefore, to avoid circular convolution, the standard zero-padding techniques [64] can be used. In fact, the inverse Fourier transform will yield the correct field in real space, if the number of cells in each dimension after zero-padding is not less than twice the number of physical cells. The latter requirement ensures that the inferred periodic boundary condition of this enlarged region with padded zeros will not affect the physical data in real space after the inverse FFT is performed. In fact, to perform the FFT which does require overall periodicity, and yet not to allow the cells in the simulated region to be affected by the fields in the extended periods, the void buffer area between a physical cell and the first image cell in the adjacent period must exceed the interaction force range given by the number of the cells $N$ [65]. Assuming that the dimensions of the zero-padded discretization grid are $2n_x,2n_y,2n_z$, along the directions $x,y,z$ respectively, and referring for instance to the $x$ component of ${\mathbf{H}_\text{m}}$, the discrete convolution (3.8) can be written as:

$$\underline{\text{H}}_{x;i,j,k}=-\sum_{i,j,k\neq i',j',k'} \sum_{... ...eta;i-i',j-j',k-k'} \underline{\text{M}}_{\eta;i',j',k'} d_x d_y d_z \quad,$$ (3.10)

where $\underline{\text{H}}_x, \underline{\text{M}}_\eta, \underline{\text{N}}_{x\eta}$ ( $\eta\in\{x,y,z\}$) are $2n_x\times 2n_y \times 2n_z$ matrices. The Discrete Fourier Transform (DFT) $\hat{\underline{\text{H}}}_x$ of $\underline{\text{H}}_x$ can be expressed as:

$$\hat{\underline{\text{H}}}_x(\omega_x,\omega_y,\omega_z)=\sum_{i=... ...{2n_x}+\frac{\omega_y j}{2n_y} + \frac{\omega_z k}{2n_z} \right)\right] \quad,$$ (3.11)

where $\iota$ is the imaginary unit $\iota=\sqrt{-1}$ and $\omega_x, \omega_y, \omega_z$ are the frequency domain variables. Analogous expressions can be written for $\hat{\underline{\text{M}}}_x$ and $\hat{\underline{\text{N}}}_{x\eta}$. Therefore, the discrete convolution theorem states that Eq. (3.10) can be written in frequency domain as sum of matrix element-by-element product:

$$\hat{\underline{\text{H}}}_x(\omega_x,\omega_y,\omega_z)= \sum_{... ...y,\omega_z) \hat{\underline{\text{M}}}_\eta(\omega_x,\omega_y,\omega_z) \quad.$$ (3.12)

The $x$ component $\underline{\text{H}}_x$ can be obtained by computing the inverse FFT of the expression (3.12). Therefore, the cost of the evaluation of the demagnetizing field can be summarized as follows:

• In the preprocessing stage of the simulation, the FFTs of the 9 demagnetizing matrices $\underline{\text{N}}_{x\eta},\underline{\text{N}}_{y\eta},\underline{\text{N}}_{z\eta}$ ( $\eta\in\{x,y,z\}$) has to be evaluated and stored in memory.

  Figure: Computational cost scaling plot. Dashed line represents $\mathcal{O}(N^2)$ scaling. Dash-dotted line represents $\mathcal{O}(N)$ scaling. Solid line represents $\mathcal{O}(N\log N)$ scaling.

  
• For each computation of the demagnetizing field, six FFTs has to be computed: three related to magnetization matrices, namely $\hat{\underline{\text{M}}}_x, \hat{\underline{\text{M}}}_y, \hat{\underline{\text{M}}}_z$ and three inverse FFTs of the components $\hat{\underline{\text{H}}}_x, \hat{\underline{\text{H}}}_y, \hat{\underline{\text{H}}}_z$ in frequency domain.

The cost of a single FFT evaluation scales according to the $\mathcal{O}(N\log N)$ behavior. Thus, by using the FFT method the computational cost for the calculation of the magnetostatic field can be considerably reduced with respect to the $\mathcal{O}(N^2)$ scaling connected with the direct evaluation of the integral (3.7) (see Fig. 3.1).