zero padding with 3/2 rule
j1f_p[1:Int64(nx/2),1:Int64(ny/2)] = j1f[1:Int64(nx/2),1:Int64(ny/2)]
j1f_p[Int64(qx):nxe,1:Int64(ny/2)] = j1f[Int64(nx/2+1):nx,1:Int64(ny/2)]
j1f_p[1:Int64(nx/2),Int64(qy):nye] = j1f[1:Int64(nx/2),Int64(ny/2+1):ny]
j1f_p[Int64(qx):nxe,Int64(qy):nye] = j1f[Int64(nx/2+1):nx,Int64(ny/2+1):ny]
j2f_p[1:Int64(nx/2),1:Int64(ny/2)] = j2f[1:Int64(nx/2),1:Int64(ny/2)]
j2f_p[Int64(qx):nxe,1:Int64(ny/2)] = j2f[Int64(nx/2+1):nx,1:Int64(ny/2)]
j2f_p[1:Int64(nx/2),Int64(qy):nye] = j2f[1:Int64(nx/2),Int64(ny/2+1):ny]
j2f_p[Int64(qx):nxe,Int64(qy):nye] = j2f[Int64(nx/2+1):nx,Int64(ny/2+1):ny]
j3f_p[1:Int64(nx/2),1:Int64(ny/2)] = j3f[1:Int64(nx/2),1:Int64(ny/2)]
j3f_p[Int64(qx):nxe,1:Int64(ny/2)] = j3f[Int64(nx/2+1):nx,1:Int64(ny/2)]
j3f_p[1:Int64(nx/2),Int64(qy):nye] = j3f[1:Int64(nx/2),Int64(ny/2+1):ny]
j3f_p[Int64(qx):nxe,Int64(qy):nye] = j3f[Int64(nx/2+1):nx,Int64(ny/2+1):ny]
j4f_p[1:Int64(nx/2),1:Int64(ny/2)] = j4f[1:Int64(nx/2),1:Int64(ny/2)]
j4f_p[Int64(qx):nxe,1:Int64(ny/2)] = j4f[Int64(nx/2+1):nx,1:Int64(ny/2)]
j4f_p[1:Int64(nx/2),Int64(qy):nye] = j4f[1:Int64(nx/2),Int64(ny/2+1):ny]
j4f_p[Int64(qx):nxe,Int64(qy):nye] = j4f[Int64(nx/2+1):nx,Int64(ny/2+1):ny]
j1f_p = j1f_p*(nxe*nye)/(nx*ny)
j2f_p = j2f_p*(nxe*nye)/(nx*ny)
j3f_p = j3f_p*(nxe*nye)/(nx*ny)
j4f_p = j4f_p*(nxe*nye)/(nx*ny)
j1 = real(ifft(j1f_p))
j2 = real(ifft(j2f_p))
j3 = real(ifft(j3f_p))
j4 = real(ifft(j4f_p))
jacp = zeros(Float64,nxe,nye)
for i = 1:nxe for j = 1:nye
jacp[i,j] = j1[i,j]*j2[i,j] - j3[i,j]*j4[i,j]
end~end
jacpf = fft(jacp)
jf = zeros(ComplexF64,nx,ny)
jf[1:Int64(nx/2),1:Int64(ny/2)] = jacpf[1:Int64(nx/2),1:Int64(ny/2)]
jf[Int64(nx/2+1):nx,1:Int64(ny/2)] = jacpf[Int64(qx):nxe,1:Int64(ny/2)]
jf[1:Int64(nx/2),Int64(ny/2+1):ny] = jacpf[1:Int64(nx/2),Int64(qy):nye]
jf[Int64(nx/2+1):nx,Int64(ny/2+1):ny] = jacpf[Int64(qx):nxe,Int64(qy):nye]
jf = jf*(nx*ny)/(nxe*nye)
return jf
```
#### **end** ✝ ✆

We compare the performance of pseudo-spectral solver for two-dimensional incompressible flow problems written in Julia with the vectorized version of Python code. The FFT operations in Julia are performed using the FFTW package. We use two different packages for FFT operations in Python: Numpy library and pyFFTW package. We measure the CPU time for two grid resolutions: *Nx* = *Ny* = 128 and *Nx* = *Ny* = 256. We test two codes at *Re* = 1000 with Δ*t* = 0.01 from time *t* = 0 to *t* = 20. The CPU time for two grid resolutions for two versions of code is provided in Table 4. It can be seen that the computational time for Julia is slightly better than Python code with the Numpy library for performing FFT operations. The pyFFTW package uses FFTW which is a very fast C library. The computational time required by vectorized Python code with pyFFTW package is better than the Julia code. This might be due to the effective interface between the FFTW library and pyFFTW package. Therefore, each programming language has some benefits over others. However, the goal of this study is not to optimize the code or find the best scientific programming language but is to present how can one implement CFD codes.

**Table 4.** Comparison of CPU time for Julia and Python version of pseudo-spectral solver at two grid resolutions. We highlight that the codes are not written for optimization purpose and the CPU time may not be the optimal time that can be achieved for each of the programming languages.


Another common technique to reduce aliasing error which is widely used in practice is to use 2/3-rule. In this technique we retain the data on 2/3 of the grid is retained, while the remaining data is assigned with zero value. For example, if we are using a spatial resolution of *Nx* = 128 and *Ny* = 128 (i.e., −64 ≤ *kx* ≤ 64 and −64 ≤ *ky* ≤ 64), we will retain the Fourier coefficients corresponding to (−42 to 42) and remaining Fourier coefficients are set zero. One of the advantages of using this method is that it is computationally efficient since the Fourier transform of the matrix having size as the power of 2 is more efficient than the matrix which has dimension non-powers of 2. If we use 2/3-rule it is equivalent to performing calculation on *Nx* = 85 and *Ny* = 85 grid. However, a data set consisted of 128<sup>2</sup> points is used in the inverse and forward FFT calls in Equation (159). The implementation of Jacobian calculation with 2/3 padding is given in Listing 35. Before we close this section, we would like to note that there is no significant aliasing error in the example presented in this section. Therefore, one might get similar results without performing any padding (1/1 rule).

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

```