Numerical PDE III: Spectral-Galerkin Method

In fluid dynamics, spectral-Galerkin methods are commonly used to solve the incompressible Navier-Stokes equations. By transforming the equations into spectral space, the nonlinear terms can be computed efficiently while maintaining conservation properties; and the viscous term can be easily solve implicitly. For incompressible flows, the vorticity-streamfunction formulation is particularly usable, as it eliminates the pressure term and reduces the computational complexity.

In this lecture, we will focus on applying spectral methods to solve the 2D incompressible hydrodynamics equations. Specifically, we will:

Derive the governing equations in spectral space.
Discuss the vorticity-streamfunction formulation and its advantages.
Implement a spectral solver to simulate the evolution of vorticity in a 2D periodic domain.

This introduction sets the stage for understanding how spectral methods leverage mathematical elegance and computational efficiency to solve complex PDEs with remarkable accuracy.

2D Incompressible Hydrodynamics¶

The hydrodynamic equations govern the conservation of mass and momentum in a fluid.

In their compressible form (that we derived in previous lectures), they are written as:

\begin{align} \frac{\partial \rho}{\partial t} + \nabla \cdot (\rho \mathbf{u}) &= 0, \quad \text{(Continuity Equation)} \\ \frac{\partial (\rho \mathbf{u})}{\partial t} + \nabla \cdot (\rho \mathbf{u} \mathbf{u}) &= -\nabla P + \mu \nabla^2 \mathbf{u} + \mathbf{F}, \quad \text{(Momentum Equation)} \end{align}

(1)

where:

$\rho$ is the density,
$\mathbf{u}$ is the velocity field,
$P$ is the pressure,
$\mu$ is the dynamic viscosity, and
$\mathbf{F}$ is an external force.

In the incompressible limit, the sound speed approaches infinite $c \rightarrow \infty$ . For simplicity, the density $\rho$ can be assumed constant, and the continuity equation reduces to the incompressibility condition:

\begin{align} \nabla \cdot \mathbf{u} = 0. \end{align}

(2)

Substituting this condition into the momentum equation simplifies the Navier-Stokes equations to:

\begin{align} \frac{\partial \mathbf{u}}{\partial t} + (\mathbf{u} \cdot \nabla) \mathbf{u} &= -\nabla p + \nu \nabla^2 \mathbf{u} + \mathbf{f}, \end{align}

(3)

where $\nu \equiv \mu / \rho$ is the kinematic viscosity and $\mathbf{f} \equiv \mathbf{F} / \rho$ is usually still called external “force”. These equations describe the flow of incompressible fluids. They are widely used in modeling low Mach number fluid from small-scale laboratory experiments to large-scale geophysical flows.

While the incompressible Navier-Stokes equations describes three-dimensional flows, many physical systems can be effectively approximated as two-dimensional. For example:

Atmospheric flows and ocean currents are largely horizontal due to their vast spatial extent compared to their depth.
Thin liquid films and confined flows are geometrically restricted to two dimensions.

To simplify the mathematical and computational treatment of 2D incompressible hydrodynamics, the governing equations are often reformulated in terms of the vorticity $w$ and the streamfunction $\psi$ .

This formulation has several advantages: it eliminates the pressure term from the equations, reduces the number of variables, and ensures incompressibility is automatically satisfied.

Vorticity-Streamfunction Formulation¶

For 2D incompressible flows, the velocity field $\mathbf{u} = (u_x, u_y)$ can be expressed in terms of a scalar function, the streamfunction $\psi(x, y, t)$ , as:

\begin{align} \mathbf{u} = \nabla \times (\psi \mathbf{\hat{z}}), \end{align}

(4)

where $\mathbf{\hat{z}}$ is the unit vector perpendicular to the 2D plane. In component form:

\begin{align} u_x = \frac{\partial \psi}{\partial y}, \quad u_y = -\frac{\partial \psi}{\partial x}. \end{align}

(5)

This representation automatically satisfies the incompressibility condition:

\begin{align} \nabla \cdot \mathbf{u} = \frac{\partial u_x}{\partial x} + \frac{\partial u_y}{\partial y} = 0. \end{align}

(6)

# Let's visualize the streamfunction

from jax import numpy as np

Lx, Ly = 2*np.pi, 2*np.pi  # domain size
Nx, Ny = 256,     256      # number of grid points

x, dx = np.linspace(-Lx/2, Lx/2, Nx, endpoint=False, retstep=True)
y, dy = np.linspace(-Ly/2, Ly/2, Ny, endpoint=False, retstep=True)
x, y  = np.meshgrid(x, y)

psi = np.exp(-(x*x + y*y)/2)
ux  = -y * psi
uy  = +x * psi

from matplotlib import pyplot as plt

fig, (ax0, ax1) = plt.subplots(1, 2, figsize=(8,4))

ax0.quiver(x[::8,::8], y[::8,::8], ux[::8,::8], uy[::8,::8])
ax0.set_aspect('equal')

ax1.contour(x, y, psi)
ax1.set_aspect('equal')

The vorticity $\omega$ is a scalar quantity in 2D, defined as the curl of the velocity field:

\begin{align} \omega = \nabla \times \mathbf{u}. \end{align}

(7)

Using the velocity components in terms of the streamfunction, the vorticity can be written as:

\begin{align} \omega = \frac{\partial u_y}{\partial x} - \frac{\partial u_x}{\partial y} = -\nabla^2 \psi. \end{align}

(8)

Thus, the vorticity and streamfunction are related by the Poisson equation:

\begin{align} \omega = -\nabla^2 \psi. \end{align}

(9)

This relationship allows the velocity field to be computed from the vorticity by solving the Poisson equation for $\psi$ , and then taking derivatives of $\psi$ to find $u_x$ and $u_y$ .

omega = -(x*x + y*y - 2) * psi

fig = plt.figure(figsize=(8,4))

ax0 = fig.add_subplot(1, 2, 1, projection='3d')
ax1 = fig.add_subplot(1, 2, 2, projection='3d')

ax0.plot_surface(x, y, omega)
ax0.set_title(r"Vorticity $\omega$")
ax1.plot_surface(x, y, psi)
ax1.set_title(r"Streamfunction $\psi$")

The vorticity transport equation is derived from the incompressible Navier-Stokes equations. Taking the curl of the momentum equation eliminates the pressure gradient term,

\begin{align} \frac{\partial\omega}{\partial t} + \mathbf{u} \cdot \nabla\omega = \nu \nabla^2\omega + \mathbf{f}_\omega, \end{align}

(10)

where:

$\mathbf{u} \cdot \nabla\omega$ represents the nonlinear advection of vorticity,
$\nu \nabla^2\omega$ accounts for viscous diffusion,
$\mathbf{f}_\omega$ is the vorticity-specific forcing term.

The term $\mathbf{u} \cdot \nabla\omega$ can be expanded using the velocity components as:

\begin{align} \mathbf{u} \cdot \nabla\omega = u_x \frac{\partial\omega}{\partial x} + u_y \frac{\partial\omega}{\partial y}. \end{align}

(11)

By substituting $u_x$ and $u_y$ in terms of $\psi$ , the nonlinear advection term is rewritten as the Jacobian determinant:

\begin{align} J(\psi, \omega) = \frac{\partial \psi}{\partial x} \frac{\partial\omega}{\partial y} - \frac{\partial \psi}{\partial y} \frac{\partial\omega}{\partial x}. \end{align}

(12)

Thus, the vorticity transport equation becomes:

\begin{align} \frac{\partial\omega}{\partial t} - J(\psi, \omega) = \nu \nabla^2\omega + \mathbf{f}_\omega. \end{align}

(13)

Ekman Damping¶

Ekman damping models frictional effects caused by the interaction of the fluid with a boundary layer. It acts as a large-scale energy sink and is proportional to the vorticity $-\mu\omega$ , where $\mu$ is the Ekman coefficient.

Including this term in the vorticity transport equation gives:

\begin{align} \frac{\partial\omega}{\partial t} - J(\psi, \omega) = \nu \nabla^2\omega - \mu\omega + \mathbf{f}_\omega. \end{align}

(14)

Ekman damping is particularly relevant in geophysical systems, where it represents energy dissipation due to the Earth’s surface or ocean floors.

$\beta$ -Plane Model¶

The $\beta$ -plane approximation models the variation of the Coriolis parameter with latitude. In the vorticity equation, this introduces a term proportional to the northward velocity component, i.e., $\beta u_y$ , where $\beta$ is the linear expansion coefficient of the Coriolis parameter.

Including this term in the vorticity transport equation gives:

\begin{align} \frac{\partial\omega}{\partial t} - J(\psi, \omega) + \beta u_y = \nu \nabla^2\omega - \mu\omega + \mathbf{f}_\omega. \end{align}

(15)

The beta-plane term is crucial for studying large-scale atmospheric and oceanic dynamics, as it leads to phenomena such as Rossby waves and geostrophic turbulence.

In addition to damping and planetary rotation effects, many fluid systems are subject to external forcing that injects energy into the flow. External forcing terms are introduced into the vorticity equation through $\mathbf{f}_\omega$ , and their form strongly influences the resulting turbulence and flow structures.

Kolmogorov Forcing¶

Kolmogorov forcing is one of the most widely used deterministic forcing schemes in studies of two-dimensional turbulence. It drives the flow through a sinusoidal body force that acts at a specific wavenumber $k_f$ . In vorticity form, the forcing can be written as

\begin{align} \mathbf{f}_\omega(x, y) = f_0 \cos(k_f y), \end{align}

(16)

where $f_0$ is the forcing amplitude and $k_f$ is the forcing wavenumber.

This forcing continuously injects energy at large scales and maintains a statistically steady turbulent state by balancing viscous dissipation. Kolmogorov forcing is particularly useful for studying cascades of energy and enstrophy in spectral simulations, since it creates a controlled, spatially periodic source of turbulence. However, note that Kolmogorov forcing does not provide constant energy input to the turbulence system.

Random Forcing¶

An alternative is to model forcing as a stochastic process. Random forcing mimics the unpredictable nature of real-world turbulence drivers, such as atmospheric convection, wind stress on the ocean surface, or small-scale instabilities.

In practice, random forcing is often applied in Fourier space, where modes in a narrow band of wavenumbers around $k_f$ are driven with random amplitudes and phases:

\begin{align} \hat{f}_\omega(\mathbf{k}, t) = \begin{cases} A(\mathbf{k}, t), & k \in [k_f - \Delta k/2, k_f + \Delta k/2], \\ 0, & \text{otherwise}, \end{cases} \end{align}

(17)

where $A(\mathbf{k}, t)$ is a stochastic process, commonly chosen as Gaussian white noise with prescribed variance.

Random forcing produces flows that are less structured than Kolmogorov forcing, and is often employed to study statistical steady states of turbulence. It provides a more realistic representation of natural geophysical systems, where energy input is irregular and multi-scale. Note that, however, random forcing provides a statistically constant energy input to the turbulence system.

Inverse Energy Cascade¶

2D turbulence shows unique features that distinguish them from 3D flows: Conservation of Enstrophy.

In 2D, the vorticity $\omega = \nabla \times \mathbf{u}$ is a scalar field. Its evolution is governed by the vorticity transport equation that we derive, which conserves both energy and enstrophy in the absence of dissipation:

\begin{align} E &= \frac{1}{2} \int |\mathbf{u}|^2 dx dy\\ Z &= \frac{1}{2} \int \omega^2 dx dy. \end{align}

(18)

Energy conservation governs the total kinetic energy of the system, while enstrophy conservation introduces a second constraint that strongly influences the flow dynamics.

A striking feature of 2D turbulence is the inverse energy cascade. In 3D turbulence, energy flows from large scales (low wavenumbers) to small scales (high wavenumbers) and is dissipated by viscosity. In 2D, however, energy flows in the opposite direction, from small scales to large scales, leading to the formation of large, coherent structures like cyclones and anticyclones. This behavior is directly tied to the dual conservation of energy and enstrophy.

Spectral Equation¶

The vorticity-streamfunction formulation simplifies the governing equations of 2D incompressible hydrodynamics, making them naturally suited to spectral methods. By representing the vorticity and streamfunction in terms of Fourier series, we can efficiently compute derivatives and nonlinear terms in spectral space.

Fourier Representation of Vorticity and Streamfunction¶

For a periodic domain of size $L_x \times L_y$ , the vorticity $\omega(x, y, t)$ and streamfunction $\psi(x, y, t)$ can be expanded as Fourier series:

\begin{align} \omega(x, y, t) &= \sum_{k_x, k_y} \hat{\omega}_{k_x, k_y}(t) e^{i (k_x x + k_y y)}, \\ \psi (x, y, t) &= \sum_{k_x, k_y} \hat{\psi }_{k_x, k_y}(t) e^{i (k_x x + k_y y)}, \end{align}

(19)

where:

$\hat{\omega}_{k_x, k_y}$ and $\hat{\psi}_{k_x, k_y}$ are the Fourier coefficients of vorticity and streamfunction, respectively.
$k_x$ and $k_y$ are the wavenumbers in the $x$ - and $y$ -directions, given by:
$\begin{align} k_x = \frac{2\pi n_x}{L_x}, \quad k_y = \frac{2\pi n_y}{L_y}, \quad n_x, n_y \in \mathbb{Z}. \end{align}$
(20)

Poisson Equation in Spectral Space¶

The Poisson equation relates the vorticity $\omega$ and the streamfunction $\psi$ :

\begin{align} \omega = -\nabla^2 \psi. \end{align}

(21)

Recalling in spectral space, derivatives with respect to $x$ and $y$ transform into simple multiplications by $ik_x$ and $ik_y$ , respectively. For example:

\begin{align} \frac{\partial \omega}{\partial x } \rightarrow ik_x \hat{\omega}_{k_x, k_y}, \quad \frac{\partial^2\omega}{\partial x^2} \rightarrow -k_x^2\hat{\omega}_{k_x, k_y}. \end{align}

(22)

This property makes spectral methods computationally efficient, as derivatives reduce to element-wise multiplications.

Hence, in Fourier space, the Laplacian $\nabla^2$ becomes multiplication by $-k^2$ , where $k^2 = k_x^2 + k_y^2$ . Thus, the Poisson equation transforms into:

\begin{align} \hat{\psi}_{k_x, k_y} = \frac{\hat{\omega}_{k_x, k_y}}{k^2}, \quad k^2 \neq 0. \end{align}

(23)

Note that the $k^2 = 0$ mode corresponds to the mean value of $\psi$ , which can be set to zero for flows with no net circulation.

Vorticity Transport Equation in Spectral Space¶

Recalling the vorticity transport equation in real space is:

\begin{align} \frac{\partial w}{\partial t} - J(\psi, w) + \beta u_y = \nu \nabla^2 w - \mu w + \mathbf{f}_w. \end{align}

(24)

In spectral space:

The Laplacian term $\nabla^2 \omega$ transforms into $-k^2 \hat{\omega}_{k_x, k_y}$ .
The Ekman damping term $\mu \omega$ becomes $-\mu\hat{\omega}_{k_x, k_y}$ .
The Coriolis term $\beta u_y$ becomes $-\beta ik_x \hat{\omega}_{k_x, k_y}/k^2$ .
The nonlinear term $J(\psi, \omega)$ is evaluated in real space and transformed back to spectral space using the Fast Fourier Transform (FFT).

The equation in spectral space becomes:

\begin{align} \frac{\partial}{\partial t}\hat{\omega}_{k_x, k_y} = \widehat{J(\psi, \omega)}_{k_x, k_y} - \left(\nu k^2 + \mu - \beta \frac{ik_x}{k^2}\right)\hat{\omega}_{k_x, k_y} + \hat{f}_{\omega k_x, k_y}. \end{align}

(25)

# Create spectral grid

kx = 2*np.pi * np.fft.fftfreq(Nx, d=Lx/Nx)
ky = 2*np.pi * np.fft.fftfreq(Ny, d=Ly/Ny)
kx, ky = np.meshgrid(kx, ky)

kk  = kx*kx + ky*ky
ikk = 1.0 / (kk + 1.2e-38)
ikk = ikk.at[0,0].set(0)  # avoid multiplied by infinity for mean mode

print(ikk[0:4,0:4])

from jax.numpy.fft import fft2

# Initialize vorticity and streamfunction in spectral space
omega = np.exp(-0.5*(x*x + y*y)*16) # example initial vorticity

# Transform to Fourier (spectral) space
Omega = fft2(omega)
Psi   = ikk * Omega

from jax.numpy.fft import ifft2

# Obtain velocity field
def vel(Psi):
    psi_x = ifft2(1j * kx * Psi).real
    psi_y = ifft2(1j * ky * Psi).real
    return psi_y, -psi_x

def plot(Omega, skip=Nx//32, cmap='bwr'):
    Psi    = ikk * Omega
    ux, uy = vel(Psi)
    plt.imshow(ifft2(Omega).real, origin='lower', extent=[-Lx/2,Lx/2,-Ly/2,Ly/2], cmap=cmap)
    plt.quiver(x[::skip,::skip], y[::skip,::skip], ux[::skip,::skip], uy[::skip,::skip])

plot(Omega, cmap='Reds')

# Compute Jacobian determinant in real space
def Jdet(Psi, Omega):
    psi_x   = ifft2(1j * kx * Psi  ).real
    psi_y   = ifft2(1j * ky * Psi  ).real
    omega_x = ifft2(1j * kx * Omega).real
    omega_y = ifft2(1j * ky * Omega).real
    return fft2(psi_x * omega_y - psi_y * omega_x)

J = Jdet(Psi, Omega)
print(np.max(abs(ifft2(J).imag)))

Handling Aliasing Errors¶

Nonlinear terms, such as $J(\psi, w)$ , involve products in real space that translate into convolutions in spectral space. Due to the finite resolution of the grid, these convolutions can introduce spurious interactions between modes known as aliasing errors.

The 2/3 rule is a widely used dealiasing technique that truncates Fourier modes beyond $2/3$ of the maximum wavenumber. For a grid with $N$ points in each dimension:

Retain modes for $|k_x|, |k_y| \leq N/3$ .
Set all other modes to zero.

The 2/3 rule ensures that spurious contributions from nonlinear interactions fall outside the resolved spectral range.

# Apply the 2/3 rule
def twothird(F):
    lx = Nx // 3
    ly = Ny // 3
    F  = F.at[:,lx:-lx+1].set(0)
    F  = F.at[ly:-ly+1,:].set(0)
    return F

# Apply de-aliasing to the Jacobian
J = twothird(J)
plt.imshow(np.log(abs(J)))
print(np.max(abs(ifft2(J).imag)))

Time Integration of the Spectral Equations¶

After reformulating the vorticity-streamfunction equations in spectral space, the next step is to integrate the vorticity transport equation in time. Time integration involves advancing the Fourier coefficients of vorticity $\hat{\omega}_{k_x, k_y}$ while accurately handling the nonlinear and linear terms. This section discusses suitable time-stepping schemes, their implementation, and how they are applied to the spectral representation of the equations.

Time-Stepping Schemes¶

We need time integration methods for the spectral vorticity transport equation must balance stability, accuracy, and computational efficiency. The commonly schemes to consider for spectral methods are:

Explicit Schemes¶

Explicit schemes, such as the Runge-Kutta (RK) family, compute the solution at the next time step based on known quantities at the current time step. They are easy to implement and efficient for problems dominated by advection or nonlinear dynamics.

Implicit Schemes¶

Implicit schemes, such as Crank-Nicolson, are unconditionally stable for linear terms but require solving a system of equations at each time step. It is uncommon to use a purely implicit scheme to integrate even the nonlinear term in spectral methods.

Semi-Implicit Methods¶

However, it makes sense to use implicit schemes for the linear terms (e.g., viscous diffusion) to allow larger time steps, while keep explicit schemes for nonlinear terms.

A common approach splits the equation into linear and nonlinear parts:

\begin{align} \frac{\partial \hat{\omega}}{\partial t} = N(\hat{\omega}) + L \hat{\omega}, \end{align}

(26)

where:

$N(\hat{\omega})$ represents nonlinear terms (e.g., Jacobian).
$L(\hat{\omega})$ represents linear terms (e.g., viscous diffusion, Ekman damping, Coriolis effects),

The linear terms are treated implicitly, while the nonlinear terms are advanced explicitly. For example, the time-discretized form is:

\begin{align} \hat{\omega}^{n+1} = (1 - \Delta t L)^{-1}[\hat{\omega}^n + \Delta t N(\hat{\omega}^n)]. \end{align}

(27)

This approach combines the stability of implicit methods with the simplicity of explicit methods.

Time Integration of the Vorticity Transport Equation¶

The spectral vorticity transport equation is:

\begin{align} \frac{\partial}{\partial t}\hat{\omega}_{k_x, k_y} = \widehat{J(\psi, \omega)}_{k_x, k_y} - \left(\nu k^2 + \mu - \beta \frac{ik_x}{k^2}\right)\hat{\omega}_{k_x, k_y} + \hat{f}_{\omega k_x, k_y}. \end{align}

(28)

Breaking this into linear and nonlinear parts:

Nonlinear Terms: Jacobian determinant $\widehat{J(\psi, w)}_{k_x, k_y}$ and external forcing $\hat{f}_{k_x, k_y}$ .
Linear Terms: Viscous diffusion $-\nu k^2 \hat{w}$ , Ekman damping $-\mu \hat{w}$ , and Coriolis effect $\beta ik_x \hat{w}/k^2$ .

Using a semi-implicit scheme:

Advance the nonlinear terms explicitly to $\hat{\omega}_{k_x, k_y}^{*n+1}$ by, e.g., Runge-Kutta methods.
Solve the linear terms implicitly using the formula:
$\begin{align} \hat{\omega}_{k_x, k_y}^{n+1} = \frac{\hat{\omega}_{k_x, k_y}^{*n+1}}{1 + \Delta t (\nu k^2 + \mu - \beta ik_x/k^2)}. \end{align}$
(29)

# Initialize vorticity and streamfunction in spectral space
from jax.random import key, normal

ux = np.where(y >= 0, 1, -1) # Example initial vorticity
uy = 0.5 * normal(key(0), shape=x.shape)

Omega = twothird(1j * kx * fft2(uy) - 1j * ky * fft2(ux))
plot(Omega)

# Define simulation parameters

N    = Nx* 10 # number of time steps
S    = Nx//10 # number of skip between outputs

dt   = 1.0/Nx # time step
nu   = 0.001  # viscosity
mu   = 0.0    # Ekman damping
beta = 0.0    # Coriolis parameter

from jax import jit

@jit
def nonlinear(Omega):
    # Obtain stream function
    Psi = ikk * Omega

    # Compute nonlinear term (Jacobian determinant)
    J = Jdet(Psi, Omega)
    J = twothird(J)  # Apply 2/3 rule

    return J

@jit
def semiimplicit(Omega):

    # Explicit RK4
    k1 = nonlinear(Omega                )
    k2 = nonlinear(Omega + 0.5 * dt * k1)
    k3 = nonlinear(Omega + 0.5 * dt * k2)
    k4 = nonlinear(Omega +       dt * k3)
    Omega = Omega + dt * (k1/6 + k2/3 + k3/3 + k4/6)

    # TODO: External forcing

    # Implicit backward Euler
    Omega = Omega / (1 + dt * (nu * kk + mu - (1j * beta) * (kx * ikk)))

    return Omega

from tqdm import tqdm

! rm -rf imgs && mkdir imgs

# Time-stepping loop
for n in tqdm(range(N//S)):
    plot(Omega)
    plt.savefig(f'imgs/{n:04d}.png', dpi=144, bbox_inches='tight')
    plt.close()

    for j in range(S):
        Omega = semiimplicit(Omega)

# Plot and save the last frame
plot(Omega)
plt.savefig(f'imgs/{n+1:04d}.png', dpi=144, bbox_inches='tight')

# Let's combine the png files into a movie:

! rm -f KH.mpg && ffmpeg -i imgs/%04d.png -loglevel panic -qmax 2 KH.mpg

To solve the full vorticity transport equation numerically, we adopted a semi-implicit time integration scheme:

The nonlinear advection term is integrated explicitly with a high-order Runge–Kutta (RK4) method, ensuring accuracy in capturing nonlinear interactions.
The linear dissipative terms (viscosity and Ekman damping) are treated implicitly using backward Euler, which provides unconditional stability against stiff linear decay.

This combination gives us a robust, efficient method capable of evolving turbulent flows for long times while maintaining both stability and accuracy.

Closing Thought¶

With this, we conclude the first half of ASTR 513 on computational methods. We have built up the foundations of numerical analysis and applied them to some of the most beautiful and challenging problems in physics: nonlinear dynamics and turbulence.

The real power of computational physics is not just in solving equations, but in enabling us to explore phenomena beyond the reach of experiment and pure theory. May the tools you now hold inspire you to look deeper, simulate boldly, and discover what lies hidden in the flow of nature.

Numerical PDE III: Spectral-Galerkin Method

2D Incompressible Hydrodynamics¶

Vorticity-Streamfunction Formulation¶

Ekman Damping¶

β\betaβ-Plane Model¶

Kolmogorov Forcing¶

Random Forcing¶

Inverse Energy Cascade¶

Spectral Equation¶

Fourier Representation of Vorticity and Streamfunction¶

Poisson Equation in Spectral Space¶

Vorticity Transport Equation in Spectral Space¶

Handling Aliasing Errors¶

Time Integration of the Spectral Equations¶

Time-Stepping Schemes¶

Explicit Schemes¶

Implicit Schemes¶

Semi-Implicit Methods¶

Time Integration of the Vorticity Transport Equation¶

Closing Thought¶

$\beta$ -Plane Model¶