First commit

# Object files

*.o

*.ko

*.obj

*.elf

# Precompiled Headers

*.gch

*.pch

# Libraries

*.lib

*.a

*.la

*.lo

# Shared objects (inc. Windows DLLs)

*.dll

*.so

*.so.*

*.dylib

# Executables

*.exe

*.out

*.app

*.i*86

*.x86_64

*.hex

# Debug files

*.dSYM/

## Config file for P-Fkp ##

REDSHIFT 0.188            // redshift for IC cooling

## Input ##

#I_TXT                  // Input from ASCII fil

I_GADGET2              // Gadget Format 2 Input

#I_NONE                  // Dummy Input 

## Output ##    

#O_TXT                   // Output to ASCII file

O_BINARY               // Output to Binary Data

## OPTIONS ##

READ_MORE_INFO			// Get more velocities from snap

VARIABLE_TIMESTEPS

#COMPRESSION             // use curve & binary spectrum compression

#COMPRESSION_INTERNAL    // use compression to store spectrum internally (saves RAM)

#COMPRESSION_RANGE_WARNING // Warn if binary compression is out of range

#TIMEVARIABLECMB

## PHYSICS MODULES ##

COMPUTE_DPP             // Make sure to enable a DPP switch as well

#DPP_BRUNETTI_07         // Reacceleration from Brunetti+ 07

#HP_BRUNETTI_07         

#Q_CASSANO_05            // use this to compare with Cassano datA          

HP_DONNERT_13

DPP_DONNERT_13

Q_SHOCK_PRIMARIES            // Primary injection at shocks

#Q_SHOCK_WITH_HADRONIC_BACKGROUND

#Q_SIMPLE_SECONDARIES   // Injection from simple cut power-law

#Q_BRUNETTI_05          // High energy approx. Secondary Brunetti+ 05

#Q_ZERO                 // Zero Injection

T_ZERO                 // Zero Escape

N0_ZERO                // Zero Initial Spectrum

## Hard Sphere Equation ##

#DPP_HARD_SPHERE         // Test Case (Park & Petrosian 1996

#HP_HARD_SPHERE 5       // set according to subcase: 0..3, 4, 5

#Q_HARD_SPHERE

#T_HARD_SPHERE

#N0_HARD_SPHERE         // Zero Initial Spectrum

## Test Case Cassano 05

#I_GADGET2              // Gadget Format 2 Input

#O_TXT

#O_BINARY               // Output to Binary Data

#DPP_BRUNETTI_07         // Reacceleration from Brunetti+ 07

#HP_BRUNETTI_07

#Q_CASSANO_05            // use this to compare with Cassano data

#T_ZERO                 // Zero Escape

#N0_ZERO                // Zero Initial Spectrum

## Test Case Expansion

#DPP_EXPANSION

#HP_EXPANSION

#Q_ZERO

#T_ZERO

#N0_POWERLAW

## P-Fkp Makefile ## 

SHELL = /bin/bash

# Target Computer

ifndef SYSTYPE

	SYSTYPE	:= $(shell hostname)

endif

#Std systype

CC		= mpicc

OPTIMIZE 	= -O3 -m64

MPI_LIBS 	=

MPI_INCL 	= 

GSL_INCL 	= -I/homes/gheller/Pleiadi/gsl-2.7.1/build/include/ -I/homes/gheller/Pleiadi/hdf5-1.12.1/build/include/

GSL_LIBS 	= -L/homes/gheller/Pleiadi/gsl-2.7.1/build/lib/

HDF5_LIBS       = -L/homes/gheller/Pleiadi/hdf5-1.12.1/build/lib/

PGPLOT_INCL = 

PGPLOT_LIBS =

ifeq ($(SYSTYPE),MSI)

CC           =  mpicc -cc=icc

OPTIMIZE     = -O0 -Wall -g -xHost

MPI_LIBS     = -lmpi -lrt -L$(LD_LIBRARY_PATH)

MPI_INCL     = -I$(INCLUDE)

GSL_INCL     =

GSL_LIBS     =

PGPLOT_INCL  =

PGPLOT_LIBS  =

endif

ifeq ($(SYSTYPE),para33.strw.leidenuniv.nl)

CC           =  mpicc

OPTIMIZE     = -O3  -m64 -g -march=native -flto -fwhole-program

MPI_LIBS     = -lmpich -lrt

MPI_INCL     =

GSL_INCL     =

GSL_LIBS     =

PGPLOT_INCL  =

PGPLOT_LIBS  =

endif

ifeq ($(SYSTYPE),mach64.ira.inaf.it)

CC           =  mpicc

OPTIMIZE     = -O3  -m64 #-g -march=bdver1 -mtune=native -mprefer-avx128 -mieee-fp -minline-all-stringops -fprefetch-loop-arrays --param prefetch-latency=300 -funroll-all-loops -flto -fwhole-program

MPI_LIBS     =  #-lmpich

MPI_INCL     =

GSL_INCL     = -I/homes/gheller/Work/gsl-2.6/build/include -I/homes/gheller/mach64/hdf5-1.12.0/build/include

GSL_LIBS     = -L/homes/gheller/Work/gsl-2.6/build/lib

HDF5_LIBS    = -L/homes/gheller/mach64/hdf5-1.12.0/build/lib

PGPLOT_INCL  =

PGPLOT_LIBS  =

endif

ifeq ($(SYSTYPE),"DARWIN")

CC           =  mpicc

OPTIMIZE     = -O3 -std=c99 -mtune=native -march=corei7

MPI_LIBS     = -lmpich -L/Users/julius/Devel/lib

MPI_INCL     = -I/Users/julius/Devel/include

GSL_INCL     = 

GSL_LIBS     = 

PGPLOT_INCL  =

PGPLOT_LIBS  =

endif

ifeq ($(SYSTYPE),"MPA")

CC           =  mpicc

OPTIMIZE     = -O3 -Wall -g

MPI_LIBS     = -mpi 

MPI_INCL     = 

GSL_INCL     = -I/afs/mpa/home/jdonnert/Libs/@sys/include

GSL_LIBS     = -L/afs/mpa/home/jdonnert/Libs/@sys/lib

PGPLOT_INCL  = -I/usr/local/include

PGPLOT_LIBS  = -L/usr/common/pdsoft/appl/pgplot/ -L/usr/X11R6/lib

endif

ifeq ($(SYSTYPE),"RZG-ODIN") 

CC       =  mpicc

OPTIMIZE = -O3  -Wall -g

MPI_LIBS = 

MPI_INCL = 

GSL_INCL = -I/u/jdonnert/Libs/include

GSL_LIBS = -L/u/jdonnert/Libs/lib

endif

# Target

EXEC 	= P-Fkp		

# Sources 

SRCDIR	= src/

OBJFILES= io/io.o \

		  		io/txt.o \

				io/gadget.o \

				io/hdf5_reader.o \

				io/binary.o \

		  		modules/zero.o	\

				modules/secondaries_brunetti.o \

				modules/simple_secondaries.o \

				modules/hard_sphere.o \

				modules/brunetti_07.o \

				modules/donnert_13.o \

				modules/expansion.o \

				modules/shock_primaries.o \

		 common.o main.o setup.o solver.o cosmo.o timestep.o sort.o  \

		 sort_particles.o timing.o compress.o unit.o print_settings.o

OBJS	= $(addprefix $(SRCDIR),$(OBJFILES))

INCLFILES	= common.h cosmo.h compress.h unit.h constants.h config.h sort.h \

			  timing.h io/io.h \

			  		io/gadget.h \

			  modules/modules.h	\

			  ../Makefile

INCL	= $(addprefix $(SRCDIR),$(INCLFILES))

LIBS = -lm -lgsl -lgslcblas  $(MPI_LIBS) $(GSL_LIBS) $(HDF5_LIBS) -lhdf5

CFLAGS = -std=c99  -fopenmp $(OPTIMIZE) $(GSL_INCL)

$(EXEC)	: $(OBJS)

	$(CC) $(CFLAGS)  $(OBJS)  $(LIBS) -o $(EXEC)

	ctags -R src/*.[ch]

$(OBJS)	: $(INCL)

$(SRCDIR)config.h : Config 

	sed '/^#/d; /^$$/d; s/^/#define /g' Config > $(SRCDIR)config.h

.ONESHELL:

$(SRCDIR)print_settings.c : Config

	echo '#include "common.h"' >  $(SRCDIR)print_settings.c

	echo 'void print_compile_time_settings(){' >> $(SRCDIR)print_settings.c

	echo 'rprintf("Compiled with : \n"' >> $(SRCDIR)print_settings.c

	sed '/^#/d; /^$$/d; s/^/"   /g; s/$$/ \\n"/g;' Config >>  $(SRCDIR)print_settings.c

	echo ');}' >> $(SRCDIR)print_settings.c

clean	: 

	rm -f  $(OBJS) $(EXEC)

#include "common.h"

struct Parameters Param;

struct ParallelInfos ThisTask;

struct SnapProperties Snap = {0};

struct Particle_Data *P=NULL, *LastP=NULL;

struct Gas_Data *SphP=NULL, *LastSphP=NULL;

struct ParallelInfos ThisTask = { 0, 0, 0, {0} };

#pragma omp threadprivate(Omp)

struct OpenMP_infos Omp = { 0 };

short *TimeBins = NULL;

float TimeBase = 0; 

float *timesteps = NULL;

size_t LowIdx = 0, HighIdx = 0;

double p[N_SPEC_BINS] = { 0 }, q[N_SPEC_BINS] = { 0 };

double dp[N_SPEC_BINS] = { 0 }, dq[N_SPEC_BINS] = { 0 }; 

#pragma omp threadprivate(Flag)

int Flag = 0, Ipart = 2343828;

/* 

 * Memory Management 

 */

void *malloc_info(size_t size, const char* file

        , const char* func, const int line)

{

    char errmsg[MAXLINELENGTH];

	void *result = malloc(size);

    sprintf(errmsg, "Allocation failed, <%zu> bytes at %s in %s line %d \n"

        ,size, file, func, line);

    Assert(result != NULL, errmsg);

	return (result);

}

void *realloc_info(void *ptr, size_t size, const char* file

        , const char* func, const int line)

{

    char errmsg[MAXLINELENGTH];

	void *result = realloc(ptr, size);

    sprintf(errmsg, "Reallocation failed, <%zu> bytes at %s in %s line %d \n"

        ,size, file, func, line);

    Assert(result != NULL && size != 0, errmsg);

	return (result);

}

void free_info(void *ptr) 

{

    Assert(ptr != NULL, "Can't free a NULL pointer");

    free(ptr);

    return;

}

/* 

 * Error Handling, we use variable arguments to be able

 * to print more informative error messages 

 */

void Assert_Info(const char *func, const char *file, int line,

		int64_t expr, const char *errmsg, ...)

{

    if (expr)

        return;

	va_list varArgList;

	va_start(varArgList, errmsg);

	/* we fucked up, tell them */

    fprintf(stderr, 

			"\nERROR Task %d: In file %s, function %s(), line %d :\n\n	", 

			ThisTask.Rank, file, func, line);

	vfprintf(stderr, errmsg, varArgList); 

	fprintf(stderr, "\n\n"); 

	fflush(stderr);

	va_end(varArgList);

    MPI_Abort(MPI_COMM_WORLD, -1); // finish him ...

    exit(EXIT_FAILURE); // ... fatality

    return;

}

/* 

 * Particle Metrics 

 */

extern void Reallocate_P(long long partTotal, 

        long long nPart[N_PART_TYPES], int sign)

{

	if (partTotal == 0)

		return;

	ThisTask.PartTotal += sign * partTotal;

	for (int type = 0; type < N_PART_TYPES; type++) {

		ThisTask.Npart[type] += sign * nPart[type];

        Assert(ThisTask.Npart[type] >= 0, 

                "Can't allocate negative particles");

	}

	P = Realloc((void *)P, sizeof(*P)*ThisTask.PartTotal);

	SphP = Realloc((void *)SphP, sizeof(*SphP)*ThisTask.Npart[0]);

	return;

}

extern void Shift_Particle_Data()

{

	rprintf("Shifting Particle data ...");

    Snap.LastSnapNum = Snap.SnapNum;

    Snap.LastTime = Snap.Time;

    Snap.LastRedshift = Snap.Redshift;

    struct Particle_Data *tmp_P = LastP;    // save the pointer, can be NULL

    LastP = P;

    P = tmp_P;

    struct Gas_Data *tmp_SphP = LastSphP;

    LastSphP = SphP;

    SphP = tmp_SphP;

    if (P != NULL)    // clear memory

        memset(P, 0, ThisTask.Npart[0] * sizeof(*P));

    if (SphP != NULL)

        memset(SphP, 0, ThisTask.Npart[0] * sizeof(*SphP));

    ThisTask.Npart[0] = 0;  // these are set at readin 

    ThisTask.PartTotal = 0;

    Snap.Npart[0] = 0;

    Snap.PartTotal = 0;

	rprintf("done\n");

    return;

}