A "metric" semi-Lagrangian Vlasov-Poisson solver

TL;DR

This paper introduces a novel semi-Lagrangian Vlasov-Poisson solver that uses metric elements to accurately track phase-space flow, enabling precise reconstruction of the distribution function in gravitational and plasma dynamics.

Contribution

The paper presents a new metric-based semi-Lagrangian approach that improves phase-space tracking and reconstruction in Vlasov-Poisson simulations, with potential for high-dimensional applications.

Findings

Accurate phase-space reconstruction using second-order geometry expansion.

Efficient remapping with higher-order splines during transport.

Algorithm successfully tested in 1D gravitational dynamics, extendable to higher dimensions.

Abstract

We propose a new semi-Lagrangian Vlasov-Poisson solver. It employs elements of metric to follow locally the flow and its deformation, allowing one to find quickly and accurately the initial phase-space position $Q(P)$ of any test particle $P$, by expanding at second order the geometry of the motion in the vicinity of the closest element. It is thus possible to reconstruct accurately the phase-space distribution function at any time $t$ and position $P$ by proper interpolation of initial conditions, following Liouville theorem. When distorsion of the elements of metric becomes too large, it is necessary to create new initial conditions along with isotropic elements and repeat the procedure again until next resampling. To speed up the process, interpolation of the phase-space distribution is performed at second order during the transport phase, while third order splines are used at the moments of remapping. We also show how to compute accurately the region of influence of each element of metric with the proper percolation scheme. The algorithm is tested here in the framework of one-dimensional gravitational dynamics but is implemented in such a way that it can be extended easily to four or six-dimensional phase-space. It can also be trivially generalised to plasmas.

Figures

Tables

Table 1. Table 1: Parameters of the simulations performed in this work, namely the grid cell size, Δ = Δ ​ x = Δ ​ y Δ Δ 𝑥 Δ 𝑦 \Delta=\Delta x=\Delta y , the initial separation between each metric element, Δ m = Δ ​ x m = Δ ​ y m subscript Δ m Δ subscript 𝑥 m Δ subscript 𝑦 m \Delta_{\rm m}=\Delta x_{\rm m}=\Delta y_{\rm m} , the number of sub-cycles, n s subscript 𝑛 s n_{\rm s} , between resamplings and the approximate total CPU time expressed in units of the time spent for simulation I. The calculation of t CPU subscript 𝑡 CPU t_{\rm CPU} was performed for the simulations with initial conditions corresponding to the random set of stationary clouds as described in §  3.2 , but the results do not change drastically for other initial conditions. A horizontal line separates the runs performed with Vlamet (I–XI) from those performed with the standard splitting algorithm (XII–XVI). Additionally, the time step Δ ​ t Δ 𝑡 \Delta t was chosen constant, given by Δ ​ t = 0.01 Δ 𝑡 0.01 \Delta t=0.01 , 0.01 0.01 0.01 and 0.0025 0.0025 0.0025 for the stationary, Gaussian and random set of halos, respectively. The two last choices come from the constraints on the time step obtained in the waterbag runs by Colombi & Touma ( 2014 ) .

Designation	$\Delta$	${\Delta }_m$	$n_s$	$t_{\mathrm{CPU}}$
I	0.002	0.02	25	1.0
II	0.002	0.04	25	1.0
III	0.002	0.01	25	1.2
IV	0.002	0.005	25	1.4
V	0.002	0.02	100	1.1
VI	0.002	0.02	50	1.1
VII	0.002	0.02	10	1.2
VIII	0.002	0.02	5	1.5
IX	0.002	0.02	1	3.2
X	0.001	0.01	25	4.0
XI	0.004	0.04	25	0.4
XII	0.0005	–	–	8.8
XIII	0.00707	–	–	4.5
XIV	0.001	–	–	2.3
XV	0.00141	–	–	1.2
XVI	0.002	–	–	0.7

Equations

$$\frac{\partial f}{\partial t}+\mathbf{v}.{\nabla }_{\mathbf{x}}f-{\nabla }_{\mathbf{x}}\varphi .{\nabla }_{\mathbf{v}}f=0,$$

$${\Delta }_{\mathbf{x}}\varphi =4\pi G\rho ,\rho (\mathbf{x},t)\equiv \int f(\mathbf{x},\mathbf{v},t)\mathrm{d}\mathbf{v},$$

$$\mathbf{P}\equiv (\mathbf{x},\mathbf{v}).$$

$$\mathbf{Q}\equiv (q_x,q_v),$$

$$\mathbf{P}(\mathbf{Q},t)\equiv [x(q_x,q_v,t),v(q_x,q_v,t)],$$

$$\mathbf{P}(\mathbf{Q},t=0)=\mathbf{Q}.$$

$$\mathbf{T}(\mathbf{Q},t)\equiv \frac{\partial \mathbf{P}}{\partial \mathbf{Q}},$$

$$\mathbf{H}(\mathbf{Q},t)\equiv \frac{{\partial }^2\mathbf{P}}{\partial {\mathbf{Q}}^2}.$$

$$\mathbf{P}\simeq {\mathbf{P}}_{\mathrm{m}}+{\mathbf{T}}_{\mathrm{m}}(\mathbf{Q}-{\mathbf{Q}}_{\mathrm{m}})+\frac{1}{2}(\mathbf{Q}-{\mathbf{Q}}_{\mathrm{m}}{)}^{\mathrm{T}}{\mathbf{H}}_{\mathrm{m}}(\mathbf{Q}-{\mathbf{Q}}_{\mathrm{m}}).$$

$$\mathbf{Q}$$

$$|\mathbf{T}|=1,$$

$$\mathbf{T}(\mathbf{Q},t=0)=\mathbf{I},$$

$$J\equiv |\mathbf{T}|=\frac{\partial x}{\partial q_x}\frac{\partial v}{\partial q_v}-\frac{\partial x}{\partial q_v}\frac{\partial v}{\partial q_x}=1,$$

$$\mathcal{I}=\frac{1}{2}[\mathrm{d}{\mathbf{P}}_1{]}^{\mathrm{T}}\mathbf{S}\mathrm{d}{\mathbf{P}}_2,$$

$${\boldsymbol{S}}\equiv\left(\begin{array}[]{cc}{\boldsymbol{0}}&-{\boldsymbol{I}}\\ {\boldsymbol{I}}&{\boldsymbol{0}}\end{array}\right)$$

$$\mathrm{d}{\mathbf{P}}_j=\frac{\partial \mathbf{P}}{\partial \mathbf{Q}}\mathrm{d}{\mathbf{Q}}_j.$$

$${\left[\frac{\partial \mathbf{v}}{\partial {\mathbf{q}}_x}\right]}^{\mathrm{T}}\frac{\partial \mathbf{x}}{\partial {\mathbf{q}}_x}-{\left[\frac{\partial \mathbf{x}}{\partial {\mathbf{q}}_x}\right]}^{\mathrm{T}}\frac{\partial \mathbf{v}}{\partial {\mathbf{q}}_x}=0,$$

$${\left[\frac{\partial \mathbf{v}}{\partial {\mathbf{q}}_v}\right]}^{\mathrm{T}}\frac{\partial \mathbf{x}}{\partial {\mathbf{q}}_v}-{\left[\frac{\partial \mathbf{x}}{\partial {\mathbf{q}}_v}\right]}^{\mathrm{T}}\frac{\partial \mathbf{v}}{\partial {\mathbf{q}}_v}=0,$$

$${\left[\frac{\partial \mathbf{x}}{\partial {\mathbf{q}}_x}\right]}^{\mathrm{T}}\frac{\partial \mathbf{v}}{\partial {\mathbf{q}}_v}-{\left[\frac{\partial \mathbf{v}}{\partial {\mathbf{q}}_x}\right]}^{\mathrm{T}}\frac{\partial \mathbf{x}}{\partial {\mathbf{q}}_v}=\mathbf{I},$$

$${\mathbf{Q}}_{\mathrm{G}}=(i\Delta x+q_{x,\mathrm{mid}},j\Delta v+q_{v,\mathrm{mid}}),$$

$$f(q_x,q_v)$$

$$(i_{\mathrm{n}},j_{\mathrm{n}})$$

$$(i_{\mathrm{c}},j_{\mathrm{c}})$$

$$\mathrm{d}{q}_x$$

$$\mathrm{d}{q}_v$$

$$|f|>f_{\mathrm{th}},$$

$$f(q_x,q_v)$$

$${\mathbf{Q}}_{\mathrm{m}}(i_{\mathrm{m}},j_{\mathrm{m}})=(\Delta x_{\mathrm{m}}{i}_{\mathrm{m}}+q_{x,\mathrm{mid}},\Delta v_{\mathrm{m}}{j}_{\mathrm{m}}+q_{v,\mathrm{mid}}),$$

$$\Delta x_{\rm m}\mathrel{\mathchoice{\vbox{\offinterlineskip\halign{\hfil$\displaystyle#$\hfil\cr>\cr\sim\cr}}}{\vbox{\offinterlineskip\halign{\hfil$\textstyle#$\hfil\cr>\cr\sim\cr}}}{\vbox{\offinterlineskip\halign{\hfil$\scriptstyle#$\hfil\cr>\cr\sim\cr}}}{\vbox{\offinterlineskip\halign{\hfil$\scriptscriptstyle#$\hfil\cr>\cr\sim\cr}}}}\Delta x,\quad\Delta v_{\rm m}\mathrel{\mathchoice{\vbox{\offinterlineskip\halign{\hfil$\displaystyle#$\hfil\cr>\cr\sim\cr}}}{\vbox{\offinterlineskip\halign{\hfil$\textstyle#$\hfil\cr>\cr\sim\cr}}}{\vbox{\offinterlineskip\halign{\hfil$\scriptstyle#$\hfil\cr>\cr\sim\cr}}}{\vbox{\offinterlineskip\halign{\hfil$\scriptscriptstyle#$\hfil\cr>\cr\sim\cr}}}}\Delta v.$$

$$(i,j)=\left(\lfloor (x_{\mathrm{m}}-x_{\mathrm{mid}})/\Delta x\rceil ,\lfloor (v_{\mathrm{m}}-v_{\mathrm{mid}})/\Delta v\rceil \right),$$

$$\delta \mathbf{Q}(i,j,i_{\mathrm{m}},j_{\mathrm{m}})$$

$$\equiv$$

$$(i_{\rm m}^{\prime},j_{\rm m}^{\prime}):\left\{\begin{array}[]{ll}i_{\rm m}^{\prime}\in[i_{\rm m}-1,i_{\rm m}+1],\\ j_{\rm m}^{\prime}\in[j_{\rm m}-1,j_{\rm m}+1],\\ (i_{\rm m}^{\prime},j_{\rm m}^{\prime})\neq(i_{\rm m},j_{\rm m}),\end{array}\right.$$

Full text

1 Introduction

The dynamics of Dark matter in the Universe and stars in galaxies is governed by Vlasov-Poisson equations

[TABLE]

where $f({\boldsymbol{x}},{\boldsymbol{v}},t)$ is the phase-space distribution function at position ${\boldsymbol{x}}$, velocity ${\boldsymbol{v}}$ and time $t$, $\phi$ the gravitational potential, $G$ the gravitational constant and $\rho$ the projected density. Usually, these equations are solved numerically with a $N$-body approach, where the phase-space distribution function is represented by an ensemble of Dirac functions interacting gravitationally with a softened potential (see, e.g., Hockney and Eastwood, 1988; Bertschinger, 1998; Colombi, 2001; Dolag et al., 2008; Dehnen and Read, 2011, and references therein). It is known that such a rough representation of the phase-space fluid is not free of biases nor systematic effects (see, e.g., Aarseth et al., 1988; Kandrup and Smith, 1991; Boily et al., 2002; Binney, 2004; Melott, 2007; Joyce et al., 2009; Colombi et al., 2015). With the computational power available today, an alternative approach is becoming possible, consisting in solving directly Vlasov dynamics in phase-space (see, e.g., Yoshikawa et al., 2013). In fact, this point of view is already commonly adopted in plasma physics (see, e.g., Grandgirard et al., 2016, and references therein), but seldom used in the gravitational case.

In this paper, we consider the case where the phase-space distribution is a smooth function of the vector

[TABLE]

This means that initial conditions have non zero initial velocity dispersion, which is relevant to the analysis of e.g. galaxies or relaxed dark matter halos. On the other hand, if the aim would be to follow accurately the dynamics of dark matter in the Cold Dark Matter paradigm, the initial velocity would be null. In this case, the problem is different since one has to describe the foldings of a $D$-dimensional hypersurface evolving in 2$D$ dimensional phase-space (see, e.g., Hahn and Angulo, 2016; Sousbie and Colombi, 2016).

In the warm case considered here, it is possible to represent $f({\boldsymbol{P}},t)$ on an Eulerian grid. Most of the numerical methods trying to resolve directly the dynamics of function $f({\boldsymbol{P}},t)$ have been developed in plasma physics and are of semi-Lagrangian nature. They usually exploit directly Liouville theorem, namely that $f$ is conserved along characteristics. In practice, a test particle is associated to each grid site where the value $f_{\rm T}$ of the phase-space distribution function has to be evaluated. This test particle is followed backwards in time to find its position at previous time step. Then, function $f$ from previous time step is interpolated at the root of the characteristic to obtain the searched value of $f_{\rm T}$, following Liouville theorem. This approach was pioneered by Cheng and Knorr (1976) in a split fashion and first applied to astrophysical systems by Fujiwara (1981), Nishida et al. (1981) and Watanabe et al. (1981). In the classical implementation, interpolation of the phase-space distribution function is performed using a third order spline.

There exist many subsequent improvements over the splitting algorithm of Cheng & Knorr (see, e.g., Shoucri and Gagne, 1978; Sonnendrücker et al., 1999; Filbet et al., 2001; Besse and Sonnendrücker, 2003; Crouseilles et al., 2010; Rossmanith and Seal, 2011; Qiu and Shu, 2011; Güçlü et al., 2014, but this list is far from comprehensive). The main problem common to all these algorithms is the unavoidable diffusion due to repeated resampling of the phase-space distribution function. Indeed, because of the finite grid resolution, repeated interpolation, as sophisticated it can be, induces information loss and augmentation of entropy. One way to fix this problem is to employ adaptive mesh refinement (see, e.g., Gutnic et al., 2004; Besse et al., 2008, 2017), i.e. to augment local resolution when the level of details asks for it. However, diffusion can also be significantly dampened if remapping of the phase-space distribution function is performed as seldom as possible. One way to achieve this is to represent $f({\boldsymbol{P}},t)$ with an ensemble of elements of which the shape is allowed to change with time according to the Lagrangian equations of motion (see, e.g., Alard and Colombi, 2005; Campos Pinto et al., 2014; Larson and Young, 2015). In the “cloudy” method that we proposed in Alard and Colombi (2005), the phase-space distribution is sampled on an ensemble of overlapping round Gaussian functions of which the contours tranform into ellipsoids following the phase-space flow. When the ellipsoids become too elongated, the phase-space distribution function is resampled with new round clouds. Because resampling is seldom performed, say at fractions of dynamical times, diffusion effects are much less prominent, even though coarse graining is still performed at the resolution scale given by the inter-cloud center spacing. The problem with this remapping approach is the computational cost: to have a smooth and accurate representation of the phase-space distribution function, the clouds need to overlap significantly: the extra cost due to this overlap becomes prohibitive in high number of dimensions.

However, one can realize in fact that it is not necessary to represent the phase-space distribution on a set of clouds: the only important piece of information is the metric intrinsically associated to each cloud of the previous approach, i.e. the information on how function $f({\boldsymbol{P}},t)$ is deformed by the motion. If one is indeed able to follow the properties of the motion up to some order in the neighbourhood of test particles, that we call metric elements, reconstructing the initial position ${\boldsymbol{Q}}({\boldsymbol{P}},t)$ of any point nearby is straightforward, and so is the calculation of $f({\boldsymbol{P}},t)$ from interpolation of initial conditions following Liouville theorem.

Our algorithm, Vlamet, described in details in § 2, can thus be summarised as follows. It consists in sampling at all times the phase-space distribution function $f({\boldsymbol{P}},t)$ on a grid of fixed resolution and following $f$ at second order in space and time using a predictor-corrector time integration scheme. The calculation of $f({\boldsymbol{P}},t)$ exploits directly Liouville theorem, namely that $f$ is conserved along the characteristics, i.e. is equal to its initial value $f_{\rm ini}\equiv f({\boldsymbol{Q}},t_{\rm ini})$ at Lagrangian position ${\boldsymbol{Q}}({\boldsymbol{P}},t)$. To follow these characteristics accurately, we define a set of elements of metric that are sparse sampling the grid. These elements of metric follow the Lagrangian equations of motion and carry during time local information about the deformation of phase-space up to second order, to account for the curvature of phase-space structures building up during the course of the dynamics. To resample $f$ at time $t$ and at a given point ${\boldsymbol{P}}$ of phase-space, we consider the closest elements of metric in Lagrangian space and reconstruct the initial position ${\boldsymbol{Q}}$ of a test particle coinciding with ${\boldsymbol{P}}$ using a second order expansion of the geometry of the flow around the metric elements. Once the initial position ${\boldsymbol{Q}}$ is found, we just need to interpolate $f_{\rm ini}$ at position ${\boldsymbol{Q}}$ and attribute it to $f({\boldsymbol{P}},t)$.

At some time $t_{\rm resample}$, the distorsion of phase-space is too strong to be simply represented at second order: $f({\boldsymbol{P}},t_{\rm resample})$ is considered as a new initial condition and new isotropic elements of metric are created. Between initial time and $t_{\rm resample}$, we use a cheap, second order interpolation to reconstruct $f({\boldsymbol{P}},t)$, while resampling at $t=t_{\rm resample}$ is more accurate, using third order splines. Indeed, although our algorithm is meant to be globally precise only at second order, a third order description of the phase-space distribution function is essential at the moment of resampling to be able to conserve quantities such total energy and total mass in the long term. Furthermore, to preserve smoothness of the resampling procedure, we perform some interpolation in Lagrangian space between the initial positions ${\boldsymbol{Q}}$ proposed by each element of metric close to the point of interest. Indeed, different neighbouring metric elements can provide different answers. The procedure is illustrated by Fig. 1.

In § 3, we thoroughly test the performances of our algorithm, by checking how a stationary solution is preserved and by simulating systems with Gaussian initial condition as well as random sets of halos. To demonstrate the potential improvements brought by our code Vlamet compared to traditional semi-Lagrangian approaches, we compare our algorithm to the traditional splitting method of Cheng and Knorr (1976) with third order spline interpolation. Results are also tested against simulations performed with the entropy conserving waterbag code of Colombi & Touma (2014). The waterbag scheme, extremely accurate but very costly, is meant to provide a supposedly “exact” solution. It is mainly used to test the convergence of Vlamet and of the splitting algorithm for the Gaussian initial condition. From the pure mathematical side, which is not really approached in this article, it is worth mentioning the recent work of Campos Pinto and Charles (2016), who discuss some convergence properties of a simplified version of our algorithm.

2 Algorithm

Now we describe the main steps of the algorithm, namely, (i) our second order representation of the metric (§ 2.1), (ii) our interpolation schemes to reconstruct $f$ (§ 2.2), (iii) our percolation algorithm to find the nearest element of metric to an Eulerian pixel ${\boldsymbol{P}}$ and the reconstruction of its Lagrangian counterpart ${\boldsymbol{Q}}$ using an interpolation between the alternatives proposed by the neighbouring elements of metric (§ 2.3), (iv) the equations of motion as well as their numerical implementation (§ 2.4) and finally, (v) the calculation of the force, its first and second derivatives, required to follow the second order representation of the local metric during motion (§ 2.5).

While the actual implementation of the algorithm is performed in 2D phase-space, generalisation to higher number of dimensions will be discussed in details. Indeed, our algorithm is meant to improve other standard semi-Lagrangian solvers in terms of diffusion while being of comparable cost, and is designed in such a way that it can be easily generalised to 4- and 6-dimensional phase-space.

2.1 Second order representation of the metric

The Lagrangian position of each fluid element in phase space is defined by

[TABLE]

where $q_{x}$ and $q_{v}$ are its initial position and velocity, respectively. Here we consider 2 dimensional phase-space, so these quantities are scalar, but in general, they are vectors. In what follows, unless specified otherwise, all the calculations apply to the general case: one just has to replace $q_{x}$ and $q_{v}$ with ${\boldsymbol{q}}_{x}$ and ${\boldsymbol{q}}_{v}$. The Eulerian position at time $t$ of an element of fluid in phase-space is given by

[TABLE]

where $x(q_{x},q_{v},t)$ and $v(q_{x},q_{v},t)=\partial x/\partial t$ are its position and velocity at time $t$. Initially

[TABLE]

In our second order approach, we need the deformation tensor,

[TABLE]

as well as the Hessian of the characteristics,

[TABLE]

With the knowledge at all times of the position ${\boldsymbol{P}}_{\rm m}$, the deformation tensor ${\boldsymbol{T}}_{\rm m}$ and the Hessian ${\boldsymbol{H}}_{\rm m}$ for each element of metric, we can predict, at second order, the position of an element of fluid initially sufficiently close to an element of metric:

[TABLE]

Of course the series expansion could be developed further, but we restrict here to an algorithm at second order in space and time. Indeed, going beyond second order would be very costly in higher number of dimensions, particularly 6-dimensional phase-space. The interesting bit, as we shall see later, is the inverse relation between Eulerian position and Lagrangian position, i.e., at second order

[TABLE]

In the two dimensional phase-space case considered here, the deformation tensor ${\boldsymbol{T}}$ is a 4 elements matrix with no special symmetry except that according to the Hamiltonian nature of the system its determinant should be conserved with time and equal to the initial value:

[TABLE]

at all times. Indeed

[TABLE]

where ${\boldsymbol{I}}$ is the identity matrix. Note that this property is general, as volume is conserved in phase-space, it is not specific to the one dimensional dynamical case. We have,

[TABLE]

hence ${\boldsymbol{T}}$ has only 3 independent elements. In practice we shall not impose equation (13) but checking to which extent it is verified can be used as a test of the good numerical behaviour of the code, or equivalently, as a test of its symplecticity.

In higher number of dimensions, phase-space volume conservation is not the only constraint. In fact it is a consequence of the symplectic nature of the system, which hence conserves Poincaré invariants, in particular symplectic forms of the kind

[TABLE]

where

[TABLE]

is the symplectic matrix, while ${\rm d}{\boldsymbol{P}}_{1}$ and ${\rm d}{\boldsymbol{P}}_{2}$ correspond to 2 sides of an elementary triangle composed of three particles following the equations of motion in phase-space. In particular,

[TABLE]

So the generalisation of equation (13) to higher number of dimensions provides the following constraints:

[TABLE]

which can be again used to double check, a posteriori, the good behaviour of the code in the general case.

In the 2-dimensional phase-space case considered here, the Hessian ${{{\boldsymbol{H}}}}$ is a $2\times 2\times 2$ tensor in which only 6 terms are independant if one just takes into account the symmetries in second derivatives, e.g. $\partial^{2}x/(\partial q_{x}\partial q_{v})=\partial^{2}x/(\partial q_{v}\partial q_{x})$. In 4- and 6- dimensional phase-space, the same symmetries imply that among its $4\times 4\times 4=64$ and $6\times 6\times 6=$ 216 elements, the Hessian includes 40 and 126 independant terms, respectively.

While exploiting the Hamiltonian nature of the system would allow us to carry less information, we chose for each element of metric to store the phase-space coordinates, the deformation tensor and the Hessian. The calculations above show that this ends up in carrying 2+4+6=12, 4+16+40=60 and 6+36+126=168 numbers for each element of metric in 2-, 4- and 6-dimensional phase-space, respectively. This might look a large number, particularly in 4- and 6-dimensional phase-space, but one has to keep in mind the fact that elements of metric are much sparser than the grid used to compute the phase-space distribution function. For instance, a factor 3 in each direction of the spacing between elements of metric compared to the grid cell size implies that there is one element of metric for $3^{6}=729$ grid sites. Another issue is the cost of applying the inverse mapping (10), but it can be reduced considerably by precomputing for each element of metric the matrixes ${\boldsymbol{T}}_{\rm m}^{-1}$ and $[{\boldsymbol{T}}_{\rm m}^{-1}]^{\rm T}{\boldsymbol{H}}_{\rm m}{\boldsymbol{T}}_{\rm m}^{-1}$ prior to any interpolation of the phase-space distribution function using the inverse mapping.

2.2 Interpolation schemes

We use two interpolation schemes. Between two coarse steps, a period of time during which the phase-space distribution function is transported along motion using the same elements of metric at each time step, we use a simple interpolation of $f$ based on local fitting of a quadratic hypersurface. Then, during remapping, i.e. when new elements of metric are set, $f$ is interpolated with a third order spline in order to have a more accurate reconstruction of phase-space. The advantage of using second order interpolation during the transport phase is that generalisation to 6-dimensional phase-space is relatively cheap, since only 28 terms contribute in this case to such an interpolation. The loss of accuracy due to the relative roughness of the reconstruction, which is not even continuous, has little consequence on the dynamics because most of the information is recovered during the resampling phase. This will be tested and demonstrated in § 3. On the other hand, because our method does not allow for splitting, $4^{6}=4096$ terms contribute to each grid site when resampling $f$ in 6 dimensions with the third order spline interpolation, if one assumes that this latter takes about 4 operation per grid site in one dimension. This rather costly step is thus roughly $4096/24\simeq 171$ times more expensive than what is expected when performing a full time step in the standard splitting algorithm.111 This calculation assumes that third order spline interpolation is used in the splitting algorithm. In a split fashion, this interpolation is performed individually on each axis, which corresponds, in 6 dimension, to $6\times 4=24$ operations per grid site for a full time step, or $9\times 4=36$ operations if synchronisation between velocities and positions is performed. It is however performed rarely, hence making the algorithm of potentially comparable cost to the standard semi-Lagrangian approach, with the advantage of being, in principle, less diffusive.

Now we provide the details for the second order interpolation we use, followed by the standard third order spline interpolation.

2.2.1 Second order interpolation

The general procedure for second order interpolation consists in understanding that a second order hypersurface requires respectively 6, 15 and 28 points in 2D, 4D and 6D phase-space. What we have access to are the values $F({\boldsymbol{G}})$ of the phase-space distribution function on a grid ${\boldsymbol{G}}$, e.g., in one dimension ${\boldsymbol{G}}=(i,j)$ where $i$ and $j$ are integers, corresponding to the following Lagrangian positions in phase space,

[TABLE]

where $\Delta x$ and $\Delta v$ are the steps chosen for the grid (Eulerian or Lagrangian) that remain, in the current version of the code, fixed during runtime and $(q_{x,{\rm mid}},q_{v,{\rm mid}})$ corresponds to the middle of the Lagrangian computational domain. The bounds of the computational domain are then related to maximum values $n_{x}$ and $n_{v}$ of $|i|$ and $|j|$.

This is the way we proceed in 2D phase-space: given a point ${\boldsymbol{Q}}=(q_{x},q_{v})$ in phase space, we first find the nearest grid point $(i_{\rm n},j_{\rm n})$ to it, then consider the cross composed of $(i_{\rm n}+1,j_{\rm n})$, $(i_{\rm n},j_{\rm n}+1)$, $(i_{\rm n}-1,j_{\rm n})$, $(i_{\rm n},j_{\rm n}-1)$ and $(i_{\rm n},j_{\rm n})$ on which the interpolation must coincide with $F$. The last point $(i,j)$ needed is taken to be the one for which a square composed of 3 points of the cross and $(i,j)$ contains $(q_{x},q_{v})$ (see Fig. 2). With this choice of the points of the grid for the interpolation, we have, at second order,

[TABLE]

with

[TABLE]

where $\lfloor.\rceil$ and $\lfloor.\rfloor$ are respectively the nearest integer and the integer part functions. Note obviously that when defining the edges of the computing domain, one makes sure that the interpolation is defined everywhere. In the present implementation, $f$ is supposed to be different from zero on a compact domain in phase-space. To avoid artificial enlarging of the computing domain during runtime due to diffusion and aliasing, the computing domain is restricted to a rectangular region containing all the points where

[TABLE]

with $f_{\rm th}$ a small number, typically $f_{\rm th}=10^{-6}$. Note the absolute value in equation (26), because our algorithm does not enforce positivity of $f$ and better total mass conservation is obtained by keeping the regions where $f<0$.

The generalisation of this second order interpolation to higher number of dimensions is straightforward: the most difficult is to chose the grid points through which the hypersurface must pass. To achieve this, similarly as in 2 dimensions, we first find the nearest grid point ${\boldsymbol{G}}_{\rm n}$ to ${\boldsymbol{Q}}$. Then take the collection of all segments centered on ${\boldsymbol{G}}_{\rm n}$ and composed of 3 points of the grid (the cross defined above in 2 dimensions): this provides us in total 9 and 13 points respectively in 4 and 6 dimensions. Then, project ${\boldsymbol{Q}}$ on the planes defined by each pair of such segments and find the corresponding grid site in this plane such that the square composed of this grid site plus 3 contained in the 2 segments contains the projection of ${\boldsymbol{Q}}$ and add it to our list of points defining the second order hypersurface. There are 6 and 15 such points in 4 and 6 dimensions, respectively, and we thus end up using this procedure with respectively 15 and 28 points as needed to interpolate our hypersurface, with straightforward generalisation of equation (LABEL:eq:int1) to perform the interpolation.

2.2.2 Standard third order spline

The third order interpolation we use employs standard cubic B-splines as very clearly described in e.g. Crouseilles et al. (2009), so we do not provide the details here, except the edge conditions, which are

[TABLE]

and this can be trivially generalised to 4- or 6-dimensional phase-space.

2.3 Tracing back the characteristics in Lagrangian space: percolation algorithm

To determine the value of $f({\boldsymbol{P}},t)$, we have to trace the trajectory of a test particle coinciding with ${\boldsymbol{P}}$ back to its Lagrangian, initial position ${\boldsymbol{Q}}$ (the position at time of remapping, in general). To do this, we need to find the closest element of metric to this test particle and use equation (10) to compute ${\boldsymbol{Q}}$. However, the closest element of metric, from the dynamical point of view, has to be found in Lagrangian space, which complicates the procedure. Furthermore, two close test particles might end up associated to different elements of metric that will propose slightly different mappings, implying that function ${\boldsymbol{Q}}({\boldsymbol{P}})$ is discontinuous if no supplementary treatment is performed. While this is not a big deal during the transport phase, it can have dramatic consequences on the long term during the remapping step which requires high level of precision when interpolating the phase-space distribution function with third order splines. Hence, each time a full remapping of $f$ is performed, one needs to find not one but several elements of metric in the neighbourhood of a test particle to be able to interpolate between the proposed Lagrangian positions provided by each of these metric elements so that function ${\boldsymbol{Q}}({\boldsymbol{P}})$ stays smooth: this is the most complex part of the algorithm, because enforcing full continuity of the interpolation is actually not trivial since it has to be performed in Lagrangian space. The percolation algorithm however allows us to define in an unambiguous way regions of influence of each element of metric which in turn will allow us to make function ${\boldsymbol{Q}}({\boldsymbol{P}})$ smooth.

We assume that elements of metric are initially set on a Lagrangian grid covering the computing domain, ${\boldsymbol{M}}=(i_{\rm m},j_{\rm m})$, corresponding to Lagrangian positions

[TABLE]

with

[TABLE]

The actual computing domain changes subsequently with time. To estimate it at each time step we use the positions of the metric elements and we define a rectangle with lower-left and upper right corners respectively given by $[{\rm min}_{i_{\rm m},j_{\rm m}}(x_{\rm m}),{\rm min}_{i_{\rm m},j_{\rm m}}(v_{\rm m})]$ and $[{\rm max}_{i_{\rm m},j_{\rm m}}(x_{\rm m}),{\rm max}_{i_{\rm m},j_{\rm m}}(v_{\rm m})]$. This means that the element of metric coverage should be wide enough to avoid missing regions where $f$ should be non equal to zero. We therefore allow for an extra layer of metric elements on all the sides of the Lagrangian computing domain. Note thus that the calculation of the computing domain is far from being optimal, which can have dramatic consequences for the cost of the algorithm in 6 dimensions. In order to fix this problem, one would need to define a more complex shape dependent computational domain by e.g. constructing a k-d tree, i.e. a spatial decomposition of phase-space with nested cubes. Such a k-d tree could also constitute the basis for an adaptive mesh refinement algorithm following in an optimal way not only the computational domain, but also all the details of the phase-space distribution function when needed. However such a sophistication is not needed in our 2-dimensional implementation so we stick to the rectangular computing domain for now.

2.3.1 Percolation

To compute the region of influence of each element of metric $(i_{\rm m},j_{\rm m})$ with coordinates $(x_{\rm m},v_{\rm m})$, we percolate around it on the Eulerian computational grid starting from its nearest grid point

[TABLE]

with $(x_{\rm mid},v_{\rm mid})$ being the time dependent position of the middle of the Eulerian computational domain. The initial Lagrangian position ${\boldsymbol{Q}}(i,j,i_{\rm m},j_{\rm m})$ of this point is estimated with equation (10) using the element of metric we started from as a reference. Note that another neighbouring element of metric $(i_{\rm m}^{\prime},j_{\rm m}^{\prime})$ would propose a slightly different value for the initial position of the test particle associated to Eulerian grid site $(i,j)$, hence the dependence on $(i_{\rm m},j_{\rm m})$ of ${\boldsymbol{Q}}$. Then one defines the function

[TABLE]

where ${\boldsymbol{Q}}_{\rm m}(i_{\rm m},j_{\rm m})$ is the Lagrangian position of the element of metric, to find out which neighbours $(i_{\rm m}^{\prime},j_{\rm m}^{\prime})$ of the element of metric are susceptible to be closer to the true initial ${\boldsymbol{Q}}(i,j)$ in Lagrangian space. Because the calculation of ${\boldsymbol{Q}}(i,j)$ is not exact we need to introduce the concept of uncertainty. Considering all the neighbours of the current element of metric under consideration

[TABLE]

we are able to chose the one which proposes the smallest value of

[TABLE]

where the dependence on $(i,j,i_{\rm m},j_{\rm m})$ of $\delta q_{x}$ and $\delta q_{v}$ has been made implicit. If there is a candidate element of metric $(i_{\rm m}^{\prime},j_{\rm m}^{\prime})$ closer than $(i_{\rm m},j_{\rm m})$ to the Eulerian grid element $(i,j)$ using this Lagrangian distance, the percolation process is stopped. Otherwise, one continues (bond) percolation i.e. proceeds with neighbouring grid sites $(i+1,j)$, $(i-1,j)$, $(i,j+1)$, $(i,j-1)$ unless they have been already processed.

Note that the percolation procedure just defined above is not always free of aliasing: it can be successful only if the separation between elements of metric is sufficiently large compared to the grid element size and if their zone of influence is not too distorted. Testing the validity of the percolation process is relatively simple, by checking if some points inside the computing domain have been left unexplored. If it is the case, percolation is performed again locally until all the unexplored points are processed, which induces a small additional cost in the algorithm.

To speed up the process, which is not really relevant here in 2D phase-space but crucial in higher number of dimensions, we reduce the range of investigations of values of $i_{\rm m}^{\prime}$ to $[i_{\rm inf},i_{\rm sup}]$ where

[TABLE]

with the dependence on $(i,j,i_{\rm m},j_{\rm m})$ of $\delta q_{x}$ still implicit and $\alpha$ is a parameter close to $1/2$. Similar conditions are set to define an interval $[j_{\rm inf},j_{\rm sup}]$. The parameter $\alpha$ is introduced to take into account the uncertainty ${\boldsymbol{E}}_{\rm map}=(E_{x},E_{v})$ on the remapping (10). This is illustrated by Fig. 3.

However, even though we can actually estimate ${\boldsymbol{E}}_{\rm map}$ as described below in § 2.3.3, we do not use it to calculate $\alpha$ which is fixed from the beginning. Indeed, our implementation favours performance over perfect accuracy, to make a future extension to 6-dimensional phase-space possible: while it is better to find in Lagrangian space the actual nearest element of metric to a test particle, this is not a must and the algorithm can still work without it, at the cost of some loss of accuracy on the long run.

To illustrate the result of the percolation process, Figure 4 displays, at $t=25$, the region of influence of each element of metric at the end of a cycle (i.e. just before resampling with new elements of metric) for the simulation II with Gaussian initial conditions described in § 3.

2.3.2 Interpolation of the Lagrangian position of a test particle

We now explain how the Lagrangian position of a test particle can be calculated in a smooth way with proper interpolation in Lagrangian space between predicates given by each of neighbouring element of metric. We define the following mapping function ${\boldsymbol{Q}}_{\rm r}({\boldsymbol{P}})$ from Eulerian to Lagrangian space

[TABLE]

where function ${\boldsymbol{Q}}({\boldsymbol{P}},i_{\rm m},j_{\rm m})$ is the obvious extension to the continuum of its discrete counterpart ${\boldsymbol{Q}}(i,j,i_{\rm m},j_{\rm m})$ defined previously, likewise for $\delta{\boldsymbol{Q}}({\boldsymbol{P}},i_{\rm m},j_{\rm m})$, and $W(x)$ is a smooth isotropic kernel function compactly supported. Clearly, function ${\boldsymbol{Q}}_{\rm r}({\boldsymbol{P}})$ as defined above is at least continuous. Our choice for $W(x)$ is the spline of Monaghan (1992):

[TABLE]

which preserves the second order nature of the Taylor expansion around each metric element position. Our choice of $W(x)$ allows us to restrict the range of investigation of values of $i_{\rm m}^{\prime}$ and $j_{\rm m}^{\prime}$ to the respective intervals $[i_{\rm m}-1,i_{\rm m}+1]$ and $[j_{\rm m}-1,j_{\rm m}+1]$ where $(i_{\rm m},j_{\rm m})$ corresponds to the closest element of metric to point ${\boldsymbol{P}}$ as computed in previous paragraph, § 2.3.1, for sites $(i,j)$ of a grid.222Obviously, we assume here that the errors in the remapping (10) are smaller than the inter-element of metric separation. Here, we indeed only need to consider points ${\boldsymbol{P}}=(\Delta x\,i+x_{\rm mid},\Delta v\,j+v_{\rm mid})$ associated to such an Eulerian grid. To reduce furthermore the number of elements of metric investigated to perform the interpolation (39) while still taking into account small variations related to the previously discussed uncertainty ${\boldsymbol{E}}_{\rm map}$ on the calculation of ${\boldsymbol{Q}}$, we restrict the search to the intervals $[i_{\rm min},i_{\rm max}]$ and $[j_{\rm min},j_{\rm max}]$ defined as follow, with the same notations as in equations (36), (37) and (38):

[TABLE]

where $\beta$ is a small parameter accounting for the error ${\boldsymbol{E}}_{\rm map}$, similarly for $j_{\rm min}$ and $j_{\rm max}$, with the additional constraints $i^{\prime}_{\rm m}j^{\prime}_{\rm m}=0$ if $i_{\rm min}i_{\rm max}=-1$ or $j_{\rm min}j_{\rm max}=-1$ to gain even more time.

Again, similarly as stated at the end of § 2.3.1, the above algorithm is not guaranteed to provide a smooth interpolation if the error on the remapping associated to the metric elements becomes larger than some threshold. However, this should have little consequence on the algorithm accuracy and the gain of speed in high dimensional phase-space is considerable with such a procedure. Indeed, if ${\boldsymbol{E}}_{\rm map}$ becomes too large, the loss of accuracy on the characteristics themselves is much more dramatic than small discontinuities introduced on their interpolation.

2.3.3 Estimates of errors on the Lagrangian position and choice of the margin parameters $\alpha$ and $\beta$

In the two subsections above, we defined two parameters $\alpha$ and $\beta$ but did not provide any numerical value. In principle these parameters could be user-defined, but we chose to fix them:

[TABLE]

The idea behind this choice is the following. Firstly, according to Figure 5, since we do not want to consider any diagonal element of metric if $i_{\rm min}i_{\rm max}=-1$ or $j_{\rm min}j_{\rm max}=1$, we must have for full consistency,

[TABLE]

Because there is some error on the reconstructed value of ${\boldsymbol{Q}}$, we take half this value, assuming that

[TABLE]

Hence $\alpha=1/2-E_{\rm max}$ is given by equation (44) to take into account this upper bound for the errors. With this setting, the algorithm is fully consistent only if the uncertainty on the remapped trajectories remains small enough, but even if $E_{x}/\Delta x_{\rm m}$ or $E_{v}/\Delta v_{\rm m}$ are larger than $E_{\rm max}$, it still works. Simply function ${\boldsymbol{Q}}_{\rm r}({\boldsymbol{P}})$ is not guaranteed to be completely smooth anymore.

Note that in our algorithm, a global error $E_{\rm map}=\sqrt{(E_{x}/\Delta x_{\rm m})^{2}+(E_{v}/\Delta v_{\rm m})^{2}}$ is estimated during the percolation process by comparing the value of ${\boldsymbol{Q}}$ proposed by a neighbouring element of metric to the value of ${\boldsymbol{Q}}$ proposed by the element of metric under consideration using the Lagrangian distance $d$ defined in equation (35). Such an error can in principle be used to force global remapping if $E_{\rm map}$ exceeds some threshold, for instance $E_{\rm max}$ to preserve full smoothness of function ${\boldsymbol{Q}}_{\rm r}(\boldsymbol{P})$. While this is an option in our code, we do not use it for the tests performed in § 3. Indeed, we found difficult to find some reliable criterion on $E_{\rm map}$ and preferred to keep the number $n_{\rm s}$ of time steps between full resampling fixed and see how the results depend on the choice of $n_{\rm s}$, being aware of the fact that this is potentially suboptimal. Note also that $E_{x}$ and $E_{v}$ could be estimated as the residuals from the third order description of the local metric elements, but this is very costly in highly dimensional phase-space.

2.4 Equations of motion

In $2D$-dimensional phase-space, the equations characterising the evolution of an element of metric with time are given by the following expressions with some obvious new notations:

[TABLE]

where ${\boldsymbol{a}}$ is the acceleration calculated by solving Poisson equation. Note thus that we need to have access at all times to the gradient and the Hessian of the gravitational force. In these equations we have dropped the “m” subscript for simplicity.

The actual time step implementation uses a splitting algorithm with a drift phase during half a time step where only spatial positions are modified, followed with a kick phase where only velocities are modified during a full time step using the acceleration computed after the drift, followed again by a drift during half a time step, where spatial positions are modified using the new velocities.

The first drift step or predictor is thus written

[TABLE]

after which Poisson equation is solved to compute the acceleration ${\boldsymbol{a}}(t+\Delta t/2)$, its gradient and its Hessian. To compute the acceleration, we first need to estimate the phase-space distribution function $f(\boldsymbol{x},\boldsymbol{v},t+\Delta t/2)$ on the Eulerian grid. To map $f$ back to initial conditions (or previous stage after full remapping), we use the elements of metric as transformed by the predictor step.

The kick phase reads

[TABLE]

Finally the second drift phase, or corrector, is expressed as follows:

[TABLE]

With our predictor-corrector scheme, it is possible to adopt a slowly variable time step. In this case, we use the following dynamical constraint (see, e.g., Alard and Colombi, 2005; Colombi & Touma, 2014),

[TABLE]

where $\rho_{\rm max}={\rm max}_{x}\rho(x)$ is the maximum value of the projected density, $\rho(x)\equiv\int f(x,v){\rm d}v$. Note that for the tests performed in § 3, we used a constant time step, but this should not have critical effects on the numerical results.

2.5 Calculation of the force and its derivatives

In practice, the force is calculated on a grid and its gradient and Hessian estimated by simple finite difference methods.

In our 2D phase-space implementation, the acceleration is calculated at discrete points $i$ with coordinate $x_{i}=x_{\rm mid}+\Delta x\,i$ corresponding to the projection of the Eulerian phase-space grid. To estimate it, we first calculate the projected mass $M(i)$ inside each pixel $i$ by summing up the phase-space density calculated on each point $(i,j)$ of the Eulerian phase-space grid. In the present work we consider isolated systems in an empty space, which means that the acceleration at a given point $x$ is simply proportional to the total mass at the right of point $x$ minus the total mass at the left of $x$. This implies, within some trivial factor depending on code units

[TABLE]

with

[TABLE]

and where $M_{\rm total}$ is the total mass of the system. Note thus that the mass in pixel $i$ does not contribute to $a(i)$ as it is supposed to be distributed symmetrically inside it.333Note that a potentially better procedure could consist in using a staggered mesh for computing the acceleration, in which the nodes would correspond to boundaries of each pixel of the grid. In larger number of dimensions, Poisson equation is more complex to solve but this can be done on a fixed grid resolution with FFT methods, whether the system is spatially periodic or isolated as supposed here.

Our estimate of the first and second derivatives of the acceleration at grid points is written

[TABLE]

which is equivalent to fit the acceleration with a second order polynomial on the three successive points corresponding to $i-1$, $i$ and $i+1$.

To compute the force at an arbitrary point of the computing domain, we perform dual TSC (triangular shaped cloud, see e.g., Hockney and Eastwood, 1988) interpolation from the grid, i.e.,

[TABLE]

and likewise for its gradient and its second derivative. Note that we also tested linear interpolation, which leads in practice to nearly the same results as TSC for all the numerical tests we discuss in next section.

Clearly, the way we estimate the force and its successive derivatives is probably too naive but actually works pretty well for all the experiments we did so we did not bother to try to improve this part of our algorithm. However, this should be investigated in a future work.

3 Numerical experiments

To test the method, we performed an ensemble of simulations of which the parameters are listed in Table 1. In § 3.2, we compare the results obtained with Vlamet to the standard splitting algorithm of Cheng and Knorr (1976) with cubic B-spline interpolation as well as the extremely accurate entropy conserving waterbag simulations of Colombi & Touma (2014, hereafter CT14). Before going into the details of these comparisons, we first provide in § 3.1 additional technical information on the performances of Vlamet and of the splitting algorithm.444Comparisons to the waterbag code in terms of performances are irrelevant here, since the future goal is to generalize Vlamet to higher number of dimensions in phase-space. The waterbag method, extremely accurate but already costly in 2-dimensional phase-space, becomes prohibitive in higher number of dimensions. For instance, in 4-dimensional phase-space, it consists in following, in an adaptive fashion, ensembles of hypersurfaces sampled with tetrahedra (see, e.g. Sousbie and Colombi, 2016). In particular, we discuss parallelisation issues as well as extension to higher number of dimensions.

3.1 Performances: parallel scaling and computing cost

Both Vlamet and the splitting algorithm are parallelised on shared memory architecture using OpenMP library, but no attempt has been made to optimize cash memory, so there is still room for improvement. Parallel scaling is displayed on Fig. 6, which shows that the two codes present a reasonable speedup, greater than 10 on 32 cores. The scaling is slightly better for the metric code than for the splitting algorithm but this comparison is irrelevant without further optimisation. The main point here is to show that OpenMP parallelisation is rather straightforward for Vlamet, while it is already well known that both distributed and shared memory parallelisations are not a problem for the splitting algorithm (see, e.g., Crouseilles et al., 2009).

As an additional note, but this is something we leave for future work, there is in principle no particular difficulty in parallelising Vlamet on a distributed memory architecture with the MPI library. Except for solving Poisson equation and for cubic B-spline interpolation that we discuss just after, all the operations performed in the algorithm, including the percolation process, are of local nature, hence trivially compatible with domain decomposition. In 4 and 6-dimensional phase-space, introducing a k-d tree structure to focus on important regions of phase space is also quite standard in parallel programming. The solution to Poisson equation is well known in many parallel implementations. For instance there exist parallel Fast Fourier Transform libraries such as FFTW to compute the gravitational potential on a fixed resolution grid. As for spline interpolation, a domain decomposition implementation is possible using the method proposed by Crouseilles et al. (2009). In this case, interpolation is localised to each domain by using Hermite boundary conditions between the domains with an ad hoc reconstruction of the derivatives. In this implementation, each domain would for instance correspond to a leaf of the k-d tree mentioned above. Of course, as long as we did not actually implement the algorithm on distributed memory architecture, this discussion stays at the speculative level.

In terms of total computational time, interesting conclusions can be driven from examination of fifth column of Table 1, which indicates, for each run, the approximate total CPU cost in units of simulation I. Neglecting the Poisson equation resolution step, both Vlamet and the splitting algorithm are expected in general to scale roughly like the number $n_{\rm cells}\propto\Delta^{-2D}$ of cells in the grid representing phase-space, where $D$ (here equal to unity) is the dimension of space. This is approximately the case in Table 1 if $\Delta$ is small enough. The cost also increases for Vlamet when resamplings of the phase-space distribution function are more frequent (in practice, if n_{\rm s}\mathrel{\mathchoice{\vbox{\offinterlineskip\halign{\hfil\displaystyle#\hfil\cr<\cr\sim\cr}}}{\vbox{\offinterlineskip\halign{\hfil\textstyle#\hfil\cr<\cr\sim\cr}}}{\vbox{\offinterlineskip\halign{\hfil\scriptstyle#\hfil\cr<\cr\sim\cr}}}{\vbox{\offinterlineskip\halign{\hfil\scriptscriptstyle#\hfil\cr<\cr\sim\cr}}}}10 for the runs we did) or when the element of metric density increases (in practice, if \Delta_{\rm m}\mathrel{\mathchoice{\vbox{\offinterlineskip\halign{\hfil\displaystyle#\hfil\cr<\cr\sim\cr}}}{\vbox{\offinterlineskip\halign{\hfil\textstyle#\hfil\cr<\cr\sim\cr}}}{\vbox{\offinterlineskip\halign{\hfil\scriptstyle#\hfil\cr<\cr\sim\cr}}}{\vbox{\offinterlineskip\halign{\hfil\scriptscriptstyle#\hfil\cr<\cr\sim\cr}}}}10\Delta for the runs we did).

When it comes to really compare Vlamet to the standard splitting algorithm, we focus on single thread runs, assuming that full parallelisation is equally optimal on both algorithms. We noticed that with the same spatial resolution $\Delta$, Vlamet is about $\alpha_{\rm slower}=1.4-2$ times slower than the splitting algorithm for n_{\rm s}\mathrel{\mathchoice{\vbox{\offinterlineskip\halign{\hfil\displaystyle#\hfil\cr>\cr\sim\cr}}}{\vbox{\offinterlineskip\halign{\hfil\textstyle#\hfil\cr>\cr\sim\cr}}}{\vbox{\offinterlineskip\halign{\hfil\scriptstyle#\hfil\cr>\cr\sim\cr}}}{\vbox{\offinterlineskip\halign{\hfil\scriptscriptstyle#\hfil\cr>\cr\sim\cr}}}}25 and \Delta_{\rm m}/\Delta\mathrel{\mathchoice{\vbox{\offinterlineskip\halign{\hfil\displaystyle#\hfil\cr>\cr\sim\cr}}}{\vbox{\offinterlineskip\halign{\hfil\textstyle#\hfil\cr>\cr\sim\cr}}}{\vbox{\offinterlineskip\halign{\hfil\scriptstyle#\hfil\cr>\cr\sim\cr}}}{\vbox{\offinterlineskip\halign{\hfil\scriptscriptstyle#\hfil\cr>\cr\sim\cr}}}}2.5. How these estimates generalize to 4 and 6 dimensions remains an open question, but we can try to infer some rough estimates. From now on we restrict to 6D phase-space. In this case, the critical step in Vlamet is the cubic spline interpolation during the resampling phase, which is expected to be about 171 times more costly than a step of the standart splitting algorithm, as discussed in introduction of § 2.2. The rest of the algorithm, on the other hand, should perform comparably well to the splitting scheme, i.e. should be $\alpha_{\rm loss}$ times slower with $\alpha_{\rm loss}$ of the order of unity. Hence we expect, in six dimensions, a factor $\alpha_{\rm slower}^{6D}\simeq(n_{\rm s}\times\alpha_{\rm loss}+171)/n_{\rm s}$ between Vlamet and the standard splitting algorithm, that is $\alpha_{\rm slower}^{6D}$ ranging from about 2 in the optimistic cases (e.g. $\alpha_{\rm loss}\simeq 1$ and $n_{\rm s}=200$) to about 20 in the pessimistic ones (e.g. $\alpha_{\rm loss}=10$ and $n_{\rm s}=15$). While a factor as large as 20 might seem potentially prohibitive, remind that our code is supposed to be less diffusive than the splitting algorithm, as we shall illustrate next, hence should require less spatial resolution. A simple factor 2 in resolution corresponding to a gain of about $2^{6}=64$ in speed and memory would make Vlamet extremely competitive.

3.2 Test runs: comparisons between codes

Before going further, it is important here to specify the units in which we are solving Vlasov-Poisson equations. Our choice is the same as in CT14, i.e. we are studying the following system of equations:

[TABLE]

To test Vlamet, three kinds of initial conditions were considered, namely the thermal equilibrium solution, Gaussian initial conditions and a random set of stationary clouds, as detailed below.

• •

The thermal equilibrium solution (Spitzer, 1942; Camm, 1950; Rybicki, 1971) reads

[TABLE]

In practice, we set

[TABLE]

with the apodizing function $g_{\eta}$ given by

[TABLE]

For the runs performed in this paper, we took $\rho_{\rm S}=4$, $\sigma_{\rm S}=0.2$ and $\eta=0.02$, which means that the maximum projected density of the system is $\rho_{\rm max}\simeq 2$. The harmonic dynamical time corresponding to such a density is, in our units,

[TABLE]

The simulation was run up to $t=100$, which corresponds to $100/T_{\rm dyn}\simeq 32$ dynamical times. Thermal equilibrium is stable and hence represents a crucial test for a code, which has to be able to maintain it accurately. Figure 7 shows how Vlamet performs for runs I, II, III and IV, when considering the difference $\delta f=f(x,v,t)-f_{\rm ini}(x,v)$ at $t=100$, i.e. checks the behaviour of the code when changing the distance $\Delta_{\rm m}$ between metric elements while keeping resolution fixed to $\Delta=0.002$. For comparison, a curve corresponding to the splitting algorithm run XVI is also shown at same resolution $\Delta$. As expected, Vlamet performs very well, as long as $\Delta_{\rm m}/\Delta$ remains sufficiently small, \Delta_{\rm m}/\Delta\mathrel{\mathchoice{\vbox{\offinterlineskip\halign{\hfil\displaystyle#\hfil\cr<\cr\sim\cr}}}{\vbox{\offinterlineskip\halign{\hfil\textstyle#\hfil\cr<\cr\sim\cr}}}{\vbox{\offinterlineskip\halign{\hfil\scriptstyle#\hfil\cr<\cr\sim\cr}}}{\vbox{\offinterlineskip\halign{\hfil\scriptscriptstyle#\hfil\cr<\cr\sim\cr}}}}10, to avoid too large errors on the reconstruction of Lagrangian positions from the metric elements. We also checked that $|\delta f|$ decreases when augmenting the spatial resolution of the simulation. Note that the splitting algorithm seems to perform even better than Vlamet, as suggested by the turquoise curve on right panel of Fig. 7. A residual error remains in all cases, since there is a clear convergence to a non zero $\delta f$. From performing the a run similar to I but with a twice larger time step and another one with a truncature parameter $\eta$ twenty times smaller, we see that this small residual can be attributed to the finite value of the time step and especially the (rather strong) truncature of $f$ with apodizing function $g_{\eta}$, which implies that the initial conditions constructed this way do not correspond exactly to thermal equilibrium. Although we do not have any rigorous mathematical argument to support the claim that follows, decreasing $\Delta t$ and the value of $\eta$ should make this residual tend to zero, both for Vlamet and the splitting algorithm, provided that resolution $1/\Delta$ becomes also arbitrarily large and the ratio $\Delta_{\rm m}/\Delta$ is kept sufficiently small.s

• •

Gaussian initial conditions are given by

[TABLE]

with

[TABLE]

In practice, to compare the output of Vlamet to the waterbag result of CT14, we take ${\cal R}=1$, $\rho_{\rm G}=4$, $\sigma_{\rm G}=0.2$ and $\eta_{\rm G}=0.02$, which makes the total mass of the system approximately equal to unity for a Gaussian truncated at 5 sigmas. The measured maximum projected density of this system, is, after relaxation, $\rho_{\rm max}=2.95$, so simulating this system up to $t=100$ corresponds to following it during approximately 39 dynamical times according to equation (76). Figure 8 shows, at various times, the phase-space distribution function obtained from run I and the waterbag run Gaussian84 of CT14. The agreement between both simulations seems visually close to perfect. A more careful examination of the details of the phase-space distribution function shows small deviations which can be reduced by augmenting the resolution of the Vlamet run. This is illustrated by Fig. 9, which displays on a small interval, the projected density for runs I, X, XI (different resolutions $\Delta$ but fixed $n_{\rm s}=25$ and fixed ratio $\Delta_{\rm m}/\Delta=10$) and Gaussian84. What is particularly interesting in this figure are the turquoise curves corresponding to runs performed with the splitting algorithm, namely XII, XIV and XVI. We notice a very good mach between Vlamet and the splitting algorithm but only if this latter is run at twice the resolution of Vlamet! This figure already demonstrates the significantly weaker level of diffusion of Vlamet compared to the splitting algorithm.

• •

The ensemble of stationary clouds corresponds to a set of halos following the stationary solution (73), created exactly as in CT14. In particular, their profile follows equation (73) with $\rho_{\rm S}=6$ and their velocity dispersion, $\sigma_{\rm S}$, ranges in the interval $[0.005,0.1]$. The components are added on the top of each other in phase-space, to obtain the desired initial distribution function $f_{\rm ini}(x,v)$, which is apodized as in equation (74) with $\eta=0.05$. The interesting fact about this system is that it was found experimentally by CT14 to be chaotic. In figure 10, the phase-space distribution function is shown at various times for simulation I and run Random of CT14. At $t=25$, the agreement between both simulations is impressive, despite the fact that the phase-space distribution function is sampled with only three waterbags in the Random run of CT14. At $t=80$, the system presents remarkable structures: some halos have been destroyed by tidal effects and a complex filamentary structure has grown. At this time, the Vlamet and the waterbag runs diverge from each other, but remain still in agreement at the macroscopic level. The fact that the simulations diverge from each other is not surprising since the run performed by CT14 samples the phase-space distribution function with only 3 waterbags. The differences are also amplified by the intrinsic chaotic nature of the system. Note that the Vlamet run is subject, at this point, to significant aliasing effects. For instance, the minimum value of the measured phase-space distribution function on lower-left panel of Fig. 10 is equal to $f_{\rm min}=-0.38$, which is of rather large magnitude given the maximum theoretical value of $f$, $f_{\rm max}=6$. These aliasing effects are inherent to the spline interpolation procedure we are using. Importantly, aliasing effects do not affect the dynamical properties of the system: increasing the resolution or decreasing it (runs X and XI) does not change the results significantly, despite the chaotic nature of the system, except at very small scales close to spatial resolution, of course. As illustrated below quantitatively, similarly as in the Gaussian case discussed above, simulations with the splitting algorithm provide comparable results to Vlamet but only if performed at better resolution.

To perform more quantitative tests, we checked for conservation of total energy, which can be conveniently written

[TABLE]

where $a(x,t)$ is the gravitational force and $M_{\rm total}$ the total mass. These integral are computed directly as sums on the grid used to sample the phase-space distribution function and the force. In addition, we check for conservation of the following Casimir

[TABLE]

which is numerically estimated the same way as total energy. This quantity would reduce to Gibbs entropy if positivity of the phase-space distribution function was preserved. While we know it is not rigorously correct, we shall still designate $S$ by entropy. We also tested for total mass conservation, as well as $L^{1}$ and $L^{2}$ norms, but do not show the results here, to avoid superfluous content. Indeed, energy and entropy were found in practice to be sufficiently representative of the system when it comes to test the accuracy of Vlamet.

Figures 11, 12 and 13 show $E$ and $S$ in the simulations we performed for the stationary, the Gaussian and the random set of halos, respectively. These figures allow us to test thoroughly the performances of Vlamet and to compare it to the splitting algorithm. They are supplemented with Figure 14, which displays the corresponding maximum mismatch between the Lagrangian positions predicates from neighbouring metric elements, as defined in § 2.3.3.

The most striking result when examining Figs. 12 and 13 is the superior behaviour of Vlamet compared to the standard splitting method with respect to entropy, in agreement with the preliminary conclusions driven from the analysis of Fig. 9. Indeed, to obtain comparable entropy conservation in Vlamet and the splitting algorithm, it is required to have a higher spatial resolution in the latter. Of course, the results depend on the choice of $n_{\rm s}$ and of the ratio $\Delta x_{\rm m}/\Delta$. For instance, Vlamet runs with $n_{\rm s}=25$ and $\Delta x_{\rm m}/\Delta=10$ perform as well as the splitting algorithm with twice the spatial resolution. But if $n_{\rm s}$ is increased to 100, it is required to increase the resolution by a factor as large as $2.8$ (middle right panel of Fig. 13) or $4$ (middle right panel of Fig. 12) in the splitting algorithm to perform as well as Vlamet. The price to pay might be energy conservation, which is in general not as good in Vlamet as in the splitting algorithm, that preserves energy at a level better than $\sim 10^{-5}$ regardless of initial conditions and spatial resolution. However, energy conservation violation in Vlamet can be kept under control with the appropriate combination of parameter $n_{\rm s}$ and ratio $\Delta_{\rm m}/\Delta$, as discussed further below. On the same ground, we notice, when examining Fig. 11, that the splitting algorithm seems to preserve better the thermal equilibrium solution than Vlamet, as already noticed on right panel of Fig. 7, unless $n_{\rm s}$ or/and $\Delta_{\rm m}/\Delta$ are chosen small enough. However, one has to be aware that the smoothness of the thermal equilibrium solution provides a very good ground for the splitting algorithm.

Further investigation of Figs. 11, 12, 13 and 14 allows us to study more in details how Vlamet behaves with the various control parameters:

• •

Increasing overall resolution, i.e. decreasing $\Delta=\Delta x=\Delta y$ and $\Delta_{\rm m}=\Delta x_{\rm m}=\Delta y_{\rm m}$ by some factor improves both energy conservation and entropy conservation, as can be seen from last line of panels of all the figures, as expected. Note that adopting $\Delta x=\Delta y$ and $\Delta x_{\rm m}=\Delta y_{\rm m}$ as we did in this work is not necessarily optimal, but is a natural choice in the code units if total mass of the system is of the order of unity.

• •

Changing the resampling frequency parameter $n_{\rm s}$ affects both entropy and energy, but to a much larger extent the former than the latter. In particular, increasing $n_{\rm s}$ up to 100 clearly improves entropy at an acceptable cost for energy conservation, of which the violation remains below the $\sim 5\times 10^{-4}$ level, as can be seen from the second line of panels of Figs. 11, 12 and 13. This is not surprising: performing less frequent resamplings of the phase-space distribution function decreases diffusion and aliasing, which is good for entropy conservation, but it worsens energy conservation, because (i) the error on the reconstruction of Lagrangian coordinates increases (middle line of panels in Fig. 14) and (ii) the phase-space distribution function is sampled more approximately during subcycles. While effect (i) can be reduced by augmenting the metric elements density as discussed below, effect (ii) always remains, which prevents arbitrarily large values of $n_{\rm s}$.

• •

Sparse sampling the elements of metric, i.e. increasing $\Delta_{\rm m}$, augments the error in the Lagrangian coordinates reconstruction (first line of panels in Fig. 14). This affects energy conservation but not so much the entropy (except for the stationary simulations, which nevertheless preserve entropy at a level better than $10^{-3}$). Entropy is indeed mainly sensitive to spatial resolution, i.e. the value of $\Delta$, and as stated just above, to the frequency of resamplings. For the experiments we did, we find that the ratio $\Delta_{\rm m}/\Delta$ should not exceed a factor 10, and the smaller it is, the better the results, obviously. For instance, using \Delta_{\rm m}/\Delta\mathrel{\mathchoice{\vbox{\offinterlineskip\halign{\hfil\displaystyle#\hfil\cr<\cr\sim\cr}}}{\vbox{\offinterlineskip\halign{\hfil\textstyle#\hfil\cr<\cr\sim\cr}}}{\vbox{\offinterlineskip\halign{\hfil\scriptstyle#\hfil\cr<\cr\sim\cr}}}{\vbox{\offinterlineskip\halign{\hfil\scriptscriptstyle#\hfil\cr<\cr\sim\cr}}}}5 leads to energy conserved better than $\sim 5\times 10^{-5}$ for $\Delta=0.002$ and $n_{\rm s}=25$. Being aware that the gain in memory is already a factor $2.5^{6}=244$ in 6D when sparse sampling the elements of metric by a factor $2.5$ in each direction compared to the computational grid, it is obviously wise to keep the elements of metric sampling sufficiently dense at the price of a small increase of the computational cost (a factor of the order of 1.4 in two-dimensional phase-space for $\Delta_{\rm m}/\Delta=2.5$, according to fifth column of Table 1). Furthermore, a higher density of metric elements allows for a larger value of $n_{\rm s}$.

4 Conclusion

We have proposed and tested a new Vlasov solver, Vlamet, which uses small elements of metric to carry information about the flow and its deformation. While slightly more costly in 2-dimensional phase-space than standard semi-Lagrangian solvers such as the splitting algorithm at the same resolution, this code is intrinsically less diffusive due to seldom actual resampling of the phase-space distribution function, as explicitly demonstrated by our numerical experiments. For the tests we performed, we indeed found that to have comparable level of details in the splitting algorithm and in Vlamet, it was needed to increase the spatial resolution by a factor 2 to 4 in the splitting algorithm runs, and this factor might turn to become even larger after additional optimisation: some work is indeed still needed to understand how to use Vlamet in the best way, while keeping errors, in particular deviations from energy conservation, under control. Indeed, despite the numerous runs we did, we did not span all the possible range of control parameters of Vlamet. In particular, we did not explore thoroughly the case when the separation between elements of metric approaches the spatial resolution, e.g 2.5\mathrel{\mathchoice{\vbox{\offinterlineskip\halign{\hfil\displaystyle#\hfil\cr<\cr\sim\cr}}}{\vbox{\offinterlineskip\halign{\hfil\textstyle#\hfil\cr<\cr\sim\cr}}}{\vbox{\offinterlineskip\halign{\hfil\scriptstyle#\hfil\cr<\cr\sim\cr}}}{\vbox{\offinterlineskip\halign{\hfil\scriptscriptstyle#\hfil\cr<\cr\sim\cr}}}}\Delta_{\rm m}/\Delta\mathrel{\mathchoice{\vbox{\offinterlineskip\halign{\hfil\displaystyle#\hfil\cr<\cr\sim\cr}}}{\vbox{\offinterlineskip\halign{\hfil\textstyle#\hfil\cr<\cr\sim\cr}}}{\vbox{\offinterlineskip\halign{\hfil\scriptstyle#\hfil\cr<\cr\sim\cr}}}{\vbox{\offinterlineskip\halign{\hfil\scriptscriptstyle#\hfil\cr<\cr\sim\cr}}}}5, which should allow for even more seldom resamplings than what we achieved.

The gain in effective spatial resolution of Vlamet compared to the splitting algorithm compensates largely for the cost of the code and makes our method potentially extremely competitive in 4 and 6 dimensional phase-space. Our current implementation treats isolated self-gravitating systems in 2 dimensional phase-space and uses parallel acceleration on shared memory architecture with the OpenMP library. The next step, quite relevant to galactic dynamics consists in generalising the code to 4 and 6 dimensional phase-space. This extension will require the introduction of adaptive mesh refinement with e.g. a k-d tree structure to focus on regions of phase-space of interest as well as parallel programming on distributed memory architecture. This should not present, in principle, any major difficulty, as discussed in § 3.1.

Note also that generalisation to plasmas is obvious and could represent a way to improve the accuracy of current codes.

Still, our implementation remains quite naive and some more work is needed for full optimisation. In particular, we use simple finite difference method to compute gradient and Hessian of the force, which calls for improvements. Also, at the time of resampling, we employ the traditional third order spline interpolation which is clearly not free of defects, in particular aliasing effects that induce strong deviations from positivity of the phase-space distribution function (while, nevertheless not modifying the good dynamical behaviour of the system). One possible way to improve over cubic B-splines could consist in employing a discontinuous Galerkin representation (see, e.g. Cockburn et al., 2000; Qiu and Shu, 2011; Besse et al., 2017), where the phase-space distribution function is described in each computational cell by some polynomial of some order. Continuity is not necessarily enforced between neighbouring cells, which provides more flexibility. However, while the splitting algorithm allows one to decompose the dynamical set up in individual one dimensional problems, this is not possible in Vlamet. We have to find a way, at the moments of remapping, to reconstruct a new discontinuous Galerkin representation by composing these local polynomials defined in a four- or six-dimensional space with the non linear displacement field provided by the metric elements, which remains a challenging task.

Acknowledgements

SC greatly acknowledges hospitality of Yukawa Institute for Theoretical Physics (YITP), where part of this work has been carried out. We thank N. Besse for fruitful discussions. We also acknowledge the support of YITP in organising the workshop “Vlasov-Poisson: towards numerical methods without particles” in Kyoto, funded by grant YITP-T-15-02, ANR grant ANR-13-MONU-0003 and by Institut Lagrange de Paris (ANR-10-LABX-63 and ANR-11-IDEX-0004-02). This work was supported in part by ANR grant ANR-13-MONU-0003.