Challenges related to Particle Regularization in SPH

1. Challenges related to Particle Regularization in SPH S. Børve (University of Oslo/Norwegian Defence Research Establishment) R. Speith (University of Tübingen) M. Omang (Norwegian Defence Estates Agency) J. Trulsen (University of Oslo)

2. Outline • RSPH overview • A closer look at particle regularization • Numerical diffusion and viscosity • Solid boundaries and free surfaces • Current application: planet-disc interaction • Conclusion

3. Background In trying to achieve accurate descriptions of shocks in compressible fluids (HD/MHD) using SPH, we addressed the following issues: • Reduced accuracy due to particle disorder (lack of reproducibility). • Low accuracy in low density regions (Lagrangian adaptivity). • Increased noise level when variable smoothing length (h) is introduced.

4. Regularized SPH (RSPH) in Short The following characterizes the current RSPH version: • Dynamic calculations performed with standard SPH equations. • Quantized smoothing length values (h-levels) are used. • All SPH particles have fixed h and only interact with particles on the same h-level. • Auxiliary particles handles the coupling of different h-levels (2nd order interpolation). • The entire particle distribution redefined at temporal intervals. • During regulariztion, h-profiles are determined according to specified functions or from calculations based on second derivative estimates.

5. WHY Particle Regularization • Can control the level of particle disorder and thereby reduce the overall numerical noise level. • Allows for tailor-made particle distributions providing optimized numerical descriptions. • Represents a useful tool in analyzing numerical errors as one might study the effects of specific particle distributions. Solar wind interaction with a 2D magnetic dipole The simulation using RSPH revealed effects similar to magnetic reconnection despite the fact that the equation were in principle ideal MHD. Numerical resistivity due to particle regularization caused it. (Børve 2001).

6. HOW TO Regularize At temporal intervals, particle regularization is performed: • For each old particle, o, determine optimal, ho. A global, quantized h-profile is made based on the optimal h-values. • Create new particles according to global h-profile by placing them regularly (cartesian/polar), with modifications due to boundaries, inside the domain. • Determine the volume occupied by each particle (old/new) and the spatial overlap between old-new pairs of particles (Voronoi scheme). • New particles inherit conservatively mass, momentum, and energy from old particles based on ratios of spatial overlap.

7. Problems with Regularization • Any modification/relaxation of the particle distribution needed at initialization, is also needed at regularization. • Conservation of linear momentum usually excludes exact conservation of angular momentum and kinetic energy. • Regularization introduces numerical diffusion and viscosity. • Currently, nothing is done to accurately reproduce free surfaces in regularization. • Regularization causes RSPH to be more complex to code than SPH.

8. Regularization & Linear (M)HD Waves Figure 1: Combination of sound and Alfvén waves with Ps/PA equal to 10 (left) and 0.1 (right). Regularization at 5x or 0.2x slower period (solid/dashed line).

9. Regularization & Compressible (M)HD Shock Waves Figure 2: Plane shock tube test with strong compression. Comparison of analytic (solid line) and RSPH (dots) solutions.

10. Regularization & Compressible (M)HD Shock Waves (cont.) Figure 3: Comparison of analytic (solid line) and SPH (dots) solutions (Price & Monaghan, 2005).

11. Simple Test of Numerical Diffusion Figure 3: A simple 1D test where a square pulse in a force free system is advected exhibits numerical diffusion similar to that found in Eulerian codes. The solid, dotted, and dashed lines indicate the pulse shape after 0, 24, and 48 regularizations, respectively. The drift velocity of the particles have been chosen so that old and new particles always are slightly misaligned (worst case scenario) when performing regualrization.

12. Viscously spreading accretion ring First proposed and solved by Lüst already in 1952 assuming small kinematic viscosity and no thermal effects. The analytic solution for the surface density is: where M is the total ring mass, τ = 12νt/R02 is the viscous time unit, x = R/R0 is dimensionless radial position, and I1/4 is the modified Bessel function of the order ¼. The viscosity model (proposed by Meglicki et al., 1993) uses the artificial viscosity with α = ν/(σcsh) and β = 0, where σ is a constant found by an integral over the specific kernel.

13. Regularization and numerical viscosity Figure 4: Comparison of 4 different runs at τ = 0.27. Runs (a-c) correspond to regularization intervals of 320, 640, and 1280 time steps. Run (d) has no regularization. All runs have h = 0.02 and an initial particle spacing of h/1.4. The estimated viscosity has been estimated to be 36% (a), 11% (b), and 1% (c) higher than the analytic viscosity.

14. Comparison with other SPH solutions Speith has recently (in prep.) studied various approaches to simulating viscosity in SPH. As part of the work he has also looked for new formulation for obtaining higher order derivatives in SPH. Figure 5: Results for the viscous ring test with the Flebbe approach (Flebbe et al. 1994) and with iSPH (Speith 2007).

15. Boundaries and free surfaces A potential problem with regularization is solid boundaries or free surfaces that cannot be accurately described by particles placed regularly on a cartesian grid. Figure 6: Some initial work was done in 2002 to apply RSPH to incompressible water waves (Børve et al. 2003). This is a close-up view of a breaking wave at 3 different stages. Notice the step-like free surface. Could we do better ?

16. Semi-regular particle distributions Figure 7: Modified particle distribution near circular boundaries. Only the distance to the surface is corrected, leading to transversal irregularities. Figure 8: Attempt to modify the distribution in the visous ring test assuming circularly shaped free surface.

17. Identifying free surfaces • The viscous ring test is the first problem where we have tried to model free surfaces in compressible fluids. Problems related to a too quickly expanding particle-filled region was observed. Two main reasons: • The Voronoi scheme do not handle free surfaces properly. • Regularly placed new particles cannot reproduce a general free surface very accurately. Figure 9: α-shapes The process of determining a convex hull of a set of particles is related to finding the Voronoi diagram. A generalization of convex hulls can be used to identify so-called α-shapes. (Edelsbrunner & Mücke, 1994)

18. Disc-planet interaction • Results from a comparative study within the EU-RTN project ”The Origin of Planetary Systems” published in 2006 (de Val-Borro et al., 2006). • Test problem: Jupiter-sized planet in a fixed circular orbit in a circumstellar disc (sun-like star). Neglecting self-gravitation. • 17 independent codes, including 2 SPH codes, tested. • Gap-formation due to tidal interaction leading to reduced resolution in SPH results near the planet orbit. • RSPH simulation performed with constant h. Only artificial (not kinematic) viscosity.

19. RSPH Figure 10: Density plotted on logarithmic scale after 100 orbits. Comparing the codes JUPITER (Pember et al., 1995), SPHTREE (Benz, 1990), PARASPH (Schäfer et al., 2004), and RSPH (Børve et al., 2001). The lower colour bar applies to the RSPH results, the upper colour bar applies to the remaining results.

20. The following can be seen from looking at radial cuts: • RSPH: High peak values and sharp shock fronts. • JUPITER: Follows the main trend among the grid codes. • PARASPH: Low resolution, but consistent with other codes. • TREESPH: Deviates more from the other codes. Figure 11: Radial cuts opposite the planet position. Results from the comparative study on the left (JUPITER: pink dotted, SPHTREE black dashed, and PARASPH red dashed), RSPH results on the right.

21. Conclusion • RSPH utilizes particle regularization to increase the overall accuracy of the SPH computations. • Particle regularization leads to artificial effects such as enhanced diffusion and viscosity. • The viscous ring test results show that an approriately chosen regularization interval can result in acceptable accuracy and noise level. • Solid boundaries and free surfaces can constitute a challenge to regularization. • Results from the planet-disc interaction test support the idea that RSPH can be an interesting alternative to standard SPH.