grafic++ Code
Brought to you by:
dougpotter
Menu
▾
▴
Tree [cb7a14] default tip / History
Doug Potter
dougpotter
Read Me
README.txt for grafic2 written by E. Bertschinger, March 18, 2001
Please send all questions, comments, bugs, etc. to edbert@mit.edu
*** This file contains important information for users of grafic2.
*** Please read and retain for reference as you run grafic2.
This package produces multiscale gaussian random fields on nested
Cartesian lattices. There are two programs: grafic1 and grafic2.
grafic1 produces the top-level periodic grid. grafic2 is used for
refinement. Refinement may be recursive to give subgrids within
subgrids like Russian matrioshka dolls.
The periodic top grid is called the level-0 grid. The first refinement
is called the level-1 grid. The level-n grid may be refined by grafic2
to produce the level-(n+1) grid.
Unavoidable numerical errors are introduced at each stage of refinement.
The errors are greatest for the velocity field and they arise from aliasing
of the tidal fields. For typical grids and a CDM-like spectrum the rms
errors are less than 5 percent at the level-2 grid.
For an example of grafic1 and grafic2, type 'make' and then 'make test'
Compare test_results/test.out and test_results/*.gif with the files in
tests_sample.
Procedure for a single refinement level.
----------------------------------------
1. Edit grafic1.inc to set the number of grid points in each dimension
for the level-0 grid (np1,np2,np3). Also set offsets for baryon and
cdm velocities (offvelb, offvelc) and the parameter sigstart that fixes
the initial redshift by setting the standard deviation of the density
fluctuation at the highest level of refinement.
2. Edit the Makefile to replace the f77 compiler name if necessary
(e.g. for Linux, F77 = g77) and to use your favorite compiler
flags for optimization, etc. Then type 'make grafic1'
3. Run grafic1 interactively once or twice to see what parameters must
be entered. The memory requirement is a little more than (np1*np2*np3)
words. Look at grafic1.inp to see an example of the input parameters.
4. The first parameters entered specify the transfer function calculation
for the baryon and CDM power spectra. (These are allowed to be different,
although they are set equal unless one uses LINGERS.)
Transfer functions for baryons and CDM may be precomputed by running
LINGERS (a modification of LINGER_SYN from the COSMICS package that
outputs multiple times, allowing interpolation of the transfer function
to an arbitrary redshift). Alternatively, transfer functions for CDM
models may be given by fitting formulae (BBKS or enter your own into
power.f), or scale-free power-laws may be used.
***BEWARE*** the BBKS transfer function significantly underestimates
the amount of large-scale power for currently favored models.
It assumes omegab=0. For realistic models one should run LINGERS,
which is provided with this package.
Scale-free power spectra are not set up in grafic to automatically
normalize multiscale initial conditions. For scale-free initial conditions
the user will have to increase the normalization k^3P(k) (a parameter
input to grafic1) at each refinement level by a factor nref**(3+n).
5. For realistic models, grafic1 can normalize the power spectrum either
with the CMB quadrupole using Q_rms-PS or with sigma_8. Whatever your
choice, grafic1 will calculate the other choice and inform you. For
accurate transfer functions computed with LINGERS, I recommend CMB
normalization with Q_rms-PS = 18 micro-K for a n=1 spectrum, following
Hinshaw et al 1997, ApJ 464, L17. Careful simulators will also run
CMBFAST to get the shape of the CMB power spectrum. Beware that CMBFAST
normalizes using the Bunn-White prescription; however, you can determine
Q_rms-PS from C_2.
6. For multiscale initial conditions, the user will need to generate
density and velocity fields for each of the grids and refinement levels
that are used. In the case of a single refinement level, there are
two sets of data files: one for the top (level-0) grid and one for the
level-1 refinement. grafic1 must be run once for each --- once to produce
output for the level-0 grid and once to produce input to grafic2 which
will produce the level-1 output. When running grafic1 (and similarly
grafic2), the user must inform the program whether to produce output
fields at the current refinement level or to produce input for further
refinement. grafic1 asks the user to specify which case is being run.
Note that output fields are smoothed with a Hanning filter (spherical
cosine window function that goes to zero at the one-dimensional Nyquist
frequency).
7. If further refinement will be performed, grafic1 prompts the user
to enter the size of the subgrid (in units of the top grid spacing dx)
to be extracted for refinement, and grid offsets to give the position
of its corner. These are all integers. The subgrid must be no larger
than half the top grid size (np1,np2,np3). The offset may be positive
or negative, from -0.5*(np1,np2,np3) to +0.5*(np1,np2,np3). Periodic
boundary conditions apply on the top grid.
8. In order to reduce numerical aliasing errors with multiple refinement
levels, grafic1 then prompts the user to enter the final subvolume size
and offsets (integers, with units given by the top grid spacing dx).
For refinement to level 1, these are the same as the "next grid"
values (see grafic1.inp).
9. The last step for grafic1 is to enter information for random numbers.
It is strongly recommended that the user save the white noise files
created and used by grafic1 (and grafic2). For multiple refinement,
the same random numbers must always be used for all refinements across
a given level of the hierarchy; saving the files and reusing them
can ensure this. Also beware that the random number generator
will not necessarily produce the same results (for a given seed)
with 32-bit and 64-bit processors. If the user wishes, s/he may
provide the random number routine (by replacing randa with a uniform
(0,1) pseudorandom generator) or the white noise files. Note that
the random fields are sampled in real space, not Fourier space.
10. grafic1 should be run once to produce final output files at the
top grid level. These will be called ic_deltab, ic_velbx, etc.
They should be put into a separate directory labeled "level_0" to
avoid confusion with the files that will be produced for refined
grids by grafic2. These output files provide the information
needed for comsological simulation codes. They are written as follows:
open(11,file='ic_deltab',form='unformatted')
rewind 11
write(11) np1,np2,np3,dx,x1o,x2o,x3o,astart,omegam,omegav,h0
do i3=1,np3
write(11) ((deltab(i1,i2,i3),i1=1,np1),i2=1,np2)
end do
close(11)
The velocity files are written the same way. The parameters are
np[1-3] = grid size, dx = grid spacing in comoving Mpc,
x[1-3]o = grid offsets in comoving Mpc (=0 unless offvelb/c is set nonzero
in grafic1.inc), astart is the starting expansion factor (with a=1 today),
omegam,omegav,h0 are self-explanatory (H0 has units of km/s/Mpc).
Note that deltab is the dimensionless density fluctuation in real space.
The velocity (for ic_velbx, etc.) is the peculiar velocity in units of
proper km/s at astart. grafic1 and grafic2 no longer output particle
positions like their predecessor grafic does. See the end of this file
for information on how to compute them from the velocities.
11. grafic1 should be run a second time to produce intermediate output
for the level-1 grid. This time the user must specify the subvolume
that will be extracted (size and offset). Both of these quantities are
given in terms of the top grid spacing. Thus, if the top grid is a
120 Mpc cube np1=np2=np3=240, and one wishes a 40 Mpc cubical subvolume
with 4 times higher resolution centered at (0,30,60) Mpc, then the
subgrid size is n1c=n2c=n3c=80 and the offsets are m1off=-40, m2off=20,
m3off=80. The refinement factor nref=4 isn't needed by grafic1.
Save it for grafic2. Note that the same random numbers must be used
for both the first and second runs of grafic1. They determine the
structure in the box.
12. The second run of grafic1 will produce an output file grafic2.top.
This should be moved to a new directory called "level_1". It is the
input data file needed by grafic2.
13. Standard output of grafic1 contains useful information about the
random numbers and the standard deviation of each field it computes.
When run in the "intermediate output" mode to produce grafic2.top,
the velocity fields are split into "inner" and "outer" parts due to
fluctuations inside and outside of the subvolume. The user may
generally ignore the standard output (except when warnings and errors
are written), although the sophisticated user will find useful
statistics.
14. To generate the level-1 subgrid, edit grafic2.inc to include the
subgrid size (n1c,n2c,n3c) (exactly as entered into grafic1) and
the refinement factor nref, i.e. the factor by which the grid spacing
dx will be decreased. The memory requirement for grafic2 is a little
more than 16*n1c*n2c*n3c*nref**3. Then make grafic2
15. Run grafic2 in the level_1 directory. Its input is similar to grafic1,
except that one must specify the grafic2.dat input filename from grafic1.
Normally this would be grafic2.top unless one changed the name. grafic2
doesn't need information about the final subvolume since these were already
given to grafic1 and saved in the grafic2.top file. To produce output
density and velocity fields at level 1, choose the correct mode.
grafic2 has 2 modes (final output, intermediate output) just like grafic1.
Finally, choose a different random number seed for grafic2. These random
numbers are used for the level-1 subgrid. grafic2 will produce output
files ic_deltab, etc exactly like grafic1. Keep these in the level_1
directory so that they won't be confused with the level_0 files.
(The headers contain the grid spacing dx for each level, allowing them
to be distinguished in case of confusion.)
16. grafic2 takes a large amount of cpu time; it performs 25-39 FFTs of
size 8*n1c*n2c*n3c*nref**3. For a size of 512**3, on the Origin 2000
the runtime is about 6 hours for the final output mode (25 FFTs).
17. grafic2 produces several temporary files (temp.*) that may be removed
after each grafic2 run. A large amount of scratch disk space is needed
for grafic2. Including the input and output files, the scratch disk
requirement can be as large as 52*n1c*n2c*n3c*nref**3 4-byte words.
Procedure for multiple refinement levels.
-----------------------------------------
1. If one wishes to produce initial conditions for more than two levels
of refinement, grafic2 may be run recursively. Producing output for
refinement level n requires running grafic1 once and grafic2 n times.
Thus, producing output at all levels (0, 1, ..., n) of an n-level
multiscale simulation requires n+1 runs of grafic1 and n(n+1)/2 runs of
grafic2. The user will have to keep track of the files at each level
and provide the correct subgrid information to grafic1 and 2 for each run.
It might be worthwhile writing a perl script to automate this.
2. For example, suppose that one has already done 2-level (level_0
and level_1) as described above, and now wishes to further refine
to produce a level_2 subgrid. The first step is to set the size and
offset of this subgrid relative to the level_0 grid. Subgrids of all
refinement levels must be aligned on level_0 cells; all subvolumes
of all refinement levels must have physical sizes that are a whole
number of top grid spacings across. Once the user has set the size
and location of the level_2 subgrid, the first step is to run grafic1.
grafic1 needs to know both the size and location of the level_1 subgrid
and the level_2 subgrid. (To produce output for level_n, grafic1 needs
to know the level_n subgrid as the "final" subgrid.) Note that the
same random numbers must be used with grafic1 as were used in the
earlier grafic runs for earlier levels of the grid hierarchy.
3. After running grafic1, the user edits grafic2.inc and sets
(n1c,n2c,n3c) equal to the level_1 subgrid size that was entered
into grafic1. The user also sets the refinement factor for level 1
and makes grafic2. grafic2 is then run with the output of grafic1
in the "intermediate output" mode, which requires specifying the size
and offsets of the level_2 subgrid that will be extracted. Be careful
here. The size is in units of the level_1 grid, and is therefore
larger by a factor of the level_1 refinement factor than the level_2
grid size entered into grafic1. Similarly, the offsets given to
grafic2 are relative to the origin of the level_1 grid and are in units
of the level_1 grid spacing. The same random numbers must be used for
this run as were used in the previous grafic2 run at level_1. grafic2
will produce grafic2.dat.
4. The user then edits grafic2.inc again to enter the level_2 subgrid size
that was entered into grafic2 in step 3 as well as the refinement factor
for level_2, i.e. the factor by which level_1 is refined to yield level_2.
grafic2 is run in final output mode, and different random numbers should
be used than were used at level_0 and level_1. (Just name the random number
files wnsub_level0, wnsub_level1, etc., and keep careful track of the
refinement level, and you should have no problem.) The result will be
ic_* for level_2.
5. To generate output for higher levels, one repeats this pattern, extending
the chain of runs of grafic2 by one for each additional refinement level.
Computing particle positions from velocities.
---------------------------------------------
Particle positions follow readily from velocities, as illustrated by
the following f77 code fragment:
c (np1,np2,np3) are the size of the level_n subgrid.
real velx(np1,np2,np3),x(np1,np2,np3)
c
print*,'Enter top grid spacing in Mpc'
read(*,*) dx0
c Change the filename for the appropriate component;
c ic_velbx will give the x-coordinate of baryon particles.
open(11,file='ic_velbx',form='unformatted')
rewind 11
read(11) np1,np2,np3,dx,x1o,x2o,x3o,astart,omegam,omegav,h0
c velocity (proper km/s) = Displacement (comoving Mpc at astart) * vfact.
c vfact = dln(D+)/dtau where tau=conformal time.
c These functions (fomega, dladt) are in time.f, so link them in.
vfact=fomega(astart,omegam,omegav)
& *h0*dladt(astart,omegam,omegav)/astart
c Offset the positions so that the (dx0/dx)**3 particles corresponding
c to one level_0 particle are centered on the level_0 particle.
offset=0.5*(dx0-dx)
c Read each slice of the velocity field.
do i3=1,np3
read(11) ((velx(i1,i2,i3),i1=1,np1),i2=1,np2)
c Unperturbed grid position.
z0=(i3-1)*dx+offset
do i2=1,np2
y0=(i2-1)*dx+offset
do i1=1,np1
x0=(i1-1)*dx+offset
c For the y-component, replace x0 with y0.
x(i1,i2,i3)=x0+velx(i1,i2,i3)/vfact
end do
end do
end do
close(11)
Additional Notes
----------------
1. The initial density fluctuations will be smaller at the coarser grid
levels. The user's multiscale evolution code may have to compute
accurate gravitational forces when the density is almost uniform over
the top grid yet highly deformed at the highest level of refinement.
The user may wish to apply the Zel'dovich approximation at the top
refinement level while using a Poisson solver only for the refinements,
until the density fluctuations become sufficiently large. You should
check how well your code reproduces the Zel'dovich approximation at
early times, using the accurate velocity and displacement fields computed
by grafic1 and grafic2.
2. Multiple subvolumes may be used at any given refinement level.
For example, if the top grid spacing is 1 Mpc, then non-overlapping
subgrids of spacing 0.25 Mpc may be put at several places in the box.
However, the small-scale fields within one will be ignored when
computing the tides felt by the other. Thus, the subgrids should be
widely separated.
3. The refinement factor may be any whole number greater than 1. Different
refinement factors may be used between different levels.
4. A utility ic2gif is provided to produce images of the ic_* files.
This can be helpful in checking. It is created by typing 'make ic2gif'
Then run it as follows: 'ic2gif ic_filename gif_filename.gif'
where ic_filename is the data file to be images and gif_filename.gif
is the output gif file. ic2gif will image one slice.
