Recent Changes - Search:


Hot !

New !

Codes and methods

Theory and analysis

Past events


edit SideBar

VlaPoly: a Vlasov-Poisson solver in 1D using the waterbag method

Warning: this page is better viewed with Firefox or Safari

Authors of the code: S. Colombi & J. Touma
Publications introducing and exploiting the code (to be cited when using VlaPoly):

Obtain VlaPoly: on request to S. Colombi -> mailto:colombi [snail] iap [period] fr

Table of Contents

The waterbag method

In the proper units, Vlasov-Poisson equations are given in one dimension by

where x is the position, v the velocity, t the time, f(x,v,t) the phase-space density distribution function, Φ(x,t) the gravitational potential and ρ(x,t) the projected density.

The waterbag method (DePackh 1962; Robert & Berk 1967; Janin 1971; Cuperman, Harten & Lecar 1971a,b) exploits directly Liouville theorem, namely that the phase-space distribution function is conserved along trajectories of test particles:

This property can be exploited in a powerful way by decomposing the initial distribution on small patches, the waterbags, where f is approximated by a constant. In one dimension, the numerical implementation is thus conceptually very simple: one just needs to follow the boundaries of the waterbag with a polygon, which is enriched with new vertices when the shape of the waterbag gets more involved (Figure 1).

Figure 1: decomposition of phase-space on triangular waterbags and subsequent evolution of the phase-space distribution function. Images: S. Colombi & J. Touma.

Practical implementation

Initial conditions

In VlaPoly, the initial phase-space distribution function can be sampled initially either with an adaptive triangulation (Figure 1) or with isocontours of f(x,v) (Figure 2). In the isocontour method, which is more optimal, the isocontours are chosen such to bound the mean square difference between the true and the sampled phase-space distribution function weighted by the waterbag thickness. To draw the iscontours, the so-called marching square algorithm is used, inspired from its famous three-dimensional alter-ego (Lorensen & Cline, 1987).

Poisson equation

The equation of motion of the polygon vertices is the same as for test particles, of which the acceleration is given in one dimension by

where Mright(x,t), Mleft(x,t) and Mtot are respectively the total mass at the right of position x, the total mass at the left of position x and the total mass.

If one approximates f by a constant with value fk within a patch, Pk, k=,1...,Npatch, then

where s is a curvilinear coordinate, the last equation making use of Green's theorem. Hence, Poisson equation resolution reduces to a (partial) circulation along the contours of each individual patch that are themselves followed with an adaptive orientated polygon. The calculation of the acceleration corresponds to the most costly part of the code, as many segments of the polygon can contribute to computing the force at a given point when mixing becomes strong.

Time integration

The vertices of the adaptive polygon are moved with a standard predictor-corrector scheme, with a slowly varying time step, of which the value is mainly determined by the condition

where ρmax is the maximum value of the projected density. This criterion on the time step is derived from harmonic oscillator dynamics considerations. Two other additional constraints are added to make the time step harmonious with refinement of the orientated polygon.


The refinement procedure used in VlaPoly is quite close in spirit to Cuperman et al. (1971a). It consists of a geometric construct using arc of circles passing through sets of three successive points of the polygon and reduces, in the small angle approximation to linearly interpolating local curvature. Refinement criterion relies on local measurement of the surface of the triangle formed by the refinement candidate point and its too neighbors. If the surface of such a triangle exceeds some threshold fixed by the user, then refinement is performed. Additionally, if the distance between two successive points of the polygon exceeds some threshold fixed by the user, a new point is added. There is as well the possibility of unrefining, but this option turns to be suboptimal due to the cumulative noise it introduces in the long run.


Diagnostics, in addition to standard conservation tests on energy, volume and mass, include output of snapshots, such as the full polygon structure. This structure can be post-treated to compute various profiles as as well as images of the phase-space distribution function in (x,v) or Action-Angle space. To generate these images, we use a modified version of the parity algorithm.

Illustrations and movies

Figure 2: "Landau damping:" evolution of an initially Gaussian phase-space distribution function in (x,v) and Action-Angle space to emphasize on the quiescent nature of mixing, as studied in Colombi S., Touma J. (2014). Note: this figure is analogous to Figure 1, but waterbag borders coincide with isocontours of the phase-space distribution function.
- A Movie (same initial conditions, but different simulation): gaussian.avi

Figure 3: Mosaic of the evolution in phase-space of systems with various initial conditions. Images: S. Colombi & J. Touma.

Figure 4: Phase-space evolution of an ensemble of interacting "halos" (Colombi S., Touma J., 2014).

Figure 5: Convergence to the cold case as studied in Colombi S., Touma J. (2014), with single waterbags of decreasing thickness 2Δp.
- Left panels: tophat_1v000.avi
- Middle panels: tophat_0v100.avi
- Right panels: tophat_0v010.avi

Figure 6: The randomly perturbed nearly cold case, as studied in Colombi S., Touma J. (2014). The zooms illustrate the exquisite accuracy of the waterbag method.
- Movie: tophat_0v003_pert0v0003.avi


  1. Cuperman S., Harten A., Lecar M., 1971a, Ap&SS, 13, 411
  2. Cuperman S., Harten A., Lecar M., 1971b, Ap&SS, 13, 425
  3. DePackh D. C., 1962, J. Electr. Contrib., 13, 417
  4. Janin G., 1971, A&A, 11, 188
  5. Lorensen W.E., Cline H.E., 1987, Comput. Graph., 21, 163
  6. Roberts K.V., Berk H.L., 1967, Phys. Rev. Lett., 19, 297