Recent Changes - Search:


New !

Codes and methods

Theory and analysis

Past events


edit SideBar

ColDICE: a parallel Vlasov-Poisson solver in 4D/6D using moving adaptive simplicial tessellation

Warning: this page is better viewed with Firefox or Safari

Author of the code: T. Sousbie

Publication introducing the code (to be cited when using ColDICE):

Filmed Presentation given at CIRM in the framework during the conference "Collisionless Boltzmann (Vlasov) Equation and Modeling of Self-Gravitating Systems and Plasmas"

ColDICE and the DICE libraries are available on a Github repository.

  • Download the package by either cloning the git repository using the command

git clone --depth 1

Table of Contents

Structure of the phase-space sheet of a two sine waves simulation performed in 4-dimensional phase-space with ColDICE. The x and y axes represent the 2D position space while the velocity ux is represented along the z-axis. The velocity uy is color coded in blue and red, with a black stripe corresponding to uy=0. From Sousbie & Colombi (2016).

The Vlasov-Poisson equations in the cold case: a D-dimensional sheet evolving in 2D-dimensional phase-space

Dark matter in the Universe can be described as a smooth self-gravitating collisionless fluid following Vlasov-Poisson equations,

where f(r,u,t) represents the phase-space density at position r, velocity u and time t, Φ is the gravitational potential and G is the gravitational constant. The ``cold'' nature of dark matter implies a virtually null initial velocity dispersion, which can be expressed as follows for the initial phase-space distribution function:

where ρi corresponds to the initial projected density and ui to the initial velocity field. In the cosmological case, the initial density field is close to constant and the initial velocity field, when subtracted from the expansion of the Universe, is very small: the large scales structures of the Universe observed today grew from small initial fluctuations in the density field through gravitational instability.

Equation (1) states that the phase-space distribution function is initially concentrated on a D-dimensional phase-space sheet in 2D-dimensional phase-space. The Hamiltonian nature of the system implies that topological properties of the phase-space distribution function are conserved during motion, in particular that the phase-space sheet remains a non self-intersecting D-dimensional hyper-surface at all times.

In ColDICE, we sample the phase-space sheet with an adaptive simplicial tessellation. Its dynamical evolution in phase-space is followed accurately in a fully parallel fashion by exploiting shared (OpenMP) and distributed (MPI) parallelisms. Both 4 (D=2) and 6-dimensional (D=3) cases are considered. Note that the two-dimensional case, D=2, corresponds to the continuum of infinite lines, for which the gravitational potential is of logarithmic nature instead of the usual 1/r form in D=3 dimensions.

Our approach is analogous to the recent scheme proposed by Hahn, Abel & Kaehler (2013) and Hahn & Angulo (2015), although completely different in the actual implementation. It also follows closely ideas discussed in Shandarin, Habib & Heitmann (2012) and Abel, Hahn & Kaehler (2012), where the concept of sampling the phase-space density with a Lagrangian tessellation was introduced in the cosmological context.

A moving, conforming, adaptive, quadratic mesh

Figure 1: Example of a possible Lagrangian tessellation of the phase-space sheet into simplices in the 2D case (left) and in the 3D case (right). From Sousbie & Colombi (2016).

Figure 2: Illustration of the refinement of four tetrahedra incident to a bisected edge. A single new vertex is inserted and uniquely defines the geometry of the eight resulting refined simplices. From Sousbie & Colombi (2016).

Figure 3: Illustration of the usage of edge tracer particles to define quadratic mesh elements. On the left, first order elements with edge tracers, on the right, second order elements. A total of 6 nodes per element (3 red tracers and 3 green vertices) is required to uniquely define the quadratic triangular elements of a 2D quadratic mesh. The 3D case is similar, with a total of 10 nodes (6 tracers and 4 vertices) required to define quadratic tetrahedra. From Sousbie & Colombi (2016).

We sample the phase-space sheet with a conforming simplicial tessellation, that is with an ensemble of joint tetrahedra that coincide exactly on facets when intersecting with each other. To start with, the phase-space sheet follows equation (1) with ρi =constant and ui=0 and is sampled with a homogeneous simplicial mesh (Figure 1). Its vertices are then perturbed by applying a smooth Lagrangian displacement field in phase-space to create the initial state. In the cosmological case, this displacement field is computed using standard Lagrangian perturbation theory (see, e.g. Bernardeau et al., 2002). Once the initial conditions are set-up, the vertices of the tessellation evolve according to the Lagrangian equations of motion like particles in a N-body simulation, except that Poisson equation is solved exploiting the continuous nature of the tessellation up to linear order, as detailed in next section.

Because of mixing, it is necessary to add sampling elements when needed, i.e. to introduce local refinement on the tessellation. To do so, we use techniques analogous to finite element methods: the phase-space sheet is refined using bisection, that is by cutting when required some simplices into two smaller simplices while preserving the conforming nature of the tessellation at all times (Figure 2). New vertices created during this procedure are placed in such a way that the local representation of the phase-space sheet remains accurate at second order, which requires the use of additional tracers (Figure 3). Indeed, because a significant amount of curvature is generated during the course of dynamics, following locally the shape of the phase-space sheet at second order greatly improves the quality of the representation of the system and is nearly a must.

In order to preserve as well as possible the Hamiltonian nature of the system, the criterion we use to decide when simplices have to be refined as well as the way they are refined relies on the measurement of local Poincaré invariants, that is for of a pair of phase-space vectors (zj, zk) coinciding with two sides of the triangle, the symplectic areas defined by

with ω the symplectic matrix:

In Hamiltonian systems, such areas should be conserved during motion: similarly to VlaPoly, our refinement criterion tries to set limits to violations of this property.

The parallel adaptive simplicial mesh structure used by ColDICE is in fact quite generic and as a part of the standalone library DICE, can be employed by other algorithms. Un-refinement is implemented as well to optimise the number of simplices in the tessellation, although it is not yet used by ColDICE.

A robust and accurate projection algorithm between a simplicial mesh and a regular rectangular mesh exact to linear order

Figure 4: Illustration of the definition of the vector triplets (P,T,N) (2D, left figure) and quadruplets (P,T,N,B) (3D, right figure) associated to each polyhedron as defined by Franklin & Kankanhalli (1993). Such multiplet exists for each possible combination of incident vertex, edge and facet in the polyhedron, and there therefore exist 12 of them for the 2D polyhedron on the left figure and 36 of them for the 3D polyhedron of the right figure. From Sousbie & Colombi (2016).

To solve Poisson equation, we project the tessellation onto a regular rectangular grid, by adapting the volume calculation method of Franklin & Kankanhalli (1993) using raytracing and generalising it to first order, that is with a hypersurface density represented at linear order inside each simplex. This elegant method consists in realising that properties of any polyhedron can be computed only as functions of its vertices and, for each vertex P, sets of vectors T, N, B (Franklin 1987) defined as follows (Figure 4):

  • T is a unit vector aligned with one edge of the polyhedron,
  • N is a unit vector orthogonal to T in the plane of a facet of the polyhedron and pointing inside it,
  • B is a unit vector orthogonal to T and N and pointing inside the polyhedron.

For instance the surface/volume of a polyhedron is given by the following expressions, respectively in two and three dimensions,

where the sums are performed over all the possible combinations of quadruplets (P,T, N, B). This algorithm can be generalised to integrate any polynomial function over a polyhedron, which we do at linear order in our projection algorithm.

The difficulty in implementing such an algorithm relies in enforcing robustness and accuracy while preserving its performances. To enforce robustness, we employ simulation of simplicity as proposed by Edelsbrunner & Mucke (1990). To warrant a given level of accuracy, we use multi-precision floating point number types, for instance double-double types, when needed. The regions where higher precision arithmetic is needed are found by a posteriori tests of accuracy on a per-voxel basis.

To speed up the process at a small cost in accuracy, we use a standard adaptive mesh refinement technique to construct the voxels/pixels grid in such a way that the size of the mesh elements is locally of the same order of that of the simplices under consideration. In this sense, our calculation of the intersections is exact at linear order but only at scales comparable to the size of simplices. Gathering the information on the final fixed resolution grid is performed by a simple donor cell procedure. Then force calculation and vertex position and velocity updates are performed just as in standard Particle-Mesh codes (see, e.g., Hockney & Eastwood 1981).

Note that our fully parallel projecting code is a component of the standalone library DICE and can be used in other algorithms needing the calculation of the projection of a simplicial grid onto a rectangular one.


  1. Abel T., Hahn O., Kaehler R., 2012, MNRAS 427, 61
  2. Bernardeau F., Colombi S., Gaztañaga E., Scoccimarro R., 2002, Physics Reports 367, 1
  3. Edelsbrunner H., Mucke E. P., 1990, Simulation of simplicity: a technique to cope with degenerate cases in geometric algorithms, ACM Transactions on Graphics (TOG) 9, 66
  4. Franklin W. R., 1987, Polygon properties calculated from the vertex neighborhoods, Proceedings of the third annual symposium on Computational geometry, ACM, p. 110
  5. Franklin W. R., Kankanhalli M. S., 1993, Volumes from overlaying 3-D triangulations in parallel, in Advances in Spatial Databases, Springer, p. 447
  6. Hahn O., Abel T. Kaehler R., 2013, MNRAS 434, 1171
  7. Hahn O., Angulo R. E., 2015, arXiv e-prints 1501.01959
  8. Hockney R. W., Eastwood J. W., 1981, Computer Simulation Using Particles, New York: McGraw-Hill
  9. Shandarin S., Habib S., Heitmann K., 2012, Physical Review D 85, 083005

Illustrations and movies

Chaotic potential in 4-dimensional phase-space

Evolution in time of the projected density produced by an initially small square patch in a chaotic potential. The initial tessellation mesh of the patch is shown plotted over the projected density as a zoom on the first figure. From Sousbie & Colombi (2016).
Movie: chaosColored.mkv

Sine waves in 4-dimensional cosmological phase-space

Projected density at various times of a cosmological two sinusoidal waves simulation. A zoom in the region [5L/16, 11L/16] where L is the simulation box size has been performed. From Sousbie & Colombi (2016).

Cosmological warm dark matter simulation in 6-dimensional phase-space

Projected density in a cosmological warm dark matter simulation, integrated along the x-axis. From Sousbie & Colombi (2016).

Phase-space structure of a halo in the warm dark matter simulation of the previous figure. The image on the bottom of each picture represents a 1 voxel thick slice perpendicular to the z axis and going through the center of the halo of the projected density field used in the simulation. The halo itself is represented in the (x,y,vx) sub-space, as the projection of a wire corresponding to a uniform grid in Lagrangian space (left panel) and Lagrangian subsets initially cubical (right panel). From Sousbie & Colombi (2016).

Movies, available as well on YouTube (Author of the movies: T. Sousbie)