Source-linked AI summary
Fundamental differences between SPH and grid methods
Oscar Agertz, Ben Moore, Joachim Stadel, Doug Potter, Francesco Miniati, Justin Read, Lucio Mayer, Artur Gawryszczak, Andrey Kravtsov, Joe Monaghan, Ake Nordlund, Frazer Pearce, Vincent Quilis, Douglas Rudd, Volker Springel, James Stone, Elizabeth Tasker, Romain Teyssier, James Wadsley, Rolf Walder
TL;DR
The paper compares grid and SPH simulations of interacting multiphase fluids, focusing on their ability to resolve astrophysically important instabilities. Grid methods resolve mixing and dynamical instabilities, whereas standard SPH poorly resolves them because density gradients generate erroneous pressure forces and a smoothing-kernel-scale gap. The methods agree during early stripping but diverge after large-scale instabilities develop.
Problem
The study addresses divergent numerical treatment of Kelvin-Helmholtz and Rayleigh-Taylor instabilities, which are important to astrophysical structure formation.
Method
The paper compares grid and SPH simulations of a cold cloud interacting with a hot ambient medium and examines the source of their discrepancies.
Results
Grid methods resolve mixing and dynamical instabilities, while standard SPH poorly resolves them; after large-scale instabilities grow, the grid cloud disrupts but the SPH cloud continues stripping.
Takeaways & Limitations
Standard SPH produces erroneous pressure forces near steep density gradients, creating a smoothing-kernel-scale gap that impedes information transfer.
Takeaways & Limitations
The methods agree during early pressure-driven stripping and diverge only after large-scale instabilities develop; resolution comparisons between grids and SPH are not straightforward.
Abstract
from arXiv · showhide
We have carried out a hydrodynamical code comparison study of interacting multiphase fluids. The two commonly used techniques of grid and smoothed particle hydrodynamics (SPH) show striking differences in their ability to model processes that are fundamentally important across many areas of astrophysics. Whilst Eulerian grid based methods are able to resolve and treat important dynamical instabilities, such as Kelvin-Helmholtz or Rayleigh-Taylor, these processes are poorly or not at all resolved by existing SPH techniques. We show that the reason for this is that SPH, at least in its standard implementation, introduces spurious pressure forces on particles in regions where there are steep density gradients. This results in a boundary gap of the size of the SPH smoothing kernel over which information is not transferred.
1 INTRODUCTION
Interacting-fluid simulations are central to astrophysics because instabilities shape structure formation. The study compares Eulerian grid and SPH methods on problems designed to test their treatment of these gas processes.
- Kelvin-Helmholtz and Rayleigh-Taylor instabilities are fundamental to astrophysical structure formation across systems from proto-planetary disks to galaxies.
- Eulerian grids follow gas through information fluxes between adjacent cells, whereas SPH follows Lagrangian particles whose properties are averaged over nearest neighbours.
- The comparison uses a dense cold cloud moving through a low-density hot medium to capture physical processes relevant to astrophysical structure formation.
- A second configuration studies shear between fluids of different densities to clarify the numerical problems revealed by the cloud test.
- Earlier studies reported divergent gas stripping outcomes for similar galaxy-cluster conditions, with SPH removing half the interstellar medium and a grid calculation removing all of it [Abadi et al. 1999; Quilis et al. 2000].
2 THE BLOB TEST
The blob test places a pressure-balanced cold cloud in a hot supersonic wind to examine shocks, stripping, and instability growth. The flow forms a bow shock, remains subsonic behind it initially, and later accelerates around the cloud.
- The blob test places a spherical cloud in a periodic wind tunnel, with ambient gas ten times hotter and ten times less dense than the cloud, maintaining pressure equilibrium.
- The setup probes ram-pressure stripping and fragmentation driven by Kelvin-Helmholtz and Rayleigh-Taylor instabilities in multiphase flows.
- A supersonic external flow creates a bow shock ahead of the cloud, while post-shock gas is initially subsonic and accelerates to supersonic speed along the cloud’s sides.
3 ANALYTICAL EXPECTATIONS
The analysis estimates cloud-crushing, Kelvin–Helmholtz, and Rayleigh–Taylor timescales for a supersonic cloud and predicts how instability growth depends on scale and acceleration. KH growth is slower than cloud crushing, while small-scale RT modes can grow rapidly at the cloud front.
- 3 ANALYTICAL EXPECTATIONS: The post-shock flow becomes subsonic before the shock disappears, approaches co-motion with the cloud, and produces turbulent boundary-layer mass loss after roughly one crushing time.The analysis also predicts compression along the flow and lateral overspilling caused by pressure differences, producing mass loss independently of instabilities.
- 3.1 The Kelvin-Helmholtz instability: Shorter-wavelength KH modes grow first, but interface broadening shifts the fastest-growing modes toward the interface thickness and ultimately the cloud scale.The cloud-destruction mode has k_cl ∼ 2π/R_cl.
- 3.1 The Kelvin-Helmholtz instability: The KH growth time exceeds cloud crushing time, with the cloud-scale mode expected to reach full growth at approximately τ = 1.6 τ_cr.This characteristic time is subsequently denoted τ_KH.
- 3.1 The Kelvin-Helmholtz instability: Compressibility is omitted from the KH estimate, and gravity, viscosity, magnetic fields, radiation, and related effects can damp the instability in more physical settings.These omissions and modifications limit the scope of the analytic estimate.
- 3.2 The Rayleigh-Taylor instability: For large modes, τ_KH < τ_RT; the largest RT mode grows slowly, whereas small-scale RT growth is expected near the stagnation point, followed later by mixed KH and RT evolution.The RT estimate uses an efficiency factor ϵ = 1, providing a lower limit on τ_RT.
4 NUMERICAL SIMULATIONS
The simulations solve the Euler equations for a perfect gas while isolating hydrodynamic-solver differences by neglecting physical viscosity, radiative processes, and gas self-gravity.
- 4 NUMERICAL SIMULATIONS: The simulations assume a perfect gas and strictly adiabatic evolution away from shocks, with heating occurring through adiabatic compression, expansion, or irreversible shock heating.Physical viscosity, radiative processes, and gas self-gravity are neglected.
4.1 Initial conditions
The blob test places a relaxed, pressure-equilibrated cloud in a periodic wind-tunnel box and uses the same particle-derived initial information for SPH and grid simulations. Initial perturbations arise from particle noise and the imposed velocity start.
- 4.1 Initial conditions: The ambient velocity is added abruptly, and particle noise provides seeds for small-scale RT and Richtmyer–Meshkov instabilities at the cloud surface.A smoother velocity ramp would be more astrophysically faithful but is not used here.
- 4.1 Initial conditions: Grid initial conditions are obtained by smoothing the SPH particle fields with a 32-neighbour kernel and mapping them onto a uniform grid, preserving the particle noise across methods.Resolution and artificial-viscosity strength are the main parameter-study variables.
4.2 The codes
The comparison includes multiple AMR grid and SPH codes that use distinct shock-handling and viscosity formulations. Grid methods employ Godunov-family solvers, while the SPH implementations rely on artificial viscosity and entropy or energy formulations.
- 4.2 The codes: About a dozen independent codes were tested, with consistent results within the grid and SPH groups motivating detailed analysis of selected representative codes.The selected implementations are summarized in Table 1.
- 4.2.1 ART (AMR): ART is an AMR code using a second-order shock-capturing Godunov solver with piecewise-linear reconstruction and shocks resolved within approximately 1–2 cells.It adds a small amount of artificial diffusion to numerical fluxes.
- 4.2.2 CHARM (AMR): CHARM is an AMR code using a higher-order Godunov method with Van Leer reconstruction and a nonlinear Riemann solver, tested here for the influence of initial conditions on cloud evolution.The method is second-order accurate in space and time.
- 4.2.3 Enzo (AMR): Enzo uses an Eulerian AMR PPM solver that is third-order in space and second-order in time, providing accurate shock treatment relative to SPH codes using artificial viscosity.PPM combines piecewise-parabolic interpolation with a nonlinear Riemann solver.
- 4.2.4 FLASH (AMR): FLASH uses an AMR PPM hydrodynamical solver with formal second-order accuracy in space and time and refinement up to the resolutions listed in Table 1.Its critical steps reach third- or fourth-order accuracy.
- 4.2.5 Gasoline (SPH): Gasoline is a Tree+SPH code using standard and shear-reduced artificial viscosity, an asymmetric energy equation, close entropy conservation, and a compact-support spline kernel.Its viscosity coefficients α and β control viscosity strength, shock capture, and particle interpenetration, with standard values α = 1 and β = 2.
- 4.2.6 GADGET-2 (SPH): GADGET-2 is an updated TreeSPH code with entropy-conserving SPH and an artificial-viscosity formulation based on a signal velocity that vanishes for approaching-free particle pairs.Its thermodynamic state is defined through specific entropy rather than specific thermal energy.
5 RESULTS OF THE SIMULATIONS
Grid simulations resolve cloud instabilities and complete fragmentation, whereas SPH simulations retain a dense cloud and lose mass mainly through gradual stripping. Resolution changes grid morphology and small-scale structure, while initial-condition changes alter instability phases but not the eventual grid-cloud disruption timescale.
- 5 RESULTS OF THE SIMULATIONS: By 2.5 τKH, grid simulations show dynamical instabilities and complete fragmentation, whereas most SPH gas remains in one cold, dense blob.The grid cloud rapidly disrupts and mixes after large-scale Kelvin–Helmholtz growth, while SPH continues stripping without comparable fragmentation.
- 5 RESULTS OF THE SIMULATIONS: SPH cloud evolution begins with bow-shock compression, lateral elongation, ablation, and downstream vorticity, but surface KH and RT instabilities later drive grid-cloud fragmentation and mixing.The two methods initially evolve similarly before their late-time behavior diverges.
- 5 RESULTS OF THE SIMULATIONS: Both methods lose similar cloud mass until roughly τKH, but after that the grid cloud rapidly mixes while SPH retains about 40% at 2.5τKH.Before τKH, vortex shedding through Bernoulli zones is the main mass-loss mechanism; afterward, dynamical instabilities dominate grid mass loss.
- 5.1 Resolution dependence: Increasing grid resolution changes instability phase, reduces diffusion, and reveals more small-scale fragments, while the overall destruction time remains similar.The dominant symmetry-axis KHI weakens from Enzo 64 to Enzo 128 and disappears at Enzo 256, as higher resolution exposes additional small-scale instabilities.
- 5.1 Resolution dependence: Grid and SPH resolution cannot be directly equated because grid cells are uniform, whereas SPH particles require neighboring particles and have effective resolution set by the smoothing kernel.The high-resolution comparison uses a 256×256×1024 grid and 10^7 SPH particles, but the numbers are not equivalent resolution elements.
- 5.1 Resolution dependence: SPH simulations do not resolve small-scale instabilities, and increasing SPH resolution does not reduce mass loss; lower resolution instead produces more rapid post-τKH loss.The latter behavior is attributed to more violent, bullet-like momentum transfer by more massive particles.
- 5.2 Initial Seeds: Initial-condition perturbations alter the phase and appearance of the destructive mode, but grid clouds still become debris and mix within a few τKH.Analytic, symmetric initial conditions with reduced perturbations produce a different instability phase, yet the cloud’s fundamental disruption outcome is unchanged.
6 WHY SO DIFFERENT?
The comparison attributes SPH–grid discrepancies primarily to artificial-viscosity effects and, more fundamentally, to SPH pressure-force errors across steep density gradients. Grid methods resolve Kelvin–Helmholtz instabilities, whereas standard SPH suppresses them when density contrasts are present.
- 6.1 Artificial viscosity: Artificial viscosity damps small-scale velocity perturbations and diffuses post-shock vorticity, suppressing instability growth and smearing turbulence.The effect is especially important for post-shock vorticity that should help destabilize the cloud surface.
- 6.1 Artificial viscosity: SPH viscosity changes cloud stability but does not recover agreement with grid-based codes, indicating a more fundamental discrepancy.Lower viscosity makes the cloud less stable, while standard viscosity produces the most stable evolution; neither matches grid results.
- 6.2 Resolving instabilities: Grid simulations resolve Kelvin–Helmholtz instabilities with growth times close to analytical expectations, whereas SPH rapidly damps velocity and density perturbations despite changes in resolution, viscosity, and initial conditions.The comparison uses two shearing fluids with different densities and imposed interface perturbations.
- 6.3 Mind the gap: SPH forms a smoothing-kernel-sized gap at steep density interfaces because density overestimation creates repulsive pressure forces that prevent physical contact between phases.The resulting vacuum layer also explains cloud mass loss through expansion from the cloud edges rather than physical stripping at the leading surface.
- 6.3 Mind the gap: When the density contrast vanishes, the gap cannot form and SPH is able to capture Kelvin–Helmholtz instability, supporting the gap as the source of suppression.The less evolved standard-viscosity run also illustrates the separate effect of viscosity.
7 SUMMARY
Grid methods resolve dynamical instabilities and mixing that standard SPH handles poorly because density gradients generate erroneous pressure forces and boundary gaps. The methods agree during early cloud stripping, but diverge once large-scale instabilities grow, with implications for astrophysical multiphase-fluid simulations.
- Grid codes resolve dynamical instabilities and mixing, whereas current SPH techniques resolve these processes poorly or not at all.
- Standard SPH produces erroneous pressure forces near steep density gradients, creating gaps between high-density regions.The gaps arise from asymmetric density within the smoothing kernel.
- These differences matter for modeling galaxy gas stripping, disk formation, star formation, feedback, turbulence, and interacting multiphase fluids.
- Grid and SPH methods agree during early pressure-driven gas stripping, but diverge after Kelvin-Helmholtz and Rayleigh-Taylor instabilities grow.
- SPH boundary gaps can prevent viscosity information from transferring between shearing fluids, obstructing Kelvin-Helmholtz instability growth.