An efficient spectral Poisson solver for the nirvana-III code: the shearing-box case with vertical vacuum boundary conditions

The stability of a differentially rotating fluid subject to its own gravity is a problem with applications across wide areas of astrophysics–from protoplanetary discs (PPDs) to entire galaxies. The shearing box formalism offers a conceptually simple framework for studying differential rotation in the local approximation. Aimed at self-gravitating, and importantly, vertically stratified PPDs, we develop two novel methods for solving Poisson’s equation in the framework of the shearing box with vertical vacuum boundary conditions (BCs). Both approaches naturally make use of multi-dimensional fast Fourier transforms for computational efficiency. While the first one exploits the linearity properties of the Poisson equation, the second, which is slightly more accurate, consists of finding the adequate discrete Green’s function (in Fourier space) adapted to the problem at hand. To this end, we have revisited the method proposed by Vico et al. (2016) and have derived an analytical Green’s function satisfying the shear-periodic BCs in the plane as well as vacuum BCs, vertically. Our spectral method demonstrates excellent accuracy, even with a modest number of grid points, and exhibits third-order convergence. It has been implemented in the NIRVANA-III code, where it exhibits good scalability up to 4096 CPU cores, consuming less than 6% of the total runtime. This was achieved through the use of P3DFFT, a fast Fourier Transform library that employs pencil decomposition, overcoming the scalability limitations inherent in libraries using slab decomposition. We have introduced two novel spectral Poisson solvers that guarantees high accuracy, performance, and intrinsically support vertical vacuum boundary conditions in the shearing-box framework. Our solvers enable high-resolution local studies involving self-gravity, such as MHD simulations of gravito-turbulence or gravitational fragmentation.

💡 Research Summary

**
The paper addresses the long‑standing challenge of solving Poisson’s equation for self‑gravity in shearing‑box simulations of vertically stratified protoplanetary discs. In the shearing‑box framework the radial (x) direction is sheared in time, breaking strict periodicity, while the azimuthal (y) direction remains periodic and the vertical (z) direction must satisfy free‑space (vacuum) boundary conditions. The authors develop two novel spectral solvers that exploit multidimensional fast Fourier transforms (FFTs) and are fully compatible with these mixed boundary conditions.

The first method, called the Superposition Analytical‑Spectral Hybrid Approach (SASHA), splits the density field into a uniform background ⟨ρ⟩ and a fluctuating part ρ′=ρ−⟨ρ⟩. The uniform component yields an analytically integrable one‑dimensional potential Φ₀(z)=4πG⟨ρ⟩(z²/2−z₀z), where z₀ is the centre‑of‑mass height. The fluctuating component is treated spectrally: after applying a shear‑advection‑by‑Fourier‑interpolation (SAFI) transformation to obtain a fully periodic (x,y) frame, the authors zero‑pad the vertical direction, perform a 3‑D FFT, and solve ˆΦ(k)=−4πG ˆρ(k)/|k|² for all non‑zero wavevectors. Subtracting the mean density guarantees that the volume integral of the source term vanishes, satisfying the mathematical requirement for a periodic Fourier transform. The final potential is the sum Φ=Φ₀+Φ′.

The second method adapts the Vico‑Greengard‑Ferrando (VGF) free‑space Green’s function to the shearing‑box geometry, yielding what the authors term the VGF‑HybridBC solver. After Fourier transforming in the periodic x‑y plane, the Poisson equation reduces to a one‑dimensional screened Helmholtz equation in z: d²Φ̃/dz²−k²Φ̃=4πG ρ̃, where k²=kₓ²+k_y². The corresponding Green’s function is Gₖ(z−z′)=½|z−z′| for k=0 and Gₖ(z−z′)=−½k⁻¹e^{−k|z−z′|} for k≠0. To embed this kernel in a finite computational box while preserving free‑space behaviour, the authors multiply it by a rectangular window (rect) of width L=αL_z (α=1 in the paper) and compute its Fourier transform analytically, obtaining a regularised kernel ˆG_L(k,k_z) that remains finite even at (k,k_z)=(0,0). The potential is then obtained by a single inverse 3‑D FFT of the product ˆΦ=4πG ˆG_L·ˆρ. This approach eliminates the singularity and automatically enforces the vacuum condition at large |z|.

Both solvers are implemented within the finite‑volume MHD code NIRVANA‑III using the P3DFFT library, which provides a pencil decomposition (two‑dimensional domain split) and thus scales efficiently on modern massively parallel architectures. The authors discuss memory considerations: SASHA’s zero‑padding doubles the vertical size but the Hermitian symmetry of real‑valued fields halves the complex storage; VGF‑HybridBC originally required a four‑fold vertical padding, but the authors adopt a two‑fold scheme with a pre‑computation step to keep memory usage manageable.

Benchmark tests include static uniform‑density and Gaussian‑density configurations, as well as a dynamic gravito‑turbulent growth test. All cases exhibit third‑order convergence, with errors reaching machine precision for modest grid resolutions (e.g., 64³). Parallel scalability tests up to 4096 CPU cores show that the Poisson solver consumes less than 6 % of the total runtime, demonstrating excellent strong scaling. The use of P3DFFT yields a factor‑of‑two speed‑up compared with slab‑decomposition FFT libraries, removing a previous bottleneck in self‑gravity calculations.

In conclusion, the paper delivers two high‑accuracy, high‑performance spectral Poisson solvers tailored to the shearing‑box with vertical vacuum boundaries. SASHA offers a simple, robust implementation based on mean‑density subtraction, while VGF‑HybridBC provides a mathematically exact Green’s‑function treatment with slightly higher accuracy. Both methods outperform traditional multigrid or screening approaches in terms of precision, computational cost, and parallel scalability. They open the door to fully resolved, large‑scale local simulations of self‑gravitating discs, including studies of gravito‑turbulence, fragmentation, and magnetic dynamo action in protoplanetary environments.