Arakawa scheme
j1 = gg*((w[2:nx+3,1:ny+2]-w[0:nx+1,1:ny+2]) \
*( s[1:nx+2,2:ny+3]-s[1:nx+2,0:ny+1]) \
*( s[2:nx+3,1:ny+2]-s[0:nx+1,1:ny+2]))
j2 = gg*( w[2:nx+3,1:ny+2]*(s[2:nx+3,2:ny+3]-s[2:nx+3,0:ny+1]) \
+ w[1:nx+2,0:ny+1]*(s[2:nx+3,0:ny+1]-s[0:nx+1,0:ny+1]))
j3 = gg*( w[2:nx+3,2:ny+3]*(s[1:nx+2,2:ny+3]-s[2:nx+3,1:ny+2]) \
+ w[2:nx+3,0:ny+1]*(s[2:nx+3,1:ny+2]-s[1:nx+2,0:ny+1]) )
jac = (j1+j2+j3)*hh
lap = aa*(w[2:nx+3,1:ny+2]-2.0*w[1:nx+2,1:ny+2]+w[0:nx+1,1:ny+2]) \
+ bb*(w[1:nx+2,2:ny+3]-2.0*w[1:nx+2,1:ny+2]+w[1:nx+2,0:ny+1])
f[1:nx+2,1:ny+2] = -jac + lap/re
return f
```
#### **8. Hybrid Arakawa-Spectral Solver**

In this section, instead of using a fully explicit method in time for solving 2D incompressible Navier-Stokes equation, we show how to design hybrid explicit and implicit scheme. The nonlinear Jacobian term in vorticity Equation (133) is treated explicitly using third-order Runge-Kutta scheme, and we treat the viscous term implicitly using the Crank-Nicolson scheme. This type of hybrid approach is useful when we design solvers for wall-bounded flows, where we cluster the grid in the boundary layer. For wall-bounded flows, we use dense mesh with a smaller grid size in the boundary layer region to capture the boundary layer flow correctly. First, we can re-write Equation (133) with nonlinear Jacobian term on the right hand side

✝ ✆

$$\frac{\partial \omega}{\partial t} = -J(\omega, \psi) + \nu \nabla^2 \omega. \tag{141}$$

The hybrid third-order Runge-Kutta/Crank-Nicolson scheme can be written as [51]

$$
\omega^{(1)} = \omega^{(n)} + \gamma\_1 \Delta t (-f^{(n)}) + \rho\_1 \Delta t (-f^{(n-1)(2)}) + \frac{\mathfrak{a}\_1 \Delta t \nu}{2} \nabla^2 (\omega^{(1)} + \omega^{(n)}),
\tag{142}
$$

$$
\omega^{(2)} = \omega^{(1)} + \gamma\_2 \Delta t (-\boldsymbol{J}^{(1)}) + \rho\_2 \Delta t (-\boldsymbol{J}^{(n)}) + \frac{a\_2 \Delta t \nu}{2} \nabla^2 (\omega^{(1)} + \omega^{(2)}),
\tag{143}
$$

$$
\omega^{(n+1)} = \omega^{(2)} + \gamma\_3 \Delta t (-J^{(2)}) + \rho\_3 \Delta t (-J^{(1)}) + \frac{\mathfrak{a}\_3 \Delta t \nu}{2} \nabla^2 (\omega^{(n+1)} + \omega^{(2)}),
\tag{144}
$$

where the term <sup>−</sup>*J*(*n*−1)(2) is the nonlinear Jacobian term at the second-step of previous time step at (*n* − <sup>1</sup>). We can choose coefficients in Equations (142)–(144) in such a way that *ρ*<sup>1</sup> vanishes and we do not have to store the solution at second-step of the previous time step. These coefficients are

$$
\alpha\_1 = \frac{8}{15}, \alpha\_2 = \frac{2}{15}, \alpha\_3 = \frac{1}{3}, \tag{145}
$$

$$
\gamma\_1 = \frac{8}{15}, \gamma\_2 = \frac{5}{12}, \gamma\_3 = \frac{3}{4} \tag{146}
$$

$$
\rho\_1 = 0, \rho\_2 = -\frac{17}{60}, \rho\_2 = -\frac{5}{12}.\tag{147}
$$

Besides, to use a hybrid scheme for time, we also design a hybrid finite difference and spectral scheme. We use Arakawa finite difference scheme for the nonlinear Jacobian term and spectral scheme for linear viscous terms. In Fourier space, Equations (142)–(144) can be written as

$$
\hat{\omega}^{(1)} = \hat{\omega}^{(n)} + \gamma\_1 \Delta t (-\J^{(n)}) + \rho\_1 \Delta t (-\J^{(n-1)(2)}) + \frac{a\_1 \Delta t \nu}{2} (-k^2) (\hat{\omega}^{(1)} + \hat{\omega}^{(n)}),
\tag{148}
$$

$$
\hat{\omega}^{(2)} = \hat{\omega}^{(1)} + \gamma\_2 \Delta t (-\vec{f}^{(1)}) + \rho\_2 \Delta t (-\vec{f}^{(n)}) + \frac{\mathfrak{a} 2 \Delta t \nu}{2} (-k^2) (\hat{\omega}^{(1)} + \hat{\omega}^{(2)}),
\tag{149}
$$

$$
\tilde{\omega}^{(n+1)} = \tilde{\omega}^{(2)} + \gamma\_3 \Delta t (-\mathbb{J}^{(2)}) + \rho\_3 \Delta t (-\mathbb{J}^{(1)}) + \frac{a\_3 \Delta t \nu}{2} (-k^2) (\tilde{\omega}^{(n+1)} + \tilde{\omega}^{(2)}),
\tag{150}
$$

where *k*<sup>2</sup> = *m*<sup>2</sup> + *n*2. Here, *m* and *n* are wavenumber in *x* and *y* directions, respectively. Rearranging above equations, we get

$$
\tilde{\omega}^{(1)} = \frac{1 - \frac{a\_1 \Delta t \upsilon k^2}{2}}{1 + \frac{a\_1 \Delta t \upsilon k^2}{2}} \tilde{\omega}^{(n)} + \frac{\gamma\_1 \Delta t}{1 + \frac{a\_1 \Delta t \upsilon k^2}{2}} (-\tilde{f}^{(n)}),
\tag{151}
$$

$$
\tilde{\omega}^{(2)} = \frac{1 - \frac{a\_2 \Delta t \text{vk}^2}{2}}{1 + \frac{a\_2 \Delta t \text{vk}^2}{2}} \tilde{\omega}^{(1)} + \frac{\gamma 2 \Delta t (-\tilde{f}^{(1)}) + \rho\_2 \Delta t (-\tilde{f}^{(n)})}{1 + \frac{a\_2 \Delta t \text{vk}^2}{2}},\tag{152}
$$

$$
\omega^{(n+1)} = \frac{1 - \frac{a\_3 \Delta t \upsilon k^2}{2}}{1 + \frac{a\_3 \Delta t \upsilon k^2}{2}} \omega^{(2)} + \frac{\gamma\_3 \Delta t (-\mathcal{J}^{(2)}) + \rho\_3 \Delta t (-\mathcal{J}^{(1)})}{1 + \frac{a\_3 \Delta t \upsilon k^2}{2}}.\tag{153}$$

We use Arakawa finite difference scheme described in Section 7 to compute the Jacobian term in physical space. The streamfunction is computed from the vorticity field using a spectral method. Once the Jacobian term is available in physical space we convert it to Fourier space using FFT. The relationship between vorticity and streamfunction is Fourier space can be written as

$$
\psi = \frac{\tilde{\omega}}{k^2}.\tag{154}
$$

The time integration using hybrid third-order Runge-Kutta/ Crank-Nicolson scheme in Julia is listed in Listing 32. The Julia implementation to compute the nonlinear Jacobian term for the Arakawa-spectral solver is given in Listing 33.

**Listing 32.** Julia implementation of hybrid third-order Runge-Kutta/ Crank-Nicolson scheme for 2D incompressible Navier-Stokes equation with Arakawa-spectral finite difference scheme. ✞ ☎

```