Arakawa finite difference scheme
gg = 1.0/(4.0*dx*dy)
hh = 1.0/3.0
for i = 2:nx+1 for j = 2:ny+1
j1 = gg*((w[i+1,j]-w[i-1,j])*(s[i,j+1]-s[i,j-1]) -
(w[i,j+1]-w[i,j-1])*(s[i+1,j]-s[i-1,j]))
j2 = gg*(w[i+1,j]*(s[i+1,j+1]-s[i+1,j-1]) -
w[i-1,j]*(s[i-1,j+1]-s[i-1,j-1]) -
w[i,j+1]*(s[i+1,j+1]-s[i-1,j+1]) +
w[i,j-1]*(s[i+1,j-1]-s[i-1,j-1]))
j3 = gg*(w[i+1,j+1]*(s[i,j+1]-s[i+1,j]) -
w[i-1,j-1]*(s[i-1,j]-s[i,j-1]) -
w[i-1,j+1]*(s[i,j+1]-s[i-1,j]) +
w[i+1,j-1]*(s[i+1,j]-s[i,j-1]))
data[i-1,j-1] = complex(-(j1+j2+j3)*hh, 0.0)
end~end
jf = fft(data) # convert Jacobian from physical to fourier space
return jf
end
✝ ✆
```
We validate our hybrid Arakawa-spectral solver for vortex-merger problem. Figure 27 shows the evolution of two point vortices merging into a single vortex computed using hybridized finite difference and spectral methodology. In the next section, we will eliminate finite difference and present a pseudo-spectral methodology.

**Figure 27.** Evolution of the vorticity field at three time steps during the co-rotating vortex interaction. The solution is computed using hybrid Arakawa-spectral solver.

#### **9. Pseudo-Spectral Solver**

In Section 8, we computed the Jacobian term using finite difference Arakawa scheme. We can also compute this term using a spectral method to develop the pseudo-spectral solver. The main motivation for the spectral solver is a unique advantage of both accuracy and efficiency considerations in simple geometries. It reduces a PDE problem to a set of ODEs without introducing a numerical discretization error. Furthermore, the code implementation will be much simpler due to the availability of FFT libraries. We refer to a paper by Mortensen and Langtangen [52] for further discussion of high performance pseudo-spectral solver implementation of 3D Navier-Stokes equations. We use the hybrid explicit and implicit Crank-Nicolson scheme for time integration similar to the previous section. The only difference is the computation of the Jacobian term. The Jacobian term involves the first derivative term and can be calculated easily using a spectral method. The first derivative of Equation (88) in *x* and *y* direction is calculated as

$$\frac{\partial}{\partial \mathbf{x}} \boldsymbol{u}\_{i,j} = \frac{1}{N\_{\mathbf{x}}} \frac{1}{N\_{\mathbf{y}}} \sum\_{m=-N\_{\mathbf{x}}/2}^{N\_{\mathbf{x}}/2-1} \sum\_{n=-N\_{\mathbf{y}}/2}^{N\_{\mathbf{y}}/2-1} (\iota m) \tilde{\boldsymbol{u}}\_{m,n} \boldsymbol{e}^{\prime 2\pi \left(\frac{\mathbf{m}\bar{\mathbf{y}}}{N\_{\mathbf{x}}} + \frac{\mathbf{n}^{j}}{N\_{\mathbf{y}}}\right)},\tag{155}$$

$$\frac{\partial}{\partial y}\mu\_{i,j} = \frac{1}{N\_x}\frac{1}{N\_y}\sum\_{m=-N\_x/2}^{N\_x/2-1}\sum\_{n=-N\_y/2}^{N\_y/2-1} \left(\mu n\right)\overline{u}\_{m,n}e^{i2\pi\left(\frac{mi}{N\_x} + \frac{nj}{N\_y}\right)}.\tag{156}$$

The Jacobian term in the Fourier space can be computed as

$$J\_{m,n} = (\iota m \vec{\psi}) \circledast (\iota m \vec{\omega}) - (\iota m \vec{\psi}) \circledast (\iota m \vec{\omega}),\tag{157}$$

where is the convolution operation. Using Equation (154), we get

$$J\_{m,n} = \left(\ln \frac{\tilde{\omega}}{k^2}\right) \circledast \left(\mu m \tilde{\omega}\right) - \left(\mu m \frac{\tilde{\omega}}{k^2}\right) \circledast \left(\mu n \tilde{\omega}\right). \tag{158}$$

The convolution operation is performed in physical space and then converted to physical space. The pseudo-code for calculating the Jacobian term is given as

$$J\_{m,n} = \text{fft}\underbrace{\text{ifft}\left(\mu\underbrace{\bar{\omega}}\_{k^2}\right)}\_{\frac{\text{\reflectbox{ $\omega$ }}{\text{\reflectbox{ $\omega$ }}}}} \otimes \underbrace{\text{ifft}\left(\mu m\bar{\omega}\right)}\_{\frac{\text{\reflectbox{ $\omega$ }}{\text{\reflectbox{ $\omega$ }}}}} - \underbrace{\text{ifft}\left(\mu\frac{\bar{\omega}}{k^2}\right)}\_{\frac{\text{\reflectbox{ $\omega$ }}{\text{\reflectbox{ $\omega$ }}}}} \otimes \underbrace{\text{ifft}\left(\mu m\bar{\omega}\right)}\_{\frac{\text{\reflectbox{ $\omega$ }}{\text{\reflectbox{ $\omega$ }}}}}\right).\tag{159}$$

*Fluids* **2019**, *4*, 159

We can observe that the Jacobian term involves the product of two functions. When we evaluate the product of two functions using Fourier transform, aliasing errors appear in the computation [12]. Aliasing errors refers to the misrepresentation of wavenumbers when the function is mapped onto a grid with not enough grid points. One way to reduce aliasing errors is to use 3/2-rule. In this method, the original grid is refined by a factor *M* = 3/2*N* in each direction. The additional grid points are assigned zero value and this is also referred to as zero-padding. For example, if we are using spatial resolution of *Nx* = 128 and *Ny* = 128 in *x* and *y* direction, the Jacobian term will be evaluated on grid with *Nx* = 192 and *Ny* = 192 (i.e., a data set consisted of 192<sup>2</sup> points is used in the inverse and forward FFT calls in Equation (159)). The main code for the pseudo-spectral solver is the same as the hybrid solver discussed in Section 8 except for the Jacobian function. The implementation of Jacobian calculation with 3/2 padding is given in Listing 34.

**Listing 34.** Julia implementation of computing the nonlinear Jacobian term in Fourier space with 3/2 padding rule. ✞ ☎

```