A semi-analytical pseudo-spectral method for 3D Boussinesq equations of rotating, stratified flows in unbounded cylindrical domains

Abstract

We present a pseudo-spectral method for solving the three-dimensional Boussinesq equations in unbounded cylindrical domains, specifically tailored for rotating, stably stratified flows subject to strong azimuthal shear. To effectively capture the global geometry without sacrificing spectral accuracy, the spatial discretization employs Fourier expansions in the azimuthal and axial directions alongside mapped associated Legendre polynomials in the radial direction. This basis spans the semi-infinite domain while analytically resolving the coordinate singularity at the origin. While this spectral framework ensures high spatial fidelity, the temporal integration of these rotating shear flows presents a formidable computational challenge due to the numerical stiffness driven by fast restorative wave forces and rapid background advection. To circumvent this, we develop an exponential time differencing (ETD) scheme that analytically integrates the fully coupled linear operator, including the radially dependent advective cross terms. By encoding the physical resonance characteristics and stability limits of the background flow directly into the integration operators, the proposed ETD formulation removes the numerical stability constraints imposed by the background shear and stratification. This permits integration time steps scaled by the slow macroscopic evolution of the physical instabilities rather than the fast background kinematics, offering significant performance gains over standard mixed implicit-explicit schemes. The method's accuracy and stability are validated through the precise conservation of energy and angular momentum, establishing a robust framework for simulating instabilities in astrophysical and geophysical vortices.

Figures (6)

• Figure 1: Mapping between the global cylindrical coordinates $(r, \theta, z)$ and the local shearing box $(\tilde{x}, \tilde{y}, \tilde{z})$ centered at $(r_0, \theta_0)$. The local coordinates are defined as $\tilde{y} = r - r_0$, $\tilde{x} = -r_0(\theta - \theta_0)$, and $\tilde{z} = z$, such that the local unit vectors correspond to $\tilde{\bm{y}} = \hat{\bm{r}}$ and $\tilde{\bm{x}} = -\hat{\bm{\theta}}$. Gravity acts in the $-z$ direction (pointing into the page). The local frame co-rotates with the background at an angular velocity $\tilde{\Omega} = \bar{\Omega}(r_0) + \Omega$, where $\bar{\Omega}$ represents the radially dependent angular velocity of the global background flow (typically following a power-law profile, such as the Keplerian $\bar{\Omega}\propto r^{-3/2}$ in accretion disks) and $\Omega$ is the constant background rotation rate of the global frame. The differential rotation of the background flow is locally approximated as a linear shear profile $\bar{u}_{\tilde{x}} = -\sigma \tilde{y}$, where the constant shear rate is defined as $\sigma = r_0 (\dd\bar{\Omega}/\dd{r})|_{r_0}$ ($\sigma < 0$ in the depicted example).
• Figure 2: Single-node strong scaling performance for a $256 \times 256 \times 128$ physical grid on dual-socket AMD EPYC 7763 processors ($128$ total cores). Legend labels denote the total step time, the pointwise nonlinear evaluation ($\mathcal{N}$), the ETD advancement ($\bm{\mathsf{E}}, \bm{\mathsf{NL}}$), the full space conversions (PPP $\leftrightarrow$ FFF), the localized radial expansions, and the solenoidal projection ($\mathbb{P}$). Note that the plotted timings for the space conversions and radial expansions represent the average duration per single operation, whereas a complete integration step encompasses multiple such transformations. The dashed line represents ideal linear scaling.
• Figure 3: Discrete unstable modes for a stratified Lamb-Oseen vortex at $\bar{N}=5$. (a) Growth rates $\lambda$ and (b) wave frequencies $\omega$ as functions of the scaled axial wavenumber $k/\bar{N}$. Open circles represent the eigenvalues computed by the present global spectral eigenvalue solver. Three points on the first mode branch ($k/\bar{N} = 1.4, 2.4, 4.8$) are distinctly marked for further structural analysis. Solid black lines denote the benchmark results of dizes_2008 obtained via a shooting method with complex contour deformation. The numerical spatial discretization employs $N_r = 252$ radial collocation points, a radial spectral truncation of $N = 244$, and a mapping parameter $L = 6.0$. Only the first three unstable branches are plotted.
• Figure 4: Frequency maps (left column) and normalized spatial eigenmode profiles $\tilde{u}_\theta$ (right column) tracking the topological evolution of the first mode at the three distinct scaled wavenumbers, as marked in figure \ref{['fig:dispersion_validation']}: (a) $k/\bar{N} = 1.4$, (b) $k/\bar{N} = 2.4$, and (c) $k/\bar{N} = 4.8$. The left panels plot the radial profiles of the local epicyclic frequencies (dashed magenta lines) and the co-rotation frequency (solid red line), alongside the mode's constant wave frequency $\omega$ (horizontal green line). The intersections of the horizontal wave frequency line with these respective curves define the radial locations of the turning points and the barotropic critical layer. At $k/\bar{N} = 4.8$, the loss of the inner turning points exposes the mode directly to the barotropic critical layer, triggering the Gibbs-like spectral ringing observed in the spatial profile.
• Figure 5: Temporal convergence of the ETD scheme. The fractional errors in global kinetic energy ($E_\mathrm{K}$) and available potential energy ($E_\mathrm{AP}$) follow the theoretical second-order scaling ($\propto \Delta t^2$) as $\Delta t$ is reduced. Errors are evaluated at $t=20$, corresponding to approximately $3.18$ core rotation periods ($T_{\text{core}} = 2\pi$) of the background Lamb-Oseen vortex, ensuring the measurement captures fully developed nonlinear dynamics.