banner

The site's hierarchy



Home page :: Private area :: Site Map :: Credits
Codes / powmes

Navigation

Mar 2009
MTWTFSS
2324252627281
2345678
9101112131415
16171819202122
23242526272829
303112345
News
  • CINES received a 150 TFlops machine
  • GENCI has acquired for the French University Supercomputing Center a 12288 core SGI Altix. It is ranked 14th in the World and first in France.
  • Site Web du CINES
  • IDRIS received a 140 Tflops BlueGene/P machine.
  • The CNRS supercomputing centre has acquired a new BlueGene machine with 40480 processors and 20 To of RAM memory.
  • See the press release.
  • Mare Nostrum simulation down to z=1.55
  • We carried out our fourth run of MareNostrum down to z=1.55 while producing 150 Million "stars". More info is given here
  • Horizon 4Pi 4096^3 down to z=0
  • We have carried a full sky 2Gpc/h 4096^3 Dark Matter simulation down to z=0. More details can be found here
  • Horizon 4Pi run 2048^3 completed
  • We have completed a full sky cone run of DM with RAMSES on platine CEA with 2048^3 particles. See a slice here
  • Mare Nostrum @ z=1.9
  • We carried out our third run of MareNostrum down to z=1.9 while producing 120 Million "stars"
  • Mare Nostrum simulation down to z=2.4
  • We carried our second run on the Mare Nostrum supercomputing center down to z=2.4, thanks to 4-fold improvement in overall performance.
  • 1024^3 500 h-1Mpc simulation completed
  • RAMSES run of the HORIZON initial condition was carried down to z=0 using 1024^3 dark matter particles in a (500 Mpc/h)^3 box with up to 16 levels of refinement (corresponding to local effective resolution of 65536).
  • Horizon 1024^3 100 Mpc/h simulation completed.
  • A RAMSES run of the HORIZON initial condition was carried down to z=0 using 1024^3 dark matter particles in a (100 Mpc/h)^3 box with up to 16 levels of refinement (corresponding to local effective resolution of 65536).
  • A 43 Tflops supercomputer at CEA
  • CCRT, the supercomputing centre of CEA, has just received a new BULL supercomputer with 6800 processors and 13.6 Tb of RAM memory. This facility is dedicated to academic and industrial research.
  • Lire le communiqué
  • Fourth Horizon Workshop 2006b
  • The 4th workshop takes place from December, the 11th to December, the 12th at Paris Observatory. Registration form available here.
  • MareNostrum simulation has started at BSC
  • After 1 week of operations, we have reached redshift 4. The simulation will proceed further more during the next quarters. We have started to post-process these initial data.
  • Horizon Grid
  • The Horizon grid is now in operation. It links the 6 quadri opteron of Meudon, Paris, Lyon, Saclay, IAP and Marseille and is open to all members of the collaboration. Voir ce lien
  • http://grille.projet-horizon.fr
  • Third Horizon Workshop 2006a
  • Horizon Workshop 2006a 10, 11 and 12 april 2006 in Lyon. It will be dedicated to progress reports. Registration form
  • Horizon is part of the DEISA "Extreme Computing Initiative"
  • Horizon will launch one extreme application for galaxy formation using the DEISA infrastructure at Mare Nostrum http://www.bsc.es
  • DEISA web site
  • Second Horizon Workshop
  • Horizon Workshop 2005 took place at Paris Observatory the 14th and 15th november. It was dedicated to internal discussions and management issues.
  • November 2005: Horizon méso-machine is fully operational at HPC1
  • October 2005: Horizon was awarded a 500 k€ grant by ANR, the French Science Foundation.
  • Horizon response to «ANR Blanche» Call for Proposals has been succesfull in the «UNIVERS» program.
  • ANR web site
  • July 2005: HP France was selected to host the Horizon "meso-computer"
  • HP company has succesfully answered the "call for proposals" issued by the Hrizon Project for a medium-size dedicated super-computer. The computer will be hosted at HPC1, the HP high performance computing center near Paris.
  • February 2005: INSU and CEA agree to fund the Horizon "meso-machine".
  • Horizon has submitted to INSU a proposal for funding a medium-size computer dedicated to the project. The proposal was accepted as a joint CEA and INSU operation.
  • September 2004: Kick-off meeting in Paris
  • During 3 days (13, 14 and 15 of september), 30 scientists have met at Observatoire de Paris to set up the basic organisation and objectives of the Horizon Project.
  • KO Meeting
  • June 2004: Project report at SF2A in Paris
  • April 2004: Horizon received official support from ASSNA.
  • The french initiative «Action Spécifique pour les Simulations Numériques en Astrophysique» have given to the Horizon project its "label" for outstanding and structuring computational project.
  • April 2004: The french Astroparticule program provides financial support for the Horizon "mini-grid"
  • The Programme Astro-Particule, a joint IN2P3 and INSU initiative, has agreed to finance 6 quad AMD64 servers dedicated to the Horizon Project.
  • January 2004: Review by PNC and PNG
  • Both PNG and PNC scientific comitees are officially supporting the Horizon Project initiative.
  • June 2003: Presentation of the project to SF2A
  • Theoretical Virtual Observatory workshop (IAP April 5-6, 2006)
  • A workshop dedicated to the Theoretical Virtual Observatory will take place at IAP on April 5-6th.

    The goal is to bring together experts of the Virtual Observatory and theoreticians who would like to make results of their simulations (e.g. databases or catalogs) or numerical codes available to the worldwild astronomical community.

  • Program and speakers

powmes

POWMES is a F90 program to measure very accurately the power spectrum in a N-body simulation.

by Pichon Christophe (Wednesday 19 November 2008)

POWMES is a F90 program to measure very accurately the power spectrum in a N-body simulation, using Taylor expansion of some order on the cosine and sine transforms. It can read GADGET format and requires FFTW2 to be installed.

The code itself is found here

GZ - 1.8 Mb
powmes tgz archive
POWMES tarball

An extended version with an example is found there

The corresponding paper is found there

!=======================================================================
!                              POWMES
!=======================================================================
! Authors : Stephane Colombi (Institut d'Astrophysique de Paris,
!                             colombiATiap.fr)
!           Dmitri Novikov   (Imperial College, London,
!                             d.novikovATimperial.ac.uk)
!
! Publication : Colombi, Jaffe, Novikov, Pichon, MNRAS, in press
!           (arXiv:0811.0313)
!
!-----------------------------------------------------------------------
!
! DISCLAIMER
! ++++++++++
!
! The POWMES software is provided `as is', and the authors (SC & DN)
! make no representations or warranties, expressed or implied, of
! merchantability or fitness for any particular purpose, and disclaim
! all liability for direct or consequential damages resulting from
! anyone's use of the program.
!
!-----------------------------------------------------------------------
!
! THE PROGRAM
! +++++++++++
!
! This programs measures the power spectrum in a 3D simulation using
! Fourier Taylor expansion of some order on the cosine and sine
! transforms.
!
! This program is still under development, and deserves further
! optimization to deal with large simulation samples.
!
! Taylor expansion :
! -------------------
!
! Approximation of zeroth order provides NGP interpolation (Nearest Grid
! Point) while higher order approximations take into account small
! displacements within cells of the mesh used to perform the Fourier
! transform increasingly accurately. Arbitrarily high order converges
! to the exact transform, namely to the sum
!
! delta_k=\Sum_i M_i exp(I.k.x_i/N), (1)
!
! where x_i is the position of the particle in [0,2.PI[, M_i its mass,
! (total mass of the system unity), k is an integer wavenumber, I^2=-1,
! and N is the size of the grid used to perform the calculations.
! The power spectrum P_k reads
!
! P_k=< |delta_k|^2 >_k,
!
! where <.>_k represents an angular average over the values of k. More
! exactly, function
!
! E_k=int(|k|+0.5)
!
! defines the size of the bin in k space (all k having the same E_k
! contribute to the same bin).
!
! Folding :
! ---------
!
! To speed up the calculation, it is possible to perform foldings of the
! particle distribution. One folding along one axis consists in replacing
! all the positions veryfing x_i >= PI with x_i-PI, and then multiplying
! all the x_i's (including the unmodified ones) by a factor 2 to span
! again the range [0,2.PI[. This operation is performed on all the
! coordinates. One can then demonstrate that the modes given by
! equation (1) remain unchanged for even values of k.
!
! Assuming now that we have a periodic simulation box and that we
! use a grid of fixed size N to perform the calculations, successive
! foldings will allow us to estimate the power spectrum for
!
! 0 folding  : k=0,1,2,3,4,...,N/2
! 1 folding  : k=2,4,6,8,...,N
! 2 foldings : k=4,8,12,16,...,2N
! 4 foldings : k=8,16,24,32,...,4N
!
! and so on.
!
! Errors :
! --------
! There are 2 possible sources of errors in the calculation.
!
! (a) the errors due to the finite number of modes sampled in a given
!     bin k. If N_k is the number of independent modes, then the relative
!     error due to the finiteness of N_k is approximately
!
!     Delta P_k/P_k=1/sqrt(N_k). (2)
!
!     Obviously Delta P_k/P_k decreases while k increases, since it
!     scales roughly like 1/k.
!
!     In the code, the number of independent modes, N_k, is output.
!     Another way of estimating the errors is also performed, by
!     simple estimation of the scatter in each bin. It turns out in
!     practice to give errors of the same order of equation (2), which
!     gives the results for an underlying random Gaussian field.
!     Actual calculation of the errors involve non Gaussian terms,
!     so the errors provided by the code are only approximations that
!     should be improved on for really accurate estimates of errorbars.
!
! (b) the biases due to the Taylor expansion and the discreteness of
!     the particle distribution, namely
!     (i) the so called 1/N_part shot noise bias where N_part is the number
!         of particles. In the Fourier Taylor method, this bias is
!         actually slightly different from 1/N_part, but this is estimated
!         accurately by the code. Note that for a slightly perturbed grid
!         pattern, such a bias does not arise and should not be corrected
!         for. Only in a relaxed, locally Poissonian stage, the correction
!         for shot noise makes really sense.
!     (ii) the bias due to the interpolation involved in the calculation.
!         It depends on the order Norder of the Taylor expansion.
!         For example, for NGP, Norder=0, the bias is given by the
!         famous sinc^2(k/2) term. This bias can be easily corrected
!         for at any order. It is performed in the code.
!     (iii) The unknown effects of aliasing of the contributions of the
!         power-spectrum at scales smaller than the sampling imposed
!         by the grid used to perform the calculations. In the limit
!         of a locally Poissonian stationary random process, these effects
!         of aliasing can be bounded. The code computes useful numbers
!         to estimate the maximum effect of aliasing (see Colombi et
!         al. paper for details). The effects of aliasing increase when
!         approaching the nyquist frequency of the sampling grid, so
!         it is wise to stay sufficiently away from the Nyquist frequency
!         to avoid too much aliasing.
!
! Obviously, by combining foldings and using (a)+(b), one can find
! the best compromise that minimizes the error on the estimate of P_k for  
! a given k at a minimum computational cost. However this is not performed
! at the present stage.
!
! In its current form, the code computes P(k) for each folding of the particle
! distribution. It computes as well shot noise contribution (b-i), and
! proposes a correction for the biase (b-ii) and usefull numbers for
! estimating (b-iii). Of course, the statistical errors are also estimated,
! (a). It outputs, for each folding, detailed data. However, for practical
! use, a global file contains all the combined information of all the
! foldings and this is what really matter for the user who does not
! necessarily want to have fully optimal results, but still fairly
! accurate ones. Right now, only that part of the code is fully explained
! in what follows.
!
!-----------------------------------------------------------------------
!
! COMPILATION
! +++++++++++
!
! Goto the directory /bin and modify the Makefile according to the
! instructions in it, type make.
!
! Notes :
! -------
! - The program needs FFTW2, that has to be installed.
! - It is possible to use openMP acceleration on shared memory computers,
!   in that case the thread part of FFTW2 has to be properly installed.
! - Possible problems of underscore while calling FFTW2, follow the
!   instructions in the Makefile.
! - Be aware of the little/big endians problems for input file reading.
!
!-----------------------------------------------------------------------
!
! EXECUTION
! +++++++++
!
! Create a config file, following the example given in powmes.config, say
! myconfigfile.config. Then type the command
! whereispowmes/powmes myconfigfile.config, where whereispowmes
! is the location of the executable.
! Typing just the command powmes assumes that the name of the config file
! is powmes.config.
!
!-----------------------------------------------------------------------
!
! PARAMETERS IN THE CONFIG FILE
! +++++++++++++++++++++++++++++
!
! verbose (.true./.false.) : verbose mode
!
! megaverbose (.true./.false.) : detailed verbose mode
!
! filein ('myfile') : name of the input file
!
! nfile (1/2/3) : type of the input file
!     1 : GADGET file
!     2 : RAMSES file
!     3 : 4 columns ascii file :
!         first line : npart where npart is the number of particles
!         next lines : x(i),y(i),z(i),mass(i) with x,y,z in [0,1[
!
! nmpi (-1/positive integer) : number of MPI files for nfile=1
!     (this parameter is thus ignored for nfile=2 or nfile=3)
!     set nmpi=-1 if there is only one file
!
! read_mass (.true./.false.) : .true. for reading the mass in the
!         file, .false. for ignoring it and giving the same mass to
!         all the particles. The parameter read_mass is ignored
!         if nfile=1 or nfile=2, where all the masses are supposed
!         to be the same for all the particles. WARNING : the code
!         has not been tested for unequal masses : this is still under
!         development.
!
! nfoldpow (integer) : folding parameter
!     - if nfoldpow >= 0, it corresponds to the number
!     of foldings to be performed. Note that the largest value of k
!     that will be probed will be kmax=k_nyquist*2^(nfoldpow-1), where
!     k_nyquist is the nyquist frequency of the grid used to perform the
!     calculations.
!     - if nfoldpow < 0, then kmax=-nfoldpow is the maximum integer
!     wavenumber up to which one wants to measure the power-spectrum.
!     The number of foldings needed to achieve that is computed
!     automatically.
!
! ngrid (even positive integer): resolution of the grid used to
!     perform the calculations.
!     ngrid must be an even number (but not necessarily a power of 2).
!
! norder (positive integer): order of the Taylor expansion.
!     Norder=3 is the recommended value.
!
! shift(float,float,float) : vector corresponding to a constant shift to
!     be applied to all the particles prior to the calculation. This
!     can be handy to improve the accuracy of the measurements for a
!     slightly perturbed grid. shift is expressed in units of the grid
!     cell size.  
!
! filepower ('powspec.dat') : the name of the main output file where all
!                the informations of the measured power-spectrum are gathered
!                in a handy (but not necessarily optimal in terms of signal
!                to noise or biases) way.
!
!     This is a 6 columns ascii file :
!
!     column 1 : the integer wave number
!     column 2 : the number of statistically independent modes C(k)
!                contributing to the bin k. This number can be useful to
!                compute errorbars : the statistical error
!                on the rough power-spectrum
!
!                            P_rough(k)=P(k)+1/Npart  
!
!                is, under the Gaussian field hypothesis,  
!
!                           DP_rough(k)/P_rough(k)=1/sqrt(C(k)) (3)
!
!     column 3 : the estimated rough power spectrum PP_rough(k) before
!                debiasing from Taylor interpolation.
!     column 4 : the debiased estimated P_rough(k). To obtain the true
!                power-spectrum, one can subtract the white noise contribution
!                (except for a slightly deformed grid, e.g. initial conditions)
!
!                           P(k)=P_rough(k)-W(k)/N_part, (4)
!
!                where W(k) is given by column 5 and is a number very
!                close to unity.
!
!     column 5 : the shot noise factor W(k) very close to unity, needed
!                to correct for shot noise contribution the power spectrum
!                given by equation (4).
!
!     column 6 : the statistical error computed self-consistently from
!                the dispersion inside each k bin.
!                The number displayed is DP_rough(k)/P_rough(k) and
!                should give in practice something very similar to
!                equation (3).
!            
! filepower_fold ('powspec'/'#junk') : root name of the output files where
!     the measurements are stored separately for each folding. If the
!     name of the file starts with a '#", these optional files are not
!     generated.
!    
!     Example : if filepower_fold='powspec', then 4 files are created :
!
!     powspec.waven  : the wavenumbers          (1st column of filepower)
!     powspec.nmodes : the number of modes file (second column of filepower)
!     powspec.power  : the rough power-spectrum, prior to debiasing
!                                               (third column of filepower)
!     powspec.powerdebiased : the debiased rough power-spectrum
!                                               (fourth column of filepower)
!     powspec.shotfac: the shot noise factor W(k), needed to correct
!                      for shot noise contribution when relevant
!                                               (5th column of filepower)
!     powspec.staterr: the statistical error    (6th column of filepower)
!          
!
!     More specifically :
!
!     powspec.waven : an ascii file with nfoldpow+1 columns, each column
!           corresponding to a given folding (first column : no fold,
!           second column : 1 fold, etc). In each column, the integer
!           wavenumber probed by the calculation is written.
!
!           Example : for ngrid=8 and nfoldpow=3 powspec.waven
!           looks as follows :
!                
!                 0    0    0    0
!                 1    2    4    8
!                 2    4    8   16
!                 3    6   12   24
!                 4    8   16   32
!
!     powspec.nmodes : similar as powspec.waven, but the number of modes
!           in the corresponding k-bin is output
!
!     powspec.power : similar as powspec.waven, but the measured value of
!           the rough power-spectrum prior to debiasing is written in
!           each column.
!
!     Similarly for powspec.powerdebiased, powspec.shotfac,powspec.staterr.
!
! filetaylor ('powspec.taylor'/'#junk') : name of the file containing
!     useful informations for the fourier-taylor calculation. This file
!     is optional and is not output if the name starts with a '#'.
!     This is a multicolumn ascii file, each line corresponds to an
!     integer value of k, starting from 0.
!
!     Column 1 : the number of independent modes C(k)
!     Column 2 : an estimate of L(k)=W(k)*U(k) where U(k) is explained
!                below
!     Column 3 : approximation for L(k) valid at small k
!     Column 4 : the estimated bias U(k) on the rough power spectrum
!                Basically, debiasing the power-spectrum (i.e. passing
!                from column 3 to column 4 of filepower) consists in
!                multiplying the rough power-spectrum by 1/U(k)
!     Column 5 : The W(k) already mentionned earlier.
!
!     For more details on filetaylor, ask to colombATiap.fr
!
!=======================================================================