Particle-Particle Particle-Mesh
Luty et al. [37] and Rajagopal et al. [40] offered an alternative approach to the Ewald
sum by extending the Particle-Particle Particle-Mesh method ( PPPM)
developed by Hockney and Eastwood [31]. The method relies on expressing the
long-range inter-particle force as the sum of two components: the short-range
force, which is only nonzero within some cutoff radius, and the ``reference''
force, that is long-ranged and smooth and can be approximated on a grid. The
analogy between PPPM and the Ewald sum is clear. In the traditional Ewald sum,
the direct sum, due to the point-charge and Gaussian counter-distributions, is
also short-ranged; while the reciprocal sum, due to the Gaussian distributions,
is a smooth function and its Fourier transform converges rapidly. The authors [37] used PPPM's standard charge distribution, a sphere with a uniform decreasing density (the S2 function), rather than the Gaussian distribution used for the Ewald sum. The S2 distribution is given by:
where a is a parameter that adjusts the S2 distribution.
The short-range potential between two particles, each with an S2 charge distribution, is evaluated with a cutoff radius,
,
given in terms of 
by:
where
The long-range potential is evaluated in Fourier space using :
where ``
''
indicates the Fourier transform of a function. The influence function is
usually given by 
but
can be optimized depending on the system size, charge shape function, and the
interpolation function. The long-range potential is computed using the
following steps, see Fig.(4): - Assign charge to a 3D grid that fills the simulation space. This
step yields the charge distribution
which is a function of both the charge distribution 
and assignment functions. Several charge assignments can be used
depending on the accuracy desired, e.g. triangle-shaped charge function
(TSC). However, for a charge assignment function to be feasible, it has to
cover a relatively small number of grid points. Moreover, the assignment
should vary smoothly with the particle's location, which requires adequate
grid spacing and in turn determines the runtime cost. - Using Fourier transform over the grid, obtain
, and calculate 
using Eq.(23). Apply the inverse Fourier Transform to
obtain 
evaluated at the grid points. - Obtain the grid-defined electrostatic forces by differentiating the
potentials numerically. Several numerical differencing techniques can be
used such as 4-point central differencing.
- Interpolate the electrostatic fields (potentials) from the grid to
particle locations using the same function of step 1.
Fast Fourier
Transforms are used in PPPM methods resulting in an 
algorithm.
It should be noted that in using this approach, the influence function 
is
system-specific and hence for each new system, and depending on system
parameters, e.g. size or charge shape, a new optimal influence function has to
be computed resulting in some loss of generalization. In addition, the PPPM
implementation of [37] implies that in order to increase the
accuracy of potential/force computations, one needs to either refine the mesh,
or use a better weighing/interpolation scheme; both choices can be
computationally expensive. In addition, since the electrostatic force
experienced at a grid point is obtained by numerically differentiating the
potential, another source of error is introduced to the results. In order to
improve the accuracy and reduce the error in the above steps, higher-order
differencing schemes have to be incorporated. There is a rich interaction
between all of the above choices ( i.e. weighing/interpolation,
differentiation, charge shape, influence function, etc.) and the user has to
experiment with the system of interest in order to utilize this algorithm
efficiently.
algorithm.
It should be noted that in using this approach, the influence function is
system-specific and hence for each new system, and depending on system
parameters, e.g. size or charge shape, a new optimal influence function has to
be computed resulting in some loss of generalization. In addition, the PPPM
implementation of [37] implies that in order to increase the
accuracy of potential/force computations, one needs to either refine the mesh,
or use a better weighing/interpolation scheme; both choices can be
computationally expensive. In addition, since the electrostatic force
experienced at a grid point is obtained by numerically differentiating the
potential, another source of error is introduced to the results. In order to
improve the accuracy and reduce the error in the above steps, higher-order
differencing schemes have to be incorporated. There is a rich interaction
between all of the above choices ( i.e. weighing/interpolation,
differentiation, charge shape, influence function, etc.) and the user has to
experiment with the system of interest in order to utilize this algorithm
efficiently. 
Figure 4: A 2D schematic of particle-mesh technique used in most Fourier-based methods (a) A system of charged particles. (b) The charges are interpolated on a 2D grid. (c) Using FFT, the potential and forces are calculated at grid points. (d) Interpolate forces back to particles and update coordinates.
Mon Jan 22 12:05:30 EST 1996
Next: Fast Fourier Poisson Up: Fourier-based Ewald Summation Previous: Particle-Particle Particle-Mesh
Particle-Mesh Ewald
The Particle-Mesh Ewald method (
PME) [14] is also inspired by Hockney and Eastwood's
particle-particle particle-mesh method (PPPM) [31,20]. Unlike PPPM, PME divides the potential
energy into Ewald's standard direct and reciprocal sums and uses the
conventional Gaussian charge distributions. The direct sum, Eq.(3), is evaluated explicitly using cutoffs while
the reciprocal sum, Eq.(4), is approximated using FFT with convolutions
on a grid where charges are interpolated to the grid points. In addition, in
contrast to particle-mesh methods, PME does not interpolate but rather
evaluates the forces by analytically differentiating the energies, thus
reducing memory requirements substantially. This method is reported to be highly efficient incurring only
overhead
over conventional truncated list-based (i.e. non-Ewald) methods at a relative
force accuracy around 
.
PME is also capable of achieving higher accuracy (
relative force error) with relatively little increase in computational cost. In computing the direct sum, the Ewald parameter
is
chosen large enough so that a fixed cutoff radius can be applied thus reducing
the complexity of the direct sum from 
to
.
To compensate for the truncation in evaluating the direct sum, the number of
reciprocal vectors is increased proportionally to N to bound errors. The reciprocal sum is computed using 3D-FFT with an overhead that grows as
.
PME is therefore an 
method.
The reciprocal sum, Eq.(4), is given by: 
where
is
defined in Eq.(15). The structure factor can be approximated by:
where
is
the 3D FFT of Q, the charge matrix. The Q matrix is a three
dimensional matrix that is obtained by interpolating the point charges to a
uniform grid of dimensions 
that
fills the simulation cell. By combining Eq.(25) with Eq.(24), the reciprocal energy can be also approximated
by: 
The above equation is rewritten, after some manipulation, as a convolution:
where
is
the reciprocal pair potential and ``
''
indicates a convolution. To evaluate the reciprocal sum, the Q matrix is
first computed over a 3D uniform grid and then transformed using inverse 3D FFT
to obtain the structure factors. The reciprocal energy is then calculated using
Eq.(27) with the aid of FFT to compute the
convolution 
.
The charge interpolation function used originally in PME was Lagrange interpolation [14]. However, an enhanced PME [15] utilizes the B-spline interpolation function, which is smoother and allows higher accuracy by simply increasing the order of interpolation. The smoothness of B-spline interpolation allows the force expressions to be evaluated analytically, with high accuracy, by differentiating the real and reciprocal energy equations rather than using finite differencing techniques.
Next: Fast Fourier Poisson Up: Fourier-based Ewald Summation Previous: Particle-Particle Particle-Mesh
Mon Jan 22 12:05:30 EST 1996
Abstract
for: Cosmological Simulations
Using
Adaptive Particle-Mesh Methods
Hugh
Couchman, University of Western Ontario
The complete manuscript is available in postscript. The code is
available from the Hydra consortium. It contains a full N-body code with potential
solver and time evolution. The final state of a Dark Matter simulation using Couchman's
code has been created using the group's visualization program: TIPSY.
It is easier
to see the three dimensional structure in a MPEG movie of a rotation about the z axis. This movie was also created using TIPSY.
The purpose
of cosmological simulations is to model the growth of structure in the
universe. The techniques described here model the mass density of the universe
with a distribution of softened particles. Evolution is simulated by following
the trajectories of these particles under their mutual gravity. Since the
region of the universe modelled is typically very much smaller than the horizon
volume, Newtonian equations of motion are adequate to describe the particle
dynamics and evolution of structure. The resulting numerical schemes are
straightforward and robust. An alternative approach is to describe the mass
distribution as a fluid and calculate the evolution from the appropriate fluid
equations (e.g., Peebles 1987). A number of other analytic and semi-analytic
approaches to gravitational evolution have been discussed. The description of
the growth of small perturbations in the matter distribution is well established,
both to first order and at higher orders (see, for example, Peebles 1980). Many
observational constraints, however, are the result of highly non-linear
gravitational evolution for which perturbation analysis is of limited use.
Numerical N-body techniques offer a simple and effective tool for investigating
non-linear cosmological gravitational evolution.
Many of the
results arising from potential theory and the theory of gravitating systems of
particles, as discussed in this volume and elsewhere, are directly applicable
to the cosmological situation and will not be addressed here in further detail.
The principal features peculiar to cosmological simulations arise from the
construction of an initial particle distribution, the necessity of describing
evolution in an expanding cosmic background and the choice of suitable boundary
conditions. The remainder of the introduction addresses these issues and
describes the popular mesh-based cosmological particle methods. Two computer
codes are discussed in this chapter: the first is used to generate an initial
distribution of particles; the second to evolve it. The algorithms are
described in Sect. 2 and the FORTRAN77 implementation in Sect. 3. The reader
interested in trying the codes immediately may skip to Sect. 4 where a test
case is presented.
__________________________________________________________________
www-hpcc@astro.washington.edu
Particle-Mesh
Method
This project is the first phase of building a complete particle mesh code. The steps in a particle mesh code are:
1. Construct the Green's function of the form 1/sqrt(r*r + eps*eps) where eps is a softening length for the potential
2. Find the Fast Fourier Transform (FFT) of the Green's function.
3. Main Loop
- Map the particle masses onto a mesh.
- Compute the FFT for mesh.
- Form a convolution of the Green's function and the particle mesh as
an element-by-element product of the FFTs.
- Take the inverse FFT of the product. This is the potential.
- Interpolate forces from the potential to the particles.
- Determine the time step.
- Update and move the particles.
- Calculate the kinetic, potential and total energies for the system.
Source Code
Example Potentials (all have eps=1.0, mass=1.0)
Particle at the center:
Two particles
Test of 256
particles with random starting positions and zero initial velocity
Plot of Energies over Time
Plot of particles for time=0
Plot of particles for time=4.71 (first
collapse)
Plot of particles for time=8.29 (first
expansion)
Plot of particles for time=13.32 (second
collapse)
This run was profiled using 'prof'. The output clearly shows that the FFT routine is responsible for the majority
(79%) of the execution time of the program.
 |
Astrophysics Source Code
Library
|
PMCODE:
Particle-Mesh code for cosmological simulations
Abstract: Particle-Mesh (PM) codes are still very useful
tools for testing predictions of cosmological models in cases when extra high
resolution is not very important. We release for public use a cosmological PM
N-body code. The code is very fast and simple. We provide a complete package of
routines needed to set initial conditions, to run the code, and to analyze the
results. The package allows you to simulate models with numerous combinations
of parameters: open/flat/closed background, with or without the cosmological
constant, different values of the Hubble constant, with or without hot
neutrinos, tilted or non-tilted initial spectra, different amount of baryons.
Subject headings: cosmology: theory - dark matter - large-scale structure of the universe
Latest Version: 1997 December 9 Subject headings: cosmology: theory - dark matter - large-scale structure of the universe
Archived: 1999 September 6
Papers: Klypin, A., Holtzman, J., Primack, J., Regos, E.; Astrophys.J.; 1993ApJ...416....1K
Klypin, A., Borgani, S., Holtzman, J., Primack, J., 1995ApJ...444....1K
Preprint: astro-ph/9712217
Language: Fortran
External Explanatory Pages: http://astro.nmsu.edu/~aklypin/PM/pmcode/index.html
Source Code: PMcode.tar.gz
NASA ADS Astronomy Abstract Service
·
Title:
|
Structure Formation
with Cold plus Hot Dark Matter
|
|
Authors:
|
||
Journal:
|
Astrophysical
Journal v.416, p.1 (ApJ Homepage)
|
|
Publication Date:
|
10/1993
|
|
Origin:
|
APJ; KNUDSEN
|
|
ApJ Keywords:
|
COSMOLOGY: DARK
MATTER, COSMOLOGY: LARGE-SCALE STRUCTURE OF UNIVERSE, COSMOLOGY: THEORY,
GALAXIES: CLUSTERING, GALAXIES: FORMATION, METHODS: NUMERICAL
|
|
Bibliographic
Code:
|
1993ApJ...416....1K
|
Abstract
Printing Options
More
Article Retrieval Options
HELP
for Article Retrieval
Astrophysics, abstract
astro-ph/0302065
From: Paul Bode <bode@astro.princeton.edu>
Date: Tue, 4 Feb 2003 21:49:40 GMT (87kb)
Tree-Particle-Mesh: an adaptive,
efficient, and parallel code for collisionless cosmological simulation
Authors: Paul Bode, Jeremiah P. Ostriker (Princeton University Observatory)
Comments: 29 pages, includes 12 figures; to appear in ApJS. Source code at this http URL
Journal-ref: Astrophys.J.Suppl. 145 (2003) 1-14
Comments: 29 pages, includes 12 figures; to appear in ApJS. Source code at this http URL
Journal-ref: Astrophys.J.Suppl. 145 (2003) 1-14
An improved
implementation of an N-body code for simulating collisionless cosmological
dynamics is presented. TPM (Tree-Particle-Mesh) combines the PM method on large
scales with a tree code to handle particle-particle interactions at small
separations. After the global PM forces are calculated, spatially distinct
regions above a given density contrast are located; the tree code calculates
the gravitational interactions inside these denser objects at higher spatial
and temporal resolution. The new implementation includes individual particle
time steps within trees, an improved treatment of tidal forces on trees, new
criteria for higher force resolution and choice of time step, and parallel
treatment of large trees. TPM is compared to P^3M and a tree code (GADGET) and
is found to give equivalent results in significantly less time. The implementation
is highly portable (requiring a Fortran compiler and MPI) and efficient on
parallel machines. The source code can be found at this http URL
Full-text: PostScript, PDF, or Other formats
References and citations for this
submission:SLAC-SPIRES HEP (refers to , cited by, arXiv reformatted);
CiteBase (autonomous citation navigation and analysis)
NASA ADS Astronomy Abstract Service
·
Title:
|
Damped lyman-alpha
systems versus cold + hot dark matter
|
|
Authors:
|
||
Affiliation:
|
AA(New Mexico State
University, Las Cruces, New Mexico, US), AB(New Mexico State University, Las
Cruces, New Mexico, US), AC(Lowell Observatory, Flagstaff, Arizona, US),
AD(Lowell Observatory, Flagstaff, Arizona, US)
|
|
Journal:
|
Astrophysical
Journal, Part 1 (ISSN 0004-637X), vol. 444, no. 1, p. 1-14 (ApJ Homepage)
|
|
Publication Date:
|
05/1995
|
|
Category:
|
Astrophysics
|
|
Origin:
|
STI
|
|
NASA/STI
Keywords:
|
COMPUTATIONAL
ASTROPHYSICS, COSMOLOGY, DAMPING, DARK MATTER, LYMAN ALPHA RADIATION,
UNIVERSE, ABSORPTION SPECTRA, APPROXIMATION, GALACTIC EVOLUTION, GALACTIC
HALOS, MANY BODY PROBLEM, NEUTRINOS, QUASARS, SIMULATION, SPACE DENSITY
|
|
Bibliographic
Code:
|
1995ApJ...444....1K
|
Abstract
Although the cold +
hot dark matter (CHDM) cosmology provides perhaps the best fit of any model to
all the available data at the current epoch (z = 0), CHDM produces structure at
relatively low redshifts and thus is very sensitive to the observed numbers of
massive objects at high redshifts. Damped lyman-alpha systems are abundant in
quasar absorption spectra and provide possibly the most significant evidence
for early structure formation, and thus a stringent constraint on CHDM. Using
the numbers of halos in N-body simulations to normalize Press-Schecter
estimates of the number densities of protogalaxies as a function of reshift, we
find that CHDM with Omegac/Omeganu/Omegab =
0.6/0.3/0.1 is compatible with the damped lyman-alpha data only at less than or
equal to 2.5, but that it is probably incompatible with the z greater than 3
damped lyman-alpha data. The situation is uncertain because there is very
little data for z greater than 3. The predictions of CHDM are quite sensitive
to hot (neutrino) fraction, and we find that Omegac/Omeganu/Omegab
= 0.725/0.20/0.075 (and possibly even Omegac/Omeganu/Omega
b = 0.675/0.25/0.075) is compatible with the z greater than 3 data.
With one massive neutrino species, using Omeganu = 0.20 instead of
0.30 corresponds to lowering the neutrino mass from 7.0 to 4.7 eV, for H0
= 50 km/s Mpc and T = 2.726 K. In CHDM, the higher redshift damped lyman-alpha
systems are predicted to have lower masses (approximately 3 x 1010solar
mass at z = 3), a prediction which can be checked by measuring the velocity
widths of the associated metal-line systems. Predictions for high-z objects
crucially depend on the effects of limited resolution and the finite box size
in N-body simulations or on the parameters of the Press-Schechter
approximation, if it is used. By analyzing our numerical simulations with
vastly different resolutions and box sizes as well as those of Ma &
Berschinger (1994), we show that for the CHDM models with Omeganu =
0.2-0.3 the Press Schechter approximation should be used with Gaussian filter
with deltac = 1.5 if halos are defined with the mean overdensity
larger than 200. If one tries to recover the total mass of a collapsed halo, a
better value for the collapse parameter is deltac = 1.40. We argue
that nonlinear effects due to waves both longer and shorter than those
considered in numerical simulations could probably result in deltac
as low as deltac = 1.3.
No hay comentarios:
Publicar un comentario