This repository contains codes to study Quasi-Geostrophic Magnetohydrodynamics (QG-MHD).
By: Jonathan Tessier and Francis J. Poulin
Nonlinear Evolution
ParallelShenfun solves the nonlinear QG-MHD equations formulated for the potential vorticity $q$, and magnetic streamfunction $A$. Assuming an optional stationary background state $(\bar q ,\bar A)$, the code evolves the pertubations $(q' ,A')$, where $q = \bar q + q'$, $A = \bar A + A'$, on a doubly-periodic domain using a pseudo-spectral method with a smooth filter and AB3 timestepping. The code is written in Python and runs in parallel using Shenfun.
The nonlinear equations read (ignoring diffusion on the background):
$$
F = \frac{L}{R_d}, \quad M = \frac{V_A}{U}, \quad R_e = \frac{\nu}{UL}, \quad R_m = \frac{\kappa}{UL},
$$
where $R_d=\sqrt{gH}/f$ is the external Rossby radius of deformation and $V_A=B_0/\sqrt{\mu\rho}$ is the Aflvén wave speed, for gravity $g$, mean depth $H$, Coriolis frequency $f$, magnetic field strength $B_0$, magnetic permeability $\mu$ and fluid density $\rho$. $R_e$ and $R_m$ are the hydrodynamic and magnetic Reynolds numbers with viscosity $\nu$ and magnetic diffusivity $\kappa$.
To turn off the Lorentz force, set $M=0$. $A$ is then a passive tracer advected around by the flow. To consider a rigid lid evolution (equivalent to 2D-MHD), set $F=0$.
Linear Stability Analysis
LinearStability solves the linear stability of parallel shear flows such as a Bickley jet. Linearizing the QG-MHD equations about a background state in geostrophic balance, we neglect the advection of the perturbation by the pertubation but still maintain the advection of the pertrubation by the mean state and the advection of the mean state by the perturbation. Assuming a normal mode decomposition on the perturbations,
$$
\psi' = \hat \psi(y) \exp[i (k x - \omega t)], \quad\quad A' = \hat A(y) \exp[i (k x - \omega t)],
$$
we can generate an eigenvalue problem where the eigenvalue $c=\omega/k$ determines the growth rate $\omega_i$ and the phase speed $\omega_r$ of the unstable mode, and the eigenfunctions $(\hat \psi, \hat A)$ show the spatial structure in the y-direction. This transforms the Laplacian operator such that $\hat\nabla^2 = \partial_{yy}-k^2$ and we find the generalized EVP to read (for background velocity $U(y)$ and field $B(y)$ )