Smoothed Particle Hydrodynamics and Magnetohydrodynamics

Daniel J. Price

Source-linked AI summary

Smoothed Particle Hydrodynamics and Magnetohydrodynamics

Daniel J. Price

arXiv:1012.1885v1Dec 8, 2010astro-ph.IMastro-ph.COastro-ph.SRphysics.flu-dynphysics.plasm-ph

TL;DR

SPH and SPMHD need formulations grounded in their density estimate, with careful treatment of misconceptions, instabilities, dissipation, and divergence control. The paper develops these foundations variationally, reviews practical remedies and current developments, and supplies a multidimensional code testbed for the numerical examples.

Problem
The paper addresses widespread misconceptions and unresolved practical issues in applying SPH and SPMHD, including neighbor-number choices, pairing and tensile instabilities, shocks, and divergence control.
Method
The review derives SPH equations of motion and energy self-consistently from the density estimate, interprets them through interpolation theory, and surveys practical SPMHD formulations and remedies.
Results
The paper provides a unified overview of SPH methodology, identifies the origins of particle-distribution and MHD instabilities, and reviews approaches for dissipation, tensile stabilization, and enforcing ∇·B = 0.
Takeaways & Limitations
The paper presents the density estimate and local conservation as central to accurate, robust SPH, while practical SPMHD requires explicit treatment of tensile instability, shocks, and magnetic divergence.
Takeaways & Limitations
Neighbor number is not an independent resolution parameter: changing it can change h relative to particle spacing, and its interpretation depends on particle-density uniformity and dimension.

Abstract

from arXiv · show

This paper presents an overview and introduction to Smoothed Particle Hydrodynamics and Magnetohydrodynamics in theory and in practice. Firstly, we give a basic grounding in the fundamentals of SPH, showing how the equations of motion and energy can be self-consistently derived from the density estimate. We then show how to interpret these equations using the basic SPH interpolation formulae and highlight the subtle difference in approach between SPH and other particle methods. In doing so, we also critique several `urban myths' regarding SPH, in particular the idea that one can simply increase the `neighbour number' more slowly than the total number of particles in order to obtain convergence. We also discuss the origin of numerical instabilities such as the pairing and tensile instabilities. Finally, we give practical advice on how to resolve three of the main issues with SPMHD: removing the tensile instability, formulating dissipative terms for MHD shocks and enforcing the divergence constraint on the particles, and we give the current status of developments in this area. Accompanying the paper is the first public release of the NDSPMHD SPH code, a 1, 2 and 3 dimensional code designed as a testbed for SPH/SPMHD algorithms that can be used to test many of the ideas and used to run all of the numerical examples contained in the paper.

1. Introduction

This pedagogical review revisits the fundamentals of SPH and develops practical guidance for robust SPH and SPMHD algorithms. It is intended to clarify common misconceptions while connecting theory with numerical practice and an accompanying testbed code.

1. Introduction: The paper presents general principles intended to apply across the astrophysical, geophysical, engineering, film, and computer-games contexts where SPH is used.The stated principles are not restricted to the topical issue of computational astrophysics.
1. Introduction: The review focuses on why SPH works, how to formulate general algorithms, and how to avoid common implementation pitfalls.It is not intended as a comprehensive survey of applications or physical-model implementations.
1. Introduction: SPH fundamentals are used as the basis for discussing robust and accurate formulations of Smoothed Particle Magnetohydrodynamics.The paper motivates this progression by the importance of MHD in astrophysical problems.
1. Introduction: The released ndspmhd code provides 1D, 2D, and 3D exercises for testing SPH/SPMHD algorithms and reproducing the paper’s numerical examples.The code is presented as a resource for advanced researchers and students.

2. The foundations of SPH: Calculating density

SPH begins with estimating a continuous density field from point particles using a weighted sum over neighbors. Kernel choice and smoothing-length control determine resolution, smoothness, sampling quality, and the interpretation of neighbor number.

2. The foundations of SPH: Calculating density: Unlike fixed meshes, adaptive local sampling addresses clustered and sparse regions by scaling the sampling volume with local particle density.A fixed-neighbor local-volume estimate can remain noisy because boundary particles switch abruptly between inclusion and exclusion.
2. The foundations of SPH: Calculating density: SPH computes density by a weighted sum over nearby particles, with kernel weight decreasing with separation according to smoothing length h.This mesh-free estimator is the paper’s foundational construction.
2. The foundations of SPH: Calculating density: A good density kernel is positive, symmetric, monotonic with distance, smoothly differentiable, and relatively flat near the center.These properties reduce sensitivity to small positional changes while preserving a stable density estimate.
2. The foundations of SPH: Calculating density: Gaussian kernels are smooth and accurate but require interactions with all particles, giving density-evaluation cost O(N^2) at particle locations.Compact-support B-splines provide a practical alternative by limiting interactions to a finite radius.
2. The foundations of SPH: Calculating density: Higher-order B-splines increase kernel smoothness and support radius while preserving the meaning of h, unlike stretching the cubic spline to change neighbor number.The M4, M5, and M6 kernels are truncated at 2h, 2.5h, and 3h, respectively.
2. The foundations of SPH: Calculating density: Increasing neighbor number alone is not a general convergence strategy because changing it through the h–ρ relation changes h relative to particle spacing.The neighbor-number parameter also varies with dimension and only tracks actual neighbors under approximately uniform particle density.
2. The foundations of SPH: Calculating density: Higher-order kernels can yield potentially negative density evaluations and require very well-sampled kernels, increasing computational cost.The paper suggests switching between high- and low-order kernels as one possible response.

3. From density to equations of motion

The paper derives SPH dynamics from a discrete Lagrangian, using variational mechanics rather than treating the particle equations as an independent prescription. This exposes both conservation structure and assumptions requiring special treatment.

3. From density to equations of motion: The SPH algorithm can be derived from the density estimate because the remaining conservative formulation follows from that estimate.This is why the review spends substantial attention on density before introducing equations of motion.
3. From density to equations of motion: A discrete Lagrangian containing kinetic and thermal energies provides the starting point for deriving particle equations of motion.The internal energy may depend on density and entropy.
3. From density to equations of motion: The discrete Hamiltonian system directly retains the symmetries and conservation properties of the discrete Lagrangian.The particle motions are constrained by those discrete-system properties without explicit reference to the continuum system.
3. From density to equations of motion: The Euler–Lagrange equations follow by applying the principle of least action to variations in particle coordinates and velocities.The derivation assumes variations vanish at the initial and final times.
3. From density to equations of motion: Exact conservation claims exclude time-integration error because the derivation does not discretize the time integral.The paper distinguishes spatial or algebraic conservation from errors introduced by the integration scheme.
3. From density to equations of motion: Because the Lagrangian is assumed differentiable, discontinuities require dissipative treatment that smooths them to a resolvable scale of a few h.This assumption excludes discontinuous solutions from the direct variational derivation.

3.3. Equations of motion

Differentiating the density estimate and accounting for smoothing-length variation yields the SPH equations of motion. Their symmetry produces exact momentum conservation, while pressure-dependent residual forces create a known numerical issue.

3.3. Equations of motion: The thermodynamic relation at constant entropy connects changes in internal energy to density changes and pressure work through the SPH volume estimate V = m/ρ.The derivation uses dV = −m/ρ^2 dρ.
3.3. Equations of motion: Differentiating the density estimate supplies the spatial derivative needed to express the Lagrangian in particle coordinates.When h depends on density, the derivative includes a smoothing-length-gradient correction.
3.3. Equations of motion: The resulting variable-smoothing-length equations reduce to the standard SPH expression when the smoothing length is constant.The general form is obtained by combining the thermodynamic and density-gradient derivatives in the Euler–Lagrange equations.
3.3. Equations of motion: The equations conserve total linear momentum exactly because the pairwise kernel-gradient contributions cancel by antisymmetry.This conservation follows from translation invariance of the Lagrangian and density estimate.
3.3. Equations of motion: The equations also conserve total angular momentum exactly through antisymmetry of pairwise terms and rotational invariance.The relevant cancellation uses the antisymmetry of the particle cross products.
3.3. Equations of motion: Even constant pressure can produce a non-zero interparticle force because the equations depend on absolute pressure, not only pressure differences.The paper identifies this as a spurious force whose importance and problems are discussed separately.

3.4. Energy equation

The SPH energy equation can be derived consistently from Hamiltonian dynamics, with equivalent choices of thermal, total, or entropy variables under a Lagrangian formulation.

Hamiltonian dynamics derives the dissipationless SPH energy equation alongside the equations of motion.The conserved total energy follows from the Hamiltonian because the Lagrangian has no explicit time dependence.
Evolving thermal energy u, total specific energy e, or entropy variable K is equivalent in SPH when equations come from the Lagrangian.Differences arise only from the time-integration algorithm, unlike Eulerian formulations where advection differencing changes the formulations more substantially.
The dissipationless SPH scheme exactly conserves total energy through antisymmetry of the particle double sum.This conservation follows from Lagrangian time-translation symmetry.
Entropy: Using K makes evolution independent of the time-integration algorithm when dissipation is absent.Its disadvantage is greater difficulty with non-ideal equations of state.
Entropy: The label ‘entropy-conserving’ refers more precisely to smoothing-length gradient accounting and a consistent energy equation than to evolving an entropy variable.Entropy per particle is also exactly conserved when the alternative energy equations include the required smoothing-length gradient terms, aside from timestepping differences.

3.5. Summary

The summary presents SPH as a discrete Hamiltonian system built from a density estimate, thermodynamics, and a Lagrangian, while identifying assumptions that become important near discontinuities.

The complete SPH system evolves density, velocity, and thermal energy, with equivalent alternatives for the energy equation.The equations are derived from the density estimate and the Lagrangian rather than introduced directly from continuum hydrodynamics.
The derivation assumes exact time derivatives, differentiable density and thermal energies, and no entropy change.It also uses the particle-volume relation dV = −m/ρ^2 dρ.
The differentiability and no-entropy-change assumptions become relevant when treating shocks and other discontinuities.The paper discusses these issues further in its treatment of shocks.

3.6. Alternative formulations

Within a Hamiltonian SPH formulation, changing the dissipationless algorithm while preserving its conservation properties offers only limited freedom.

Two consistent changes are available: alter the density calculation or introduce additional physical terms with associated constraints.This restriction applies apart from changing the first law of thermodynamics.

4. Kernel interpolation theory and SPH derivatives

The review uses kernel interpolation primarily to interpret Hamiltonian-SPH equations and analyze derivative errors, while contrasting linear accuracy with nonlinear conservation properties.

The paper deliberately introduces kernel interpolation to interpret derived SPH equations and guide additional physics, rather than defining SPH through linear errors.The authors argue that linear-error focus can miss nonlinear behavior important to practical performance.
The summation interpolant discretizes kernel smoothing over particles and becomes the SPH density estimate when A = ρ.It also provides a general interpolation of particle-defined quantities at arbitrary spatial points.
Differentiating the density summation yields a discrete continuity equation that is an exact, time-independent solution because particle masses remain fixed.This connects the density estimate directly to mass conservation.
The Lagrangian-derived force and energy equations provide valid discrete hydrodynamic equations, but they differ from basic kernel-interpolation operators.Their construction follows from the discrete Hamiltonian system rather than directly from interpolation formulae.
Errors: Interpolation accuracy depends on discrete normalisation conditions, particle regularity, and the ratio h/∆p.Regular particles and an appropriate smoothing-length-to-spacing ratio are important for accurate density estimates.
Errors: Gradient interpolation has O(h^2) integral error, and explicit error terms can construct more accurate gradient operators.The corresponding discrete accuracy depends on how well particle summations reproduce the required integral conditions.
Errors: Matrix-based gradient normalisation requires nine extra stored quantities per particle in three dimensions and prior knowledge of density.This makes the correction cumbersome in practice.
Alternative formulations: Linear-error formulations can improve accuracy for regular particles in linear or weakly nonlinear problems, but may sacrifice exact momentum and energy conservation.The Hamiltonian formulation instead emphasizes nonlinear symmetries, conservation, and global constraints relevant to long-term stability.

5. Why a bad derivative leads to good derivatives: The importance of local conservation

SPH deliberately uses locally conservative, pairwise forces whose lower formal gradient accuracy helps maintain regular particle arrangements and therefore improves practical gradient estimates. The same conservation-driven behavior explains settling, remeshing, and instabilities that constrain kernel and stress choices.

Local conservation: The Lagrangian formulation prioritizes symmetry and conservation over the most accurate linear-error gradient estimate.The resulting acceleration operator is formally poor but is the only one preserving the Lagrangian’s symmetry and conservation properties.
Local conservation: Pairwise momentum conservation makes particles respond to irregular arrangements, producing an intrinsic remeshing effect that maintains local order.Non-conservative accurate-gradient formulations lack this regularisation and may require explicit remeshing.
Example 1: Settling: In the settling experiment, conservative SPH rapidly regularises random particles and obtains good gradients, even without artificial viscosity.Artificial viscosity accelerates settling, while the non-conservative relative-pressure formulation leaves particles randomised and causes exponential total-energy growth.
Example 2: Shock tube: In a multidimensional shock, compression creates anisotropic particle arrangements that subsequently remesh through transverse, resolution-scale motions.These motions arise from SPH regularisation rather than the physical shock itself.
Tensile instability: Negative stresses reverse pairwise repulsion into attraction, causing the tensile instability and unphysical particle clumping in MHD and elastic dynamics.The instability occurs when the employed stress tensor permits physically negative stresses.

6. Second derivatives and dissipation terms in SPH and SPMHD

The paper develops second-derivative operators and dissipation terms for SPH, emphasizing kernel forms that improve accuracy and conductivity terms that correct contact-discontinuity errors. It also interprets artificial viscosity in terms of physical transport and demonstrates improved shock behavior with higher-order kernels.

6.1. The SPH Laplacian: Second-derivative estimates based directly on compact bell-shaped kernel second derivatives are poorly sampled because those derivatives change sign and may be discontinuous.The resulting normalization conditions require extremely good sampling, motivating a separate second-derivative kernel.
6.1. The SPH Laplacian: The Brookshaw Laplacian can be interpreted as using the first-derivative kernel divided by particle spacing, equivalently applying a new kernel Yab.This interpretation extends to formulations of thermal conductivity and other dissipative terms, although vector second derivatives require a more involved treatment.
6.1. The SPH Laplacian: Constructing Y′′(q) = −2w′(q)/q from kernel gradients produces monotonically decreasing, positive functions that estimate second derivatives more accurately.The paper compares these constructed functions for the M4 cubic and M6 quintic kernels against standard second derivatives and normalization conditions.
6.3.1. Interpretation of SPH artificial viscosity terms: Artificial viscosity is Galilean invariant, vanishes for rigid rotation, conserves linear and angular momentum, and behaves like physical shear and bulk viscosity proportional to resolution length.For shock capturing it acts only during particle approach, so the effective coefficients in uniform shear are approximately half the nominal values.
6.3.2. General formulation of dissipative terms in SPH and SPMHD: Artificial viscosity smooths shocks but leaves a pressure blip at contact discontinuities, suppressing Kelvin–Helmholtz mixing through an artificial surface-tension effect.Adding artificial conductivity removes the pressure blip and restores mixing in the two-dimensional Kelvin–Helmholtz problem.
6.3.4. Examples 4 and 5: One and two dimensional shock tubes, and Kelvin-Helmholtz instabilities: Using the M6 quintic kernel instead of the cubic spline substantially reduces particle-resettling noise behind a two-dimensional shock, producing results comparable to the one-dimensional case.The example illustrates how higher-order kernels can be used to obtain convergence in SPH.

7. Smoothed Particle Magnetohydrodynamics from a Lagrangian

The paper derives SPMHD equations by combining an MHD Lagrangian with the induction equation and SPH derivative formulations. The resulting discretization conserves linear momentum and total energy, and current-loop advection remains unchanged after 1000 domain crossings without explicit resistivity.

7.1. The MHD Lagrangian: SPMHD is constructed by writing an MHD Lagrangian, imposing physical constraints, and consistently deriving the equations of motion.The MHD Lagrangian differs from the hydrodynamic one by subtracting a magnetic-energy term.
7.2. The induction equation: Because the magnetic field is not directly a function of particle coordinates, deriving the equations requires a variational formulation together with an SPH representation of the induction equation.The induction equation supplies the relation between magnetic-field and particle-coordinate perturbations.
7.3. Equations of motion: The resulting SPMHD force uses the MHD stress tensor and conserves linear momentum exactly through pairwise force symmetry.The anisotropic magnetic force derives from the numerical induction equation, while the isotropic term derives from magnetic energy and density perturbations.
7.4. Energy equation: The SPMHD formulation also conserves total energy exactly, yielding an evolution equation for the specific energy.The energy result follows from the discrete acceleration, continuity, and induction equations.
7.5. Examples: The current-loop test shows no magnetic-field change after 1000 computational-domain crossings without explicit resistivity, independent of propagation direction, resolution, or advection velocity.The result is described as exact advection in SPMHD for this test.

8. The tensile instability in MHD

The paper explains that conservative SPMHD becomes vulnerable when magnetic stress turns attractive, producing tensile instability along field lines. It reviews corrections that trade exact conservation for stability and shows that Morris and Børve-type approaches stabilize an Alfvén-wave test.

8. The tensile instability in MHD: The conservative SPMHD force becomes unstable when total stress is negative, because particle interactions become attractive rather than repulsive.In the multidimensional stability criterion, instability occurs when magnetic pressure exceeds gas pressure: 1/2 B^2/μ0 > P.
8. The tensile instability in MHD: The instability originates from a monopole force proportional to numerically non-zero magnetic-field divergence and directed along the magnetic field.This source term can drive attraction when pressure cannot counteract it, although the physical divergence should vanish.
8. The tensile instability in MHD: Exact conservation and exact perpendicularity to the magnetic field cannot apparently be achieved simultaneously in SPMHD by enforcing the divergence constraint.The paper notes that the tensile instability can occur even in one dimension despite trivial enforcement of Bx = const.
8.1. Fix 1: subtract a constant from the stress: Subtracting the maximum negative stress removes the instability while conserving momentum, but it does not conserve total energy and can produce unphysical effects if the correction is large.The method has relatively few side effects when the correction is small and is useful under large constant external stresses.
8.2. Fix 2: the Morris approach: The Morris approach retains conservative isotropic forces but uses a non-conservative anisotropic derivative, sacrificing exact momentum and energy conservation while preserving correct MHD-shock jump conditions.The paper reports that its non-conservation errors are small in practice, including on strong shock-tube problems.
8.3. Fix 3: subtract the unphysical source term: The Børve correction subtracts the monopole source term using the same divergence discretization as the force equation and can be switched off where the conservative formulation is stable.The correction violates exact momentum and energy conservation only insofar as the discrete divergence is non-zero.
8.4. Example: In the circularly polarized Alfvén-wave test, the conservative formulation destroys the particle arrangement, whereas Morris and Børve formulations keep the wave stable and propagating correctly.The conservative calculation shows particle attraction along the magnetic field lines before failing at later times.

9. Dissipation terms and shocks in SPMHD

SPMHD shock capturing requires dissipation tailored to jumps in velocity, thermal energy, and magnetic field. The paper illustrates these choices on standard MHD shock-tube problems, including a 2D result comparable to the 1D solution.

Dissipation for MHD shocks: Artificial viscosity, thermal conductivity, and resistivity treat jumps in v, u, and B, respectively, enabling capture of all seven MHD shock-tube discontinuities.The seven structures are slow, Alfvén, and fast waves propagating in both directions, plus the contact discontinuity.
Artificial resistivity: The resistivity signal velocity can differ from viscosity’s, with ndspmhd using an averaged Alfvén speed as a simple choice.
Artificial resistivity: Artificial resistivity can use magnetic-field components perpendicular to the line of sight, although the total-energy-jump formulation performed better in practice.
Dissipation for MHD shocks: Artificial thermal conductivity and resistivity should generally be applied whether particle pairs are approaching or receding.
Shock-tube examples: 800 × 24 particles with the M6 quintic kernel produced a 2D Brio–Wu shock-tube solution comparable to the 1D solution.The comparison uses the numerical solution from Balsara (1998).

10. The divergence constraint in SPMHD

SPMHD must control numerical divergence errors because ∇·B=0 is imposed initially but is not maintained indefinitely under truncation error. The paper compares ignoring, cleaning, and preventing these errors, while documenting important practical limitations of each approach.

The divergence problem: Truncation error causes the numerical ∇·B=0 condition to deteriorate, motivating approaches that ignore, clean, or prevent monopole growth.
Monopole-conserving formulations: The 8-wave formulation uses a monopole-conserving induction equation and, in SPMHD, requires momentum-source subtraction to stabilize the tensile instability.
Ignoring divergence errors: Monitoring divergence without correction can work for some test problems, but other problems require more active control.The monitored dimensionless error is typically hoped to remain below a few percent.
Divergence cleaning: Elliptic cleaning is only approximate with a Green’s-function projection and is impractical or inefficient for 3D simulations with individual particle timesteps.
Divergence cleaning: Hyperbolic/parabolic cleaning can reduce maximum divergence by at most a factor of ~2 in real problems, while poor parameters may increase errors through constructive interference.
Preventing divergence errors: Constrained transport cannot be directly applied to SPH because it requires surface rather than volume integrals.
Euler potentials: Euler potentials enforce ∇·B=0 by construction and have enabled magnetized star-formation calculations without divergence-driven numerical blow-up.
Euler potentials: Euler potentials omit full MHD topology because A·B is identically zero, limiting them to initially simple fields and preventing multiple complete windings.

11. Summary

The paper grounds SPH in density estimation and derives its motion and energy equations self-consistently through a variational principle. It also interprets the method through interpolation theory and discusses accuracy and general derivative construction.

Summary: SPH formulations begin with a density estimate, from which the equations of motion and energy can be self-consistently derived using a variational principle.
Summary: Kernel interpolation theory provides an interpretation of the SPH equations and supports construction of more accurate, general derivatives through linear error analysis.