A comparison between grid and particle methods on the statistics of driven, supersonic, isothermal turbulence

We compare the statistics of driven, supersonic turbulence at a high Mach number using flash, a widely used Eulerian grid-based code, and phantom, a Lagrangian smoothed particle hydrodynamics (SPH) code at resolutions of up to 512³ in both grid cells and SPH particles. We find excellent agreement between codes on the basic statistical properties: a slope of k^-1.95 in the velocity power spectrum for hydrodynamic, Mach 10 turbulence, evidence in both codes for a Kolmogorov-like slope of k^-5/3 in the variable ρ^1/3 v as suggested by Kritsuk et al. and a lognormal probability distribution function (PDF) with a width that scales with the Mach number and proportionality constant b = 0.33 -0.5 in the density variance-Mach number relation. The measured structure function slopes are not converged in either code at 512³ elements. We find that for measuring volumetric statistics such as the power spectrum slope and structure function scaling, SPH and grid codes give roughly comparable results when the number of SPH particles is approximately equal to the number of grid cells. In particular, to accurately measure the power spectrum slope in the inertial range, in the absence of sub-grid turbulence models, requires at least 512³ computational elements in either code. On the other hand the SPH code was found to be better at resolving dense structures, giving maximum densities at a resolution of 128³ particles that were similar to the maximum densities resolved in the grid code at 512³ cells, reflected also in the high density tail of the PDF. We find SPH to be more dissipative at comparable numbers of computational elements in statistics of the velocity field, but correspondingly less dissipative than the grid code in the statistics of density-weighted quantities such as ρ^1/3 v. For SPH simulations of high Mach number turbulence, we find it important to use sufficient non-linear β-viscosity in order to prevent particle interpenetration in shocks (we require β_visc = 4 instead of the widely used default value, β_visc = 2).