Recent Changes - Search:


Codes and methods

Theory and analysis

Past events


edit SideBar

VlaSolve: an improved traditional semi-Lagrangian code in spherical symmetry

Warning: this page is better viewed with Firefox

Author of the code: T. Sousbie
Publication introducing and exploiting the code (to be cited when using VlaSolve):

Download VlaSolve: vlasolve_v0.9.0.tar.gz

Table of Contents

The traditional semi-Lagrangian splitting algorithm

In spherical symmetry, Vlasov equation reads

where v is the radial component of the velocity, j is the angular momentum, Mr=M(r'<r ) the mass inside a sphere of radius r and G the gravitational constant.

To solve equation (1) numerically, we wrote an improved implementation of the splitting algorithm of Cheng & Knorr (1976), following closely Fujiwara (1983). This algorithm exploits directly the Liouville theorem, namely that the phase-space density is conserved along motion: the equations of the dynamics during a time step are separated in "drift" and "kick" parts in terms of Hamiltonian dynamics and are resolved backwards,

In the drift part, which includes, in addition to the trivial ballistic part, the inertial contribution from the conserved angular momentumr(t-dt/2) and v(t-dt/2) can be estimated analytically (see, e.g., Colombi & Touma 2008), while the kick part accounts for the gravitational contribution, computed from f*.

The phase-space distribution function is sampled on a mesh, and each step is performed by using test particles coinciding with mesh sites and following the equations of motion split as above. Resampling of f*, f**, and finally the phase-space distribution function at next time step is performed by using third order spline interpolation on the grid.

In the current implementation of VlaSolve, time step is a constant parameter chosen by the user. Phase space is sampled into a rectangular mesh of size (Nr,Nv,Nj) which spans the domain of interest as follows: the position RminrRmax is sampled in a logarithmic fashion, the velocity -vmaxvvmax linearly, and the angular momentum 0 ≤ jJmax is such that the position of the kth angular momentum slice is jk=Jmax(k-1/2)2/Nj2, corresponding formally to the interval [Jmax (k-1)2/Nj2, Jmax k2/Nj2].

Reflecting boundaries with time delay



Figure 1: Comparison between the reflecting central sphere method (upper panel) and our improved delayed central sphere implementation (lower panel), in (r,j) space. The artificial features introduced by instantaneous reflection are suppressed by introducing time delay. From Colombi et al. (2015).

Using a logarithmic scale along the radial direction proves particularly advantageous for tracing tiny features close to the center of the system that subsequently expand to larger scales as they get away from it. One downside of this procedure is that a small sphere of radius Rmin is necessarily missing from the computing domain. This shortcoming is usually dealt with by considering the central part as a small reflecting sphere (see, e.g, Gott 1973, Fujiwara 1983). Although simple, this method presents the major drawback of introducing a systematic lag between orbits: particles reaching the reflective kernel boundary instantly travel the 2Rmin distance through the central region, while they should actually take a finite time depending on their radial velocity and angular momentum.

In VlaSolve, we improve the reflecting sphere method by taking into account the actual time spent spent by matter elements traveling inside the region rRmin, which is made easily possible by neglecting the gravitational force. This assumption, combined with the fact that we use a constant time step, allows us to precompute the trajectories of these matter elements once and for all by creating link lists associated to each grid site whose corresponding test particle radial position r half a time step backward in time is such that rRmin. A comparison of the results obtained with the reflective central sphere to our improved delayed central sphere is shown on figure 1. The improvements are unquestionable.


Thanks to its hydrid parallel nature, VlaSolve can be used to run large simulations. Shared memory parallelism with OpenMP library is achieved in j space by taking advantage of the fact that the angular momentum is a conserved quantity. Distributed memory parallelism with MPI library is implemented by performing domain decomposition in (r,v) space, following the approach of Crouseilles et al. (2009), who propose to localize the cubic spline interpolation to each domain by using Hermite boundary conditions between the domains with an ad hoc reconstruction of the derivatives.

Illustration: simulation of a Hénon sphere and comparison to the N-body result

The Hénon sphere (Hénon 1964) corresponds to a very simple initial set up for the phase-space distribution function:

with (4π/3) ρ0 RH3=M, the total mass of the system. The main parameter defining the properties of this system is its virial ratio R=|2T/W|=5 RH σv2/M. Figure 2 shows the phase-space distribution function obtained at various times for a Hénon sphere with total mass unity, RH=2 and R≈0.5. The results obtained with VlaSolve are compared to a N-body simulations performed with the public code Gadget-2 (Springel, Yoshida & White 2001; Springel 2005). In order to avoid exaggerate aliasing in the Vlasov code, the Hénon sphere was apodized as follows:

with Δ=1/2. The agreement between both codes is striking. One notices that the semi-Lagrangian code is more diffusive but less noisy than the N-body simulation, as expected. More detailed comparisons are performed in Colombi et al. (2015).

f(r,v) in a thin angular momentum slice

f(r,v) integrated over angular momentum

Figure 2: Comparison between the Vlasov code VlaSolve and the N-body code Gadget-2 for a Hénon sphere. The two top lines of panels display, from left to right, the evolution of the phase-space distribution function in a thin angular momentum slice. The two bottom lines of panels are similar, but the phase-space distribution function has been integrated over all dynamic range in j. The VlaSolve simulation (top and third line of panels) has a resolution given by (Nr,Nv,Nj)=(1024,1024,512) while the Gadget-2 simulation (second and bottom line of panels) involves 10 millions particles. From Colombi et al. (2015).
Movie (created by S. Peirani & T. Sousbie): Vlasov_versus_Nbody.mpg


  1. Cheng C. Z., Knorr G., 1976, JCoPh, 22, 330
  2. Colombi S., Touma J., 2008, CNSNS, 13, 46
  3. Crouseilles N., Latu G., Sonnendrücker E., 2009, JCoPh, 228, 1429
  4. Fujiwara T., 1983, PASJ, 35, 547
  5. Gott J.R. III, 1973, ApJ, 186, 481
  6. Hénon M., 1964, AnAp, 27, 83
  7. Springel V., 2005, MNRAS, 364, 1105
  8. Springel V., Yoshida N., White S.D.M., 2001, NewA, 6, 51