WVTICs -- SPH initial conditions for everyone
Alexander Arth Julius Donnert Ulrich Steinwandel Ludwig Böss Timo Halbesma Martin Pütz David Hubber Klaus Dolag
astro-ph.IM cs.CE

	

	
	usm,excluster,mpe]Arth, A.
	arth@usm.uni-muenchen.de
	led]Donnert, J.
	[url]https://github.com/jdonnert/WVTICs
	usm,mpa]Steinwandel, U. P.
	usm]Böss, L.
	mpa]Halbesma, T.
	konstanz]Pütz, M.
	usm,excluster]Hubber, D.
	usm,mpa]Dolag, K.
	[usm]University Observatory Munich, Ludwig-Maximilians-University Munich, Scheinerstr.1, D-81679 Munich, Germany
	[mpe]Max-Planck-Institute for Extraterrestrial Physics, Giessenbachstrasse 1, 85748 Garching, Germany
	[led]Leiden Observatory, PO Box 9513, NL-2300 RA Leiden, The Netherlands
	[mpa]Max Planck Institute for Astrophysics, Karl-Schwarzschild-Str. 1, D-85741 Garching, Germany
	[konstanz]Universität Konstanz, Fachbereich für Physik, 78457 Konstanz, Germany
	[excluster]Excellence Cluster Universe, Boltzmannstr. 2, D-85748 Garching, Germany
	
	
		We present a novel and fast application to generate glass-like initial conditions for Lagrangian hydrodynamic schemes (e.g. Smoothed Particle Hydrodynamics (SPH)) following arbitrary density models based on weighted Voronoi tessellations and combine it with improved initial configurations and an additional particle reshuffling scheme. We show our application's ability to sample different kinds of density features and to converge properly towards the given model density as well as a glass-like particle configuration. We analyse convergence with iterations as well as with varying particle number. Additionally, we demonstrate the versatility of the implemented algorithms by providing an extensive test suite for standard (magneto-) hydrodynamic test cases as well as a few common astrophysical applications. We indicate the potential to bridge further between observational astronomy and simulations as well as applicability to other fields of science by advanced features such as describing a density model using gridded data for exampling from an image file instead of an analytic model.
	
	
	
		methods: numerical; smoothed particle hydrodynamics; initial conditions; open source
	




§ INTRODUCTION

	

	With rising computing power, simulations have become an integral part of modern astrophysics over the last few decades. In the past years people have traversed more and more from theoretical calculations with pen and paper to complex computer driven computations. These range from high precision simulations of idealised systems on various scales, as for example <cit.>, to large cosmological boxes which model a significant portion of the visible universe <cit.>. Several numerical techniques have been developed, refined and compared in order to improve the tools to widen our understanding of the physical processes which take place in the universe.

	What they all have in common is the prerequisite of some sort of initial conditions for their simulations. These may be physically motivated as in the examples given above or simply pose the means to test the behaviour of one's code with an analytically understood problem. Not only a proper definition of the physical key quantities in these initial conditions is important, but also their numerical representation. Often simulations analyse physical processes in the highly non-linear regime, where noise produced from badly sampled initial conditions or even the round off error introduced by finite accuracy floating point numbers becomes important (e.g. <cit.>).

	Depending on the numerical scheme, obtaining optimal initial conditions can be a nearly trivial task or rather complicated. The former being the case in finite difference and finite volume discretisation schemes and the latter, as in the focus of this paper, for  mass discretisation schemes such as Smoothed Particle Hydrodynamics (SPH) <cit.> or Meshless Finite Mass / Volume <cit.>. In this paper, we present an open source code which can generate relaxed initial conditions for SPH given any arbitrary physical description of initial conditions based on the work of <cit.>. We start by defining the actual task in section <ref>, followed by a description of the presented algorithms in section <ref>. In section <ref> we show the performance of our implementation measured by way of several test problems and finally present a few applications in section <ref> before closing with describing the actual usage of the application in section <ref>.



§ OVERVIEW IC GENERATION

	
	
	

 §.§ Requirements and degrees of freedom

		
	
		Initial conditions contain a complete physical description of all physical quantities as a function of position inside the simulation volume, such as the density ρ( r⃗), internal energy u ( r⃗) and the velocity v⃗( r⃗). This mathematical description can be translated into a volume discretisation by volume averaging the functions, or in case of finite difference schemes, evaluating the functions at the discretisation points. For a mass discretisation scheme this step is not straight forward, because mass sampling particles have to be placed in a way that the physical quantities calculated from the ensemble are reproduced with minimal error. To this end, obtaining the SPH smoothing length & density already requires solving the continuity equation <cit.>. The other quantities then fall into place and need not to be discussed further.


		In the SPH formalism the density is calculated as the kernel weighted sum over the nearest neighbours:
		
    ρ_i = ∑_j^#NGBs m_j W ( | r⃗_i - r⃗_j |, h_i)

		The kernel function depends on the distance between both particles i and j and the smoothing length h_i. As standard kernel in our implementation and all tests presented in this paper we chose the Wendland C6 kernel function with 295 neighbours in three dimensions, following <cit.>.
		Since this is the only SPH specific type of equation we need for this paper, we refer the interested reader to <cit.> for a recent review of the SPH method.

		To reach a certain model density one can choose to use either constant or variable mass particles. While the latter makes it possible to choose a universally ideal particle placement by adjusting the masses to fit the density, it is very disadvantageous. A significant contrast in particle masses leads to either a large variance in the amount of neighbour inside a kernel and to high computational cost in low density regions or even to numerical instabilities. <cit.> have shown, that one can not reach a low energy state with variable mass SPH particles with the classical SPH method.[There exist modern hybrid methods for which this is not necessarily the case, like the MFV method presented in <cit.>.] Therefore, we use constant mass SPH particles.
	
	

 §.§ Particle placement

		
		
		
		Even in the simplest case of a constant density several common approaches exist for particle placement. These include for example random positions, lattice configurations and so called glass files.

		The former is easy to implement, but the poison noise translates directly into a large relative error in the density estimate.

		In the left panel of Fig. <ref> we present a two dimensional, random particle distribution, which reveals an irregular structure full of holes and clumps. Furthermore, particles are more likely to be closer to each other than the mean particle distance, than further away, which introduces an additional anisotropy onto the density noise.

		Lattice configurations on the other hand come in different flavours, for example based on a cubic unit cell (cP, bcc or fcc) or in hexagonal configuration (hcp, <cit.>). The advantage of these configurations is their regularity which results in an easy equilibrium configuration, i.e. a global minimal energy configuration. However, effects known from grid codes are introduced, like the introduction of principle axes which may lead to orientation-dependent results.

		Finally, in a glass file particles are distributed in a similar fashion as molecules in a physical glass. It defines a quasi hydrostatic equilibrium condition for the system and therefore poses (at least a local) a low energy state. Particles are distributed non-regularly, however, with a very confined distribution of inter particle distances which helps to reduce noise <cit.>. This can be seen in the right panel of Fig. <ref>.
		
	

 §.§ Variable density modifications

		
		
		The initial random sampling can be easily adjusted to follow a different density distribution than the uniform one. We will later exploit this in section <ref>. However, the sampling error remains. Therefore, we rule out this method if it is used without additional effort to reduce the sampling error.[Random sampling still poses as a valid starting state for our approach.] When considering lattices and glass files, we remark that these can not be easily translated to a variable density. The homogeneous setup has to be stretched and compressed to follow the underlying density model. This process is not straight forward and it results in a loss of regularity for the lattice, while the glass transforms to a less relaxed ensemble. Nevertheless, stretching of uniform distributions is still commonly applied in the community <cit.>. We refer to <cit.> for a comparison of different approaches.

		<cit.> present the novel SPH/N-body code Gandalf that uses an alternative method for the generation of initial conditions. This method relaxes a set of particles towards a density model and a glass distribution simultaneously by applying two forces between particle pairs, one for each target property. The algorithm we present in this paper works similarly at it's core and we proceed by describing our approach in the next section.



§ CODE DESCRIPTION

	
	
	Our implementation is based on a Weighted Voronoi Tesselation (WVT) as presented in <cit.> and aims to produce relaxed, low energy glass files for any arbitrary given physical model. In the following subsections we will first describe the core parts of the algorithm, then clarify issues that typically arise and finally present a quick study on how to improve the convergence by choosing the best possible starting configuration.

    

 §.§ WVT relaxation

	    
	    
	    [style=pseudobase,basicstyle=,frame=single,float,floatplacement=H,label=lst:pseudocode,caption=Algorithm in pseudo code]res/algorithm.pseudo
	    We do not need to implement the actual Voronoi tessellation, although the original algorithm is based on a Voronoi Tesselation as shown by <cit.> (hence the name WVT). A short pseudo code describing the main routine is given in listing <ref>. The relevant parts of the relaxation algorithm are wrapped inside a generic SPH code which calculates densities and smoothing lengths for the current particle configuration. In each relaxation step a neighbour sum in classical SPH fashion is executed and for each particle pair a force is calculated to push the particles apart. The force is weighted with both the target smoothing length and the distance of both particles in order to converge to the given density model as well as a glass-like structure. In addition, a scaling factor for the movement of particles is set depending on the mean inter-particle distance. The goal is to find the optimum for the magnitude of particle pushes between wasting computing time through small steps and the possibility of overshooting. The net particle displacement of particle i due to particle j is then given by
	    
    Δr⃗_i = const· h_ij^m · W ( | r⃗_ij|, h_ij^m ) ·r⃗_ij/|r⃗_ij|

	    with the expected model smoothing length
	    
    h_i^m = h^m ( r⃗_i ) = ( N_NGB· m/4/3πρ^m ( r⃗_i ))^1/3,

	    which is normalised such that it equals to
	    
    h_i^m = ( 1/4/3π·ρ^m ( r⃗_i )/∑_j ρ^m ( r⃗_j ))^1/3.

	    Furthermore, it enters in symmetrised form depending on both particles i and j:
	    
    h_ij^m = ( h_i + h_j )/2.

	    There are two main changes made compared to the algorithm proposed by <cit.>: First, instead of an arbitrary weighting function which declines with distance we employ the most natural choice and use the kernel function from the surrounding SPH code, following <cit.>. Second, we switch out the multiplicative factor of h_i^m and use a symmetrised version of it. Consider an ideal distribution of particles that exactly follows the model density. The algorithm should then not move particles, meaning Δr⃗_i = 0  ∀ i. In a simplified picture the push from a neighbour on one side and one on the opposing end should cancel each other out. Although, that is not the case any more with the symmetrised version, the use case for our code is to start with an imperfect particle distribution and iterate towards the desired solution. Tests show that we never reach the exact model solution globally and that the symmetrised version gives much better convergence towards it, which results in lower density errors. Whether we employ the arithmetic or geometric mean does not make any notable difference, so we keep the arithmetic one.

	    Furthermore, the algorithm can in principle get stuck in a local energy minimum, because it is formulated locally, on neighbouring particle pairs. However, under the condition that the density structures are well resolved (meaning kernels are not large enough to smooth them out), convergence is generally not an issue. We present an analysis of the convergence in section <ref>.

        Optimal detection of convergence poses a challenging task as well. Typical criteria formulated on the average and the maximum density error have proven to waste computation time without the further improving the sampled density profile. Therefore, we employ a more pragmatic criterion and stop the relaxation as soon as the most particles are moved less than a small fraction of inter-particle distance. In addition, we put an upper limit onto the maximum number of iterations to further reduce computational cost.

    

 §.§ Initial state

	    
	    
	    The rate of convergence is drastically influenced by the initial state of the particle distribution. Thus, we compare a few different approaches:
	    
	    	
  * Uniform random distribution
	    	
  * Von Neumann rejection sampling
	    	
  * Statistical model along a space filling curve
	    
		The simplest approach for an initial setup is just a random arrangement of particles with positions drawn from a uniform distribution. As already discussed in section <ref> this itself does not pose a very good density estimate but can still serve as a valid starting point for relaxation. Due to the distance dependence in our algorithm the most crucial flaw of a random distribution, namely holes and clumps of particles, are dissolved quite rapidly.

		In cases where the density model varies clearly from a uniform density, one can improve the initial state by changing the random number distribution to the density function itself. This ensures that the particle density already follows the model much better without changing other properties. This is known as Von Neumann rejection sampling <cit.> and is basically a Monte Carlo integration method. While the rejection can sometimes be costly for extremely heterogeneous distributions, the computational effort is still negligible compared to the relaxation afterwards.

		Despite the Poisson noise due to particle placement, random distributions can still be useful when the noise is energetically unimportant compared to a physical process, e.g. in driven supersonic turbulence simulations.

		We experimented with a different approach, formulating a statistical particle placement model long a space filling curve (e.g. a Peano curve <cit.>) which gives very similar results while being more complicated to formulate properly. Therefore, we do not go into further details.
		
	

 §.§ Additional redistribution

		
		
		[style=pseudobase,basicstyle=,frame=single,float,floatplacement=H,mathescape=true,label=lst:redistributionpseudo,caption=Additional redistribution in pseudo code]res/redistribution.pseudo
		As we will show, the WVT algorithm converges asymptotically towards the given model. It generates a glass distribution of particles settled in a local minimum of the potential energy defined by the pairwise inter-particle force. Since the algorithm operates only locally using an SPH kernel weighting, a global force term can help to avoid local energetic minima and improve global convergence. To this end, we include the possibility for an additional global particle redistribution step in our code, which is carried out every few relaxation steps. 

		The idea behind this approach is that particles from regions with ρ > ρ_model are taken and directly moved into regions with ρ < ρ_model. We design an algorithm which is similar to the well known Metropolis algorithm <cit.> with the relative density error serving as the energy analogue. On top of that we put an upper limit to the amount of particles sampled, so that less particles are redistributed the closer we come to the model solution. We present corresponding pseudo code in listing <ref>. We randomly select particles to be moved and particles which we can use as tracers of an under dense region. Then we place the respective particles into the close proximity of these regions by drawing a random position with maximum distance of 0.3 of the kernel support radius. To accept a particle i for redistribution, we check if a random number r_i fulfils
		
    r_i ∈ [0,1]  < ( ρ_i - ρ^m ( x⃗_i )/ρ^m ( x⃗_i )).

		Similarly a particle j is viable as destination if another random number r_j fulfils
		
    r_j ∈ [0,1]  < ρ^m ( x⃗_j ) - ρ_j/ρ^m ( x⃗_j ).

		This ensures that particles with large density errors are in favour of being redistributed and that more particles are probed the closer the density comes to the model. We limit ourselves to a percentage of the particles in order to preserve the overall density profile and allow us to omit density calculations in between this process. Otherwise the computational cost would be too high for the method to be feasible. The maximum number of particles to be redistributed and to be probed is given by runtime parameters. In order to assist convergence we also let both percentages decay exponentially over time. For example we start at 1 % of particles and decay down to 0.1 % of particles until a certain iteration step. Between each redistribution we carry out several relaxation iterations to smooth out the distribution. At some point we stop the redistribution and only the relaxation step is carried out. The impact of this additional algorithm is investigated in section <ref>.
		


§ TEST PROBLEMS

	
	
	In this section we present algorithmic behaviour in typical use cases. While it is fairly straight forward to measure the accuracy of the resulting density distribution, the quality of the glass is more difficult to judge on only by using properties of the particle distribution. We start with a quick analysis of generated constant density distributions and then go one step further to our actual application of variable density models. All test cases are carried out in 3D with periodic boundaries if not mentioned otherwise. In order to save computational time we generally set up non-cubic boxes. We run all tests for a maximum of 1024 iterations but terminate in general earlier when reaching a quasi steady state. The additional redistribution of particles is switched off. We discuss the effect of the redistribution scheme later.

	

 §.§ Constant density

		
		
		
		We show how the density evolves for a constant density model during relaxation in Fig. <ref>. The top plot shows the density values of all SPH particles after different number of relaxation steps while the lower plot presents the according distribution of density values zoomed onto a ± 10% region of density values. Initially we sample the typical Poisson distribution, with noise of the order of ± 50% around the target value, visible in the top plot of Fig. <ref>, and quickly reduce to a reasonable agreement with the density model. One can clearly see the asymptotic behaviour of the convergence until the algorithm stops after 652 iterations.

		To assess the quality of the particle distribution, an analysis similar to Fig. <ref> can reveal the emerging patterns only in an obscured way due to the integrated third dimension. Furthermore, it lacks the power of a proper quantitative analysis with a parameter to measure the quality. A better judgement can be made by looking at the distribution of nearest neighbour distances (Fig. <ref>) and the radial distance autocorrelation function[Or just radial distribution function in the field of molecular dynamics.] (Fig. <ref>).

		
		We compare these two measurements for a random distribution and a relaxed glass distribution for a 3D constant density model with only 10^4 particles in a periodic box to prevent running into memory issues in our analysis script. We scale the distance by the expected distance calculated from a perfect distribution of the volume in the periodic box using
		
    d = ( V_p)^1.0/3.0 = ( V / N_p)^1.0/3.0

		with V the volume of the whole box and N_p the particle number. Note that while this is the optimal distribution of volume onto all particles it does not present a feasible configuration. Even in best close packing of spheres, there is still about 22.2% of the volume lost <cit.>. This is reached by a proper hcp lattice and depicts the asymptotic target for a well relaxed glass distribution. The distributions in Fig. <ref> can be fitted by a Gaussian. The deviation of the mean value from the expected mean and the width of the fit can be taken as a measurement of quality of the distribution. Since particles tend to clump in a random distribution the whole distribution shifts to smaller distances. The width is dominated by the Poisson noise seen in the density distribution of Fig. <ref>. While smaller distances are in principal desirable, this figure completely obscures that when particles clump some also spread too far apart from each other. Therefore, we have to not only consider the next neighbours but better analyse the whole distribution of particles using the autocorrelation function. 
		
		It gives us even more insight into how well our particle distribution resembles a physical glass. We bin the distances of all particle pairs and divide each resulting value by
		
    Norm = V_bin<ρ> · N_p/m,

		with V_bin the bin's spherical shell volume, <ρ> the mean density in the bin, N_p the particle number and m the particle mass as presented in <cit.>.

		We analyse the data by creating a histogram of the distribution (Fig. <ref>) which we do not further normalise or stretch in any way. The results from the different approaches can be compared to a Lennard-Jones fluid model as presented for example in <cit.>. Two criteria of quality are the height of the first peak and the behaviour on large scales. First of all a glass is characterised by a main peak height of 2.85 according to the Hansen-Verlet criterion <cit.> which characterises the transition from a liquid to a frozen state, demonstrating that our glass is indeed well relaxed. Also the peak's location is important. Since the first peak is mainly dominated by next neighbours we did this comparison already in Fig. <ref>. While the random sample damps down pretty quickly to a fairly constant distribution of values, we see dampened wave-like behaviour for our glass distribution around a central value of about one, as expected. The zoom onto large scales in the lower panel of Fig. <ref> especially shows, that this extends out properly until the maximum distances given by half of the box' size, showing that the glass distribution we produce is of high quality.

	
	

 §.§ Variable density

		
		
		
		Transitioning to variable density models, the next test consists of two plateaus connected with linear slopes. Here we examine the behaviour of our implementation with respect to linear gradients and kinks in the density function. In Fig. <ref> we present the evolution of the resulting density profile along with the analytic model. In addition to the overall profile, we also show zooms onto one upper and one lower kink of the distribution. The analytic model is given by
		
    ρ(x) = {[               0.5,       for  x < 0.1;  0.5 + (x-0.1)/0.3 for  0.1 < x < 0.4;                1.5 for  0.4 < x < 0.6;  1.5 - (x-0.6)/0.3 for  0.6 < x < 0.9;                0.5      for  x > 0.9. ].

		Again the scatter reduces quite rapidly and we converge towards the model density. As the scatter decreases, the density flattens marginally with respect to the actual maxima and minima of the given density function. Therefore, a low energy state and following the proper density model seem to be slightly exclusive. Since the extremal values are not fully resolved, the gradients are marginally shallower than imposed. Also the constant density parts are not totally flat but show a slight curvature which is expected due to the influence of their surroundings. SPH can not produce any sharp edges by construction, only up to the resolution level given by the kernel size. We over plot two circles in both the high and low density region in order to give an idea of the size of the respective kernels in these regions to show that the model is well resolved with 10^6 particles. In the end, after 1024 iterations steps, we reach a steady solution with a very low degree of scatter which deviates from the analytic solution by a few percent and shows slight curvature in the plateau and smoothed out kinks which, however, sit at the expected x-values.

		
		Next, we look at a smoother density model to investigate the error in maxima / minima further without the introduction of sharp edges or linear gradients. We set the density model according to
		
    ρ(x) = 1 + 1/2·sin (2 π x)

		and plot the resulting density evolution in Fig. <ref>. Due to the smooth nature of this test it already converges after 865 iterations. When reducing the scatter to a very low degree the density function again slightly deviates from the model at the maxima and minima, proving that this is not an effect of non-smooth density definitions. The lower panel shows the resulting relative density error compared to the imposed model for all plotted steps. Although slightly asymmetric the error being of the order of a few percent is definitely satisfactory.

		
		Finally we  investigate our algorithm's behaviour when confronted with sudden density jumps. We set the target density to
		
    ρ(x) = 1/2 + 2 ·( x 1/2)

		and plot the result in Fig. <ref>. While we can observe slightly more scatter here than in the previous tests, indicating that the final state (1009 iterations) is probably less converged due to the presence of the prominent density discontinuity, the jump itself is well resolved. Already after 32 iterations hardly any particles populate intermediate densities between the high and the low state. Also the resulting density error is in the same regime as before. We conclude that our implementation is able to reproduce all kinds of density models within the boundaries of SPH itself and the imposed resolution. Only density models which contain regions with ρ = 0 are challenging to model and require some fine tuning in our current neighbour finding routine.

		
		Considering the quality of the produced glasses,  we plot again the autocorrelation function for 10^4 particles distributed according to the sine wave model density in Fig. <ref>. While the mean density is still at the same level as in section <ref>, this introduces higher and lower density regions. High density regions contain more particles than the average. This leads to a shift in the distribution to smaller distances, especially noticeable at the main peak. In addition, we see a more strongly damped wave pattern compared to the previous tests. Furthermore, we notice that the distribution declines slightly when going to large scales. We have to truncate this plot at smaller scales since we divide the distances by the box size in x direction while the y and z direction are slightly smaller here to save computing time.

		We conclude that our 3D particle distributions actually form satisfactory glass configurations.
	
	

 §.§ Including Particle Redistribution

		
		
		Now we investigate the effects of our particle redistribution scheme (section <ref>). We start to probe ten percent and to redistribute at maximum one percent of all particles. We carry out this algorithm every 5 iterations until we reach the 512th iteration step and let the percentages decay exponentially to go down a factor of ten until the end. We carry out all three presented variable density tests with this setup and compare to the results presented in section <ref>.

		The sine wave density shows a bit larger scatter until the redistribution is switched off and the end results improves only slightly. However, the algorithm now requires about twice the iterations as before. While this can be mitigated by adjusting parameters, the improvement due to redistribution seems hardly worth it. The sawtooth setup shows minimal improvement, only the discontinuity becomes a bit sharper and less round.
		
		The plateau density setup shows the biggest effect, we plot the result in Fig. <ref>. The colours indicate the same iteration steps as before in Fig. <ref> for better comparison and we include the final result from before as a black line. As can be clearly seen, the distribution settles closer to the actual model density. In this case the amount of iterations the algorithm  takes to converge is actually about the same in both cases. However, due to the additional disturbance from redistribution, the density plateau comes out less straight than before. 
		
	    
	    To investigate the iteration process further, we sample the density plateau and linear gradient tests with 10^6 particles and particle redistribution. Each iteration we write out on the global density error, the movement of all particles in the last step and the degree of particle displacements. We plot the results in Fig. <ref>. In the beginning we observe large particle displacements due to large density errors. These damp down rather quickly. At about 100 iterations, the maximum density error has gone down to a few percent while the mean error is even about a factor 10 smaller.

	    As the errors decline, a periodic pattern emerges which is the imprint of our particle redistribution scheme. As it removes particles from local over-densities and puts them into under-dense regions, injecting potential energy in the particle distribution. The algorithm damps this perturbation very quickly. Nevertheless, redistribution counteracts convergence and we can not see the improvements in the final density correctness at this resolution. The amount of particles marked for redistribution quickly reaches the allowed maximum and the amount of particles actually touched declines drastically down to 10^-4 and then follows the imposed decay. We mark where this part of the code is shut off in the plots by a dotted vertical line. We do not improve on the density error any more but rather sacrifice a bit of density correctness for a better glass distribution and therefore less scatter in the density. The net particle kicks decrease as the distribution relaxes and movement damps down quickly. One can clearly see, that the net particle movement is well suited for an abortion criterion since it experiences the quickest change and directly indicates a steady state. 

		We conclude that the improvement due to particle redistribution is highly problem dependent and not always worth the extra computational effort.
	    
	

 §.§ Convergence with Particle Number

		 

	    
	    To investigate convergence on the algorithm with particle number, we plot the density for the plateau test case with 10^4, 10^5 and 10^6 particles next to each other in Fig. <ref>, zooming onto the top plateau. In all three runs we adjust the parameters such, that the initial degree of particle movement is of the same order to have a fair comparison. While the random sampling produces very similar results we see that the particle distribution develops differently in all three runs. The lower resolution runs converge much quicker than the higher resolution, as the mean particle distance is larger in this case. The plateau approaches the imposed model with increasing resolution, but shows more curvature as observed before. We can distinguish effects of the whole algorithm and of the inherent kernel smoothing. The latter is most apparent in the lowest resolution, where we see stronger curvature at the position of kinks in the density model. Finally, the scatter around the mean in the particle distribution is of the same order in all three runs.

	    
	    To investigate the behaviour of convergence with iterations and changing amount of particles quantitatively, we present the L1-error of all of our three standard tests in three different resolution steps in Fig. <ref>. The L1 error is defined in e.g. <cit.> as
	    
    L1 = ∑_i | ρ_i - ρ_m |/N_p

	    All tests have been carried out with a similar degree of movement in the initial setup, but without the redistribution algorithm. Since the L1 error is a global measurement, we only see the globally most prominent evolution in the beginning when the distribution quickly approaches the imposed density model. After about 100 steps, we observe hardly any change. Nevertheless, we observe a clear decrease of the L1 error with increasing particle number. According to the argumentation of <cit.>, the L1-error can not become better than 10^-2 - 10^-3 in classical SPH due to the E_0 gradient error of particle disorder <cit.>. This poses a hard lower bound on convergence with particle number.
	    


 §.§ Artificial bias correction

		
To account for the over- and undershooting of the resulting density we employ an optional bias correction. 
For that we calculate the density deviation following the equation:

    δ = ρ - ρ_m/ρ_m - ρ̅ + ρ_m - ρ̅/|ρ_m - ρ̅|ρ̅,

with ρ, ρ_m and ρ̅ being achieved density, model density and mean density respectively. 
The resulting bias is calculated at the end of the run and prompted to the user.
This value can be passed in the parameter-file for the next run where the model density is then increased or decreased in the areas where we over- and undershoot:

    ρ = ρ_m + (ρ_m - ρ̅) ·δ,

where δ denotes the bias factor from Eq. <ref>.
This approximates the target density significantly more accurate, as can be seen in Fig. <ref> where the plateau test was rerun with the bias correction switched on. Likewise, with this bias correction in place the L1 error (Fig. <ref>) drops almost one order of magnitude and thereby approaches the theoretical hard lower bound of 10^-2 - 10^-3 referenced in the previous section. 



    


§ APPLICATIONS

	
	
	In this section we demonstrate typical use cases in astrophysics from standard code tests to studies of astrophysical objects.
	
	

 §.§ Code testing: (Magneto-) Hydrodynamics

		
		
		
		
		Our implementation can be used to provide high quality initial conditions for standard (magneto-) hydrodynamical benchmark tests (<cit.>, <cit.>, <cit.>, <cit.>, <cit.>, <cit.> <cit.>, <cit.>). In particular, the WVT technique is useful to generate a well defined initial particle distribution with a low degree of noise. Since the method we present in this paper is capable of sampling arbitrary density profiles, we provide implementations of common benchmark tests listed in table <ref>, alongside references for the implemented models. 

		Here we demonstrate a more complex example the Zeldovich Pancake <cit.>. We use the same initial configuration, which has been presented in <cit.>:
		
    x(q,z) = q - 1+z_c/1+zsin(kq)/k,
 
		with the wave number k=2π /λ and the unperturbed coordinate q, as well as the collapse redshift z_c=1. The density profile is then given via
		
    ρ(x,z) = ρ_0/1- 1+z_c/1+zcos(kq),

		with the initial density ρ_0. We use a Box with 64 Mpc side length and sample the density profile with 10^6 gas particles. The Zeldovich-Pancake is an important test problem for both numerical hydrodynamics and numerical integration in a comoving frame where time t is replaced with the Hubble-function H(t). The resulting initial conditions are shown in Fig. <ref>. The quality of the initial conditions are very well in line with our previous results.
		
		
	

 §.§ Binary star formation

		
		
		In this section, we present the initial conditions for the Boss and Bodenheimer test case <cit.> to simulate the formation of a binary system of stars. We follow the implementation presented in <cit.> with the density profile
		
    ρ = ρ_0[1+0.1cos(2φ)],

		using a central density ρ_0 = 3.82 · 10^-18 g/cm^3. Furthermore, we set up a sphere of radius R=5.16 · 10^16 cm, a total mass of 1M_⊙, speed of sound c_s=1.66 · 10^4 cm/s and a solid body rotation of ω = 7.2 · 10^-13 1/s, as of <cit.>. In Fig. <ref> we show the rendered density profile in x-y plane with the perturbation in direction of φ.
		
		
	

 §.§ Isolated galaxy cluster

		
	
		The baryonic content of clusters of galaxies, in general, follows a beta-model profile <cit.>, given by 
		
    ρ_gas(r)    = ρ_0[ 1 + ( r / r_c )^2 ]^-3 β / 2   ,

		where the density within the core radius r_c is constant and decays as a power-law outside the core. <cit.> set up initial conditions for idealised binary merger simulations in the now open source code Toycluster[https://github.com(jdonnert/toycluster)], further improved by <cit.> by adding WVT relaxation. For simplicity, here we ignore the dark matter halo of the cluster as well as the temperature structure. The `cluster' is set up with general parameters ρ_0 = 10^-26 g/cm^-3, r_c=20 kpc, and β=2/3 which could resemble a cool-core cluster of galaxies. We sample the density profile using 10^6 SPH particles.

				
		We show the resulting density in Fig. <ref> at 32, 64, 128 and 1024 iterations, where we hit the maximum number of iterations and the algorithm terminated. This happens due to the large range of more than three orders of magnitudes in density which we try to sample with only 10^6 particles. In the low density region we can not properly resolve the gradient any more because the SPH kernels become very large, leading to a locally large spread in the density distribution early on which it relaxes away nicely. Furthermore, this results in quite big density errors in the low density regions of the resulting distribution. Judging these errors relatively to the mean density they are, however, quite small. Consequently, since the particle masses are calculated from the mean density, this over-density directly leads to an under-density in the central regime. Nevertheless our results are consistent with <cit.>, who successfully uses similar initial conditions for actual simulations of galaxy cluster mergers.
		
	

 §.§ Reading in an initial setup: Image processing

		
		
		
		Finally, we present a special feature of our application which allows us to read in a given set of gridded data for the definition of a model density, instead of providing a hard coded analytic solution. An example where we directly used the brightness scale of an image utilising our built in PNG reading routine as proxy for the density, is shown in Fig. <ref>. The goal of astrophysical simulations is always the comparison to observations in order to learn more about the processes in the universe. Mostly only statistical comparisons can be made because it is extremely difficult to simulate an observed object exactly in all it's details. This feature might actually help here, taking in an observational image in order to model the underlying matter distribution. The implementation is currently being expanded to take in proper 3D data which might be for example de-projected from observational images. Other applications outside of astrophysics can be found, where one might for example build a CAD model of a 3D distribution <cit.>.
		


§ CODE USAGE

	
	
	In this final section we talk briefly about the actual usage of our application. It can be downloaded from github <cit.> and will be made publicly available under the MIT license with the release of this paper. The current code version described here is v2.3. This section is split into two parts: First, we briefly state the compilation requirements. Second, we show the currently required parameters and describe how to actually run the application successfully.
	
	

 §.§ Requirements and compilation

		
		
		The code is written in OpenMP parallel C with the C99 standard and requires the gsl and gslcblas libraries. If the PNG reader shall be used, libpng is necessary in addition. There are no special requirements regarding the choice of compiler or any known issues with optimisation flags.

		
		For compilation we provide two identical approaches, a direct Makefile or, as required by some IDEs, a CMake file which then generates an own Makefile. Both approaches provide the possibility to adjust the compilation by setting some options which are translated into preprocessor defines. An overview is given in table <ref>.
				
	

 §.§ Parameters and running the application

		
		
		
		In order to run the application, one has to provide a parameter file by path. The main code directory contains already a sample parameter file called ics.par. It lists all provided runtime configuration options to adjust the main algorithms of the code and chose the respective initial conditions model. A list of parameters including standard values and short explanations are given in table <ref>.

		Besides these parameters one can change the behaviour of the OpenMP parallelisation by setting the usual environment variables to appropriate values, like
OMP_NUM_THREADS.

		The output is written in standard GADGET type 2 binary format.
    


§ SUMMARY AND CONCLUSIONS

	
	
	In this paper we presented a novel application to generate arbitrary initial conditions for SPH with a glass-like particle structure. The applied algorithms are based on a combination of weighted Voronoi Tesselation <cit.>, density aware random sampling <cit.> and an additional reshuffling of particles <cit.>. We showed an overview of common methods to solve this task and illustrated the challenges that come along. We ran several test problems with different density features, like linear and non-linear gradients, kinks and jumps, and investigated thoroughly how our implementation performs when confronted with them. We found that our algorithms converge reasonably well towards the given density models while also iterating towards a glass structure, which we analysed in 3D using a radial distance distribution function. We showed, that converging towards the correct density and a glass distribution at the same time is extremely difficult, since both tasks slightly contradict each other. Particle have to be pushed apart from each other in order to remove clumps and holes in the particle distribution, but this results in particles being moved out of high density regions resulting in lower density maxima and consequently also higher minima. Therefore, one has to settle for an application dependent trade-off. We also analysed the convergence of our algorithms with variable particle number and see desired behaviour.

	Our application also inhibits a suite of various sets of initial conditions, including many standard test problems for (magneto-) hydrodynamics and the ability to produce idealised setups for a galaxy cluster <cit.>. These tests are frequently required for code development of SPH codes and often the ability to judge the quality of a code was limited by pour initial conditions.

	We also presented a special feature which allows the code to sample a density model according to a PNG image instead of an analytic model. This might be useful to bridge observations and simulations further, since it can lead to the ability to simulate observed objects better. For that more development has to be done, though.

	Finally, we also presented the actual usage of our application, including a short description of all options and parameters which can be set at the moment (v2.1). The source code  will be made publicly available with the release of this paper. Further development of the project, may include support for more read in / write out formats, extending the possible choices of SPH kernels, expanding the 2D image reader eventually to 3D and even improvements on the core algorithms themselves. WVTICs is a flexible and versatile tool applicable for all kinds of problems. We welcome feedback and contributions via github.
    


§ ACKNOWLEDGEMENTS


	The authors would like to thank their colleague's Martin Zintl and Stefan Heigl for useful discussions.
	UPS is funded  by the  Deutsche  Forschungsgemeinschaft  (DFG,  German 
Research Foundation) with the project number MO 2979/1-1.
	This product includes software developed by Greg Roelofs and contributors for the book, PNG: The Definitive Guide, published by O'Reilly and Associated.





