Frontiers in Computational Astrophysics: The Origin of Stars, Planets and Galaxies

An ESF exploratory workshop

Here we describe several tests that we would like participants to try in readiness for the meeting. Please attempt what you can. At the meeting we will spend some time discussing results and to refine the tests that we think highlight the positive and negative characteristics of different codes. We may also come up with more appropriate tests at the meeting. For now you should look at standard diagnostics, such as temperature, pressure or density maps through central slices etc. We can decide on detailed diagnostic comparisons later. We hope to eventually publish a joint paper summarising our efforts.

We decided to set up the tests using particles and provide a code for transforming particle distributions to an accurate realisation on a grid. This tool currently outputs binary files that can be read using the HDF tools. We have given particle IC's at several different resolutions which you can run directly. If you are transforming to a grid, please use the highest resolution particle IC's as an input to the conversion tool and you can specify the size of the output grid. If you run the simulations at several resolutions we can try to determine at which resolution different codes agree - please run these at the highest resolution you can in the time we have left.

The particle initial conditions are provided in both standard tipsy RPC-binary (no endianness problems) format and in ascii format. TIPSY is a nice particle visualization tool written by Neal Katz and Thomas Quinn. There are a few common conversions to other formats to be found following the above link. If you need to write your own conversion macro, the structure of the tipsy ascii file is described here. One odditiy is the hsmooth field for the gas particles which, by its name, should contain the smoothing length for SPH actually contains the gravitational softening for the particle (the same as the eps field for dark and star particles). The reason for this is that the gravitational softening is a more fundamental parameter, while the smoothing length is typically determined from the 32 nearest neighbors in our case. TIPSY gas particles should really have both fields, but there are always historical constraints with data formats.

tipgrid & smooth

# tipgrid –

Will give you a brief description of the cmdline arguments, and I will try to put together a man page for this tool in the next week. Also, there will likely by several updates to this code to strengthen the types of visualizations it can produce (right now just projected density) adding slices, max-in-column and other variables (temperature, entropy, pressure).

The code TIPGRID creates HDF5-ready files in binary, but can also produce ascii files if requested. This conversion to grid is broken up into 3 steps:

Assume the provided input particle file (in standard TIPSY binary format) is called test3.bin in what follows.

• Step 1: Calculate the SPH smoothing lengths for the particles in the file (a quantity which is not normally stored within the tipsy binary format).

# smooth -s 32g -px 2000 -py 2000 -pz 8000 hsmooth -o test3 < test3.bin

• Replace gravitational softenings in the file with the smoothing lengths calculated by the above (note: don't use the resulting output binary for a simulation with gravity).

# hsmtoeps test3.bin test3.hsm > test3.eps

• Step 2: Create hdf5-ready files for a grid of given dimensions from the particle file test3.eps. This step can take quite a while for large grids so try a small one first with a larger projected density plane (-m 4) to see that all is ok before going for a big grid (like -n 256 256 1024). Type tipgrid – to see all the options.

# tipgrid -v -n 64 64 256 -l 2000 -xpd 1 -ypd 2 -p -dmin 0 -tmin 0 -msol 2.3262e5 -m 4 -o test3 < test3.eps

• Step 2a: Alternatively, create ascii output files of the 5 arrays by adding the -ascii flag to the tipgrid command line. This produces an additional 5 files as well as a gzipped tar file containing these 5 ascii files. We provide this .tgz (gzipped tar file) for the grid intial conditions below. The x-dimension in these ascii files is again the most rapidly changing. To extract the individual ascii files from the tar files given below just type the following.

# tar xfvz test1-128.gasc.tgz

• Step 3a: Create hdf-files from the output of tipgrid using the h5import utility from the hdf5 package. The tipgrid code produces one binary file of floats and a corresponding configuration file for each of the 5 Euler variables. In addition it produces another binary file which is a 2-D projection plane (column density) of the box and associated configuration file which should also be converted to hdf5 using h5import. This projected density plane is produced to make pretty pictures and to quickly assess the results of tipgrid (see step 3b).

# h5import test3.grid.d -c test3.cfg.d test3.grid.t -c test3.cfg.t test3.grid.vx -c test3.cfg.vx test3.grid.vy -c test3.cfg.vy test3.grid.vz -c test3.cfg.vz -o test3.h5

• This is a bit long and tedious, but simple. At this point you have an hdf5 file, a grid realization of the particle distribution. To projected density plane outputs are handled in the same way.

# h5import test3.plane.d -c test3.pcfg.d -o test3.proj.h5

• Step 3b: To visualize the projected density plane you will need the h5topng utility in the h5-utils package. Then you can use xv or any other image viewer to see the results. The -c parameter specifies the desired colormap; there are quite a few.

# h5topng -c jet test3.proj.h5 # eog/xv test3.proj.h5.png

For each of the four tests we provide this sequence of commands to generate the grid initial condition.

Units used in the simulations are based on a G=1 system, which provides a convenient system when the length and velocity units are in kpc and km/s respectively; namely that the time unit is very nearly Gyrs. This system of units is thus very practical in astrophysical applications. For converting to other unit systems the following should be helpful: the time unit is given by the gravitational timescale t_unit = 1/sqrt(G_physical_value * density_unit); the energy per unit mass unit is given by energy_per_unit_mass_unit = G_physical * mass_unit / length_unit. Thus the gas constant in system units is given by gas_const = (length_unit/(G_physical_value*mass_unit))*k_boltzman/mass_hydrogen(or amu).

To convert temperature to thermal energy U = gas_const/(mean_molar_mass)/(gamma - 1)*Temperature. The mass_unit = 2.3262e5 solar masses, the length_unit = 1 kpc, the velocity_unit = 1 km/s, the time_unit = 0.9777 Gyrs. For SPH codes one can convert masses, positions, velocities, softening radii, and temperatures to other system units if this is required using the above relations. For GRID codes we provide the 5 Euler equation quantities in these units (density, velocity in x, y and z, and thermal energy) as 5 separate one-dimension arrays which can be mapped to 3-dimensional arrays using A3D(ix,iy,iz) = A1D((iz-1)*(nyDim*nxDim) + (iy-1)*(nxDim) + ix), i.e., x-dimension changes most rapidly. We will also specify the grid dimensions in these cases.

As mentioned, for very high resolution grid realizations one should use the provided tool to produce an accurate grid realization from the SPH particle initial conditions. Plots produced in the above system of units are preferred, as this will make comparisons more straight forward.

We are meant to solve only the Euler equations so we will neglect physical viscosity in all the tests.In addition all tests should be performed as adiabatic. What we mean by adiabatic is what is the standard practice in the field of structure/galaxy formation. Therefore for adiabatic we mean the gas can undergo heating and cooling by adiabatic compression/decompression + irreversible heating by shocks (achieved through an artificial viscosity term in the thermal energy equation for the case SPH codes like ours). In other words we allow conversion of kinetic energy into thermal energy. Entropy will be generated in shocks so adiabatic here is not equivalent to isentropic.

Our code employs a spline kernel for the gravitational softening (the same type of function used in the calculation of the SPH smoothing length). The gravitational force between two particles becomes keplerian at twice the softening length (see for example the Appendix in Hernquist & Katz 1989 (ApJ, 70, 419) for a detailed description of the function). Therefore softening lengths quoted in this page are always given assuming a spline kernel. Slightly different lengths will have to be specified with different softening functions in order for the gravitational force resolution to be comparable (for example for a Plummer softening function the length should be about 30% lower than the spline).