Combination Finite Difference & Spectral Solution to Advection-Diffusion PDE

Vorticity $\omega$ is a measure of local rotation (in a fluid). It is related to the angular velocity of a solid. Vorticity is generally a vector (valued function). However, it can be viewed as a scalar in 2D because it is always perpendicular to the plane. If fluid velocity $\vec{v}$ (in 2D) has components $u$ and $v$ $$\begin{equation} \vec{v}(t) = \begin{pmatrix} u(t) \\ v(t) \end{pmatrix} \end{equation}$$ then these vector components are related to a scalar quantity called the \textbf{stream function} $\psi$ as follows $$\begin{align} u = -\psi_y \\ v = \psi_x \end{align}$$ so that $$\begin{equation} \vec{v}(t) = \begin{pmatrix} u(t) \\ v(t) \end{pmatrix} = \begin{pmatrix} -\psi_y \\ \psi_x \end{pmatrix} \end{equation}$$ So $\psi$ is a scalar valued function while the components of the fluid velocity vector can be derived from it. $$\begin{equation} \label{eq:hw4_laplacian_eqn} \nabla^2 \psi = \omega \end{equation}$$ Vorticity-Streamfunction PDE $$\begin{align} \omega_t + [\psi,\omega] &= \eta \nabla^2 \omega \\ \omega_t + \psi_x \omega_y - \psi_y \omega_x &= \eta ( \omega_{xx} + \omega_{yy} ) \\ \omega_t + \psi_x \omega_y - \psi_y \omega_x &= \eta D_L \omega \\ \omega_t &= -(\psi_x \omega_y - \psi_y \omega_x) + \eta D_L \omega \end{align}$$ where $D_L$ is the Laplacian operator $\frac{\partial^2}{\partial x^2} + \frac{\partial^2}{\partial y^2}$. With discretization in space, we can write a finite difference approximation for the linear operator $D_L$. This allows us to solve the Poisson equation by solving the linear system: $$\begin{equation} D_L \vec{\psi} = \vec{\omega} \end{equation}$$ This obtains $\psi$ from $\omega$ at step n. If we use a finite difference approximation in time too we can write $$\begin{align} \omega_t &= -(\psi_x \omega_y - \psi_y \omega_x) + \eta D_L \omega \\ \frac{\omega^{k+1} - \omega^k }{ \Delta t} &= -(\psi_x \omega_y - \psi_y \omega_x) + \eta D_L \omega \end{align}$$ with RHS $$f = -(\psi_x \omega_y - \psi_y \omega_x) + \eta D_L \omega$$ The following video is made in MATLAB after integrating the system with ODE45. At each time step, the Poisson equation is solved in Fourier space (to speed things up)