35.3 Gravity Test Problems

35.3.1 Jeans Instability

The linear instability of self-gravitating fluids was first explored by Jeans (1902) in connection with the problem of star formation. The nonlinear phase of the instability is currently of great astrophysical interest, but the linear instability still provides a very useful test of the coupling of gravity to hydrodynamics in FLASH.

The Jeans problem allows one to examine the behavior of sinusoidal, adiabatic density perturbations in both the pressure-dominated and gravity-dominated limits. This problem uses periodic boundary conditions. The equation of state is that of a perfect gas. The initial conditions at $t=0$ are

$$\rho({\bf x}) = \rho_0\left[ 1 + \delta {\rm cos}({\bf k}\cdot{\bf x})\right]$$

$$p({\bf x}) = p_0\left[ 1 + \gamma\delta {\rm cos}({\bf k}\cdot{\bf x})\right]$$ (35.30)

$${\bf v}({\bf x}) = 0,$$

where the perturbation amplitude $\delta\ll 1$. The stability of the perturbation is determined by the relationship between the wavenumber $k\equiv \vert{\bf k}\vert$ and the Jeans wavenumber $k_J$, where $k_J$ is given by

$$k_J \equiv {\sqrt{4\pi G\rho_0}\over c_0} ,$$ (35.31)

and where $c_0$ is the unperturbed adiabatic sound speed

$$c_0 = \sqrt{\gamma p_0\over\rho_0}$$ (35.32)

(Chandrasekhar 1961). If $k > k_J$, the perturbation is stable and oscillates with frequency

$$\omega = \sqrt{c_0^2k^2 - 4\pi G\rho_0} ;$$ (35.33)

otherwise, it grows exponentially, with a characteristic timescale given by $\tau = (i\omega)^{-1}$.

We checked the dispersion relation (35.33) for stable perturbations with $\gamma=5/3$ by fixing $\rho_0$ and $p_0$ and performing several runs with different $k$. We followed each case for roughly five oscillation periods using a uniform grid in the box $[0,L]^2$. We used $\rho_0 = 1.5\times10^7$ g cm$^{-3}$ and $p_0 = 1.5\times10^7$ dyn cm$^{-2}$, yielding $k_J = 2.747$ cm$^{-1}$. The perturbation amplitude $\delta$ was fixed at $10^{-3}$. The box size $L$ is chosen so that $k_J$ is smaller than the smallest nonzero wavenumber that can be resolved on the grid

$$L = {1\over 2}\sqrt{\pi\gamma p_0\over G\rho_0^2} .$$ (35.34)

This prevents roundoff errors at wavenumbers less than $k_J$ from being amplified by the physical Jeans instability. We used wavevectors ${\bf k}$ parallel to and at 45 degrees to the $x$-axis. Each test calculation used the multigrid Poisson solver together with its default settings.

The resulting kinetic, thermal, and potential energies as functions of time for one choice of ${\bf k}$ are shown in Figure 35.58 together with the analytic solution, which is given in two dimensions by

$$T(t) = {\rho_0\delta^2\vert\omega\vert^2L^2\over 8k^2}\left[1-{\rm cos}(2\omega t) \right]$$

$$U(t)-U(0) = -{1\over 8}\rho_0c_0^2\delta^2L^2\left[1-{\rm cos}(2\omega t) \right]$$ (35.35)

$$W(t) = -{\pi G\rho_0^2\delta^2L^2\over 2k^2}\left[1+{\rm cos}(2\omega t) \right] .$$

The figure shows that FLASH obtains the correct amplitude and frequency of oscillation. We computed the average oscillation frequency for each run by measuring the time interval required for the kinetic energy to undergo exactly ten oscillations. Figure 35.59 compares the resulting dispersion relation to (35.33). It can be seen from this plot that FLASH correctly reproduces (35.33). At the highest wave number ($k = 100$), each wavelength is resolved using only about 14 cells on a six-level uniform grid, and the average timestep (which depends on $c_0$, $\Delta x$, and $\Delta y$, and has nothing to do with $k$) turns out to be comparable to the oscillation period. Hence the frequency determined from the numerical solution for this value of $k$ is somewhat more poorly determined than for the other runs. At lower wavenumbers, however, the frequencies are correct to less than 1%.

Figure: Kinetic, internal, and potential energy versus time for a stable Jeans mode with $k=10.984$. Points indicate numerical values found using FLASH 3.0 with a fixed four-level adaptive grid. The analytic solution for each form of energy is shown using a solid line.

Figure 35.59: Computed versus expected Jeans dispersion relation (for stable modes) found using FLASH 1.62 with a six-level uniform grid.

Table 35.10: Runtime parameters used with the Jeans test problem.

Variable	Type	Default	Description
rho0	real	$1.5\times 10^{7}$	Initial unperturbed density ($\rho_0$)
p0	real	$1.5\times 10^{7}$	Initial unperturbed pressure ($p_0$)
amplitude	real	0.001	Perturbation amplitude ($\delta$)
lambdax	real	0.572055	Perturbation wavelength in $x$ direction ( $\lambda_x = 2\pi/k_x$)
lambday	real	$1.0\times 10^{10}$	Perturbation wavelength in $y$ direction ( $\lambda_y = 2\pi/k_y$)
lambdaz	real	$1.0\times 10^{10}$	Perturbation wavelength in $z$ direction ( $\lambda_z = 2\pi/k_z$)
delta_ref	real	0.01	Refine a block if the maximum density contrast relative to $\rho_{\rm ref}$ is greater than this
delta_deref	real	-0.01	Derefine a block if the maximum density contrast relative to $\rho_{\rm ref}$ is less than this
reference_density	real	$1.5\times 10^{7}$	Reference density for grid refinement ( $\rho_{\rm ref}$). Density contrast is used to determine which blocks to refine; it is defined as $$\max_{\rm block}\left\{\left\vert{\rho_{ijk}\over\rho_{\rm ref}}-1\right\vert\right\}$$

The additional runtime parameters supplied with the Jeans problem are listed in Table 35.10. This problem is configured to use the perfect-gas equation of state (gamma) with gamma set to 1.67 and is run in a two-dimensional unit box. The refinement marking routine (Grid_markRefineDerefine.F90) supplied with this problem refines blocks whose mean density exceeds a given threshold. Since the problem is not spherically symmetric, the multigrid Poisson solver should be used.

35.3.2 Homologous Dust Collapse

The homologous dust collapse problem DustCollapse is used to test the ability of the code to solve self-gravitating problems in which the flow geometry is spherical and gas pressure is negligible. The problem was first described by Colgate and White (1966) and has been used by Mönchmeyer and Müller (1989) to test hydrodynamical schemes in curvilinear coordinates. We solve this problem using a 3D Cartesian grid.

The initial conditions consist of a uniform sphere of radius $r_0$ and density $\rho_0$ at rest. The pressure $p_0$ is taken to be constant and very small

$$p_0 \ll {4\pi G\over\gamma}\rho_0^2 r_0^2 .$$ (35.36)

We refer to such a nearly pressureless fluid as `dust'. A perfect-gas equation of state is used, but the value of $\gamma$ is not significant. Outflow boundary conditions are used for the gas, while isolated boundary conditions are used for the gravitational field.

The collapse of the dust sphere is self-similar; the cloud should remain spherical with uniform density as it collapses. The radius of the cloud, $r(t)$, should satisfy

$$\left({8\pi G\over 3}\rho_0\right)^{1/2} t = \left(1-{r(t)\over r... ...ft({r(t)\over r_0}\right)^{1/2} + \sin^{-1}\left(1-{r(t)\over r_0}\right)^{1/2}$$ (35.37)

(Colgate & White 1966). Thus. we expect to test three things with this problem: the ability of the code to maintain spherical symmetry during an implosion (in particular, no block boundary effects should be evident); the ability of the code to keep the density profile constant within the cloud; and the ability of the code to obtain the correct collapse factor. The second of these is particularly difficult, because the edge of the cloud is very sharp and because the Cartesian grid breaks spherical symmetry most dramatically at the center of the cloud, which is where all of the matter ultimately ends up.

Results of a DustCollapse run using FLASH 3.0 appear in Figure 35.60, which shows plots of density and the X component of velocity in menacing color scheme. The values are plotted at the end of the run from an X-Y plane in the center of the physical domain; density is in logarithmic scale. This run used a resolution of $128^3$, and the results were compared against a similar run using FLASH 2.5. We have also included figures from an earlier higher resolution run using FLASH2 which used $4^3$ top-level blocks and seven levels of refinement, for an effective resolution of $2048^3$. In both the runs, the multipole Poisson solver was used with a maximum multipole moment $\ell = 0$. The initial conditions used $\rho_0 = 10^9$ g cm$^{-3}$ and $r_0 = 6.5\times10^8$ cm. In Figure 35.61a, the density, pressure, and velocity are scaled by $2.43\times10^9$ g cm$^{-3}$, $2.08\times10^{17}$ dyn cm$^{-2}$, and $7.30\times10^9$ cm s$^{-1}$, respectively. In Figure 35.61b they are scaled by $1.96\times10^{11}$ g cm$^{-3}$, $2.08\times10^{17}$ dyn cm$^{-2}$, and $2.90\times10^{10}$ cm s$^{-1}$. Note that within the cloud, the profiles are very isotropic, as indicated by the small dispersion in each profile. Significant anisotropy is only present for low-density material flowing in through the Cartesian boundaries. In particular, it is encouraging that the velocity field remains isotropic all the way into the center of the grid; this shows the usefulness of refining spherically symmetric problems near $r=0$. However, as material flows inward past refinement boundaries, small ripples develop in the density profile due to interpolation errors. These remain spherically symmetric but increase in amplitude as they are compressed. Nevertheless, they are still only a few percent in relative magnitude by the second frame. The other numerical effect of note is a slight spreading at the edge of the cloud. This does not appear to worsen significantly with time. If one takes the radius at which the density drops to one-half its central value as the radius of the cloud, then the observed collapse factor agrees with our expectation from (35.37). Overall our results, including the numerical effects, agree well with those of Mönchmeyer and Müller (1989).

This problem is configured to use the perfect-gas equation of state (gamma) with gamma set to 1.67 and is run in a three-dimensional box. The problem uses the specialized refinement marking routine supplied under the Grid interface of Grid_markRefineSpecialized which refines blocks containing the center of the cloud.

Figure 35.60: XY plane of Density (a) and X component of Velocity (b) are shown at the center of the domain for the DustCollapse problem. The velocity is in normal scale, while density is logscale.

[] []

Figure: Density (black), pressure (red), and velocity (blue) profiles in the homologous dust collapse problem at (a) $t=0.0368$ sec and (b) $t=0.0637$ sec. The density, pressure, and velocity are scaled as discussed in the text.

[] []

35.3.3 Huang-Greengard Poisson Test

The PoisTest problem tests the convergence properties of the multigrid Poisson solver on a multidimensional, highly (locally) refined grid. This problem is described by Huang and Greengard (2000). The source function consists of a sum of thirteen two-dimensional Gaussians

$$\rho(x,y) = \sum_{i=1}^{13} e^{-\sigma_i[(x-x_i)^2+(y-y_i)^2]} ,$$ (35.38)

where the constants $\sigma_i$, $x_i$, and $y_i$ are given in Table 35.11. The very large range of widths and ellipticities of these peaks forces the mesh structure to be highly refined in some places. The density field and block structure are shown for a 14-level mesh in Figure 35.62.

Figure: Density field and block structure for a 14-level mesh applied to the Huang-Greengard PoisTest problem. The effective resolution of the mesh is $65,536^2$.

Table 35.11: Constants used in the poistest problem.

$i$	1	2	3	4	5	6	7
$x_i$	0	-1	-1	0.28125	0.5	0.3046875	0.3046875
$y_i$	0	0.09375	1	0.53125	0.53125	0.1875	0.125
$\sigma_i$	0.01	4000	20000	80000	16	360000	400000
$i$	8	9	10	11	12	13
$x_i$	0.375	0.5625	-0.5	-0.125	0.296875	0.5234375
$y_i$	0.15625	-0.125	-0.703125	-0.703125	-0.609375	-0.78125
$\sigma_i$	2000	18200	128	49000	37000	18900

The PoisTest problem uses one additional runtime parameters sim_smlRho, the smallest allowed value of density. Runtime parameters from the Multigrid unit (both Gravity and GridSolvers) are relevant; see Sec:GridSolversMultigrid.

35.3.4 MacLaurin

The gravitational potential at the surface of, and inside a homogeneous spheroid called a “MacLaurin spheroid” is expressible in terms of analytical functions. This handy result was first determined by MacLaurin (1801), and later summarized by, amongst others, Chandrasekhar (1989). These properties allow validation of the FLASH4 gravitational solvers against the analytical solutions.

As a test case, an oblate ( $a_1 = a_2 > a_3$) Maclaurin spheroid, of a constant density $\rho = 1$ in the interior, and $\rho = \epsilon \rightarrow 0$ outside (in FLASH4 smlrho is used). The spheroid is motionless and in hydrostatic equilibrium. The gravitational potential of such object is analytically calculable, and is:

$$\phi({\bf x}) = \pi G \rho \left[ 2A_1 a_1^2 - A_1(x^2+y^2) + A_3(a_3^2-z^2) \right] ,$$ (35.39)

for a point inside the spheroid. Here

$$A_1 = \frac{\sqrt{1-e^2}}{e^3} \sin^{-1}e - \frac{1-e^2}{e^2} ,$$ (35.40)

$$A_3 = \frac{2}{e^2} - \frac{2\sqrt{1-e^2}}{e^3} \sin^{-1}e ,$$ (35.41)

where $e$ is the ellipticity of a spheroid:

$$e = \sqrt{1 - \left( \frac{a_3}{a_1} \right) ^2} .$$ (35.42)

For a point outside the spheroid, potential is:

$$\phi({\bf x}) = \frac{2a_3}{e^2} \pi G \rho \left[a_1 e \tan^{-1}... ...( \tan^{-1}h - \frac{h}{1+h^2} \right) + 2z^2 (h-\tan^{-1}h) \right) \right] ,$$ (35.43)

where

$$h = \frac{a_1 e}{\sqrt{a_3^2 + \lambda}} ,$$ (35.44)

and $\lambda$ is the positive root of the equation

$$\frac{x^2}{a_1^2 + \lambda} + \frac{y^2}{a_2^2 + \lambda} + \frac{z^2}{a_3^2 + \lambda} = 1 .$$ (35.45)

This test is also useful because the spheroid has spherical symmetry in the X-Y plane, but also lack of such symmetry in X-Z and Y-Z planes. The density distribution of the spheroid is shown in Figure 35.3.4. Spherical symmetry is simple to reproduce with a solution using multipole expansion. However, the non-symmetric solution requires an infinite number of multipole moments, while the code calculates solution up to a certain $l_{max}$, specified by the user as runtime parameter mpole_lmax. The error is thus expected to be dominated by the first non-zero term in the remainder of expansion. Also, the solution for any point inside the spheroid is the sum of monopole and dipole moments.

Figure: Density of the MacLaurin spheroid (left X-Y plane, right Y-Z plane) with ellipticity $e=0.9$. The FLASH4 block structure is shown on top.

The simulation is calculated on a MacLaurin spheroid with eccentricity $e=0.9$; several other values for eccentricity were tried with results qualitatively the same. All tests used 3D Cartesian coordinates. The gravitational potential is calculated on an adaptive mesh, and the relative error is investigated:

$$\epsilon = \left\vert \frac{\phi_{\rm analytical} - \phi_{\rm FLASH}}{\phi_{\rm analytical}} \right\vert$$ (35.46)

from zone to zone.

Table 35.12: Minimal and maximal relative error in all zones of the simulation, calculated using (35.46). Last row is approximate time for one full timestep (gravity only).

$l_{max}$	min $(\epsilon)$	max $(\epsilon)$	relative $L^2_{in}$ norm	relative $L^2_{out}$ norm	approx. time $[s]$
0	$4.5 10^{-6}$	$0.182$	$7.1 10^{-2}$	$6.8 10^{-2}$	$9.8$
1	$4.5 10^{-6}$	$0.182$	$7.1 10^{-2}$	$6.8 10^{-2}$	$14.5$
2	$9.8 10^{-6}$	$0.062$	$1.4 10^{-2}$	$1.7 10^{-2}$	$34.7$
4	$1.0 10^{-8}$	$0.026$	$4.0 10^{-3}$	$5.0 10^{-3}$	$55.4$
6	$6.1 10^{-9}$	$0.013$	$1.6 10^{-3}$	$2.5 10^{-3}$	$134.9$
8	$7.8 10^{-9}$	$0.007$	$8.7 10^{-4}$	$1.2 10^{-3}$	$210.2$
10	$6.7 10^{-9}$	$0.004$	$5.5 10^{-4}$	$7.0 10^{-4}$	$609.7$

As expected, increasing spatial resolution improves the solution quality, but here we focus on how the solution depends on the choice of $l_max$, the cutoff $\ell$ in (8.15). In Figure 35.64-35.64 the gravitational potential for the Maclaurin spheroid, the FLASH4 solution, and relative errors for several $l_{max}$'s are shown. A similar figure produced for $l_{max}=1$ shows no difference from Figure 35.64, indicating that the last dipole term in the multipole expansion does not contribute to the accuracy of the solution but does increase computational cost. Because gravity sources are all of the same sign, and the symmetry of the problem, all odd-$l$ moments are zero: reasonable, physically motivated values for setting mpole_lmax should be an even number.

In the X-Y plane, where the solution is radially symmetric, the first monopole term is enough to qualitatively capture the correct potential. As expected, the error is the biggest on the spheroid boundary, and decreases both outwards and inwards. Increasing the maximum included moment reduces errors. However, in other non-symmetric planes, truncating the potential to certain $l_{max}$ leads to an error whose leading term will be the spherical harmonic of order $l_{max}+2$, as can be nicely seen in the lower right sections of Figure 35.64 - 35.66. Increasing $l_{max}$ reduces the error, but also increases the required time for computation. This computational increase is not linear because of the double sum in (8.17). Luckily, convergence is rather fast, and already for $l_{max} = 4$, there are only few zones with relative error bigger than 1%, while for the most of the computational domain the error is several orders of magnitude less.

Figure: Maclaurin spheroid: $l_{max} = 0$, 6 refinement levels. Left column is X-Y plane, cut through z=0.5, right column is Y-Z plane cut through x=0.5 . From top to bottom: analytical solution for the gravitational potential introduced on FLASH grid; solution of FLASH multipole solver; relative error.

Figure: Maclaurin spheroid: $l_{max} = 2$, 6 refinement levels. Left column is X-Y plane, cut through z=0.5, right column is Y-Z plane cut through x=0.5 . From top to bottom: analytical solution for the gravitational potential introduced on FLASH grid; solution of FLASH multipole solver; relative error.

Figure: Maclaurin spheroid: $l_{max} = 10$, 6 refinement levels. Left column is X-Y plane, cut through z=0.5, right column is Y-Z plane cut through x=0.5 . From top to bottom: analytical solution for the gravitational potential introduced on FLASH grid; solution of FLASH multipole solver; relative error.