© The Authors 2026

1 Introduction

Many astrophysical phenomena require accounting for the gravity of a dilute gas or plasma, often referred to as self-gravity (SG). This includes processes such as core collapse of molecular clouds (Shu 1977; Norman et al. 1980; He & Ricotti 2023), episodic outbursts in FU Orionis stars (Armitage et al. 2001; Vorobyov & Basu 2006), and accretion and angular momentum transport in young protoplanetary discs, which can lead to fragmentation for efficient cooling (Gammie 2001; Kratter & Lodato 2016; Baehr & Zhu 2021). Recent theoretical and numerical studies in the local frame, also known as the shearing box approximation (see, e.g. Gressel & Ziegler 2007; Stone & Gardiner 2010), have demonstrated the potential of gravitational instability to amplify weak magnetic fields through a process called the gravitoturbulent dynamo (Riols & Latter 2019; Deng et al. 2020; Löhnert & Peeters 2023).

The numerical computation of SG is crucial for these processes, which necessitates to solve the Poisson equation $${\mathrm{\nabla }}^2\mathrm{\Phi }=4\pi G\rho ,$$(1)

where Φ and ρ represent the gravitational potential and volume density, respectively. A standard approach in Cartesian uniform grids involves discretizing this equation, leading to a linear system typically solved using iterative multigrid methods, as demonstrated by Ziegler (2005) in the finite-volume magnetohydrodynamic (MHD) code NIRVANA-III. However, this requires an accurate evaluation of the SG potential at the numerical domain boundaries, which generally demands resorting to a multipole expansion. The screening method, which enables computation of the exact potential at the domain boundaries, is more accurate (James 1977; Moon et al. 2019; Gressel & Ziegler 2024).

The numerical computation of the gravitational potential can also be performed using spectral methods. These techniques often involve solving the Poisson equation in Fourier space, $$\hat{\mathrm{\Phi }}(\mathit{k})=-\frac{4\pi G\hat{\rho }}{\parallel \mathit{k}{\parallel }^2}(\text{\hspace{0.17em}for\hspace{0.17em}}\mathit{k}\ne \mathbf{0}),$$(2)

followed by an inverse transform in order to obtain the potential in real space. Their primary advantages are accuracy and efficiency, and – with an optimum numerical complexity of N log (N) – they are well suited for large-scale problems. The spectral approach, however, relies on the discrete Fourier transform (DFT), or, more precisely, the fast Fourier transform (FFT); the crucial advantage of the latter is that it overcomes the quadratic complexity order of naive DFT algorithms in favour of the above-mentioned scaling. However, a central assumption in such discrete methods is that the input signal represents one complete period of a periodic signal. This is a significant limitation for gravitational problems, as it implies the existence of infinite copies of the numerical window in all directions, which is often unrealistic or impractical. Additionally, for consistency, the volume integral of the source term of the Poisson equation must vanish (Binney & Tremaine 2008; Mandal et al. 2023). Therefore, in Cartesian grids it becomes desirable to develop a full spectral solver that not only ensures spectral accuracy and efficiency but also captures unbound – or ‘free-space’ – BCs for the potential.

The Hockney-Eastwood method, which is second-order accurate (Hockney & Eastwood 1981), meets most of the above constraints, except for spectral accuracy. Another method, proposed by Vico et al. (2016, hereafter VGF), which is widely known in plasma physics and condensed matter, involves modifying the Green’s function to account for the unbound nature of the potential outside the numerical domain. However, this method appears to have been overlooked by the astrophysics community. Most importantly, under the VGF formalism the Green’s function in Fourier space has an analytical form and is regularized at the singularity, making it suitable for FFT methods. Benchmarks have shown that the relative errors of this method already reach machine accuracy (i.e. for uniform grids) at a modest number of points, outperforming any fixed order of convergence algorithm (Zou et al. 2021; Mayani et al. 2024). However, the VGF method has not yet been adapted to the shearing box approximation (Goldreich & Lynden-Bell 1965), that is, for the case where the vertical structure has to be included (but see Koyama & Ostriker 2009, for a complementary approach). In our case, the box features two periodic dimensions (i.e. in the horizontal plane), as well as one free-space BC (i.e. in the vertical dimension).

In this short article, we present two new spectral methods for the Poisson equation intended to be used in Cartesian geometry. We demonstrate convergence and performance for our reference implementation within the finite-volume code NIRVANA-III. Specifically, we detail the derivation and implementation of a new analytical Green’s function meant for 3D shearing box simulations. We describe our new spectral solver and its numerical specifics in Sect. 2. In Sect. 3, we benchmark the accuracy of our solver with one dynamic and two static tests. Finally, in Sect. 4 we quantify the performance of our solver, employing distributed-memory parallelism with standard message passing interface (MPI) routines.

2 General description

In the following, we discuss some basic considerations that provide the foundations for our spectral Poisson solver. We then introduce two methods, which to our knowledge are new, for studying SG in stratified discs.

2.1 Mapping to fully periodic points

The methods discussed in this paper require periodicity in both the x- and y-directions. In the shearing box, where the x boundaries are permanently in motion, strict periodicity is only satisfied in the y-direction, and during special points in time for the x-direction. To address this, one first needs to perform a coordinate transformation along the y-direction (Gammie 2001), using either a simple linear interpolation (Gressel & Ziegler 2007) or the Shear Advection by Fourier Interpolation (SAFI) scheme (Johansen et al. 2009). It is important to note that this mapping introduces a coordinate transformation that must be accounted for in the Poisson equation (Zier & Springel 2023), modifying the wave-vector in the periodic frame as follows: $$\mathit{k}(t)=\left(\begin{array}{c}{k}_x\\ k_y\\ k_z\end{array}\right)=\left(\begin{array}{c}{K}_x+\frac{\mathrm{\Delta }{y}_0(t)}{L_x}{K}_y\\ K_y\\ K_z\end{array}\right),$$(3)

where $\mathrm{\Delta }{y}_0(t)=mod(\frac{3}{2}\mathrm{\Omega }{L}_{x}t,L_y)$, and K is the wave-vector in the initial frame. We refer only to the periodic frame and its associated wave-vector k, until the end of the paper. Since Δy₀/L_x periodically varies between 0 and L_y/L_x, the wave-vector space is a rectangular cuboid whose x length varies continuously between K_x and $K_x+\frac{L_y}{L_x}{K}_y$. Second, the Poisson equation is then solved in this new, fully-periodic (i.e. in the horizontal plane) reference frame. Finally, the obtained potential is subsequently mapped back to the original coordinate system, via reverse mapping.

In the following sections, we focus primarily on the second step of the previously outlined scheme: solving the Poisson equation with two periodic BCs, and a vacuum BC in the vertical direction. Henceforth, any reference to wave-vectors will pertain to k. We now present our two techniques, based on the spectral approach to solving Poisson’s equation.

2.2 The Superposition Analytical-Spectral Hybrid Approach (SASHA)

In this section, we aim to solve Eq. (1), subject to mixed BCs. By invoking the principle of superposition, we can decompose Poisson’s equation into a constant term and a deviation from this zeroth-order term, $${\mathrm{\nabla }}^2\mathrm{\Phi }=4\pi G⟨\rho ⟩+4\pi G(\rho -⟨\rho ⟩),$$(4)

where ⟨ρ⟩ is the volume-averaged mass density. First, we solve ∇²Φ₀ = 4πG⟨ρ⟩, which can be integrated subject to the BCs and symmetry requirements. As can be checked, a straightforward solution is $${\mathrm{\Phi }}_0(z)=4\pi G⟨\rho ⟩z(\frac{1}{2}z-z_0),$$(5)

where z₀ = ⟨zρ⟩/⟨ρ⟩ is the vertical coordinate of the centre of mass. In the second step, we utilize the numerically efficient spectral method, as presented in Eq. (2), to solve for ∇²Φ_P = 4πG (ρ − ⟨ρ⟩). Specifically, DFTs require periodic BCs in all directions, which is only satisfied in the x and y directions. Therefore, we use an aperiodic convolution in the vertical direction, which simply requires us to artificially make the domain periodic in the vertical direction using the zero-padding (ZP) technique (Hockney & Eastwood 1981). This technique consists of doubling the input signal with zeros, not allowing periodic boxes in the vertical direction to affect the computational domain, and improving the frequency resolution in Fourier space. As we see in Sect. 3, this does not alter the solution. We stress that the requirement of periodic BCs in all directions implies that the volume integral of the left-hand side of the Poisson equation must vanish: $$\iiint \mathbf{\nabla }\cdot \mathbf{\nabla }{\mathrm{\Phi }}_P{d}^3\mathit{r}=\iint \frac{\mathrm{\partial }{\mathrm{\Phi }}_P}{\mathrm{\partial }\mathit{n}}{d}^{2}S=0.$$(6)

Therefore, the subtraction of the mean density from the source term for this subproblem ensures mathematical consistency (Binney & Tremaine 2008; Mandal et al. 2023). In the final step, the full gravitational potential in the shearing box is obtained as the sum of the two aforementioned potentials, Φ = Φ₀ + Φ_P, and is consistent with the treatment of gravity.

2.3 The Vico–Greengard–Ferrando with hybrid boundary conditions (VGF-HybridBC) method

Before discussing our second method, it is instructive to study the one-dimensional case first. Here, the mass density, ρ(z), is assumed to be constant along each (infinite) horizontal slab and therefore only depends on the vertical coordinate, z. In this case, the Laplacian becomes ${\mathrm{\nabla }}^2=\frac{{\mathrm{d}}^2}{\mathrm{d}{z}^2}$, and the partial differential equation (PDE) is transformed into a simple ordinary differential equation (ODE), $$\frac{{\mathrm{d}}^2\mathrm{\Phi }}{\mathrm{d}{z}^2}=4\pi G\rho (z),$$(7)

which can be integrated using a convolution approach. This procedure employs the inhomogeneous 1D Green’s function $${\mathcal{G}}_{1D}(z,z^{\mathrm{\prime }})=\frac{1}{2}|z-z^{\mathrm{\prime }}|,$$(8)

which produces the solution Φ(z) = 4πG ∫ 𝒢_1D(z, z′)ρ(z′) dz′. This procedure yields the correct symmetry and asymptotic behaviour, contrary to methods using the homogeneous Green’s function, which require the potential to be augmented by a linear function in z. Indeed, the force norm must become symmetric when extending far above and below the midplane¹. Specifically, in the asymptotic limit (i.e. at great height above or below the midplane), the potential, Φ(z), becomes a linear function, resulting in constant gravitational acceleration. This is in contrast to the 3D case, where the acceleration drops with distance. The explanation is in the assumed infinite extent of each slab of material; the further one departs from the midplane, the larger the part of the density distribution that can be ‘seen’ (under a constant opening angle) at any given location. Since the seen area grows quadratic with height, the quadratic dependence on separation – the very hallmark of Newtonian gravity – gets entirely compensated in that case.

Extending this approach, we now introduce a Fourier-based convolution method that is inherently three-dimensional and fully compatible with vertical vacuum BCs. We aim to derive the appropriate Green’s function for the Poisson equation, adapting the VGF method by Vico et al. (2016) to a scenario with two periodic BCs and one vacuum BC. Given the periodic nature of the problem in the x and y directions, we apply a Fourier transform in these dimensions, resulting in a 1D Helmholtz equation (or screened Poisson equation) in the z-direction: $$(\frac{{\mathrm{d}}^2}{\mathrm{d}{z}^2}-k^2)\tilde{\mathrm{\Phi }}(z)=4\pi G\tilde{\rho }(z),$$(9)

where $k^2=k_x^2+k_y^2$, and $\tilde{\mathrm{\Phi }}$ and $\tilde{\rho }$ represent the partial Fourier transforms of the potential and density in the x and y directions, respectively. The associated Green’s function with the above linear differential equation is: $${\mathcal{G}}_k(z-z^{\mathrm{\prime }})=\{\begin{array}{ccc}\frac{1}{2}|z-z^{\mathrm{\prime }}|& \text{\hspace{0.17em}if\hspace{0.17em}}& k=0\\ -\frac{1}{2k}{e}^{-k|z-z^{\mathrm{\prime }}|}& \text{\hspace{0.17em}if\hspace{0.17em}}& k\ne 0.\end{array}$$(10)

We observe a discontinuity in the Green’s function as k → 0. Indeed, for k = 0, the Green’s function diverges linearly at infinity, indicating a net monopole (uncancelled zero-mode) and reflecting a constant gravitational field (see previous paragraph). Conversely, the k ≠ 0 Green’s function mitigates this divergence by introducing screening. Both cases combined permit us to satisfy physical BCs at infinity (in the vertical direction) for a stratified disc. This Green’s function enables us to obtain the potential via convolution: $$\tilde{\mathrm{\Phi }}(z)=4\pi G{\int }_{-L_z/2}^{L_z/2}{\mathcal{G}}_k(z-z^{\mathrm{\prime }})\tilde{\rho }\left(z^{\mathrm{\prime }}\right)\mathrm{d}{z}^{\mathrm{\prime }}.$$(11)

The first trick of the VGF method involves recognizing that the integral form for the potential remains unchanged in a finite-sized box, if the Green’s function is replaced by: $${\mathcal{G}}^L(z-z^{\mathrm{\prime }})=\mathrm{rect}\left(\frac{z-z^{\mathrm{\prime }}}{2L}\right){\mathcal{G}}_k(z-z^{\mathrm{\prime }}),$$(12)

where L = αL_z is a suitable enclosure (choosing α > 1), L_z is the vertical extent of the actual numerical box, and ‘rect’ is the rectangular function. For the rest of the paper we adopt a value α = 1.1 (Mayani et al. 2024).

The second trick involves computing the Fourier transform of the above Green’s function in the vertical direction $${\hat{\mathcal{G}}}^L(k,k_z)={\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}{\mathcal{G}}^L(z)e^{i{k}_{z}z}\mathrm{d}z,$$(13)

where three cases must be distinguished (see Appendix A): $$=\{\begin{array}{ll}{L}^2/2& \text{\hspace{0.17em}if\hspace{0.17em}}k=0,k_z=0,\\ \frac{[k_{z}L\sin \left(k_{z}L\right)+\cos \left(k_{z}L\right)-1]}{k_z^2}& \text{\hspace{0.17em}if\hspace{0.17em}}k=0,k_z\ne 0,\\ -\frac{e^{-kL}(\frac{k_z}{k}\sin \left(k_{z}L\right)-\cos \left(k_{z}L\right))+1}{k^2+k_z^2}& \text{\hspace{0.17em}otherwise.\hspace{0.17em}}\end{array}$$(14)

Finally, the potential in real space is simply obtained by an inverse Fourier transform of the product between the Fourier-transformed Green’s function and density: $$\mathrm{\Phi }(x,y,z)=4\pi G{\mathcal{F}}^{-1}({\hat{G}}^L\cdot \hat{\rho }),$$(15)

where $\hat{\rho }$ is the full Fourier transform of the density. It is worth noting that our calculation naturally resolves the Green’s function singularity for (k, k_z) = (0, 0).

2.4 Implementation details

It is customary to use the FFTW3 library (Frigo & Johnson 2005) for spectral methods; however, it only allows parallelization in one direction, affecting globally the performance of the hydrodynamical solver of NIRVANA-III. Indeed, the granularity of the hydrodynamical solver is optimal for a 3D parallel decomposition. For this reason, we opted for the P3DFFT library (Pekurovsky 2012), which uses a ‘pencil’ decomposition (i.e. along two out of three spatial dimensions, leaving the third local in memory). This is a good compromise, which balances the performance of the hydro module and the spectral gravitational solver of NIRVANA-III. The P3DFFT version uses stride-1 decomposition, transforming (Z, X, Y) arrays directly into (Y, X, Z) complex arrays – eliminating MPI_ALL_TO_ALL calls and transpositions, resulting in significant time and memory savings. In this context, we note that Koyama & Ostriker (2009) presented an FFT method applicable to shearing boxes. However, it involves four steps: (1) in-plane Fourier transform, (2) modification of the density in the mirror regions and vertical FFT, (3) convolution, and (4) backward transform. This particular procedure, specifically the decomposition of the 3D transform in two steps, precludes the use of parallel three-dimensional FFT, hindering the use of P3DFFT and therefore limiting performance. The advantage of our method is that the Fourier transform in the three directions is done in a single step.

The SASHA method for linear decomposition, discussed in Sect. 2.2, employs a standard ZP technique due to the vacuum BCs in the vertical direction. Although ZP (applied along one dimension only) doubles the memory usage, the Hermitian property of the Fourier transforms of real quantities reduces the number of complex elements by half, compensating for this effect. Additionally, most FFT algorithms are optimized for sizes that are multiples of two. The VGF-HybridBC method, highlighted in Sect. 2.3, also necessitates padding, but requires quadrupling the domain in the vertical direction. This is due to the handling of the aperiodic convolution in the vertical direction and the oscillatory nature of 𝒢^L (Vico et al. 2016; Mayani et al. 2024). To address this issue, Vico et al. (2016) proposed a method to reduce the grid size from quadruple to double by using a pre-computation step. However, this approach is not feasible in our specific case, due to the time-dependent nature of the Green’s function making any pre-computation impractical. We note that ZP can, in principle, be combined with so-called ‘pruned transforms’, which make the FFT library aware of the regions filled with zeros, resulting in memory savings and avoiding certain unnecessary communications between processors. However, in P3DFFT, this filtering is only possible for the wave-spectrum data and not the real data. Indeed, we speculate that the pruned transform module of P3DFFT was likely created with hydrodynamical turbulence applications in mind, where the two-thirds filtering rule is used to avoid aliasing.

For production runs, a crucial consideration is determining the potential value in ghost cells in order to compute forces. Initially, it was tempting to use the gravitational potential from the zero-padded region as the buffer zone. However, these values lack physical meaning and cannot be related to actual data. Therefore, we decided to ensure that the gradient of the full potential, Φ, in ghost cells matched that of the last active cell in the vertical direction.

3 Tests and numerical convergence

In this section, we benchmark the accuracy of our method and check our numerical implementation with tests, two static and one dynamic, respectively. We start by retrieving the potential associated with a Gaussian vertical distribution.

3.1 1D vertical test

Through this test, we aim to demonstrate that both our methods, coupled with the ZP technique, provide correct and accurate results in a setup where vertical vacuum BCs are considered. Therefore, we initiated the density profile as a shifted Gaussian, $$\rho (z)=\{\begin{array}{ll}{e}^{-\frac{1}{2}{(z-z_0)}^2}& \text{\hspace{0.17em}if\hspace{0.17em}}\\ 0& \text{\hspace{0.17em}else\hspace{0.17em}}\end{array}z\in [-L_z/2,L_z/2],$$(16)

where the extent of the numerical window in the vertical direction is [−L_z/2, L_z/2]. Here, it is particularly important to also solve the Poisson equation in the vacuum regions by ensuring the continuity of the potential and its derivative at the vertical boundaries of the box, and ensuring that the force associated with the potential is symmetric with respect to z far away from the midplane, that is, dΦ/dz(−∞) = −dΦ/dz (+∞). With these considerations in mind, the gravitational potential in the numerical window (i.e. [−L_z/2, L_z/2]) is: $$\frac{\mathrm{\Phi }(z)}{4\pi G}=\sqrt{\frac{\pi }{2}}(z-z_0)\mathrm{erf}\left(\frac{z-z_0}{\sqrt{2}}\right)+e^{-\frac{1}{2}{(z-z_0)}^2}+az,$$(17)

with $$a=-\frac{1}{2}\sqrt{\frac{\pi }{2}}[\mathrm{erf}\left(\frac{-L_z/2-z_0}{\sqrt{2}}\right)+\mathrm{erf}\left(\frac{L_z/2-z_0}{\sqrt{2}}\right)],$$(18)

and where we took the integration constant equal to zero.

For this test, we chose the fiducial parameters L_x = L_y = L_z = 6, as well as N_x = N_y = N_z = 128, and z₀ = [0, 2.0]. The choice of a uniform cubic grid, although not necessary, allowed us to infer the limit of optimal accuracy. Similarly, the offset parameter, z₀, allowed us to compare a symmetric and an off-centre density distribution, respectively.

In the top panel of Fig. 1, we show the analytic solution for the gravitational potential and compare it with the numerical results provided by both our spectral methods. For the sake of comparison, in these graphics, we subtracted a constant such that all potentials have a minimum equal to 1. The analytical prediction and numerical estimates for the centred and shifted distributions are indistinguishable. More specifically, in the bottom panel of Fig. 1 we depicted the relative error, $$ϵ=\left|\frac{{\mathrm{\Phi }}_{num}-{\mathrm{\Phi }}_a}{{\mathrm{\Phi }}_a}\right|,$$(19)

which ranges between 10⁻⁶ and 10⁻⁴ for the shifted distribution. For the centred distribution, the numerical and analytical potential are identical, except in four points close to the midplane where the relative error is only about 10⁻⁷.

Fig. 1 Potential associated with a one-dimensional vertical Gaussian distribution. We considered two cases: (i) when centred around the midplane, and (ii) with a shifted Gaussian profile.

3.2 3D test with mixed boundary conditions

A more challenging test consists of benchmarking our solver in a setup with two periodic BCs and one unbound direction. Therefore, we chose a density distribution (defined on z ∈ [−L_z/2, L_z/2]), $$\rho (x,y,z)=\cos \left(k_{x}x\right)\cos \left(k_{y}y\right)e^{-\frac{1}{2}(z/H{)}^2},$$(20)

whose associated potential (applicable in the same domain) is: $$\mathrm{\Phi }(x,y,z)=4\pi G\cos \left(k_{x}x\right)\cos \left(k_{y}y\right)[2a\cosh (kz)+F_{+}+F_{-}],$$(21)

with $k=\sqrt{k_x^2+k_y^2}$, as well as functions $$\begin{array}{r}{F}_{\pm }(z)=\sqrt{\frac{\pi }{2}}\frac{H}{2k}{e}^{\frac{kH}{2}(kH\pm 2z/H)}\mathrm{erf}\left(\frac{kH\pm z/H}{\sqrt{2}}\right)\\ a=-\sqrt{\frac{\pi }{2}}\frac{H}{2k}{e}^{\frac{(kH{)}^2}{2}}\mathrm{erf}\left(\frac{kH+L_z/(2H)}{\sqrt{2}}\right).\end{array}$$(22)

As in Sect. 3.1, we also obtained this solution by solving the Poisson equation in the vacuum regions and enforcing the continuity of both the potential and its derivative at the domain boundaries. However, unlike in Sect. 3.1, we additionally imposed the condition that the potential vanishes at large distances from the midplane (see Sect. 2.3).

For this study, we set N_x = N_y = N_z = 128, as well as L_x = L_y = L_z = 6, and chose the scale height a) three times smaller than the vertical extent of the box, and b) as large as the latter – in other words H ∈ [1.0, 3.0]. In addition, we set k_x = k_y = 2π/L_x.

In Fig. 2, we display the gravitational potential, for x = y = 0, as a function of z in the top panel and the associated relative error in the bottom panel. As previously, we subtracted a constant from all potentials to ensure that their minima were equal to 1. For both configurations, the relative error of the VGF-HybridBC approach reached its maximum at the vertical boundaries. The upper limits on the error are about 10⁻⁶ and 10⁻⁴, for H = 1.0 and H = 3.0, respectively. Note that the error vanishes exactly in the symmetry plane (i.e. at z = 0). In comparison, the SASHA approach is slightly less accurate (by a factor of a few in terms of relative error) than the VGF-HybridBC method revisited.

For completeness, in Fig. 3 we show the colour-coded relative error in the z = 0 plane (top) and the y = 3 plane (bottom), respectively. The overall error was less than 10⁻⁵ (with the maximum reached near the vacuum boundaries), including larger regions where the error was exactly zero (i.e. below machine accuracy).

Fig. 2 Potential along the z direction associated with a 3D density distribution: periodic in x and y, Gaussian in the vertical direction.

3.3 Numerical convergence

For the numerical convergence test, we conducted 1D and 3D tests as outlined in Sects. 3.1 and 3.2, varying the number of cells from 16³ to 256³. Figure 4 displays the L₂ norm of the relative error, $$\parallel ϵ{\parallel }_2=\frac{1}{N}\sqrt{\sum _{i,j,k}{ϵ}_{i,j,k}^2}$$(23)

for both 1D (top panel) and 3D (bottom panel) cases. We begin by analysing the outcomes of the SASHA method. In the 1D test, errors exhibit third-order convergence, with the highest error observed for the shifted distribution. Conversely, in the 3D test, this method demonstrates second-order convergence. Moving on to the VGF-HybridBC method, we find that in both scenarios, errors converge at a third order. The worst-case error is approximately 10⁻⁵, with the best-case error reaching 10⁻⁷ for 16³ cells. For 256³ cells, the error for all tests is consistently below 10⁻⁹ and reaches as low as 10⁻¹¹ for the 3D test with H = 1.

3.4 Jeans instability

This 1D dynamical test consists of triggering gravitational instability of a uniform density distribution at equilibrium, and is particularly tailored to solving the Poisson equation with periodic BC. Therefore, we consider a 1D box with 128 cells is the x direction and consider a uniform density distribution, augmented by a small sinusoidal perturbation: $$\rho (x)={\rho }_0(1+A\cos \left(\frac{2\pi x}{\lambda }\right))$$(24)

with A = 0.01. This perturbation can lead to either oscillating short-wavelength perturbations or to collapsing long-wavelength modes (Hubber et al. 2006; Mandal et al. 2023), whose periods and time, respectively, are $$T_{\mathrm{o}\mathrm{s}\mathrm{c}}={\left(\frac{\pi }{G{\rho }_0}\right)}^{1/2}\frac{\lambda }{{({\lambda }_J^2-{\lambda }^2)}^{1/2}};T_{\mathrm{c}\mathrm{o}\mathrm{l}\mathrm{l}}={\left(\frac{1}{4\pi G{\rho }_0}\right)}^{1/2}\frac{\lambda }{{({\lambda }^2-{\lambda }_J^2)}^{1/2}},$$(25)

where ${\lambda }_J=\sqrt{\frac{\pi c_s^2}{G{\rho }_0}}$ is the Jeans length. Specifically, the collapse time is the time it takes for the peak density to grow from Aρ₀ to Aρ₀ cosh (1). As in Mandal et al. (2023), we fix λ = 2 and vary the ratio λ/λ_J, which is equivalent to changing the sound speed. This choice makes it possible to have one perturbation normal mode in the initial state for all setups.

In Fig. 5, we depict the normalised time (i.e. with respect to the free fall time) with respect to the normalised wavelength. The star and circles correspond to the fitted numerical points in the oscillatory and collapsing regimes, respectively, while the solid lines are the theoretical estimates. We observe that the error between the theoretical prediction and the numerical estimate is smaller than 0.3% in the oscillatory regime, and about 2% in the collapsing regime, which validates this test.

Fig. 3 Log-scale relative error cuts for the 3D test, employing the VGF-HybridBC method, with mixed BCs. Top: z = 0 plane. Bottom: y = 3 plane.

Fig. 4 Convergence test: 1D (top) and 3D (bottom) results. The SASHA technique achieves second-order convergence, while the VGF-HybridBC method reaches slightly better than third-order convergence.

4 Performance estimates

We benchmarked the real-life performance of our algorithm on the Romeo cluster at the HPC Centre at the TU Dresden². Specifically, that machine is a general-purpose NEC cluster featuring 188 nodes, each equipped with dual AMD EPYC 7702 processors (64 cores @ 2.0 GHz, multithreading enabled), 512 GB of DDR4-3200 RAM, and 200 GB of local SSD storage. We employed the exclusive option of SLURM to mitigate interference with other jobs.

For the presented weak scaling test, each MPI task handled a fixed workload of 32³ cells. We emphasize that for a fixed number of cores, P3DFFT performs optimally for a processor grid configuration satisfying P_y ≥ P_x, where P_x and P_y denote the number of processors in the x and y directions, respectively. Additionally, P_y should be set to the maximum number of logical cores per node, which is 128 for the ROMEO cluster. In contrast, the MHD solver achieves optimal performance with a grid decomposition that minimizes the surface-to-volume ratio. For our setup, this implies that P_x should be as close as possible to P_y. To reconcile these conflicting requirements, we tested multiple processor grid decompositions for a fixed number of cores, adopting the configuration with the minimal runtime as the baseline for each core count.

In Fig. 6, we present the scaling of the VGF-HybridBC method, implemented in the finite volume code nirvana-iii, from 128 to 4096 processors. For completeness, we compared the timings of our spectral solver with those of the orbital advection, mesh synchronization, and MHD solver. We found that the time spent in a) the interpolation to the nearest periodic point, and b) the spectral solver, increases with the total number of cell points, but its share in the total computation time remains almost constant and below 6%.

We emphasize that the timings of our solver are problem-independent: spectral methods do not require iterations that depend on the problem, and the number of operations is fixed for a given grid size, regardless of the density distribution. Therefore we expect to maintain the obtained performance in any kind of production runs.

Fig. 5 Left: characteristic times of oscillatory and collapse modes with respect to wavelength. Right: relative error.

Fig. 6 Weak scaling study for a workload of 323 cells per MPI task and using pencil decomposition (P3DFFT). The time spent in the spectral solver is less than 6% of the whole runtime.

5 Conclusion

We have presented two novel spectral methods for solving the Poisson equation, tailored for Cartesian shearing box simulations with vertical vacuum BCs. Our first method, SASHA, combines an analytical solution with a spectral solution via the superposition principle. Our second approach, VGF-HybridBC, leverages the free-space nature of the Green’s function, yielding an analytical and regularized expression in Fourier space (Eq. (14)). This last approach allows the gravitational potential to be evaluated in a single step as a 3D convolution in Fourier space. This enables the use of fast Fourier algorithms, achieving high performance and accuracy.

We implemented these new methods in the Finite Volume code NIRVANA-III. To optimize performance, we integrated both methods with the P3DFFT library. This allows for pencil decomposition, balancing the computational load between the hydrodynamic solver and the spectral solver. Our results demonstrate excellent accuracy, with a relative error of at least 10⁻⁵ using only 16³ grid cells. The methods, SASHA and VGF-HybridBC, demonstrate second- and third-order convergence, respectively. In a standard production run using orbital advection, our spectral solver accounts for less than 6% of the total runtime, even when scaled up to 4096 processors.

Acknowledgements

This work was supported by the European Union (ERC-CoG, EPOCH-OF-TAURUS, No. 101043302). Views and opinions expressed are however those of the author(s) only and do not necessarily reflect those of the European Union or the European Research Council. Neither the European Union nor the granting authority can be held responsible for them. The authors gratefully acknowledge the computing time made available to them on the high-performance computer at the NHR Center of TU Dresden. This center is jointly supported by the Federal Ministry of Education and Research and the state governments participating in the NHR (www.nhr-verein.de/unsere-partner).

References

Appendix A Green’s function in Fourier space

In this section we detail the derivation of the Fourier transform of the Green’s function associated with the VGF-HybridBC method: $${\hat{\mathcal{G}}}^L(k,k_z)={\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}{\mathcal{G}}^L(z)e^{i{k}_{z}z}\mathrm{d}z.$$(A.1)

Following the definition of the Green’s function in real space (Eqs. 10 and 12), two cases must be distinguished:

• The case k = 0: $$\begin{array}{rl}{\hat{\mathcal{G}}}^L(k=0,k_z)& ={\int }_{-L}^L\frac{|z|}{2}{e}^{i{k}_{z}z}\mathrm{d}z,\\ & ={\int }_0^{L}z\cos \left(k_{z}z\right)\mathrm{d}z,.\\ & =\frac{k_{z}L\sin \left(k_{z}L\right)+\cos \left(k_{z}L\right)-1}{k_z^2}\end{array}$$(A.2)

  Expanding above expression to second order in k_z and taking the limit k_z → 0, we find $${\hat{\mathcal{G}}}^L(k=0,k_z=0)=\frac{L^2}{2}.$$(A.3)

• The case k ≠ 0: $$\begin{array}{rl}{\hat{\mathcal{G}}}^L(k=0,k_z)& =-\frac{1}{2k}{\int }_{-L}^L{e}^{-k|z|}{e}^{i{k}_{z}z}\mathrm{d}z,\\ & =-\frac{1}{k}{\int }_0^L\cos \left(k_{z}z\right)e^{-kz}\mathrm{d}z,\\ & =-\frac{e^{-kL}(\frac{k_z}{k}\sin \left(k_{z}L\right)-\cos \left(k_{z}L\right))+1}{k^2+k_z^2}.\end{array}$$(A.4)

All Figures

	Fig. 1 Potential associated with a one-dimensional vertical Gaussian distribution. We considered two cases: (i) when centred around the midplane, and (ii) with a shifted Gaussian profile.
In the text

	Fig. 2 Potential along the z direction associated with a 3D density distribution: periodic in x and y, Gaussian in the vertical direction.
In the text

	Fig. 3 Log-scale relative error cuts for the 3D test, employing the VGF-HybridBC method, with mixed BCs. Top: z = 0 plane. Bottom: y = 3 plane.
In the text

	Fig. 4 Convergence test: 1D (top) and 3D (bottom) results. The SASHA technique achieves second-order convergence, while the VGF-HybridBC method reaches slightly better than third-order convergence.
In the text

	Fig. 5 Left: characteristic times of oscillatory and collapse modes with respect to wavelength. Right: relative error.
In the text

	Fig. 6 Weak scaling study for a workload of 323 cells per MPI task and using pencil decomposition (P3DFFT). The time spent in the spectral solver is less than 6% of the whole runtime.
In the text