GADGET-2 testing simulation

GADGET-2 simulation requires some libraries in order to work.

HDF5: library to manage complex dataset of “all size”.

https://support.hdfgroup.org/ftp/HDF5/releases/hdf5-1.6/hdf5-1.6.9/src/hdf5-1.6.9.tar.gz

Download and extract HDF5 library 1.6.9 (avoid latest release) then go inside the folder.

./configure
make
make install

Then HDF5 needs two dependency libraries. If you haven’t homebrew yet on your system just type this on Terminal:

/usr/bin/ruby -e "$(curl -fsSL https://raw.githubusercontent.com/Homebrew/install/master/install)"

ZLIB

brew install zlib

SZIP

brew install szip

FFTW: subrutine library for computin discrete Fourier transform.

http://www.fftw.org/fftw-2.1.5.tar.gz

Download and extract FFTW library 2.1.5 (avoid latest release) then go inside the folder.

./configure
make
make install

GSL: scientific library.

brew install gsl

MPICH: MPI implementation - Message Passing Interface for parallel computing architectures.

brew install mpich

Now you can download GADGET-2: “A code for cosmological simulations of structure formation”.

http://wwwmpa.mpa-garching.mpg.de/gadget/gadget-2.0.7.tar.gz

Extract the .tar file and go inside the folder. You have a folder called “parameterfiles” where you find some parameter file including Makefile example depending on which simulation you want to run into. Keep copy of these file before you start to modify them. Let say you want to run into cluster simulation. You need to copy cluster.Makefile in parent folder and replace existing Makefile renaming it. Then you need to modify Makefile in some points.

Comment periodic boundary conditions “-DPERIODIC” in “Basic operation mode of code”

#OPT   +=  -DPERIODIC

In “Select target computer” you need to create you own SYSTYPE and comment with # the other.

SYSTYPE="mycomputer"
#SYSTYPE="MPA"

Then after “Adjust settings for target computer” you should copy a configuration from another SYSTYPE (let say MPA) and change it according to your configuration:

ifeq ($(SYSTYPE),"mycomputer")
CC       =  mpicc
OPTIMIZE =  -O3 -Wall
GSL_INCL =  -I/usr/local/Cellar/gsl/2.3/include
GSL_LIBS =  -L/usr/local/Cellar/gsl/2.3/lib
FFTW_INCL=  -I/Users/yourname/Documents/Gadget-2/fftw-2.1.5/include
FFTW_LIBS=  -L/Users/yourname/Documents/Gadget-2/fftw-2.1.5/lib
MPICHLIB =  -L/usr/local/Cellar/mpich/3.2_3/lib
HDF5INCL =  -I/Users/yourname/Documents/Gadget-2/hdf5-1.6.9/hdf5/include
HDF5LIB  =  -L/Users/yourname/Documents/Gadget-2/hdf5-1.6.9/hdf5/lib -lhdf5 -lz
endif

Libaray installed with brew should be the same path above. Please check only the version name. FFTW and HDF5 should be where you want. I put it inside Gadget-2 parent folder. Change yourname with your computer name. After these changes you can save the Makefile and

make

Probably you should see some warnings but you can skip them.

Next step is to start the simulation. Before you need to move ICs folder inside Gadegt2 folder and create “cluster” folder where the simulation puts the results.

mpirun -np 32 ./Gadget2 parameterfiles/cluster.param

“-np 32” is the number of processors you want to use. “./Gadget2” is the simulation and “parameterfiles/cluster.param” is the parameter file you want to use.

The simulation begins to calculate every step and periodically save a snapshot inside “cluster” folder. Snapshots are a binary file where you could extract many parameters of each particle simulated. The period of the step can be changed inside “cluster.param” file as well as the snapshot timing (if you want to specify manually this timing you need to set “OutputListFilename” with the name of your timing file:

OutputListFilename parameterfiles/outputs_cluster.txt

And put the time of each snampshot, one for each row:

0.166667
0.250000
0.333333 
0.5
1.0

And set in “cluster.param” the parameter “OutputListOn” to 1:

OutputListOn             1

For more information and to customize your simulation, check these links:

http://wwwmpa.mpa-garching.mpg.de/gadget/users-guide.pdf

Now you have your snapshots in binary files and you can use it converting them through the following code. “read_snapshot_to_txt.c”

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>



struct io_header_1
{
    int npart[6];
    double mass[6];
    double time;
    double redshift;
    int flag_sfr;
    int flag_feedback;
    int npartTotal[6];
    int flag_cooling;
    int num_files;
    double BoxSize;
    double Omega0;
    double OmegaLambda;
    double HubbleParam;
    char fill[256 - 6 * 4 - 6 * 8 - 2 * 8 - 2 * 4 - 6 * 4 - 2 * 4 - 4 * 8];    /* fills to 256 Bytes */
} header1;



int NumPart, Ngas;

struct particle_data
{
    float Pos[3];
    float Vel[3];
    float Mass;
    int Type;
    
    float Rho, U, Temp, Ne;

} *P;

int *Id;

double Time, Redshift;


/* Here we load a snapshot file. It can be distributed
 * onto several files (for files>1).
 * The particles are brought back into the order
 * implied by their ID's.
 * A unit conversion routine is called to do unit
 * conversion, and to evaluate the gas temperature.
 */
int main(int argc, char **argv)
{
    char path[200], input_fname[200], basename[200];
    int type, snapshot_number, files;
    
    
    sprintf(path, "./");
    sprintf(basename, "snapshot");
    snapshot_number = 0;        /* number of snapshot */
    files = 1;            /* number of files per snapshot */
    
    printf("read snapshot total: %i\n\n",atoi(argv[1]));

    for(int i = 0; i < atoi(argv[1]); i++){
        printf("read snapshot%d\n",i);
        sprintf(input_fname, "%s/%s_%03d", path, basename, i);
        load_snapshot(input_fname, files);
        reordering();            /* call this routine only if your ID's are set properly */
        unit_conversion();        /* optional stuff */
        saveSnapshotInTxtFile(i);
    }
}





/* here the particle data is at your disposal
 */
int saveSnapshotInTxtFile(int suffix)
{
    printf("do_what_you_want.\n");

    FILE *saved = stdout;
    
    char filename[200];
    sprintf(filename,"export%i.txt",suffix);

//    printf("writing on file: %c",filename);
    
    stdout = fopen(filename, "a");
    
    for (int i = 0; i < NumPart; i++){
        printf("[%i] %f %f %f %f %f %f %f %fs\n", Id[i], Time, Redshift, P[i].Pos[0], P[i].Pos[1], P[i].Pos[2], P[i].Vel[0], P[i].Vel[1], P[i].Vel[2]);
    }
    
    fclose(stdout);
    stdout = saved;
    printf("number paricles %i; time: %f and redshift: %f\n", NumPart, Time, Redshift);
    
    
}


/* this template shows how one may convert from Gadget's units
 * to cgs units.
 * In this example, the temperate of the gas is computed.
 * (assuming that the electron density in units of the hydrogen density
 * was computed by the code. This is done if cooling is enabled.)
 */
int unit_conversion(void)
{
    double GRAVITY, BOLTZMANN, PROTONMASS;
    double UnitLength_in_cm, UnitMass_in_g, UnitVelocity_in_cm_per_s;
    double UnitTime_in_s, UnitDensity_in_cgs, UnitPressure_in_cgs, UnitEnergy_in_cgs;
    double G, Xh, HubbleParam;
    
    int i;
    double MeanWeight, u, gamma;
    
    /* physical constants in cgs units */
    GRAVITY = 6.672e-8;
    BOLTZMANN = 1.3806e-16;
    PROTONMASS = 1.6726e-24;
    
    /* internal unit system of the code */
    UnitLength_in_cm = 3.085678e21;    /*  code length unit in cm/h */
    UnitMass_in_g = 1.989e43;    /*  code mass unit in g/h */
    UnitVelocity_in_cm_per_s = 1.0e5;
    
    UnitTime_in_s = UnitLength_in_cm / UnitVelocity_in_cm_per_s;
    UnitDensity_in_cgs = UnitMass_in_g / pow(UnitLength_in_cm, 3);
    UnitPressure_in_cgs = UnitMass_in_g / UnitLength_in_cm / pow(UnitTime_in_s, 2);
    UnitEnergy_in_cgs = UnitMass_in_g * pow(UnitLength_in_cm, 2) / pow(UnitTime_in_s, 2);
    
    G = GRAVITY / pow(UnitLength_in_cm, 3) * UnitMass_in_g * pow(UnitTime_in_s, 2);
    
    
    Xh = 0.76;            /* mass fraction of hydrogen */
    HubbleParam = 0.65;
    
    
    for(i = 1; i <= NumPart; i++)
    {
        if(P[i].Type == 0)    /* gas particle */
        {
            MeanWeight = 4.0 / (3 * Xh + 1 + 4 * Xh * P[i].Ne) * PROTONMASS;
            
            /* convert internal energy to cgs units */
            
            u = P[i].U * UnitEnergy_in_cgs / UnitMass_in_g;
            
            gamma = 5.0 / 3;
            
            /* get temperature in Kelvin */
            
            P[i].Temp = MeanWeight / BOLTZMANN * (gamma - 1) * u;
        }
    }
}





/* this routine loads particle data from Gadget's default
 * binary file format. (A snapshot may be distributed
 * into multiple files.
 */
int load_snapshot(char *fname, int files)
{
    FILE *fd;
    char buf[200];
    int i, j, k, dummy, ntot_withmasses;
    int t, n, off, pc, pc_new, pc_sph;
    
#define SKIP fread(&dummy, sizeof(dummy), 1, fd);
    
    for(i = 0, pc = 1; i < files; i++, pc = pc_new)
    {
        if(files > 1)
            sprintf(buf, "%s.%d", fname, i);
        else
            sprintf(buf, "%s", fname);
        
        if(!(fd = fopen(buf, "r")))
        {
            printf("can't open file `%s`\n", buf);
            exit(0);
        }
        
        printf("reading `%s' ...\n", buf);
        fflush(stdout);
        
        fread(&dummy, sizeof(dummy), 1, fd);
        fread(&header1, sizeof(header1), 1, fd);
        fread(&dummy, sizeof(dummy), 1, fd);
        
        if(files == 1)
        {
            for(k = 0, NumPart = 0, ntot_withmasses = 0; k < 6; k++)
                NumPart += header1.npart[k];
            Ngas = header1.npart[0];
        }
        else
        {
            for(k = 0, NumPart = 0, ntot_withmasses = 0; k < 6; k++)
                NumPart += header1.npartTotal[k];
            Ngas = header1.npartTotal[0];
        }
        
        for(k = 0, ntot_withmasses = 0; k < 6; k++)
        {
            if(header1.mass[k] == 0)
                ntot_withmasses += header1.npart[k];
        }
        
        if(i == 0)
            allocate_memory();
        
        SKIP;
        for(k = 0, pc_new = pc; k < 6; k++)
        {
            for(n = 0; n < header1.npart[k]; n++)
            {
                printf("read pos %i\n",pc_new);
                fread(&P[pc_new].Pos[0], sizeof(float), 3, fd);
                pc_new++;
            }
        }
        SKIP;
        
        SKIP;
        for(k = 0, pc_new = pc; k < 6; k++)
        {
            for(n = 0; n < header1.npart[k]; n++)
            {
                printf("read vel %i\n",pc_new);
                
                fread(&P[pc_new].Vel[0], sizeof(float), 3, fd);
                pc_new++;
            }
        }
        SKIP;
        
        
        SKIP;
        for(k = 0, pc_new = pc; k < 6; k++)
        {
            for(n = 0; n < header1.npart[k]; n++)
            {
                fread(&Id[pc_new], sizeof(int), 1, fd);
                pc_new++;
            }
        }
        SKIP;
        
        
        if(ntot_withmasses > 0)
            SKIP;
        for(k = 0, pc_new = pc; k < 6; k++)
        {
            for(n = 0; n < header1.npart[k]; n++)
            {
                P[pc_new].Type = k;
                
                if(header1.mass[k] == 0)
                    fread(&P[pc_new].Mass, sizeof(float), 1, fd);
                else
                    P[pc_new].Mass = header1.mass[k];
                pc_new++;
            }
        }
        if(ntot_withmasses > 0)
            SKIP;
        
        
        if(header1.npart[1] > 0)
        {
            SKIP;
            for(n = 0, pc_sph = pc; n < header1.npart[0]; n++)
            {
                fread(&P[pc_sph].U, sizeof(float), 1, fd);
                pc_sph++;
            }
            SKIP;
            
            SKIP;
            for(n = 0, pc_sph = pc; n < header1.npart[0]; n++)
            {
                fread(&P[pc_sph].Rho, sizeof(float), 1, fd);
                pc_sph++;
            }
            SKIP;
            
            if(header1.flag_cooling)
            {
                SKIP;
                for(n = 0, pc_sph = pc; n < header1.npart[0]; n++)
                {
                    fread(&P[pc_sph].Ne, sizeof(float), 1, fd);
                    pc_sph++;
                }
                SKIP;
            }
            else
                for(n = 0, pc_sph = pc; n < header1.npart[0]; n++)
                {
                    P[pc_sph].Ne = 1.0;
                    pc_sph++;
                }
        }
        
        fclose(fd);
    }
    
    Time = header1.time;
    Redshift = header1.redshift;
    
}




/* this routine allocates the memory for the
 * particle data.
 */
int allocate_memory(void)
{
    printf("allocating memory...\n");
    
    if(!(P = malloc(NumPart * sizeof(struct particle_data))))
    {
        fprintf(stderr, "failed to allocate memory.\n");
        exit(0);
    }
    
    P--;                /* start with offset 1 */
    
    
    if(!(Id = malloc(NumPart * sizeof(int))))
    {
        fprintf(stderr, "failed to allocate memory.\n");
        exit(0);
    }
    
    Id--;                /* start with offset 1 */
    
    printf("allocating memory...done\n");
}




/* This routine brings the particles back into
 * the order of their ID's.
 * NOTE: The routine only works if the ID's cover
 * the range from 1 to NumPart !
 * In other cases, one has to use more general
 * sorting routines.
 */
int reordering(void)
{
    int i, j;
    int idsource, idsave, dest;
    struct particle_data psave, psource;
    
    
    printf("reordering....\n");
    
    for(i = 1; i <= NumPart; i++)
    {
        if(Id[i] != i)
        {
            psource = P[i];
            idsource = Id[i];
            dest = Id[i];
            
            do
            {
                psave = P[dest];
                idsave = Id[dest];
                
                P[dest] = psource;
                Id[dest] = idsource;
                
                if(dest == i)
                    break;
                
                psource = psave;
                idsource = idsave;
                
                dest = idsource;
            }
            while(1);
        }
    }
    
    printf("done.\n");
    
    Id++;
    free(Id);
    
    printf("space for particle ID freed\n");
}

You need to compile it going inside the folder and type:

gcc -o read_snapshot_to_txt read_snapshot_to_txt.c

Then you can run your compiled file to extract data from binary and create your own txts files:

./read_snapshot_to_txt 2

The parameter specifies the number of snapshots you want to convert.